Merge pull request #1149 from gridap/polytopal

Minor changes to PolytopalQuadrature constructors

Author: JordiManyer

.github/workflows/ci.yml

@@ -44,9 +44,9 @@ jobs:
       - uses: julia-actions/julia-processcoverage@v1
       - uses: codecov/codecov-action@v5
         with:
-          file: lcov.info
           verbose: true
           token: ${{ secrets.CODECOV_TOKEN }}
+          flags: tests
 
   drivers:
     name: Drivers ${{ matrix.version }} - ${{ matrix.os }} - ${{ matrix.arch }} - ${{ github.event_name }}
@@ -69,8 +69,12 @@ jobs:
       - uses: julia-actions/cache@v2
       - uses: julia-actions/julia-buildpkg@v1
       - run: |
-          julia --color=yes --project=. --check-bounds=yes --depwarn=error -e '
-            using Pkg; Pkg.instantiate()'
+          julia --color=yes --project=. --depwarn=yes -e 'using Pkg; Pkg.instantiate()'
       - run: |
-          julia --color=yes --project=. --check-bounds=yes --depwarn=error -e '
-            include("test/GridapTests/runtests.jl")'
+          julia --color=yes --project=. --depwarn=yes -e 'include("test/GridapTests/runtests.jl")'
+      - uses: julia-actions/julia-processcoverage@v1
+      - uses: codecov/codecov-action@v5
+        with:
+          verbose: true
+          token: ${{ secrets.CODECOV_TOKEN }}
+          flags: drivers

NEWS.md

