Gridap 0.20

.github/workflows/CompatHelper.yml

@@ -17,7 +17,8 @@ jobs:
         run: julia -e 'import Pkg; Pkg.add("CompatHelper")'
       - name: Run CompatHelper
         # Opens a pull request when a dependency releases a version outside [compat].
-        # Gridap is pinned (=0.19.9) on purpose; its PRs can be closed until HyperFEM upgrades.
+        # Gridap is restricted to one minor series (0.20) on purpose: the engine reads Gridap's internals,
+        # so a PR for a new Gridap series needs the full test suite before it is merged.
         run: julia -e 'import CompatHelper; CompatHelper.main()'
         env:
           GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }}

CHANGELOG.md

@@ -1,5 +1,11 @@
 # Changelog
 
+## Unreleased
+
+- **Gridap 0.20** (breaking: Gridap 0.19 is no longer supported; `[compat] Gridap = "0.20"`). H(div)
+  (Raviart–Thomas) spaces are read from the cell basis of Gridap 0.20, where the sign flips are a diagonal
+  change of basis per cell. The rest of the engine needed no change.
+
 ## 0.2.3
 
 - Second derivatives (`∇∇`, `Δ`) on skeleton integrals: `jump(∇∇(u))`, `mean(∇∇(u))`, `mean(Δ(u))`, … of

Project.toml

@@ -14,7 +14,7 @@ StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 
 [compat]
 ForwardDiff = "1.3"
-Gridap = "=0.19.9"
+Gridap = "0.20"
 Krylov = "0.10"
 LinearAlgebra = "1.12"
 RuntimeGeneratedFunctions = "0.5"

README.md

