Added div, curl, and trace operators

Author: fverdugo

src/FESpaces/FEBases.jl

@@ -7,6 +7,9 @@ export FEBasis
 import Gridap: CellBasis
 import Gridap: gradient
 import Gridap: symmetric_gradient
+import Base: div
+import Gridap: trace
+import Gridap: curl
 import Gridap: inner
 import Base: +, -, *
 import Gridap: restrict
@@ -23,7 +26,7 @@ function FEBasis(fespace::FESpace)
   FEBasis(b,trian)
 end
 
-for op in (:+, :-, :(gradient),:(symmetric_gradient))
+for op in (:+,:-,:(gradient),:(symmetric_gradient),:(div),:(trace),:(curl))
   @eval begin
     function ($op)(a::FEBasis)
       FEBasis($op(a.cellbasis),a.trian)

src/FESpaces/FEFunctions.jl

@@ -12,6 +12,9 @@ import Gridap: value_type
 
 import Gridap: evaluate, gradient, return_size
 import Gridap: symmetric_gradient
+import Base: div
+import Gridap: trace
+import Gridap: curl
 import Gridap: reindex
 import Gridap: restrict
 import Gridap: Triangulation
@@ -85,16 +88,14 @@ end
 
 evaluate(f::FEFunction{D,Z},q::CellPoints{Z}) where {D,Z} = evaluate(f.cellfield,q)
 
-function gradient(f::FEFunction)
-  g = gradient(f.cellfield)
-  trian = Triangulation(f)
-  _attach_triangulation(g,trian)
-end
-
-function symmetric_gradient(f::FEFunction)
-  g = symmetric_gradient(f.cellfield)
-  trian = Triangulation(f)
-  _attach_triangulation(g,trian)
+for op in (:+,:-,:(gradient),:(symmetric_gradient),:(div),:(trace),:(curl))
+  @eval begin
+    function ($op)(a::FEFunction)
+      g = $op(a.cellfield)
+      trian = Triangulation(a)
+      _attach_triangulation(g,trian)
+    end
+  end
 end
 
 for op in (:+, :-, :*)

src/Fields/CellFieldsOperations.jl

@@ -17,7 +17,10 @@ export lincomb
 export attachgeomap
 export compose
 export symmetric_gradient
+export curl
 export ε
+import Base: div
+import Gridap: trace
 import Gridap: evaluate
 import Gridap: gradient
 import Gridap: reindex
@@ -72,6 +75,31 @@ end
 
 const ε = symmetric_gradient
 
+function trace(f::CellFieldLike{D,T}) where {D,T<:TensorValue}
+  apply(trace,f,broadcast=true)
+end
+
+function div(f::CellFieldLike{D,T}) where {D,T<:VectorValue}
+  g = gradient(f)
+  trace(g)
+end
+
+function curl(f::CellFieldLike{D,T}) where {D,T<:VectorValue}
+  g = gradient(f)
+  apply(_curl_kernel,g,broadcast=true)
+end
+
+function _curl_kernel(∇u::TensorValue{2})
+  ∇u[1,2] - ∇u[2,1]
+end
+
+function _curl_kernel(∇u::TensorValue{3})
+  c1 = ∇u[2,3] - ∇u[3,2]
+  c2 = ∇u[3,1] - ∇u[1,3]
+  c3 = ∇u[1,2] - ∇u[2,1]
+  VectorValue(c1,c2,c3)
+end
+
 function _compose(::Val{true},f,u...)
   fgrad = gradient(f)
   v = apply(f,u...,broadcast=true)

src/Integration/Triangulations.jl

@@ -18,10 +18,14 @@ import Gridap: CellGeomap
 import Gridap: CellField
 import Gridap: evaluate, gradient, return_size
 import Gridap: symmetric_gradient
+import Base: div
+import Gridap: trace
+import Gridap: curl
 import Gridap: reindex
 import Base: IndexStyle
 import Base: size
 import Base: getindex
