Added macro @law to simplify the definition of constitutive laws

src/FESpaces/FEBases.jl

@@ -10,6 +10,7 @@ import Gridap: symmetric_gradient
 import Gridap: inner
 import Base: +, -, *
 import Gridap: restrict
+import Gridap: Triangulation
 
 struct FEBasis{B<:CellBasis,T<:Triangulation}
   cellbasis::B
@@ -54,6 +55,8 @@ function inner(a::FEBasis,b::FEBasis)
   varinner(a.cellbasis,b.cellbasis)
 end
 
+Triangulation(a::FEBasis) = a.trian
+
 function CellBasis(
   trian::Triangulation{D,Z},
   fun::Function,

src/FESpaces/FEFunctions.jl

@@ -13,8 +13,10 @@ import Gridap: evaluate, gradient, return_size
 import Gridap: symmetric_gradient
 import Gridap: reindex
 import Gridap: restrict
+import Gridap: Triangulation
 import Base: size
 import Base: getindex
+import Base: +, -, *
 
 import Base: zero
 
@@ -41,6 +43,8 @@ diri_dofs(f::FEFunction) = f.diri_dofs
 
 FESpace(f::FEFunction) = f.fespace
 
+Triangulation(f::FEFunction) = Triangulation(f.fespace)
+
 value_type(::FEFunction{D,Z,T}) where {D,Z,T} = T
 
 """
@@ -80,9 +84,26 @@ end
 
 evaluate(f::FEFunction{D,Z},q::CellPoints{Z}) where {D,Z} = evaluate(f.cellfield,q)
 
-gradient(f::FEFunction) = gradient(f.cellfield)
+function gradient(f::FEFunction)
+  g = gradient(f.cellfield)
+  IndexCellFieldWithFEFunction(g,f)
+end
 
-symmetric_gradient(f::FEFunction) = symmetric_gradient(f.cellfield)
+function symmetric_gradient(f::FEFunction)
+  g = symmetric_gradient(f.cellfield)
+  IndexCellFieldWithFEFunction(g,f)
+end
+
+for op in (:+, :-, :*)
+  @eval begin
+    function ($op)(a::FEFunction,b::CellField)
+      IndexCellFieldWithFEFunction($op(a.cellfield,b),a)
+    end
+    function ($op)(a::CellField,b::FEFunction)
+      IndexCellFieldWithFEFunction($op(a,b.cellfield),b)
+    end
+  end
+end
 
 return_size(f::FEFunction,s::Tuple{Int}) = return_size(f.cellfield,s)
 
@@ -114,4 +135,46 @@ restrict(uh::FEFunction,trian::BoundaryTriangulation) = restrict(uh.cellfield,tr
 
 restrict(uh::FEFunction,trian::SkeletonTriangulation) = restrict(uh.cellfield,trian)
 
+"""
+Type used to represent the result of operations involving FEfunction objects.
+It keeps track of the corresponding FEFunction.
+"""
+struct IndexCellFieldWithFEFunction{
+  Z,C<:IndexCellField,F<:FEFunction,R} <: IndexCellValue{R,1}
+  cellfield::C
+  fefunction::F
+end
+
+function IndexCellFieldWithFEFunction(
+  cellfield::IndexCellField,
+  fefunction::FEFunction{D,Z}) where {D,Z}
+
+  C = typeof(cellfield)
+  F = typeof(fefunction)
+  R = eltype(cellfield)
+  IndexCellFieldWithFEFunction{Z,C,F,R}(cellfield,fefunction)
+end
+
+function evaluate(f::IndexCellFieldWithFEFunction{Z},q::CellPoints{Z}) where Z
+  evaluate(f.cellfield,q)
+end
+
+function gradient(f::IndexCellFieldWithFEFunction)
+  g = gradient(f.cellfield)
+  IndexCellFieldWithFEFunction(g,f.fefunction)
+end
+
+function symmetric_gradient(f::IndexCellFieldWithFEFunction)
+  g = symmetric_gradient(f.cellfield)
+  IndexCellFieldWithFEFunction(g,f.fefunction)
+end
+
+return_size(f::IndexCellFieldWithFEFunction,s::Tuple{Int}) = return_size(f.cellfield,s)
+
+getindex(f::IndexCellFieldWithFEFunction,i::Integer) = f.cellfield[i]
+
+size(f::IndexCellFieldWithFEFunction) = (length(f.cellfield),)
+
+Triangulation(f::IndexCellFieldWithFEFunction) = Triangulation(f.fefunction)
+
 end # module

src/Integration/Triangulations.jl

@@ -10,6 +10,7 @@ export CellRefFEs
 export CartesianTriangulation
 export test_triangulation
 export ncells
+export @law
 import Gridap: CellQuadrature
 import Gridap: CellPoints
 import Gridap: CellBasis
@@ -103,6 +104,29 @@ end
 
 _setup_cell_field(f) = cellnewaxis(f,dim=1)
 
+macro law(fundef)
+  s = "The @law macro is only allowed in function definitions"
+  @assert isa(fundef,Expr) s
+  @assert fundef.head in (:(=), :function) s
+  funname = fundef.args[1].args[1]
+  nargs = length(fundef.args[1].args)-1
+  if nargs == 0
+    x =  Symbol[]
+  else
+    x = fundef.args[1].args[3:end]
+  end
+  fundef2 = quote
+    function $(funname)($(x...))
+      trian = Triangulation($(x[1]))
+      CellBasis(trian,$(funname),$(x...))
+    end
+  end
+  quote
+    $(esc(fundef))
+    $(esc(fundef2))
+  end
+end
+
 """
 Factory function to create CellQuadrature objects in a convenient way
 """

src/MultiField/MultiFEBases.jl

@@ -10,12 +10,15 @@ import Base: length, getindex
 import Gridap.FESpaces: FEBasis
 import Gridap: CellBasis
 import Gridap: restrict
+import Gridap: Triangulation
 
 struct FEBasisWithFieldId{B<:FEBasis}
   febasis::B
   fieldid::Int
 end
 
+Triangulation(a::FEBasisWithFieldId) = Triangulation(a.febasis)
+
 function CellBasis(
   trian::Triangulation{D,Z},
   fun::Function,

test/FESpacesTests/FEBasesTests.jl

@@ -10,6 +10,8 @@ fespace = ConformingFESpace(Float64,model,order,diritag)
 
 bh = FEBasis(fespace)
 
+@test isa(Triangulation(bh),Triangulation)
+
 ufun(x) = x[1]
 
 uh = interpolate(fespace,ufun)

test/FESpacesTests/FEFunctionsTests.jl

@@ -5,6 +5,8 @@ using Gridap
 using Test
 using Gridap
 
+using Gridap.FEFunctions: IndexCellFieldWithFEFunction
+
 model = CartesianDiscreteModel(domain=(0.0,1.0,0.0,1.0), partition=(4,4))
 order = 1
 diritag = "boundary"
@@ -14,14 +16,36 @@ ufun(x) = x[1]
 
 uh = interpolate(fespace,ufun)
 
+@test isa(Triangulation(uh),Triangulation)
+
 tags = [7,6]
 btrian = BoundaryTriangulation(model,tags)
-quad = CellQuadrature(btrian,order=0)
+bquad = CellQuadrature(btrian,order=0)
 
 buh = restrict(uh,btrian)
 
-cn = integrate(buh,btrian,quad)
+cn = integrate(buh,btrian,bquad)
 @test isa(cn,CellNumber)
 _ = collect(cn)
 
+cf = IndexCellFieldWithFEFunction(uh.cellfield,uh)
+
+trian = Triangulation(model)
+quad = CellQuadrature(trian,order=2)
+q = coordinates(quad)
+v = collect(evaluate(cf,q))
+g = collect(evaluate(gradient(cf),q))
+
+@test isa(∇(cf),IndexCellFieldWithFEFunction)
+
+test_index_cell_field(cf,q,v,g)
+
+@test isa(∇(uh),IndexCellFieldWithFEFunction)
+@test isa(Triangulation(∇(uh)),Triangulation)
+
+u = CellField(trian,ufun)
+
+e = u - uh
+@test isa(e,IndexCellFieldWithFEFunction)
+
 end # module

test/FESpacesTests/VectorValuedFEOperatorsTests.jl

@@ -81,7 +81,7 @@ const λ = (E*ν)/((1+ν)*(1-2*ν))
 const μ = E/(2*(1+ν))
 
 # Constitutive law
-σfun(x,ε) = λ*tr(ε)*one(ε) + 2*μ*ε
+@law σ(x,ε) = λ*tr(ε)*one(ε) + 2*μ*ε
 
 # Define manufactured functions
 ufun(x) = VectorValue(x[1] + x[2],x[1])
@@ -108,7 +108,6 @@ quad = CellQuadrature(trian,order=2)
 
 # Terms in the volume
 bfield = CellField(trian,bfun)
-σ(ε) = CellBasis(trian,σfun,ε)
 a_elast(v,u) = inner( ε(v), σ(ε(u)) )
 b(v) = inner(v,bfield)
 t_Ω = AffineFETerm(a_elast,b,trian,quad)

test/MultiFieldTests/MultiFEBasesTests.jl

@@ -42,6 +42,8 @@ V = MultiFESpace([V1,V2])
 u = FEBasis(U)
 v = FEBasis(V)
 
+@test isa(Triangulation(u[1]),Triangulation)
+
 a(v,u) = inner(∇(v[1]),∇(u[1])) + inner(v[1],u[2]) + inner(∇(v[2]),∇(u[2]))
 
 b(v) = inner(v[1],b1field) + inner(v[2],b2field)