@@ -5,12 +5,17 @@ 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]
+## [0.19.4] - 2025-08-09
 
 ### Added
 
 - Added low-level API `Arrays.invalidate_cache!(cache)` that forces the re-computation of values in `LazyArray` caches. Since PR[#1137](https://github.com/gridap/Gridap.jl/pull/1137).
 - Added support for star-patch integration, i.e patch integration with masked patch boundaries. Since PR[#1138](https://github.com/gridap/Gridap.jl/pull/1138).
+- Added compatibility for `DataStructures.jl` v0.19.0. Since PR[#1149](https://github.com/gridap/Gridap.jl/pull/1149).
+
+### Changed
+
+- Small modifications to `PolytopalQuadrature` to make it more flexible. Since PR[#1149](https://github.com/gridap/Gridap.jl/pull/1149).
 
 ### Fixed
 

Project.toml

@@ -1,7 +1,7 @@
 name = "Gridap"
 uuid = "56d4f2e9-7ea1-5844-9cf6-b9c51ca7ce8e"
 authors = ["Santiago Badia <santiago.badia@monash.edu>", "Francesc Verdugo <f.verdugo.rojano@vu.nl>", "Alberto F. Martin <alberto.f.martin@anu.edu.au>"]
-version = "0.19.3"
+version = "0.19.4"
 
 [deps]
 AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c"
@@ -37,7 +37,7 @@ Aqua = "0.8"
 BSON = "0.3.4"
 BlockArrays = "1"
 Combinatorics = "1"
-DataStructures = "0.18.13"
+DataStructures = "0.18.13, 0.19"
 DocStringExtensions = "0.8.1, 0.9"
 FastGaussQuadrature = "0.4.5, 1"
 FileIO = "1.5"

src/Adaptivity/EdgeBasedRefinement.jl

@@ -421,9 +421,9 @@ function setup_edge_based_rrules(
   # Propagate red/green flags
   # Queue invariant: Cells in queue are RED, every RED cell is only touched once
   q = Queue{Int}()
-  map(c->enqueue!(q,c),cells_to_refine)
+  map(c->push!(q,c),cells_to_refine)
   while !isempty(q)
-    c = dequeue!(q)
+    c = popfirst!(q)
     c_edges = c2e_map[c]
 
     # For each non-red neighboring cell
@@ -439,7 +439,7 @@ function setup_edge_based_rrules(
         cell_color[nbor] = RED + Int8(cell_types[nbor]-1)
         is_red[nbor] = true
         is_refined[nbor_edges] .= true
-        enqueue!(q,nbor)
+        push!(q,nbor)
       end
     end
   end

src/Adaptivity/PolytopalRefinement.jl

@@ -1,4 +1,12 @@
 
+# Notes: 
+# We need to restrict the input to PolytopalDiscreteModel. This would not be strictly 
+# necessary (we can still aggregate other models), but we are currently limited 
+# because we need that all polytopes are positively-oriented rotation systems (like GeneralPolytopes). 
+# Unfortunately our ExtrusionPolytopes are not.
+# To make it work for any model, we would have to reindex the face nodes to ensure the faces
+# are correctly oriented.
+
 function coarsen(model::Geometry.PolytopalDiscreteModel,ptopo::Geometry.PatchTopology)
   new_polys, new_connectivity = generate_patch_polytopes(model,ptopo)
 

src/FESpaces/Assemblers.jl

@@ -571,7 +571,7 @@ function collect_and_merge_cell_matrix_and_vector(contributions...)
 end
 
 function merge_assembly_data(data...)
-  @assert allequal(length,data)
+  @assert allequal(map(length,data))
   map(vcat,data...)
 end
 

src/ReferenceFEs/ExtrusionPolytopes.jl

@@ -305,25 +305,25 @@ function positive_simplexify_hypercube(dim::Integer)
 end
 
 function simplexify_hypercube(dim::Int)
-    @assert dim ≥ 0
-    # Determine simplices recursively
-    simplices = Vector{Int}[]
-    corner = 0
-    vertices = [corner]
-    next_corner!(simplices, dim, vertices, corner)
-    # Correct vertex numbering
-    for s in simplices
-        s .+= 1
+  @assert dim ≥ 0
+  # Determine simplices recursively
+  simplices = Vector{Int}[]
+  corner = 0
+  vertices = [corner]
+  next_corner!(simplices, dim, vertices, corner)
+  # Correct vertex numbering
+  for s in simplices
+    s .+= 1
+  end
+  # Check output
+  @assert length(simplices) == factorial(dim)
+  for s in simplices
+    @assert length(s) == dim+1
+    for v in s
+      @assert 1 ≤ v ≤ 2^dim
     end
-    # Check output
-    @assert length(simplices) == factorial(dim)
-    for s in simplices
-        @assert length(s) == dim+1
-        for v in s
-            @assert 1 ≤ v ≤ 2^dim
-        end
-    end
-    return simplices
+  end
+  return simplices
 end
 
 """
@@ -346,22 +346,22 @@ the top. The algorithm below finds all possible paths.
 """
 function next_corner!(simplices::Vector{Vector{Int}}, dim::Int,
                       vertices::Vector{Int}, corner::Int)
-    @assert count_ones(corner) == length(vertices) - 1
-    if length(vertices) == dim + 1
-        # We have all vertices; save the simplex
-        push!(simplices, vertices)
-        return
-    end
-    # Loop over all neighbouring corners
-    for d in 1:dim
-        bit = 1 << (d-1)
-        if (corner & bit) == 0
-            new_corner = corner | bit
-            new_vertices = [vertices; new_corner]
-            next_corner!(simplices, dim, new_vertices, new_corner)
-        end
-    end
+  @assert count_ones(corner) == length(vertices) - 1
+  if length(vertices) == dim + 1
+    # We have all vertices; save the simplex
+    push!(simplices, vertices)
     return
+  end
+  # Loop over all neighbouring corners
+  for d in 1:dim
+    bit = 1 << (d-1)
+    if (corner & bit) == 0
+      new_corner = corner | bit
+      new_vertices = [vertices; new_corner]
+      next_corner!(simplices, dim, new_vertices, new_corner)
+    end
+  end
+  return
 end
 
 function Base.show(io::IO,p::ExtrusionPolytope)

src/ReferenceFEs/PolytopalQuadratures.jl

@@ -1,26 +1,36 @@
 
-struct PolytopalQuadrature{D,T,P<:GeneralPolytope{D},Q<:Quadrature{D,T}} <: Quadrature{D,T}
-  poly :: P
-  quad :: Q
-  conn :: Vector{Vector{Int32}}
+struct PolytopalQuadrature{D,Dp,T,Q<:Quadrature{D,T}} <: Quadrature{D,T}
+  quad   :: Q
+  conn   :: Vector{Vector{Int32}}
+  coords :: Vector{Point{Dp,T}}
+  function PolytopalQuadrature(
+    quad::Quadrature{D,T}, 
+    conn::Vector{Vector{Int32}}, 
+    coords::Vector{Point{Dp,T}}
+  ) where {D,Dp,T}
+    Q = typeof(quad)
+    new{D,Dp,T,Q}(quad, conn, coords)
+  end
 end
 
 function Quadrature(poly::GeneralPolytope{D},args...;kwargs...) where D
   conn, simplex = simplexify(poly)
   quad = Quadrature(simplex,args...;kwargs...)
-  PolytopalQuadrature(poly,quad,conn)
+  coords = get_vertex_coordinates(poly)
+  PolytopalQuadrature(quad,conn,coords)
 end
 
 function Quadrature(poly::GeneralPolytope{D},quad::Quadrature{D}) where D
   conn, simplex = simplexify(poly)
+  coords = get_vertex_coordinates(poly)
   @check get_polytope(quad) == simplex
-  PolytopalQuadrature(poly,quad,conn)
+  PolytopalQuadrature(quad,conn,coords)
 end
 
 # Quadrature API
 
-function get_name(q::PolytopalQuadrature{D,T,P,Q}) where {D,T,P,Q}
-  "PolytopalQuadrature{$(D),$(T),$(P),$(Q)}"
+function get_name(q::PolytopalQuadrature{D,T,Q}) where {D,T,Q}
+  "PolytopalQuadrature{$(D),$(T),$(Q)}"
 end
 
 get_coordinates(q::PolytopalQuadrature) = evaluate(PQuadCoordsMap(),q)
@@ -40,8 +50,7 @@ end
 struct PQuadCoordsMap <: Map end
 
 function Arrays.return_cache(::PQuadCoordsMap, q::PolytopalQuadrature)
-  conn = q.conn
-  coords = get_vertex_coordinates(q.poly)
+  conn, coords = q.conn, q.coords
   cmap = affine_map(Tuple(coords[first(conn)]))
 
   x = get_coordinates(q.quad)
@@ -54,8 +63,7 @@ end
 
 function Arrays.evaluate!(cache,::PQuadCoordsMap, q::PolytopalQuadrature)
   y_cache, cmap_cache = cache
-  coords = get_vertex_coordinates(q.poly)
-  conn = q.conn
+  conn, coords = q.conn, q.coords
 
   x = get_coordinates(q.quad)
   setsize!(y_cache,(length(x)*length(conn),))
@@ -78,8 +86,7 @@ function Arrays.return_cache(::PQuadWeightsMap, q::PolytopalQuadrature)
 end
 
 function Arrays.evaluate!(cache,::PQuadWeightsMap, q::PolytopalQuadrature)
-  conn = q.conn
-  coords = get_vertex_coordinates(q.poly)
+  conn, coords = q.conn, q.coords
 
   w = get_weights(q.quad)
   setsize!(cache,(length(w)*length(conn),))

src/ReferenceFEs/TensorProductQuadratures.jl

@@ -19,16 +19,16 @@ end
 # Low level constructor
 
 function _tensor_product_legendre(degrees;T::Type{<:AbstractFloat}=Float64)
-    D = length(degrees)
-    npoints = [ ceil(Int,(degrees[i]+1.0)/2.0) for i in 1:D ]
-    quads = [ gauss(T, npoints[i]) for i in 1:D ]
-    for i in 1:D
-      quads[i][1] .+= 1;
-      quads[i][1] .*= 1.0/2.0
-      quads[i][2] .*= 1.0/2.0
-    end
-    (coords, weights) = _tensor_product(Point{D,T},quads,npoints)
-    GenericQuadrature(coords,weights,"Tensor product of 1d Gauss-Legendre quadratures of degrees $degrees")
+  D = length(degrees)
+  npoints = [ ceil(Int,(degrees[i]+1.0)/2.0) for i in 1:D ]
+  quads = [ gauss(T, npoints[i]) for i in 1:D ]
+  for i in 1:D
+    quads[i][1] .+= 1;
+    quads[i][1] .*= 1.0/2.0
+    quads[i][2] .*= 1.0/2.0
+  end
+  (coords, weights) = _tensor_product(Point{D,T},quads,npoints)
+  GenericQuadrature(coords,weights,"Tensor product of 1d Gauss-Legendre quadratures of degrees $degrees")
 end
 
 # Helpers