Merge pull request gridap#1067 from gridap/pullbacks

Pullbacks

NEWS.md

@@ -11,6 +11,11 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 - Added AMR-related methods `mark` and `estimate` to `Adaptivity` module. Implemented Dorfler marking strategy. Since PR[#1063](https://github.com/gridap/Gridap.jl/pull/1063).
 
+### Changed
+
+- Low level optimisations to reduce allocations. `AffineMap` renamed to `AffineField`. New `AffineMap <: Map`, doing the same as `AffineField` without struct allocation. New `ConstantMap <: Map`, doing the same as `ConstantField` without struct allocation. Since PR[#1043](https://github.com/gridap/Gridap.jl/pull/1043).
+- `ConstantFESpaces` can now be built on triangulations. Since PR[#1069](https://github.com/gridap/Gridap.jl/pull/1069)
+
 ## [0.18.8] - 2024-12-2
 
 ### Added
@@ -25,6 +30,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 - Fixed issue where barycentric refinement rule in 3D would not produce oriented meshes. Since PR[#1055](https://github.com/gridap/Gridap.jl/pull/1055).
 
 ### Changed
+
 - Optimized MonomialBasis low-level functions. Since PR[#1059](https://github.com/gridap/Gridap.jl/pull/1059).
 
 ## [0.18.7] - 2024-10-8

docs/make.jl

@@ -23,6 +23,7 @@ pages = [
   "Gridap.Adaptivity" => "Adaptivity.md",
   "Developper notes" => Any[
     "dev-notes/block-assemblers.md",
+    "dev-notes/pullbacks.md",
   ],
 ]
 

docs/src/Adaptivity.md

@@ -156,6 +156,4 @@ When a cell is refined, we need to be able to evaluate the fields defined on the
 
 ```@docs
 FineToCoarseField
-FineToCoarseDofBasis
-FineToCoarseRefFE
 ```

docs/src/ReferenceFEs.md

@@ -73,8 +73,18 @@ Pages   = ["LagrangianRefFEs.jl","LagrangianDofBases.jl","SerendipityRefFEs.jl",
 
 ### Moment-Based ReferenceFEs
 
+#### Framework
+
+```@autodocs
+Modules = [ReferenceFEs,]
+Order   = [:type, :constant, :macro, :function]
+Pages   = ["MomentBasedReferenceFEs.jl","Pullbacks.jl"]
+```
+
+#### Available Moment-Based ReferenceFEs
+
 ```@autodocs
 Modules = [ReferenceFEs,]
 Order   = [:type, :constant, :macro, :function]
-Pages   = ["MomentBasedReferenceFEs.jl","RaviartThomasRefFEs.jl","NedelecRefFEs.jl","BDMRefFEs.jl","CRRefFEs.jl"]
+Pages   = ["RaviartThomasRefFEs.jl","NedelecRefFEs.jl","BDMRefFEs.jl","CRRefFEs.jl"]
 ```

docs/src/dev-notes/pullbacks.md

@@ -0,0 +1,82 @@
+
+# FE basis transformations
+
+## Notation
+
+Consider a reference polytope ``\hat{K}``, mapped to the physical space by a **geometrical map** ``F``, i.e. ``K = F(\hat{K})``. Consider also a linear function space on the reference polytope ``\hat{V}``, and a set of unisolvent degrees of freedom represented by moments in the dual space ``\hat{V}^*``.
+
+Throughout this document, we will use the following notation:
+
+- ``\varphi \in V`` is a **physical field** ``\varphi : K \rightarrow \mathbb{R}^k``. A basis of ``V`` is denoted by ``\Phi = \{\varphi\}``.
+- ``\hat{\varphi} \in \hat{V}`` is a **reference field** ``\hat{\varphi} : \hat{K} \rightarrow \mathbb{R}^k``. A basis of ``\hat{V}`` is denoted by ``\hat{\Phi} = \{\hat{\varphi}\}``.
+- ``\sigma \in V^*`` is a **physical moment** ``\sigma : V \rightarrow \mathbb{R}``. A basis of ``V^*`` is denoted by ``\Sigma = \{\sigma\}``.
+- ``\hat{\sigma} \in \hat{V}^*`` is a **reference moment** ``\hat{\sigma} : \hat{V} \rightarrow \mathbb{R}``. A basis of ``\hat{V}^*`` is denoted by ``\hat{\Sigma} = \{\hat{\sigma}\}``.
+
+## Pullbacks and Pushforwards
+
+We define a **pushforward** map as ``F^* : \hat{V} \rightarrow V``, mapping reference fields to physical fields. Given a pushforward ``F^*``, we define:
+
+- The **pullback** ``F_* : V^* \rightarrow \hat{V}^*``, mapping physical moments to reference moments. Its action on physical dofs is defined in terms of the pushforward map ``F^*`` as ``\hat{\sigma} = F_*(\sigma) := \sigma \circ F^*``.
+- The **inverse pushforward** ``(F^*)^{-1} : V \rightarrow \hat{V}``, mapping physical fields to reference fields.
+- The **inverse pullback** ``(F_*)^{-1} : \hat{V}^* \rightarrow V^*``, mapping reference moments to physical moments. Its action on reference dofs is defined in terms of the inverse pushforward map ``(F^*)^{-1}`` as ``\sigma = (F_*)^{-1}(\hat{\sigma}) := \hat{\sigma} \circ (F^*)^{-1}``.
+
+## Change of basis
+
+In many occasions, we will have that (as a basis)
+
+```math
+\hat{\Sigma} \neq F_*(\Sigma), \quad \text{and} \quad \Phi \neq F^*(\hat{\Phi})
+```
+
+To maintain conformity and proper scaling in these cases, we define cell-dependent invertible changes of basis ``P`` and ``M``, such that
+
+```math
+\hat{\Sigma} = P F_*(\Sigma), \quad \text{and} \quad \Phi = M F^*(\hat{\Phi})
+```
+
+An important result from [1, Theorem 3.1] is that ``P = M^T``.
+
+!!! details
+    [1, Lemma 2.6]: A key ingredient is that given ``M`` a matrix we have ``\Sigma (M \Phi) = \Sigma (\Phi) M^T`` since
+    ```math
+      [\Sigma (M \Phi)]_{ij} = \sigma_i (M_{jk} \varphi_k) = M_{jk} \sigma_i (\varphi_k) = [\Sigma (\Phi) M^T]_{ij}
+    ```
+    where we have used that moments are linear.
+
+We then have the following diagram:
+
+```math
+\hat{V}^* \xleftarrow{P} \hat{V}^* \xleftarrow{F_*} V^* \\
+\hat{V} \xrightarrow{F^*} V \xrightarrow{P^T} V
+```
+
+!!! details
+    The above diagram is well defined, since we have
+    ```math
+    \hat{\Sigma}(\hat{\Phi}) = P F_* (\Sigma)(F^{-*} (P^{-T} \Phi)) = P \Sigma (F^* (F^{-*} P^{-T} \Phi)) = P \Sigma (P^{-T} \Phi) = P \Sigma (\Phi) P^{-1} = Id \\
+    \Sigma(\Phi) = F_*^{-1}(P^{-1}\hat{\Sigma})(P^T F^*(\hat{\Phi})) = P^{-1} \hat{\Sigma} (F^{-*}(P^T F^*(\hat{\Phi}))) = P^{-1} \hat{\Sigma} (P^T \hat{\Phi}) = P^{-1} \hat{\Sigma}(\hat{\Phi}) P = Id
+    ```
+
+From an implementation point of view, it is more natural to build ``P^{-1}``  and then retrieve all other matrices by transposition/inversion.
+
+## Interpolation
+
+In each cell ``K`` and for ``C_b^k(K)`` the space of functions defined on ``K`` with at least ``k`` bounded derivatives, we define the interpolation operator ``I_K : C_b^k(K) \rightarrow V`` as
+
+```math
+I_K(g) = \Sigma(g) \Phi \quad, \quad \Sigma(g) = P^{-1} \hat{\Sigma}(F^{-*}(g))
+```
+
+## Implementation notes
+
+!!! 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 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.
+
+## References
+
+[1] [Kirby 2017, A general approach to transforming finite elements.](https://arxiv.org/abs/1706.09017)
+
+[2] [Aznaran et al. 2021, Transformations for Piola-mapped elements.](https://arxiv.org/abs/2110.13224)

src/Adaptivity/RefinementRules.jl

@@ -25,7 +25,7 @@ end
 # - The reason why we are saving both the cell maps and the inverse cell maps is to avoid recomputing
 #   them when needed. This is needed for performance when the RefinementRule is used for MacroFEs.
 #   Also, in the case the ref_grid comes from a CartesianGrid, we save the cell maps as 
-#   AffineMaps, which are more efficient than the default linear_combinations.
+#   AffineFields, which are more efficient than the default linear_combinations.
 # - We cannot parametrise the RefinementRule by all it's fields, because we will have different types of
 #   RefinementRules in a single mesh. It's the same reason why we don't parametrise the ReferenceFE type.
 

src/Adaptivity/deprecated/FineToCoarseReferenceFEs.jl

@@ -0,0 +1,132 @@
+"""
+"""
+struct FineToCoarseDofBasis{T,A,B,C} <: AbstractVector{T}
+  dof_basis :: A
+  rrule     :: B
+  child_ids :: C
+
+  function FineToCoarseDofBasis(dof_basis::AbstractVector{T},rrule::RefinementRule) where {T<:Dof}
+    nodes = get_nodes(dof_basis)
+    child_ids = map(x -> x_to_cell(rrule,x),nodes)
+
+    A = typeof(dof_basis)
+    B = typeof(rrule)
+    C = typeof(child_ids)
+    new{T,A,B,C}(dof_basis,rrule,child_ids)
+  end
+end
+
+Base.size(a::FineToCoarseDofBasis) = size(a.dof_basis)
+Base.axes(a::FineToCoarseDofBasis) = axes(a.dof_basis)
+Base.getindex(a::FineToCoarseDofBasis,i::Integer) = getindex(a.dof_basis,i)
+Base.IndexStyle(a::FineToCoarseDofBasis) = IndexStyle(a.dof_basis)
+
+ReferenceFEs.get_nodes(a::FineToCoarseDofBasis) = get_nodes(a.dof_basis)
+
+# Default behaviour
+Arrays.return_cache(b::FineToCoarseDofBasis,field) = return_cache(b.dof_basis,field)
+Arrays.evaluate!(cache,b::FineToCoarseDofBasis,field) = evaluate!(cache,b.dof_basis,field)
+
+# Spetialized behaviour
+function Arrays.return_cache(s::FineToCoarseDofBasis{T,<:LagrangianDofBasis},field::FineToCoarseField) where T
+  b     = s.dof_basis
+  cf    = return_cache(field,b.nodes,s.child_ids)
+  vals  = evaluate!(cf,field,b.nodes,s.child_ids)
+  ndofs = length(b.dof_to_node)
+  r = ReferenceFEs._lagr_dof_cache(vals,ndofs)
+  c = CachedArray(r)
+  return (c, cf)
+end
+
+function Arrays.evaluate!(cache,s::FineToCoarseDofBasis{T,<:LagrangianDofBasis},field::FineToCoarseField) where T
+  c, cf = cache
+  b     = s.dof_basis
+  vals  = evaluate!(cf,field,b.nodes,s.child_ids)
+  ndofs = length(b.dof_to_node)
+  T2    = eltype(vals)
+  ncomps = num_indep_components(T2)
+  @check ncomps == num_indep_components(eltype(b.node_and_comp_to_dof)) """\n
+  Unable to evaluate LagrangianDofBasis. The number of components of the
+  given Field does not match with the LagrangianDofBasis.
+  If you are trying to interpolate a function on a FESpace make sure that
+  both objects have the same value type.
+  For instance, trying to interpolate a vector-valued function on a scalar-valued FE space
+  would raise this error.
+  """
+  ReferenceFEs._evaluate_lagr_dof!(c,vals,b.node_and_comp_to_dof,ndofs,ncomps)
+end
+
+function Arrays.return_cache(s::FineToCoarseDofBasis{T,<:MomentBasedDofBasis},field::FineToCoarseField) where T
+  b     = s.dof_basis
+  cf    = return_cache(field,b.nodes,s.child_ids)
+  vals  = evaluate!(cf,field,b.nodes,s.child_ids)
+  ndofs = num_dofs(b)
+  r = ReferenceFEs._moment_dof_basis_cache(vals,ndofs)
+  c = CachedArray(r)
+  return (c, cf)
+end
+
+function Arrays.evaluate!(cache,s::FineToCoarseDofBasis{T,<:MomentBasedDofBasis},field::FineToCoarseField) where T
+  c, cf = cache
+  b     = s.dof_basis
+  vals  = evaluate!(cf,field,b.nodes,s.child_ids)
+  dofs  = c.array
+  ReferenceFEs._eval_moment_dof_basis!(dofs,vals,b)
+  dofs
+end
+
+
+"""
+  Wrapper for a ReferenceFE which is specialised for
+  efficiently evaluating FineToCoarseFields.
+"""
+struct FineToCoarseRefFE{T,D,A} <: ReferenceFE{D}
+  reffe     :: T
+  dof_basis :: A
+
+  function FineToCoarseRefFE(reffe::ReferenceFE{D},dof_basis::FineToCoarseDofBasis) where D
+    T = typeof(reffe)
+    A = typeof(dof_basis)
+    new{T,D,A}(reffe,dof_basis)
+  end
+end
+
+ReferenceFEs.num_dofs(reffe::FineToCoarseRefFE)      = num_dofs(reffe.reffe)
+ReferenceFEs.get_polytope(reffe::FineToCoarseRefFE)  = get_polytope(reffe.reffe)
+ReferenceFEs.get_prebasis(reffe::FineToCoarseRefFE)  = get_prebasis(reffe.reffe)
+ReferenceFEs.get_dof_basis(reffe::FineToCoarseRefFE) = reffe.dof_basis
+ReferenceFEs.Conformity(reffe::FineToCoarseRefFE)    = Conformity(reffe.reffe)
+ReferenceFEs.get_face_dofs(reffe::FineToCoarseRefFE) = get_face_dofs(reffe.reffe)
+ReferenceFEs.get_shapefuns(reffe::FineToCoarseRefFE) = get_shapefuns(reffe.reffe)
+ReferenceFEs.get_metadata(reffe::FineToCoarseRefFE)  = get_metadata(reffe.reffe)
+ReferenceFEs.get_orders(reffe::FineToCoarseRefFE)    = get_orders(reffe.reffe)
+ReferenceFEs.get_order(reffe::FineToCoarseRefFE)     = get_order(reffe.reffe)
+
+ReferenceFEs.Conformity(reffe::FineToCoarseRefFE,sym::Symbol) = Conformity(reffe.reffe,sym)
+ReferenceFEs.get_face_own_dofs(reffe::FineToCoarseRefFE,conf::Conformity) = get_face_own_dofs(reffe.reffe,conf)
+
+
+function ReferenceFEs.ReferenceFE(p::Polytope,rrule::RefinementRule,name::ReferenceFEName,order)
+  FineToCoarseRefFE(p,rrule,name,Float64,order)
+end
+
+function ReferenceFEs.ReferenceFE(p::Polytope,rrule::RefinementRule,name::ReferenceFEName,::Type{T},order) where T
+  FineToCoarseRefFE(p,rrule,name,T,order)
+end
+
+function FineToCoarseRefFE(p::Polytope,rrule::RefinementRule,name::ReferenceFEName,::Type{T},order) where T
+  @check p == get_polytope(rrule)
+  reffe = ReferenceFE(p,name,T,order)
+  dof_basis = FineToCoarseDofBasis(get_dof_basis(reffe),rrule)
+  return FineToCoarseRefFE(reffe,dof_basis)
+end
+
+# FESpaces constructors
+
+function FESpaces.TestFESpace(model::DiscreteModel,rrules::AbstractVector{<:RefinementRule},reffe::Tuple{<:ReferenceFEName,Any,Any};kwargs...)
+  @check num_cells(model) == length(rrules)
+  @check all(CompressedArray(get_polytopes(model),get_cell_type(model)) .== lazy_map(get_polytope,rrules))
+  basis, reffe_args, reffe_kwargs = reffe
+  reffes = lazy_map(rr -> ReferenceFE(get_polytope(rr),rr,basis,reffe_args...;reffe_kwargs...),rrules)
+  return TestFESpace(model,reffes;kwargs...)
+end

src/Arrays/Arrays.jl

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

src/Arrays/Maps.jl

@@ -296,3 +296,16 @@ function inverse_map(f)
   Function inverse_map is not implemented yet for objects of type $(typeof(f))
   """
 end
+
+struct InverseMap{F} <: Map
+  original::F
+end
+
+function evaluate!(cache,k::InverseMap,args...)
+  @notimplemented """\n
+  The inverse evaluation is not implemented yet for maps of type $(typeof(k.original))
+  """
+end
+
+inverse_map(k::Map) = InverseMap(k)
+inverse_map(k::InverseMap) = k.original

src/FESpaces/ConstantFESpaces.jl

@@ -1,49 +1,63 @@
 
 """
-struct ConstantFESpace <: SingleFieldFESpace
-  # private fields
-end
+    struct ConstantFESpace <: SingleFieldFESpace
+      # private fields
+    end
+
+    ConstantFESpace(model::DiscreteModel; vector_type=Vector{Float64}, field_type=Float64)
+    ConstantFESpace(trian::Triangulation; vector_type=Vector{Float64}, field_type=Float64)
+
+FESpace that is constant over the provided model/triangulation. Typically used as  
+lagrange multipliers. The kwargs `vector_type` and `field_type` are used to specify the
+types of the dof-vector and dof-value respectively.
 """
 struct ConstantFESpace{V,T,A,B,C} <: SingleFieldFESpace
-  model::DiscreteModel
+  trian::Triangulation
   cell_basis::A
   cell_dof_basis::B
   cell_dof_ids::C
 
-  function ConstantFESpace(model;
-                         vector_type::Type{V}=Vector{Float64},
-                         field_type::Type{T}=Float64) where {V,T}
-    function setup_cell_reffe(model::DiscreteModel,
-        reffe::Tuple{<:ReferenceFEName,Any,Any}; kwargs...)
-      basis, reffe_args,reffe_kwargs = reffe
-      cell_reffe = ReferenceFE(model,basis,reffe_args...;reffe_kwargs...)
-    end
+  function ConstantFESpace(
+    model::DiscreteModel{Dc},
+    trian::Triangulation{Dc};
+    vector_type::Type{V}=Vector{Float64},
+    field_type::Type{T}=Float64
+  ) where {Dc,V,T}
+    @assert num_cells(model) == num_cells(trian)
 
-    reffe = ReferenceFE(lagrangian,T,0)
-    cell_reffe = setup_cell_reffe(model,reffe)
+    basis, reffe_args, reffe_kwargs = ReferenceFE(lagrangian,T,0)
+    cell_reffe = ReferenceFE(model,basis,reffe_args...;reffe_kwargs...)
     cell_basis_array = lazy_map(get_shapefuns,cell_reffe)
 
     cell_basis = SingleFieldFEBasis(
-      cell_basis_array,
-      Triangulation(model),
-      TestBasis(),
-      ReferenceDomain())
+      cell_basis_array, trian, TestBasis(), ReferenceDomain()
+    )
+    cell_dof_basis = CellDof(
+      lazy_map(get_dof_basis,cell_reffe),trian,ReferenceDomain()
+    )
+    cell_dof_ids = Fill(Int32(1):Int32(num_indep_components(field_type)),num_cells(trian))
 
-    cell_dof_basis_array = lazy_map(get_dof_basis,cell_reffe)
-    cell_dof_basis = CellDof(cell_dof_basis_array,Triangulation(model),ReferenceDomain())
-
-    cell_dof_ids = Fill(Int32(1):Int32(num_indep_components(field_type)),num_cells(model))
     A = typeof(cell_basis)
     B = typeof(cell_dof_basis)
     C = typeof(cell_dof_ids)
-    new{V,T,A,B,C}(model, cell_basis, cell_dof_basis, cell_dof_ids)
+    new{V,T,A,B,C}(trian, cell_basis, cell_dof_basis, cell_dof_ids)
   end
 end
 
+function ConstantFESpace(model::DiscreteModel; kwargs...)
+  trian = Triangulation(model)
+  ConstantFESpace(model,trian; kwargs...)
+end
+
+function ConstantFESpace(trian::Triangulation; kwargs...)
+  model = get_active_model(trian)
+  ConstantFESpace(model,trian; kwargs...)
+end
+
 TrialFESpace(f::ConstantFESpace) = f
 
 # Delegated functions
-get_triangulation(f::ConstantFESpace) = Triangulation(f.model)
+get_triangulation(f::ConstantFESpace) = f.trian
 
 ConstraintStyle(::Type{<:ConstantFESpace}) = UnConstrained()