Merge branch 'moment-based-reffes' of github.com:gridap/Gridap.jl int…

…o moment-based-reffes

NEWS.md

@@ -24,6 +24,9 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ### 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)
+
 - Existing Jacobi polynomial bases/spaces were renamed to Legendre (which they were).
 - `Monomial` is now subtype of the new abstract type`Polynomial <: Field`
 - `MonomialBasis` is now an alias for `UniformPolyBasis{...,Monomial}`
@@ -56,6 +59,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

@@ -4,7 +4,87 @@ CurrentModule = Gridap.ReferenceFEs
 
 # Gridap.ReferenceFEs
 
+```@contents
+Pages = ["ReferenceFEs.md"]
+Depth = 2:3
+```
+
+## Polytopes
+
+### Abstract API
+
+```@autodocs
+Modules = [ReferenceFEs,]
+Order   = [:type, :constant, :macro, :function]
+Pages   = ["/Polytopes.jl"]
+```
+
+### Extrusion Polytopes
+
+```@autodocs
+Modules = [ReferenceFEs,]
+Order   = [:type, :constant, :macro, :function]
+Pages   = ["ExtrusionPolytopes.jl"]
+```
+
+### General Polytopes
+
+```@autodocs
+Modules = [ReferenceFEs,]
+Order   = [:type, :constant, :macro, :function]
+Pages   = ["GeneralPolytopes.jl"]
+```
+
+## Quadratures
+
+### Abstract API
+
+```@autodocs
+Modules = [ReferenceFEs,]
+Order   = [:type, :constant, :macro, :function]
+Pages   = ["/Quadratures.jl"]
+```
+
+### Available Quadratures
+
+```@autodocs
+Modules = [ReferenceFEs,]
+Order   = [:type, :constant, :macro, :function]
+Pages   = ["TensorProductQuadratures.jl","DuffyQuadratures.jl","StrangeQuadratures.jl","XiaoGimbutasQuadratures.jl"]
+```
+
+## ReferenceFEs
+
+### Abstract API
+
 ```@autodocs
 Modules = [ReferenceFEs,]
+Order   = [:type, :constant, :macro, :function]
+Pages   = ["ReferenceFEInterfaces.jl","Dofs.jl","LinearCombinationDofVectors.jl"]
 ```
 
+### Nodal ReferenceFEs
+
+```@autodocs
+Modules = [ReferenceFEs,]
+Order   = [:type, :constant, :macro, :function]
+Pages   = ["LagrangianRefFEs.jl","LagrangianDofBases.jl","SerendipityRefFEs.jl","BezierRefFEs.jl","ModalC0RefFEs.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   = ["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/Adaptivity.jl

@@ -45,7 +45,6 @@ export DorflerMarking, mark, estimate
 include("RefinementRules.jl")
 include("FineToCoarseFields.jl")
 include("OldToNewFields.jl")
-include("FineToCoarseReferenceFEs.jl")
 include("AdaptivityGlues.jl")
 include("AdaptedDiscreteModels.jl")
 include("AdaptedTriangulations.jl")

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/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()
 

src/FESpaces/FESpaces.jl

@@ -219,9 +219,7 @@ include("UnconstrainedFESpaces.jl")
 
 include("ConformingFESpaces.jl")
 
-include("DivConformingFESpaces.jl")
-
-include("CurlConformingFESpaces.jl")
+include("Pullbacks.jl")
 
 include("FESpaceFactories.jl")
 
@@ -253,8 +251,6 @@ include("CLagrangianFESpaces.jl")
 
 include("DirichletFESpaces.jl")
 
-#include("ExtendedFESpaces.jl")
-
 include("FESpacesWithLinearConstraints.jl")
 
 include("DiscreteModelWithFEMaps.jl")