+import Base: +, -, *
 
 """
 Minimal interface for a mesh used for numerical integration
@@ -275,14 +279,13 @@ function evaluate(f::IndexCellFieldWithTriangulation{Z},q::CellPoints{Z}) where
   evaluate(f.cellfield,q)
 end
 
-function gradient(f::IndexCellFieldWithTriangulation)
-  g = gradient(f.cellfield)
-  IndexCellFieldWithTriangulation(g,f.trian)
-end
-
-function symmetric_gradient(f::IndexCellFieldWithTriangulation)
-  g = symmetric_gradient(f.cellfield)
-  IndexCellFieldWithTriangulation(g,f.trian)
+for op in (:+,:-,:(gradient),:(symmetric_gradient),:(div),:(trace),:(curl))
+  @eval begin
+    function ($op)(a::IndexCellFieldWithTriangulation)
+      g = $op(a.cellfield)
+      IndexCellFieldWithTriangulation(g,a.trian)
+    end
+  end
 end
 
 return_size(f::IndexCellFieldWithTriangulation,s::Tuple{Int}) = return_size(f.cellfield,s)

src/MultiField/MultiFEBases.jl

@@ -4,6 +4,10 @@ using Gridap
 using Gridap.Helpers
 
 import Gridap: gradient
+import Gridap: symmetric_gradient
+import Base: div
+import Gridap: trace
+import Gridap: curl
 import Gridap: inner
 import Base: +, -, *
 import Base: length, getindex
@@ -29,7 +33,7 @@ function CellBasis(
   FEBasisWithFieldId(febasis,b.fieldid)
 end
 
-for op in (:+, :-, :(gradient))
+for op in (:+,:-,:(gradient),:(symmetric_gradient),:(div),:(trace),:(curl))
   @eval begin
     function ($op)(a::FEBasisWithFieldId)
       FEBasisWithFieldId($op(a.febasis),a.fieldid)

test/FESpacesTests/FEBasesTests.jl

@@ -60,6 +60,8 @@ ufun(x) = VectorValue(x[2],x[1])
 uh = interpolate(fespace,ufun)
 bh = FEBasis(fespace)
 @test isa(ε(bh),FEBasis)
+@test isa(div(bh),FEBasis)
+@test isa(curl(bh),FEBasis)
 @test isa(inner(ε(bh),ε(bh)),CellMap{Point{2},1,Float64,3})
 @test isa(inner(ε(bh),ε(uh)),CellMap{Point{2},1,Float64,2})
 

test/FESpacesTests/FEFunctionsTests.jl

@@ -51,4 +51,14 @@ u = CellField(trian,ufun)
 e = u - uh
 @test isa(e,IndexCellFieldWithTriangulation)
 
+T = VectorValue{2,Float64}
+fespace = ConformingFESpace(T,model,order,diritag)
+ufun(x) = VectorValue(x[2],x[1])
+uh = interpolate(fespace,ufun)
+
+@test isa(div(uh),IndexCellFieldWithTriangulation)
+@test isa(curl(uh),IndexCellFieldWithTriangulation)
+
+
+
 end # module

test/FieldsTests/CellFieldsOperationsTests.jl

@@ -189,5 +189,23 @@ test_iter_cell_field_without_grad(sg,cp,r)
 
 @test ε == symmetric_gradient
 
+gr = ∇(cf)
+sg = trace(gr)
+
+ri = fill(6.0,4)
+r = fill(ri,l)
+
+test_iter_cell_field_without_grad(sg,cp,r)
+
+sg = div(cf)
+
+test_iter_cell_field_without_grad(sg,cp,r)
+
+sg = curl(cf)
+
+ri = fill(0.0,4)
+r = fill(ri,l)
+test_iter_cell_field_without_grad(sg,cp,r)
+
 end # module
 

test/MultiFieldTests/MultiFEBasesTests.jl

@@ -2,6 +2,7 @@ module MultiFEBasesTests
 
 using Test
 using Gridap
+using Gridap.MultiFEBases: FEBasisWithFieldId
 
 import Gridap: ∇
 
@@ -78,4 +79,21 @@ mca = integrate(mcm,trian,quad)
 
 @test isa(mca,MultiCellArray{Float64,1})
 
+T = VectorValue{2,Float64}
+order = 1
+diritag = "boundary"
+fespace = ConformingFESpace(T,model,order,diritag)
+
+V1 = TestFESpace(fespace)
+V2 = V1
+
+V = MultiFESpace([V1,V2])
+
+bh = FEBasis(V)
+
+@test isa(∇(bh[1]),FEBasisWithFieldId)
+@test isa(ε(bh[1]),FEBasisWithFieldId)
+@test isa(div(bh[1]),FEBasisWithFieldId)
+@test isa(curl(bh[1]),FEBasisWithFieldId)
+
 end # module MultiFEBasesTests