Misc changes in order to support RT FEs in general manifolds (i.e., num_cell_dims(model) < num_point_dims(model))

src/FESpaces/DivConformingFESpaces.jl

@@ -137,7 +137,14 @@ function return_cache(::typeof(_transform_rt_dof_basis),
                                  get_face_moments(dofs)
   db = MomentBasedDofBasis(nodes,nf_moments,nf_nodes)
 
-  face_moments = deepcopy(nf_moments)
+  # Determine the type of the elements
+  # in the range of Field. This is the
+  # element type of the face_moments array.
+  D = num_dims(reffe)
+  pt = Point{D,et}
+  T=return_type(phi,zero(pt))
+  face_moments = [ similar(i,T)  for i in nf_moments ]
+
   Jt_q_cache = return_cache(∇(phi),db.nodes)
   cache = (db.nodes, db.face_nodes, nf_moments, face_moments, Jt_q_cache)
   cache
@@ -160,9 +167,9 @@ function evaluate!(cache,
       num_qpoints, num_moments = size(moments)
       for i in 1:num_qpoints
         Jt_q_i = Jt_q[nf_nodes[face][i]]
-        change = sign * meas(Jt_q_i) * inv(transpose(Jt_q_i))
+        change = sign * meas(Jt_q_i) * pinvJt(Jt_q_i)
         for j in 1:num_moments
-          face_moments[face][i,j] = nf_moments[face][i,j] ⋅ change
+          face_moments[face][i,j] = change ⋅ nf_moments[face][i,j]
         end
       end
     end

src/Fields/Fields.jl

@@ -28,7 +28,7 @@ import LinearAlgebra: det, inv, transpose, tr, cross
 import LinearAlgebra: ⋅, dot
 
 import Base: +, -, *, /
-import Gridap.TensorValues: ⊗, ⊙, symmetric_part, outer
+import Gridap.TensorValues: ⊗, ⊙, symmetric_part, outer, meas
 
 import Gridap.Arrays: IndexStyle
 import Gridap.Arrays: return_cache

src/Fields/FieldsInterfaces.jl

@@ -348,7 +348,7 @@ return_value(op::Broadcasting{<:Operation},x::Field...) = OperationField(op.f.op
 
 # Define some well known operations
 
-for op in (:+,:-,:*,:⋅,:⊙,:⊗,:inv,:det,:tr,:grad2curl,:symmetric_part,:transpose)
+for op in (:+,:-,:*,:⋅,:⊙,:⊗,:inv,:det,:meas,:pinvJt,:tr,:grad2curl,:symmetric_part,:transpose)
   @eval ($op)(a::Field...) = Operation($op)(a...)
 end
 

src/ReferenceFEs/RaviartThomasRefFEs.jl

@@ -149,8 +149,8 @@ function _RT_face_moments(p, fshfs, c_fips, fcips, fwips,phi)
   # Hence, J = transpose(grad(phi))
 
   Jt = fill(∇(phi),nc)
-  Jt_inv = lazy_map(Operation(inv),Jt)
-  det_Jt = lazy_map(Operation(det),Jt)
+  Jt_inv = lazy_map(Operation(pinvJt),Jt)
+  det_Jt = lazy_map(Operation(meas),Jt)
   change = lazy_map(*,det_Jt,Jt_inv)
   change_ips = lazy_map(evaluate,change,fcips)
 
@@ -222,10 +222,10 @@ function _RT_cell_values(p,et,order,phi)
   cell_moments = _RT_cell_moments(p, cbasis, ccips, cwips )
 
   # Must scale weights using phi map to get the correct integrals
-  # scaling = det(grad(phi))
+  # scaling = meas(grad(phi))
   Jt = ∇(phi)
-  Jt_inv = inv(Jt)
-  det_Jt = det(Jt)
+  Jt_inv = pinvJt(Jt)
+  det_Jt = meas(Jt)
   change = det_Jt*Jt_inv
   change_ips = evaluate(change,ccips)
 
@@ -399,7 +399,7 @@ function evaluate!(cache,::ContraVariantPiolaMap,
                    Jt::Number,
                    detJ::Number,
                    sign_flip::Bool)
-  ((-1)^sign_flip*v)⋅((1/abs(detJ))*Jt)
+  ((-1)^sign_flip*v)⋅((1/detJ)*Jt)
 end
 
 function evaluate!(cache,
@@ -408,7 +408,7 @@ function evaluate!(cache,
                    phi::Field,
                    sign_flip::AbstractVector{<:Field})
   Jt = ∇(phi)
-  detJ = Operation(det)(Jt)
+  detJ = Operation(meas)(Jt)
   Broadcasting(Operation(k))(v,Jt,detJ,sign_flip)
 end
 
@@ -419,7 +419,7 @@ function lazy_map(
   sign_flip::AbstractArray{<:AbstractArray{<:Field}})
 
   cell_Jt = lazy_map(∇,cell_map)
-  cell_detJ = lazy_map(Operation(det),cell_Jt)
+  cell_detJ = lazy_map(Operation(meas),cell_Jt)
 
   lazy_map(Broadcasting(Operation(k)),cell_ref_shapefuns,cell_Jt,cell_detJ,sign_flip)
 end

test/FESpacesTests/DivConformingFESpacesTests.jl

@@ -66,4 +66,30 @@ el2 = sqrt(sum( ∫( e⋅e )*dΩ ))
 #
 #writevtk(trian,"test",order=3,cellfields=["vh"=>vh])
 
+# Tests on manifold
+
+# Create domain
+domain = (0,1,0,1,0,1)
+cells  = (2,2,2)
+model  = CartesianDiscreteModel(domain,cells)
+
+# Restrict model to cube surface (using new BoundaryDiscreteModel)
+labels = get_face_labeling(model)
+bgface_to_mask = get_face_mask(labels,"boundary",2)
+Γface_to_bgface = findall(bgface_to_mask)
+Dc2Dp3model = BoundaryDiscreteModel(Polytope{2},model,Γface_to_bgface)
+
+order  = 0
+degree = 1
+
+reffe_rt = ReferenceFE(raviart_thomas,Float64,order)
+V  = FESpace(Dc2Dp3model, reffe_rt ; conformity=:HDiv)
+uh = FEFunction(V,rand(num_free_dofs(V)))
+vh = interpolate_everywhere(uh,V)
+
+Tₕ_K = Triangulation(Dc2Dp3model)
+Qₕ_K = CellQuadrature(Tₕ_K,degree)
+e=sqrt(sum(∫((uh-vh)⋅(uh-vh))Qₕ_K))
+@test e < 1.0e-12
+
 end # module