@@ -90,7 +90,7 @@ Until the package appears in the General registry, add it from the repository in
 pkg> add https://github.com/MultiSimOLab/HyperFEMKernels.jl
 ```
 
-It requires Julia ≥ 1.12 and Gridap 0.19.9 (the same version as HyperFEM).
+It requires Julia ≥ 1.12 and Gridap 0.20.
 
 ## Usage
 
@@ -326,7 +326,7 @@ run with both backends as part of the test suite.
   are not covered yet (use `LaggedAssembled` with HyperFEM's block preconditioners).
 - Shared-memory parallelism (threads). Distributed runs with GridapDistributed and GPU execution
   are design targets: kernels are allocation-free and states use a GPU-friendly layout.
-- Gridap is pinned to 0.19.9, because the engine reads Gridap's internal expression tree.
+- Gridap is restricted to the 0.20 series, because the engine reads Gridap's internal expression tree.
 
 ## Testing
 

benchmark/solvers/Project.toml

@@ -18,4 +18,4 @@ HyperFEM = {path = "../../../HyperFEM.jl"}
 HyperFEMKernels = {path = "../.."}
 
 [compat]
-Gridap = "=0.19.9"
+Gridap = "0.20"

benchmark/solvers/setup.jl

@@ -9,6 +9,6 @@ root = normpath(joinpath(@__DIR__, "..", ".."))
 Pkg.develop([PackageSpec(path=root), PackageSpec(path=joinpath(root, "..", "HyperFEM.jl"))])
 Pkg.add(["Gridap", "GridapGmsh", "GridapPETSc", "GridapSolvers", "SparseMatricesCSR", "SparseArrays", "LinearAlgebra",
          "Printf", "Dates", "DelimitedFiles"])
-Pkg.compat("Gridap", "=0.19.9")
+Pkg.compat("Gridap", "0.20")
 Pkg.instantiate()
 Pkg.status()

docs/README.md

@@ -20,5 +20,5 @@ with [Registrator](https://github.com/JuliaRegistries/Registrator.jl); GitHub re
 Requirements on the repository: the JuliaRegistrator and TagBot GitHub apps installed
 (organisation settings → GitHub Apps), the repository public, `LICENSE` in the root.
 `CompatHelper.yml` opens pull requests when dependencies release versions outside `[compat]`;
-Gridap is pinned on purpose (`=0.19.9`, HyperFEM's version), so its pull requests are closed until
-HyperFEM upgrades.
+Gridap is restricted to one minor series on purpose (`0.20`: the engine reads Gridap's internals), so a pull
+request for a new Gridap series needs the full test suite before it is merged.

src/tables.jl

@@ -185,14 +185,22 @@ wcols(::FieldTable{D}) where {D} = 1 + D + nhess(D)
 # Gridap builds the cell basis of a div-conforming space as the contravariant Piola map of the
 # reference shape functions with a sign flip on the dofs of the facets whose cell is the slave:
 #   φ_j(x) = σ_j / |det Jt| · φ̂_j(ξ) ⋅ Jt,   ∇φ_j = Jt⁻¹ ⋅ (σ_j / |det Jt| · ∇̂φ̂_j ⋅ Jt)
-# (Gridap 0.19 takes the Jacobian as constant in the gradient; the formula above is reproduced
+# (Gridap takes the Jacobian as constant in the gradient; the formula above is reproduced
 # exactly, also on non-affine cells, so that kernel and Gridap agree to round-off). The reference
 # functions are the same in every cell; the signs are a dense ±1 per (dof, cell).
+#
+# The lazy array that holds this basis (Gridap ≥ 0.20) is
+#   lazy_map(linear_combination, M, lazy_map(Broadcasting(Operation(ContraVariantPiolaMap())), φ̂, Jt))
+# with the change of basis `M` a `Diagonal` of ±1 per cell (`NormalSignMap`).
+
+"A lazy array of contravariant Piola maps of reference fields."
+is_piola_map(a) = a isa LazyArray && a.maps isa Fill && a.maps.value isa Broadcasting{<:Operation{ContraVariantPiolaMap}}
 
-"The cell basis of a div-conforming space (the lazy Piola-map array), or `nothing`."
+"The cell basis of a div-conforming space (the signed lazy Piola-map array), or `nothing`."
 function piola_basis(space::SingleFieldFESpace)
     a = get_data(get_fe_basis(space))
-    (a isa LazyArray && a.maps isa Fill && a.maps.value isa Broadcasting{<:Operation{ContraVariantPiolaMap}}) ? a : nothing
+    (a isa LazyArray && a.maps isa Fill && a.maps.value === linear_combination &&
+     length(a.args) == 2 && is_piola_map(a.args[2])) ? a : nothing
 end
 is_piola(space::SingleFieldFESpace) = piola_basis(space) !== nothing
 
@@ -214,7 +222,7 @@ function PiolaTable(space::SingleFieldFESpace, xq::AbstractVector, D::Int)
     a === nothing && error("HyperFEMKernels: not a Piola-mapped space")
     length(first(xq)) == D ||
         error("HyperFEMKernels: H(div) spaces on a manifold triangulation (cells of lower dimension than the space) are not supported")
-    ref = uniform_value(a.args[1], "the reference H(div) basis")
+    ref = uniform_value(a.args[2].args[1], "the reference H(div) basis")
     vals = evaluate(ref, xq)                         # np × nd, VectorValue{D}
     grads = evaluate(Broadcasting(∇)(ref), xq)       # np × nd, TensorValue{D,D}
     np, nd = size(vals)
@@ -234,20 +242,22 @@ field_table(space::SingleFieldFESpace, xq::AbstractVector, D::Int) =
     cell_signs(space, nd) -> Matrix{Int8} (nd × ncells)
 
 Sign of every cell dof of an H(div) space (−1 on the facet dofs of the slave cell, read from the
-sign flips of Gridap's cell basis); a 0 × 0 matrix for other spaces (no signs).
+sign changes of Gridap's cell basis); a 0 × 0 matrix for other spaces (no signs).
 """
 function cell_signs(space::SingleFieldFESpace, nd::Int)
     a = piola_basis(space)
     a === nothing && return Matrix{Int8}(undef, 0, 0)
-    flips = a.args[4].args[1]                        # Bool per cell dof (lazy, per cell)
-    nc = length(flips)
+    changes = a.args[1]                              # Diagonal of ±1 per cell (lazy)
+    nc = length(changes)
     s = Matrix{Int8}(undef, nd, nc)
-    cache = array_cache(flips)
+    cache = array_cache(changes)
     for cell in 1:nc
-        f = getindex!(cache, flips, cell)
-        length(f) == nd || error("HyperFEMKernels: cells with different numbers of dofs are not supported yet")
+        M = getindex!(cache, changes, cell)
+        M isa Diagonal || error("HyperFEMKernels: unexpected change of basis $(typeof(M)) in an H(div) space")
+        length(M.diag) == nd || error("HyperFEMKernels: cells with different numbers of dofs are not supported yet")
         for j in 1:nd
-            s[j, cell] = f[j] ? Int8(-1) : Int8(1)
+            abs(M.diag[j]) == 1 || error("HyperFEMKernels: H(div) spaces with scaled dofs (`scale_dof`) are not supported")
+            s[j, cell] = Int8(M.diag[j])
         end
     end
     s