Merge branch 'adaptivity' into adaptivity

NEWS.md

@@ -5,6 +5,30 @@ All notable changes to this project will be documented in this file.
 The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/),
 and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 
+## [0.18.6] - 2024-08-29
+
+### Fixed
+
+- Improved performance of PR[#967](https://github.com/gridap/Gridap.jl/pull/967). Along the way, opened the door to Triangulations of different type in SkeletonTriangulation. Since PR[#1026](https://github.com/gridap/Gridap.jl/pull/1026).
+
+## [0.18.5] - 2024-08-28 
+
+### Changed
+
+- Misc changes required to support facet integration on non-conforming meshes. These changes do not involve methods of the public API. Since PR[#967](https://github.com/gridap/Gridap.jl/pull/967)
+
+## [0.18.4] - 2024-08-09
+
+### Changed
+
+- Added WriteVTK kwargs to control the output encoding for vtk files. Since PR[#1016](https://github.com/gridap/Gridap.jl/pull/1016).
+
+### Fixed
+
+- Passing `kwargs` from `refine` to `simplexify` functions in Adaptivity. Since PR[#1015](https://github.com/gridap/Gridap.jl/pull/1015).
+- Fixed `interpolate` for `ZeroMeanFESpace`. Since PR[#1020](https://github.com/gridap/Gridap.jl/pull/1020).
+- Fixed `gather_free_and_dirichlet_values!` for `FESpaceWithConstantFixed`. Since PR[#1020](https://github.com/gridap/Gridap.jl/pull/1020).
+
 ## [0.18.3] - 2024-07-11
 
 ### Added

Project.toml

@@ -1,7 +1,7 @@
 name = "Gridap"
 uuid = "56d4f2e9-7ea1-5844-9cf6-b9c51ca7ce8e"
 authors = ["Santiago Badia <[email protected]>", "Francesc Verdugo <[email protected]>", "Alberto F. Martin <[email protected]>"]
-version = "0.18.3"
+version = "0.18.6"
 
 [deps]
 AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c"

docs/src/Adaptivity.md

@@ -37,6 +37,8 @@ that determines the refinement strategy to be used. The following strategies are
 
 - `"red_green"` :: Red-Green refinement, default.
 - `"nvb"` :: Longest-edge bisection (only for meshes of TRIangles)
+- `"barycentric"` :: Barycentric refinement (only for meshes of TRIangles)
+- `"simplexify"` :: Simplexify refinement. Same resulting mesh as the `simplexify` method, but keeps track of the parent-child relationships.
 
 Additionally, the method takes a kwarg `cells_to_refine` that determines which cells will be refined.
 Possible input types are:
@@ -45,8 +47,7 @@ Possible input types are:
 - `AbstractArray{<:Bool}` of size `num_cells(model)` :: Only cells such that `cells_to_refine[iC] == true` get refined.
 - `AbstractArray{<:Integer}` :: Cells for which `gid ∈ cells_to_refine` get refined
 
-The algorithms try to respect the `cells_to_refine` input as much as possible, but some additional cells
-might get refined in order to guarantee that the mesh remains conforming.
+The algorithms try to respect the `cells_to_refine` input as much as possible, but some additional cells might get refined in order to guarantee that the mesh remains conforming.
 
 ```julia
   function refine(model::UnstructuredDiscreteModel;refinement_method="red_green",kwargs...)
@@ -56,17 +57,29 @@ might get refined in order to guarantee that the mesh remains conforming.
 
 ## CartesianDiscreteModel refining
 
-The module provides a `refine` method for `CartesianDiscreteModel`. The method takes a `Tuple` of size `Dc`
-(the dimension of the model cells) that will determine how many times cells will be refined in
-each direction. For example, for a 2D model, `refine(model,(2,3))` will refine each QUAD cell into
-a 2x3 grid of cells.
+The module provides a `refine` method for `CartesianDiscreteModel`. The method takes a `Tuple` of size `Dc` (the dimension of the model cells) that will determine how many times cells will be refined in each direction. For example, for a 2D model, `refine(model,(2,3))` will refine each QUAD cell into a 2x3 grid of cells.
 
 ```julia
   function refine(model::CartesianDiscreteModel{Dc}, cell_partition::Tuple) where Dc
     [...]
   end
 ```
 
+## Macro Finite-Elements
+
+The module also provides support for macro finite-elements. From an abstract point of view, a macro finite-element is a finite-element defined on a refined polytope, where polynomial basis are defined on each of the subcells (creating a broken piece-wise polynomial space on the original polytope). From Gridap's point of view, a macro finite-element is a `ReferenceFE` defined on a `RefinementRule` from an array of `ReferenceFE`s defined on the subcells.
+
+Although there are countless combinations, here are two possible applications:
+
+- Linearized High-Order Lagrangian FESpaces: These are spaces which have the same DoFs as a high-order Lagrangian space, but where the basis functions are linear on each subcell.
+- Barycentrically-refined elements for Stokes-like problems: These are spaces where the basis functions for velocity are defined on the barycentrically-refined mesh, whereas the basis functions for pressure are defined on the original cells. This allows for exact so-called Stokes sequences (see [here](https://arxiv.org/abs/2002.02051)).
+
+The API is given by the following methods:
+
+```@docs
+  MacroReferenceFE
+```
+
 ## Notes for users
 
 Most of the tools provided by this module are showcased in the tests of the module itself, as well as the following tutorial (coming soon).
@@ -82,11 +95,35 @@ However, we want to stress a couple of key performance-critical points:
 
 ### RefinementRule API
 
-Given a `RefinementRule`, the library provides a set of methods to compute the mappings between parent (coarse) face ids and child (fine) face ids (and vice-versa). The ids are local to the `RefinementRule`.
+Given a `RefinementRule`, the library provides a set of methods to compute the mappings between parent (coarse) face ids and child (fine) face ids (and vice-versa).
+
+The most basic information (that can directly be hardcoded in the RefinementRule for performance) are the mappings between parent face ids and child face ids. These are provided by:
 
 ```@docs
 get_d_to_face_to_child_faces
 get_d_to_face_to_parent_face
+```
+
+On top of these two basic mappings, a whole plethora of additional topological mappings can be computed.
+These first set of routines extend the ReferenceFEs API to provide information on the face-to-node mappings and permutations:
+
+```@docs
+ReferenceFEs.get_face_vertices
+ReferenceFEs.get_face_coordinates
+ReferenceFEs.get_vertex_permutations
+ReferenceFEs.get_face_vertex_permutations
+```
+
+We also provide face-to-face maps:
+
+```@docs
+get_cface_to_num_own_ffaces
+get_cface_to_own_ffaces
+get_cface_to_ffaces
+get_cface_to_own_ffaces_to_lnodes
+get_cface_to_ffaces_to_lnodes
+get_cface_to_fface_permutations
+aggregate_cface_to_own_fface_data
 get_face_subface_ldof_to_cell_ldof
 ```
 

src/Adaptivity/AdaptedDiscreteModels.jl

@@ -2,9 +2,9 @@
 """
 
   `DiscreteModel` created by refining/coarsening another `DiscreteModel`.
-  
-  The refinement/coarsening hierarchy can be traced backwards by following the 
-  `parent` pointer chain. This allows the transfer of dofs 
+
+  The refinement/coarsening hierarchy can be traced backwards by following the
+  `parent` pointer chain. This allows the transfer of dofs
   between `FESpaces` defined on this model and its ancestors.
 
 """
@@ -77,7 +77,7 @@ end
 """
   function adapt(model::DiscreteModel,args...;kwargs...) :: AdaptedDiscreteModel
 
-  Returns an `AdaptedDiscreteModel` that is the result of adapting (mixed coarsening and refining) 
+  Returns an `AdaptedDiscreteModel` that is the result of adapting (mixed coarsening and refining)
   the given `DiscreteModel`.
 """
 function adapt(model::DiscreteModel,args...;kwargs...) :: AdaptedDiscreteModel
@@ -153,22 +153,22 @@ function _get_cartesian_domain(desc::CartesianDescriptor{D}) where D
   return Tuple(domain)
 end
 
-@generated function _c2v(idx::Union{NTuple{N,T},CartesianIndex{N}},sizes::NTuple{N,T}) where {N,T}    
+@generated function _c2v(idx::Union{NTuple{N,T},CartesianIndex{N}},sizes::NTuple{N,T}) where {N,T}
   res = :(idx[1])
   for d in 1:N-1
     ik = :((idx[$(d+1)]-1))
     for k in 1:d
         ik = :($ik * sizes[$k])
     end
-    res = :($res + $ik) 
+    res = :($res + $ik)
   end
   return res
 end
 
 @generated function _create_cartesian_f2c_maps(nC::NTuple{N,T},ref::NTuple{N,T}) where {N,T}
   J_f2c   = Meta.parse(prod(["(",["1+(I[$k]-1)÷ref[$k]," for k in 1:N]...,")"]))
   J_child = Meta.parse(prod(["(",["1+(I[$k]-1)%ref[$k]," for k in 1:N]...,")"]))
-  
+
   return :(begin
     nF = nC .* ref
     f2c_map   = Vector{Int}(undef,prod(nF))
@@ -180,7 +180,7 @@ end
       f2c_map[i] = _c2v(J_f2c,nC)
       child_map[i] = _c2v(J_child,ref)
     end
-        
+
     return f2c_map, child_map
   end)
 end

src/Adaptivity/AdaptivityGlues.jl

@@ -124,7 +124,7 @@ function get_o2n_faces_map(ncell_to_ocell::Vector{T}) where {T<:Integer}
   end
   Arrays.length_to_ptrs!(ptrs)
 
-  data = fill(zero(T),ptrs[end])
+  data = fill(zero(T),ptrs[end]-1)
   for iF in 1:nF
     iC = ncell_to_ocell[iF]
     data[ptrs[iC]] = iF

src/Adaptivity/CompositeQuadratures.jl

@@ -21,6 +21,12 @@ end
 
 struct CompositeQuadrature <: QuadratureName end
 
+"""
+    Quadrature(rr::RefinementRule,degree::Integer;kwargs...)
+
+Creates a CompositeQuadrature on the RefinementRule `rr` by concatenating 
+quadratures of degree `degree` on each subcell of the RefinementRule.
+"""
 function ReferenceFEs.Quadrature(rr::RefinementRule,degree::Integer;kwargs...)
   return ReferenceFEs.Quadrature(ReferenceFEs.get_polytope(rr),CompositeQuadrature(),rr,degree;kwargs...)
 end
@@ -52,6 +58,29 @@ function ReferenceFEs.Quadrature(
   end
   fpoints = map(get_coordinates,quads)
   ids = FineToCoarseIndices(conn)
-  coordinates = FineToCoarseArray(rr,cpoints,fpoints,ids)
+  coordinates = FineToCoarseArray{eltype(cpoints)}(rr,cpoints,fpoints,ids)
+  return GenericQuadrature(coordinates,weights,"Composite quadrature")
+end
+
+"""
+    CompositeQuadrature(quad::Quadrature,rr::RefinementRule)
+
+Creates a CompositeQuadrature on the RefinementRule `rr` by splitting
+the quadrature `quad` into the subcells of the RefinementRule.
+"""
+function CompositeQuadrature(
+  quad::Quadrature,rr::RefinementRule
+)
+  @check ReferenceFEs.get_polytope(quad) === ReferenceFEs.get_polytope(rr)
+
+  weights = get_weights(quad)
+  cpoints = get_coordinates(quad)
+
+  fpoints, conn = evaluate(CoarseToFinePointMap(),rr,cpoints)
+  fpoints = collect(Vector{eltype(cpoints)},fpoints)
+  conn = collect(Vector{Int32},conn)
+
+  ids = FineToCoarseIndices(conn)
+  coordinates = FineToCoarseArray{eltype(cpoints)}(rr,cpoints,fpoints,ids)
   return GenericQuadrature(coordinates,weights,"Composite quadrature")
 end