Started adding change of basis

docs/src/dev-notes/pullbacks.md

@@ -71,7 +71,7 @@ I_K(g) = \Sigma(g) \Phi \quad, \quad \Sigma(g) = P^{-1} \hat{\Sigma}(F^{-*}(g))
 
 !!! note
     In [2], Covariant and Contravariant Piola maps preserve exactly (without any sign change) the normal and tangential components of a vector field.
-    I am quite sure that the discrepancy is coming from the fact that the geometrical information in the reference polytope is globaly oriented.
+    I am quite sure that the discrepancy is coming from the fact that the geometrical information in the reference polytope is globally oriented.
     For instance, the normals ``n`` and ``\hat{n}`` both have the same orientation, i.e ``n = (||\hat{e}||/||e||) (det J) J^{-T} \hat{n}``. Therefore ``\hat{n}`` is not fully local. See [2, Equation 2.11].
     In our case, we will be including the sign change in the transformation matrices, which will include all cell-and-dof-dependent information.
 

src/Arrays/Arrays.jl

@@ -47,9 +47,8 @@ export testargs
 export inverse_map
 
 export Broadcasting
-
 export Operation
-
+export InverseMap
 
 # LazyArray
 

src/FESpaces/Pullbacks.jl

@@ -1,51 +1,56 @@
 
 # TODO: We probably want to export these from Gridap.ReferenceFEs
-using Gridap.ReferenceFEs: PushforwardRefFE, Pushforward, InversePullback
+using Gridap.ReferenceFEs: PushforwardRefFE, Pushforward, Pullback
 using Gridap.ReferenceFEs: ContraVariantPiolaMap, CoVariantPiolaMap
 
 function get_cell_dof_basis(
   model::DiscreteModel, cell_reffe::AbstractArray{<:GenericRefFE{T}}, conformity::Conformity
 ) where T <: PushforwardRefFE
-  pushforward = Pushforward(T)
-  cell_args = get_cell_pushforward_arguments(pushforward, model, cell_reffe, conformity)
-  cell_dofs = lazy_map(get_dof_basis, cell_reffe)
-  return lazy_map(InversePullback(pushforward), cell_dofs, cell_args...)
+  pushforward, cell_change, cell_args = get_cell_pushforward(
+    Pushforward(T), model, cell_reffe, conformity
+  )
+  cell_ref_dofs = lazy_map(get_dof_basis, cell_reffe)
+  cell_phy_dofs = lazy_map(inverse_map(Pullback(pushforward)), cell_ref_dofs, cell_args...)
+  return lazy_map(linear_combination, cell_change, cell_phy_dofs) # TODO: Inverse
 end
 
 function get_cell_shapefuns(
   model::DiscreteModel, cell_reffe::AbstractArray{<:GenericRefFE{T}}, conformity::Conformity
 ) where T <: PushforwardRefFE
-  pushforward = Pushforward(T)
-  cell_args = get_cell_pushforward_arguments(pushforward, model, cell_reffe, conformity)
-  cell_dofs = lazy_map(get_shapefuns, cell_reffe)
-  return lazy_map(pushforward, cell_dofs, cell_args...)
+  pushforward, cell_change, cell_args = get_cell_pushforward(
+    Pushforward(T), model, cell_reffe, conformity
+  )
+  cell_ref_fields = lazy_map(get_shapefuns, cell_reffe)
+  cell_phy_fields = lazy_map(pushforward, cell_ref_fields, cell_args...)
+  return lazy_map(linear_combination, cell_change, cell_phy_fields)
 end
 
