Merge branch 'master' of github.com:gridap/Gridap.jl into benchmarks

NEWS.md

@@ -5,13 +5,20 @@ 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).
 
+## [Unreleased]
+
+### Added
+
+- Added MacroFElements. These are defined as having the basis/dof-basis of a FESpace created on top of a RefinementRule. Since PR[#1024](https://github.com/gridap/Gridap.jl/pull/1024).
+- Added Barycentric refinement rule in 2D and 3D. Added Simplexify refinement rule. Since PR[#1024](https://github.com/gridap/Gridap.jl/pull/1024).
+
 ## [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 
+## [0.18.5] - 2024-08-28
 
 ### Changed
 

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/Adaptivity.jl

@@ -47,6 +47,7 @@ include("FineToCoarseReferenceFEs.jl")
 include("AdaptivityGlues.jl")
 include("AdaptedDiscreteModels.jl")
 include("AdaptedTriangulations.jl")
+include("MacroFEs.jl")
 include("CompositeQuadratures.jl")
 include("EdgeBasedRefinement.jl")
 include("SimplexifyRefinement.jl")

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

@@ -1,37 +1,39 @@
 
-function CellData.CellQuadrature(trian::Triangulation,
-                                 rrules::AbstractVector{<:RefinementRule},
-                                 degree::Integer;
-                                 kwargs...)
+function CellData.CellQuadrature(
+  trian::Triangulation,
+  rrules::AbstractVector{<:RefinementRule},
+  degree::Integer;
+  kwargs...
+)
   return CellData.CellQuadrature(trian,rrules,degree,PhysicalDomain();kwargs...)
 end
 
-function CellData.CellQuadrature(trian::Triangulation,
-                                 rrules::AbstractVector{<:RefinementRule},
-                                 degree::Integer,
-                                 ids::DomainStyle;
-                                 kwargs...)
+function CellData.CellQuadrature(
+  trian::Triangulation,
+  rrules::AbstractVector{<:RefinementRule},
+  degree::Integer,
+  ids::DomainStyle;
+  kwargs...
+)
   cell_quad = lazy_map(rr -> Quadrature(rr,degree;kwargs...),rrules)
   return CellData.CellQuadrature(trian,cell_quad,integration_domain_style=ids)
 end
 
-struct BundleQuadrature{D,T,C <: Table{Point{D,T}},W <: AbstractVector{T}} <: Quadrature{D,T}
-  coordinates::C
-  weights::W
-  name::String
-end
-
-ReferenceFEs.get_coordinates(q::BundleQuadrature) = q.coordinates
-ReferenceFEs.get_weights(q::BundleQuadrature) = q.weights
-ReferenceFEs.get_name(q::BundleQuadrature) = q.name
-
 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
 
-function ReferenceFEs.Quadrature(p::Polytope,::CompositeQuadrature,rr::RefinementRule,degree::Integer;bundle_points::Bool=false,kwargs...)
+function ReferenceFEs.Quadrature(
+  p::Polytope,::CompositeQuadrature,rr::RefinementRule,degree::Integer;kwargs...
+)
   @check p === ReferenceFEs.get_polytope(rr)
   subgrid  = get_ref_grid(rr)
   cellmap  = get_cell_map(rr)
@@ -41,25 +43,42 @@ function ReferenceFEs.Quadrature(p::Polytope,::CompositeQuadrature,rr::Refinemen
   npts = sum(map(num_points,quads))
   WT = eltype(get_weights(first(quads)))
   CT = eltype(get_coordinates(first(quads)))
-  weights     = Vector{WT}(undef,npts)
-  coordinates = Vector{CT}(undef,npts)
-  ptrs = Vector{Int}(undef,length(quads)+1); ptrs[1] = 1;
+  weights = Vector{WT}(undef,npts)
+  cpoints = Vector{CT}(undef,npts)
+  conn = Vector{Vector{Int32}}(undef,length(quads))
   k = 1
   for (iq,q) in enumerate(quads)
     n = num_points(q)
     w = get_weights(q)
     c = get_coordinates(q)
     weights[k:k+n-1] .= w .* measures[iq]
-    coordinates[k:k+n-1] .= (cellmap[iq])(c)
-    ptrs[iq+1] = ptrs[iq] + n
+    cpoints[k:k+n-1] .= (cellmap[iq])(c)
+    conn[iq] = collect(k:k+n-1)
     k = k+n
   end
+  fpoints = map(get_coordinates,quads)
+  ids = FineToCoarseIndices(conn)
+  coordinates = FineToCoarseArray{eltype(cpoints)}(rr,cpoints,fpoints,ids)
+  return GenericQuadrature(coordinates,weights,"Composite quadrature")
+end
 
-  name = "Composite quadrature"
-  if bundle_points
-    coordinates = Table(coordinates,ptrs)
-    return BundleQuadrature(coordinates,weights,name)
-  else
-    return GenericQuadrature(coordinates,weights,name)
-  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
+)
+  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