Remove the derivable Geometry.subspaces field

The single-level Geometry.subspaces Dict was consumed only by amg, for the
fine-level :full / :uniform embeddings — which are exactly I(n_doubled) and
ones(n_doubled, 1). Derive them in _assemble_amg_dicts instead and drop the
field plus its M_sub type parameter.

Updates every Geometry{...} signature (the M_sub slot, always just before the
discretization parameter, is removed), all construction sites (dropping both
the type parameter and the subspaces argument), the now-dead single-level
subspace dicts in the mesh constructors, and the CUDA-extension Geometry
conversions. MultiGrid.subspaces (the per-level hierarchy) is unchanged;
spectral amg never read the field.

Author: Sebastien Loisel

ext/MultiGridBarrierCUDAExt/block_ops.jl

@@ -1264,8 +1264,8 @@ function MultiGridBarrier._structurize_multigrid(
     XT = typeof(geom.x); WT = typeof(geom.w)
     M_op_type = valtype(operators_new)
     DiscT = typeof(geom.discretization)
-    new_geom = MultiGridBarrier.Geometry{T,XT,WT,M_op_type,M_sub,DiscT}(
-        geom.discretization, geom.x, geom.w, geom.subspaces, operators_new)
+    new_geom = MultiGridBarrier.Geometry{T,XT,WT,M_op_type,DiscT}(
+        geom.discretization, geom.x, geom.w, operators_new)
 
     return MultiGridBarrier.MultiGrid(new_geom, subspaces_new, refine_new, coarsen_new)
 end

ext/MultiGridBarrierCUDAExt/conversion.jl

@@ -36,7 +36,7 @@ Convert a native (CPU) `Geometry` to CUDA GPU types.
 - `SparseMatrixCSC{T,Int}` → `CuSparseMatrixCSR{T,Int32}`
 - Dense `Matrix{T}` operators → `CuMatrix{T}` (spectral)
 """
-function MultiGridBarrier.native_to_cuda(g::Geometry{T, Array{T,3}, Vector{T}, SparseMatrixCSC{T,Int}, SparseMatrixCSC{T,Int}, Discretization};
+function MultiGridBarrier.native_to_cuda(g::Geometry{T, Array{T,3}, Vector{T}, SparseMatrixCSC{T,Int}, Discretization};
                                          Ti::Type{<:Integer}=Int32) where {T, Discretization}
     x_cuda = CuArray{T,3}(g.x)
     w_cuda = CuVector{T}(g.w)
@@ -51,21 +51,15 @@ function MultiGridBarrier.native_to_cuda(g::Geometry{T, Array{T,3}, Vector{T}, S
         operators_cuda[key] = convert_sparse(g.operators[key])
     end
 
-    subspaces_cuda = Dict{Symbol, MType}()
-    for key in sort(collect(keys(g.subspaces)))
-        subspaces_cuda[key] = convert_sparse(g.subspaces[key])
-    end
-
-    Geometry{T, CuArray{T,3}, CuVector{T}, MType, MType, Discretization}(
-        g.discretization, x_cuda, w_cuda, subspaces_cuda, operators_cuda)
+    Geometry{T, CuArray{T,3}, CuVector{T}, MType, Discretization}(
+        g.discretization, x_cuda, w_cuda, operators_cuda)
 end
 
 # BlockDiag-operator variant: FEM Geometry whose `operators` came in via
 # `subdivide` (or `geometric_mg(...; structured=true).geometry`). Carry the
 # BlockDiag to GPU as a CuArray-backed BlockDiag; subspaces stay sparse CSR.
 function MultiGridBarrier.native_to_cuda(g::Geometry{T, Array{T,3}, Vector{T},
                                                       MultiGridBarrier.BlockDiag{T, A3},
-                                                      SparseMatrixCSC{T,Int},
                                                       Discretization};
                                          Ti::Type{<:Integer}=Int32) where {T, A3<:Array{T,3}, Discretization}
     x_cuda = CuArray{T,3}(g.x)
@@ -76,24 +70,18 @@ function MultiGridBarrier.native_to_cuda(g::Geometry{T, Array{T,3}, Vector{T},
         MultiGridBarrier.BlockDiag{T, typeof(gpu_data)}(op.p, op.q, op.N, gpu_data)
     end
 
-    sparse_to_gpu = op -> CuSparseMatrixCSR(
-        SparseMatrixCSC{T,Ti}(op.m, op.n, Ti.(op.colptr), Ti.(op.rowval), op.nzval))
-
     sample_op = op_to_gpu(first(values(g.operators)))
     OpType = typeof(sample_op)
-    SubType = CuSparseMatrixCSR{T,Ti}
 
     operators_cuda = Dict{Symbol, OpType}(
         key => op_to_gpu(g.operators[key]) for key in keys(g.operators))
-    subspaces_cuda = Dict{Symbol, SubType}(
-        key => sparse_to_gpu(g.subspaces[key]) for key in keys(g.subspaces))
 
-    Geometry{T, CuArray{T,3}, CuVector{T}, OpType, SubType, Discretization}(
-        g.discretization, x_cuda, w_cuda, subspaces_cuda, operators_cuda)
+    Geometry{T, CuArray{T,3}, CuVector{T}, OpType, Discretization}(
+        g.discretization, x_cuda, w_cuda, operators_cuda)
 end
 
 # Dense spectral variant.
-function MultiGridBarrier.native_to_cuda(g::Geometry{T, Array{T,3}, Vector{T}, Matrix{T}, Matrix{T}, Discretization};
+function MultiGridBarrier.native_to_cuda(g::Geometry{T, Array{T,3}, Vector{T}, Matrix{T}, Discretization};
                                          kwargs...) where {T, Discretization}
     x_cuda = CuArray{T,3}(g.x)
     w_cuda = CuVector{T}(g.w)
@@ -103,13 +91,8 @@ function MultiGridBarrier.native_to_cuda(g::Geometry{T, Array{T,3}, Vector{T}, M
         operators_cuda[key] = CuMatrix{T}(g.operators[key])
     end
 
-    subspaces_cuda = Dict{Symbol, CuMatrix{T}}()
-    for key in sort(collect(keys(g.subspaces)))
-        subspaces_cuda[key] = CuMatrix{T}(g.subspaces[key])
-    end
-
-    Geometry{T, CuArray{T,3}, CuVector{T}, CuMatrix{T}, CuMatrix{T}, Discretization}(
-        g.discretization, x_cuda, w_cuda, subspaces_cuda, operators_cuda)
+    Geometry{T, CuArray{T,3}, CuVector{T}, CuMatrix{T}, Discretization}(
+        g.discretization, x_cuda, w_cuda, operators_cuda)
 end
 
 # ============================================================================
@@ -204,7 +187,7 @@ end
 
 # Sparse FEM Geometry (CSR operators).
 function MultiGridBarrier.cuda_to_native(g::Geometry{T, <:CuArray{T,3}, <:CuVector{T},
-                                                      <:CuSparseMatrixCSR{T}, <:CuSparseMatrixCSR{T},
+                                                      <:CuSparseMatrixCSR{T},
                                                       Discretization}) where {T, Discretization}
     x_native = Array{T,3}(Array(g.x))
     w_native = Vector{T}(Array(g.w))
@@ -219,18 +202,13 @@ function MultiGridBarrier.cuda_to_native(g::Geometry{T, <:CuArray{T,3}, <:CuVect
         operators_native[key] = convert_back(g.operators[key])
     end
 
-    subspaces_native = Dict{Symbol, SparseMatrixCSC{T,Ti}}()
-    for key in sort(collect(keys(g.subspaces)))
-        subspaces_native[key] = convert_back(g.subspaces[key])
-    end
-
-    Geometry{T, Array{T,3}, Vector{T}, SparseMatrixCSC{T,Ti}, SparseMatrixCSC{T,Ti}, Discretization}(
-        g.discretization, x_native, w_native, subspaces_native, operators_native)
+    Geometry{T, Array{T,3}, Vector{T}, SparseMatrixCSC{T,Ti}, Discretization}(
+        g.discretization, x_native, w_native, operators_native)
 end
 
 # Dense spectral Geometry.
 function MultiGridBarrier.cuda_to_native(g::Geometry{T, <:CuArray{T,3}, <:CuVector{T},
-                                                      <:CuMatrix{T}, <:CuMatrix{T},
+                                                      <:CuMatrix{T},
                                                       Discretization}) where {T, Discretization}
     x_native = Array{T,3}(Array(g.x))
     w_native = Vector{T}(Array(g.w))
@@ -240,18 +218,13 @@ function MultiGridBarrier.cuda_to_native(g::Geometry{T, <:CuArray{T,3}, <:CuVect
         operators_native[key] = Matrix{T}(Array(g.operators[key]))
     end
 
-    subspaces_native = Dict{Symbol, Matrix{T}}()
-    for key in sort(collect(keys(g.subspaces)))
-        subspaces_native[key] = Matrix{T}(Array(g.subspaces[key]))
-    end
-
-    Geometry{T, Array{T,3}, Vector{T}, Matrix{T}, Matrix{T}, Discretization}(
-        g.discretization, x_native, w_native, subspaces_native, operators_native)
+    Geometry{T, Array{T,3}, Vector{T}, Matrix{T}, Discretization}(
+        g.discretization, x_native, w_native, operators_native)
 end
 
 # Structured FEM Geometry (block ops).
 function MultiGridBarrier.cuda_to_native(g::Geometry{T, <:CuArray{T,3}, <:CuVector{T},
-                                                      <:Any, <:Any, Discretization}) where {T, Discretization}
+                                                      <:Any, Discretization}) where {T, Discretization}
     x_native = Array{T,3}(Array(g.x))
     w_native = Vector{T}(Array(g.w))
 
@@ -267,13 +240,8 @@ function MultiGridBarrier.cuda_to_native(g::Geometry{T, <:CuArray{T,3}, <:CuVect
         operators_native[key] = convert_to_native(g.operators[key])
     end
 
-    subspaces_native = Dict{Symbol, SparseMatrixCSC{T,Ti}}()
-    for key in sort(collect(keys(g.subspaces)))
-        subspaces_native[key] = convert_to_native(g.subspaces[key])
-    end
-
-    Geometry{T, Array{T,3}, Vector{T}, SparseMatrixCSC{T,Ti}, SparseMatrixCSC{T,Ti}, Discretization}(
-        g.discretization, x_native, w_native, subspaces_native, operators_native)
+    Geometry{T, Array{T,3}, Vector{T}, SparseMatrixCSC{T,Ti}, Discretization}(
+        g.discretization, x_native, w_native, operators_native)
 end
 
 # MultiGrid cuda → native (sparse FEM).

src/AlgebraicMultiGridBarrier.jl

@@ -203,7 +203,7 @@ Base.showerror(io::IO, e::MGBConvergenceFailure) = print(io, "MGBConvergenceFail
 end
 
 """
-    Geometry{T,X<:AbstractArray{T,3},W,M_op,M_sub,Discretization}
+    Geometry{T,X<:AbstractArray{T,3},W,M_op,Discretization}
 
 Single-level container for discretization geometry. Holds only the fine-level mesh and
 operators — no multigrid hierarchy. Use `amg(geom)` to attach an algebraic-multigrid
@@ -215,7 +215,6 @@ Type parameters
 - `X<:AbstractArray{T,3}`: type of the mesh tensor `x` (typically `Array{T,3}`).
 - `W`: type of the weight storage `w` (typically `Vector{T}`).
 - `M_op`: matrix type for operators (e.g. `SparseMatrixCSC{T,Int}`, `BlockDiag{T}`).
-- `M_sub`: matrix type for subspace embeddings (e.g. `SparseMatrixCSC{T,Int}`).
 - `Discretization`: front-end descriptor (e.g. `FEM1D{T}`, `FEM2D_P2{T}`, `SPECTRAL1D{T}`).
 
 Fields
@@ -225,14 +224,12 @@ Fields
   flat layout via `reshape(x, V*N, D)` (zero-copy). Spectral discretizations
   use `N = 1` (a single notional "element" comprising every Chebyshev node).
 - `w::W`: quadrature weights matching the flattened node order (length `V*N`).
-- `subspaces::Dict{Symbol,M_sub}`: fine-level selection/embedding matrices (one per symbol).
 - `operators::Dict{Symbol,M_op}`: fine-level discrete operators (e.g. `:id`, `:dx`, `:dy`).
 """
-struct Geometry{T,X<:AbstractArray{T,3},W,M_op,M_sub,Discretization}
+struct Geometry{T,X<:AbstractArray{T,3},W,M_op,Discretization}
     discretization::Discretization
     x::X
     w::W
-    subspaces::Dict{Symbol,M_sub}
     operators::Dict{Symbol,M_op}
 end
 
@@ -532,7 +529,7 @@ function amg end
 # nominally rides is immaterial. Each `(sym, nodes)` entry of `dirichlet_nodes`
 # adds one zero-trace continuous subspace built by the discretization-specific
 # `build_dirichlet(nodes) -> (refine, coarsen, sub)` closure.
-function _assemble_amg_dicts(::Type{T}, geom,
+function _assemble_amg_dicts(::Type{T}, geom, n_doubled::Int,
         dirichlet_nodes::Dict{Symbol,Vector{Tuple{Int,Int}}},
         refine_full::Vector{SparseMatrixCSC{T,Int}},
         coarsen_full::Vector{SparseMatrixCSC{T,Int}},
@@ -544,8 +541,11 @@ function _assemble_amg_dicts(::Type{T}, geom,
         sub_full[kk]    = sparse(one(T) * I, sizes_full[kk], sizes_full[kk])
         sub_uniform[kk] = sparse(ones(T, sizes_full[kk], 1))
     end
-    sub_full[L_full]    = SparseMatrixCSC{T,Int}(geom.subspaces[:full])
-    sub_uniform[L_full] = SparseMatrixCSC{T,Int}(geom.subspaces[:uniform])
+    # Fine-level (level-L) embeddings: `:full` is the entire broken space
+    # (identity) and `:uniform` is the constant column. Both depend only on the
+    # broken-DOF count `n_doubled`.
+    sub_full[L_full]    = sparse(one(T) * I, n_doubled, n_doubled)
+    sub_uniform[L_full] = sparse(ones(T, n_doubled, 1))
 
     subspaces = Dict{Symbol,Vector{SparseMatrixCSC{T,Int}}}(:full => sub_full, :uniform => sub_uniform)
     refine_d  = Dict{Symbol,Vector{SparseMatrixCSC{T,Int}}}(:full => refine_full, :uniform => refine_full)
@@ -2608,7 +2608,7 @@ struct MGBSOL{T,X,W,Discretization,G}
     log::String
     geometry::G
 end
-function MGBSOL(z::X, sf, sm, log::String, geometry::Geometry{T,<:Any,W,<:Any,<:Any,Discretization}) where {T,X,W,Discretization}
+function MGBSOL(z::X, sf, sm, log::String, geometry::Geometry{T,<:Any,W,<:Any,Discretization}) where {T,X,W,Discretization}
     MGBSOL{T,X,W,Discretization,typeof(geometry)}(z, sf, sm, log, geometry)
 end
 plot(sol::MGBSOL,k::Int=1;kwargs...) = plot(sol.geometry,sol.z[:,k];kwargs...)

src/BlockMatrices.jl

@@ -1110,10 +1110,10 @@ end
 function _default_block_size end
 
 # Rebuild a Geometry with a new (structured) operators dict, preserving everything else.
-function _replace_operators(geom::Geometry{T,X,W,<:Any,M_sub,Disc},
-                            operators_new::Dict{Symbol,M_op_new}) where {T,X,W,M_sub,Disc,M_op_new}
-    Geometry{T,X,W,M_op_new,M_sub,Disc}(
-        geom.discretization, geom.x, geom.w, geom.subspaces, operators_new)
+function _replace_operators(geom::Geometry{T,X,W,<:Any,Disc},
+                            operators_new::Dict{Symbol,M_op_new}) where {T,X,W,Disc,M_op_new}
+    Geometry{T,X,W,M_op_new,Disc}(
+        geom.discretization, geom.x, geom.w, operators_new)
 end
 
 # Structurize a MultiGrid: convert geometry's operators to BlockDiag,

src/Mesh3d/Geometry.jl

@@ -126,14 +126,10 @@ function _geometric_fem3d_mg(::Type{T}=Float64; L::Int=2,
     end
 
     x_fine = reshape(meshes[L], sk3, N_fine, 3)
-    geom = Geometry{T, Array{T,3}, Vector{T}, SparseMatrixCSC{T,Int}, SparseMatrixCSC{T,Int}, FEM3D{T}}(
+    geom = Geometry{T, Array{T,3}, Vector{T}, SparseMatrixCSC{T,Int}, FEM3D{T}}(
         disc,
         x_fine,
         weights[L],
-        Dict{Symbol,SparseMatrixCSC{T,Int}}(
-            :full      => subspaces[:full][end],
-            :uniform   => subspaces[:uniform][end],
-            :dirichlet => subspaces[:dirichlet][end]),
         ops,
     )
     return MultiGrid(geom, subspaces, refine_ops, coarsen_ops)
@@ -270,14 +266,8 @@ function _fem3d_structured(disc::FEM3D{T}, meshes, weights, L, k, ref_el) where
     s_   = k + 1
     sk3_ = s_^3
     x_fine = reshape(meshes[L], sk3_, div(size(meshes[L], 1), sk3_), 3)
-    geom = Geometry{T, Array{T,3}, Vector{T}, BlockDiag{T,Array{T,3}},
-                    SparseMatrixCSC{T,Int}, FEM3D{T}}(
-        disc, x_fine, weights[L],
-        Dict{Symbol,SparseMatrixCSC{T,Int}}(
-            :full      => subspaces[:full][end],
-            :uniform   => subspaces[:uniform][end],
-            :dirichlet => subspaces[:dirichlet][end]),
-        operators)
+    geom = Geometry{T, Array{T,3}, Vector{T}, BlockDiag{T,Array{T,3}}, FEM3D{T}}(
+        disc, x_fine, weights[L], operators)
     return MultiGrid(geom, subspaces, refine, coarsen)
 end
 
@@ -388,7 +378,7 @@ end
 Evaluate the finite element field `u` at point `x_eval`. Returns the value and a flag
 indicating whether the point was found.
 """
-function evaluate_field(g::Geometry{T,X,W,<:Any,<:Any,FEM3D{T}}, u::Vector{T}, x_eval::Vector{T}) where {T,X,W}
+function evaluate_field(g::Geometry{T,X,W,<:Any,FEM3D{T}}, u::Vector{T}, x_eval::Vector{T}) where {T,X,W}
     k = g.discretization.k
     n_nodes_per_elem = (k+1)^3
     x_flat = _xflat(g.x)

src/Mesh3d/Plotting.jl

@@ -97,7 +97,7 @@ plot(sol; slice=(normal=[1,1,1], origin=[0,0,0]))
 plot(sol; plotter=(window_size=(1600, 1200),))
 ```
 """
-function plot(geo::Geometry{T,X,W,<:Any,<:Any,FEM3D{T}}, u::Vector{T};
+function plot(geo::Geometry{T,X,W,<:Any,FEM3D{T}}, u::Vector{T};
                        plotter::NamedTuple=(window_size=(800, 600),),
                        volume=(;),
                        scalar_bar_args=(title="",position_x=0.6,position_y=0.0,width=0.35,height=0.05,use_opacity=false),
@@ -269,7 +269,7 @@ function create_vtk_cells(k::Int, n_total_nodes::Int)
 end
 
 """
-    plot(M::Geometry{T,X,W,<:Any,<:Any,FEM3D{T}}, ts::AbstractVector, U::Matrix{T}; kwargs...)
+    plot(M::Geometry{T,X,W,<:Any,FEM3D{T}}, ts::AbstractVector, U::Matrix{T}; kwargs...)
 
 Create an animated 3D visualization from a time series of solutions.
 
@@ -293,7 +293,7 @@ base64-encoded data in an HTML5 video tag.
 Color limits (`clim`) and `isosurfaces` are automatically computed from the global
 range of U across all frames to ensure consistent visualization throughout the animation.
 """
-function plot(M::Geometry{T,X,W,<:Any,<:Any,FEM3D{T}}, ts::AbstractVector, U::Matrix{T};
+function plot(M::Geometry{T,X,W,<:Any,FEM3D{T}}, ts::AbstractVector, U::Matrix{T};
               frame_time::Real = max(0.001, minimum(diff(ts))),
               kwargs...) where {T,X,W}
 
@@ -390,7 +390,7 @@ Plot an animated 3D visualization of a parabolic solution (FEM3D).
 # Returns
 - `HTML5anim`: An HTML5 video that displays in Jupyter notebooks.
 """
-function plot(sol::ParabolicSOL{T,X,W,<:Any,<:Geometry{T,<:Any,<:Any,<:Any,<:Any,FEM3D{T}}}, k::Int=1; kwargs...) where {T,X,W}
+function plot(sol::ParabolicSOL{T,X,W,<:Any,<:Geometry{T,<:Any,<:Any,<:Any,FEM3D{T}}}, k::Int=1; kwargs...) where {T,X,W}
     return plot(sol.geometry, collect(sol.ts), hcat([sol.u[j][:, k] for j=1:length(sol.ts)]...); kwargs...)
 end
 

src/Parabolic.jl

@@ -49,7 +49,7 @@ struct ParabolicSOL{T,X,W,Discretization,G}
     ts::Vector{T}
     u::Vector{X}
 end
-function ParabolicSOL(geometry::Geometry{T,<:Any,W,<:Any,<:Any,Discretization}, ts, u::Vector{X}) where {T,X,W,Discretization}
+function ParabolicSOL(geometry::Geometry{T,<:Any,W,<:Any,Discretization}, ts, u::Vector{X}) where {T,X,W,Discretization}
     ParabolicSOL{T,X,W,Discretization,typeof(geometry)}(geometry, collect(T, ts), u)
 end
 
@@ -75,7 +75,7 @@ function Base.show(io::IO, ::MIME"text/html", A::HTML5anim)
     print(io, A.anim)
 end
 
-function plot(M::Geometry{T,X,W,<:Any,<:Any,Discretization}, ts::AbstractVector{T}, U::AbstractMatrix{T};
+function plot(M::Geometry{T,X,W,<:Any,Discretization}, ts::AbstractVector{T}, U::AbstractMatrix{T};
         frame_time::Real = max(0.001, minimum(diff(ts))),
         embed_limit=200.0,
         printer=(animation)->nothing