-function get_cell_pushforward_arguments(
+function get_cell_pushforward(
   ::Pushforward, model::DiscreteModel, cell_reffe, conformity
 )
   @abstractmethod
 end
 
 # ContraVariantPiolaMap
 
-function get_cell_pushforward_arguments(
+function get_cell_pushforward(
   ::ContraVariantPiolaMap, model::DiscreteModel, cell_reffe, conformity
 )
   cell_map  = get_cell_map(get_grid(model))
   Jt = lazy_map(Broadcasting(∇),cell_map)
-  sign_flip = lazy_map(Broadcasting(constant_field),get_sign_flip(model, cell_reffe))
-  return Jt, sign_flip
+  change = get_sign_flip(model, cell_reffe)
+  return ContraVariantPiolaMap(), change, (Jt,)
 end
 
 # CoVariantPiolaMap
 
-function get_cell_pushforward_arguments(
+function get_cell_pushforward(
   ::CoVariantPiolaMap, model::DiscreteModel, cell_reffe, conformity
 )
   cell_map = get_cell_map(get_grid(model))
   Jt = lazy_map(Broadcasting(∇),cell_map)
-  return Jt
+  change = lazy_map(r -> Diagonal(ones(num_dofs(r))), cell_reffe) # TODO: Replace by edge-signs
+  return CoVariantPiolaMap(), change, (Jt,)
 end
 
 # SignFlipMap
@@ -72,7 +77,7 @@ function return_cache(k::SignFlipMap,reffe,facet_own_dofs,cell)
   cell_facets = get_faces(topo, Dc, Dc-1)
   cell_facets_cache = array_cache(cell_facets)
 
-  return cell_facets, cell_facets_cache, CachedVector(Bool)
+  return cell_facets, cell_facets_cache, CachedVector(Float64)
 end
 
 function evaluate!(cache,k::SignFlipMap,reffe,facet_own_dofs,cell)
@@ -81,19 +86,19 @@ function evaluate!(cache,k::SignFlipMap,reffe,facet_own_dofs,cell)
 
   setsize!(sign_flip_cache, (num_dofs(reffe),))
   sign_flip = sign_flip_cache.array
-  sign_flip .= false
+  fill!(sign_flip, one(eltype(sign_flip)))
 
   facets = getindex!(cell_facets_cache,cell_facets,cell)
   for (lfacet,facet) in enumerate(facets)
     owner = facet_owners[facet]
     if owner != cell
       for dof in facet_own_dofs[lfacet]
-        sign_flip[dof] = true
+        sign_flip[dof] = -1.0
       end
     end
   end
 
-  return sign_flip
+  return Diagonal(sign_flip)
 end
 
 function get_sign_flip(model::DiscreteModel{Dc}, cell_reffe) where Dc

src/ReferenceFEs/Dofs.jl

@@ -1,5 +1,58 @@
 abstract type Dof <: Map end
 
+"""
+    struct LinearCombinationDofVector{T} <: AbstractVector{Dof}
+      values::Matrix{T}
+      dofs::AbstractVector{<:Dof}
+    end
+
+Type that implements a dof basis (a) as the linear combination of a dof basis
+(b). The dofs are first evaluated at dof basis (b) (field `dofs`) and the
+dof values are next mapped to dof basis (a) applying a change of basis (field
+`values`).
+
+Fields:
+
+- `values::Matrix{T}` the matrix of the change from dof basis (b) to (a)
+- `dofs::AbstractVector{<:Dof}` A type representing dof basis (b)
+"""
+struct LinearCombinationDofVector{T,V,F} <: AbstractVector{T}
+  values::V
+  dofs::F
+  function LinearCombinationDofVector(
+    values::AbstractMatrix{<:Number},
+    dofs::AbstractVector{<:Dof}
+  )
+    T = eltype(dofs)
+    V = typeof(values)
+    F = typeof(dofs)
+    new{T,V,F}(values,dofs)
+  end
+end
+
+@inline Base.size(a::LinearCombinationDofVector) = size(a.values,2)
+@inline Base.IndexStyle(::LinearCombinationDofVector) = IndexLinear()
+@inline Base.getindex(::LinearCombinationDofVector{T},::Integer) where T = T()
+
+function linear_combination(a::AbstractMatrix{<:Number}, b::AbstractVector{<:Dof})
+  LinearCombinationDofVector(a,b)
+end
+
+function return_cache(b::LinearCombinationDofVector,field)
+  k = LinearCombinationMap(:)
+  cf = return_cache(b.dofs,field)
+  fx = evaluate!(cf,b.dofs,field)
+  ck = return_cache(k,b.values,fx)
+  return cf, ck
+end
+
+function evaluate!(cache,b::LinearCombinationDofVector,field)
+  cf, ck = cache
+  k = LinearCombinationMap(:)
+  fx = evaluate!(cf,b.dofs,field)
+  return evaluate!(ck,k,b.values,fx)
+end
+
 """
     test_dof(dof,field,v;cmp::Function=(==))
 

src/ReferenceFEs/Pullbacks.jl

@@ -46,10 +46,9 @@ function Arrays.evaluate!(
   ::Broadcasting{typeof(∇)},
   a::Fields.BroadcastOpFieldArray{<:Pushforward}
 )
-  v, Jt, sign_flip = a.args
-  ∇v = Broadcasting(∇)(v)
-  k = ContraVariantPiolaMap()
-  Broadcasting(Operation(k))(∇v,Jt,sign_flip)
+  v, pf_args... = a.args
+  grad_v = Broadcasting(∇)(v)
+  Broadcasting(Operation(a.op))(grad_v,pf_args...)
 end
 
 # InversePushforward
@@ -106,6 +105,12 @@ function evaluate!(
   Broadcasting(Operation(k))(f_phys,args...)
 end
 
+function evaluate!(
+  cache, k::InversePushforward, f_phys::Field, args...
+)
+  Operation(k)(f_phys,args...)
+end
+
 # MappedDofBasis
 
 """
@@ -208,7 +213,6 @@ const InversePullback{PB} = InverseMap{PB} where PB <: Pullback
 function Arrays.lazy_map(
   ::typeof(evaluate), k::LazyArray{<:Fill{<:InversePullback}}, phys_cell_fields::AbstractArray
 )
-  println("InversePullback dispatch")
   pb = inverse_map(k.maps.value)
   ref_cell_dofs, pf_args... = k.args
   ref_cell_fields = lazy_map(inverse_map(pb.pushforward), phys_cell_fields, pf_args...)
@@ -227,25 +231,23 @@ end
 struct ContraVariantPiolaMap <: Pushforward end
 
 function evaluate!(
-  cache, ::ContraVariantPiolaMap, v_ref::Number, Jt::Number, sign_flip::Bool
+  cache, ::ContraVariantPiolaMap, v_ref::Number, Jt::Number
 )
-  sign  = (-1)^sign_flip
   idetJ = 1. / meas(Jt)
-  return (sign * v_ref) ⋅(idetJ * Jt)
+  return v_ref ⋅ (idetJ * Jt)
 end
 
 function return_cache(
-  ::InversePushforward{ContraVariantPiolaMap}, v_phys::Number, Jt::Number, sign_flip::Bool
+  ::InversePushforward{ContraVariantPiolaMap}, v_phys::Number, Jt::Number
 )
   nothing
 end
 
 function evaluate!(
-  cache, ::InversePushforward{ContraVariantPiolaMap}, v_phys::Number, Jt::Number, sign_flip::Bool
+  cache, ::InversePushforward{ContraVariantPiolaMap}, v_phys::Number, Jt::Number
 )
-  sign = (-1)^sign_flip
   detJ = meas(Jt)
-  return (sign * v_phys) ⋅ (detJ * pinvJt(Jt))
+  return v_phys ⋅ (detJ * pinvJt(Jt))
 end
 
 # TODO: Should this be here? Probably not...
@@ -263,11 +265,13 @@ function Fields.DIV(a::LazyArray{<:Fill{typeof(linear_combination)}})
 end
 
 function Fields.DIV(f::LazyArray{<:Fill{Broadcasting{Operation{ContraVariantPiolaMap}}}})
-  ϕrgₖ       = f.args[1]
-  fsign_flip = f.args[3]
-  div_ϕrgₖ   = lazy_map(Broadcasting(divergence),ϕrgₖ)
-  fsign_flip = lazy_map(Broadcasting(Operation(x->(-1)^x)), fsign_flip)
-  lazy_map(Broadcasting(Operation(*)),fsign_flip,div_ϕrgₖ)
+  ϕrgₖ = f.args[1]
+  return lazy_map(Broadcasting(divergence),ϕrgₖ)
+end
+
+function Fields.DIV(f::Fill{<:Fields.BroadcastOpFieldArray{ContraVariantPiolaMap}})
+  ϕrgₖ = f.value.args[1]
+  return Fill(Broadcasting(divergence)(ϕrgₖ),length(f))
 end
 
 # CoVariantPiolaMap

src/ReferenceFEs/ReferenceFEs.jl

@@ -265,6 +265,4 @@ include("BezierRefFEs.jl")
 
 include("ModalC0RefFEs.jl")
 
-include("LinearCombinationDofVectors.jl")
-
 end # module