diff --git a/NEWS.md b/NEWS.md index 09d3ff686..20e9ca313 100644 --- a/NEWS.md +++ b/NEWS.md @@ -5,7 +5,15 @@ 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.11] - 2026-04-09 + +### Fixed +- Fixed a bug in the outer constructor of UnstructuredGrid for the particular case in which the value of the kw-arg has_affine_map is nothing. Since PR[#1273](https://github.com/gridap/Gridap.jl/pull/1273) + +## [0.19.10] - 2026-04-02 + +### Added +- Changed inner constructor signature of UnstructuredDiscreteModel with flexibility in mind (as per required by GridapGeosciences.jl). Added a new outer constructor as the old inner contructor with backward compatibility in mind. Since PR[#1264](https://github.com/gridap/Gridap.jl/pull/1264) ## [0.20.6] - 2026-05-05 diff --git a/docs/src/dev-notes/constraints.md b/docs/src/dev-notes/constraints.md new file mode 100644 index 000000000..d472762e3 --- /dev/null +++ b/docs/src/dev-notes/constraints.md @@ -0,0 +1,31 @@ + +# FESpaceWithLinearConstraints + +This is a space built from an initial FESpace plus a set of linear constraints. +The initial FESpace may also be a `FESpaceWithLinearConstraints`. + +## Assumptions + +- A constrained dof depends only on master dofs, i.e., only one level of constraints. This can be easily relaxed in the future. +- Free dofs can be constrained both by free and Dirichlet dofs. Dirichlet dofs can be constrained by Dirichlet dofs only. +- All constrained dofs belong to the original space. Master dofs will generally also belong to the original space, but they may be external. This allows us to have, e.g. + - External master Dirichlet dofs, to implement affine constraints. + - In distributed context, masters may belong to other processors and therefore NOT be local (and thus not belong to the original local space). + +## Notation + +We will setup some notation to differentiate between the different dof numberings that appear throughout the implementation: + +We have two ways of numbering dofs: + +- We will refer to `dof` in non-capital letters to refer to signed dof indices. + The sign will indicate if the dof is free (positive) or Dirichlet (negative). +- We will refer to `DOF` in capital letters to refer to unsigned dof indices. In this case, Dirichlet dofs are also represented with a positive id "pas the end" of free dof ids. I.e., we have `DOF = dof` if `dof > 0` and `DOF = -dof + n_fdofs` if `dof < 0`, where `n_fdofs` is the number of free dofs in the original space. + +We will differentiate between three different sets of dofs: + +- First, DoFs in the original space will be denoted by `dof`/`DOF`. They can either be free (`fdof`) or Dirichlet (`ddof`). +- Master DoFs are the non-constrained degrees of freedom. They will be denoted by `mdof`/`mDOF`. They can either be free (`fmdof`) or Dirichlet (`dmdof`). +- Slave DoFs are the constrained degrees of freedom. They will be denoted by `sdof`/`sDOF`. + +Moreover, we will sometimes need to refer to local dofs in a cell, for which we will use `ldof` for the original space and `lmdof` for the constrained space. diff --git a/src/FESpaces/ConstraintHandlers.jl b/src/FESpaces/ConstraintHandlers.jl new file mode 100644 index 000000000..d41c80cd3 --- /dev/null +++ b/src/FESpaces/ConstraintHandlers.jl @@ -0,0 +1,486 @@ + +# --------------------------------------------------------------------------- +# ConstraintLine +# --------------------------------------------------------------------------- + +# One constrained (slave) DoF and its linear relation to master DoFs. +struct ConstraintLine{T<:Number} + dof::Int + masters::Vector{Int} + coeffs::Vector{T} + offset::T +end + +function Base.show(io::IO, l::ConstraintLine) + terms = join(["$(c)*x[$(d)]" for (d, c) in zip(l.masters, l.coeffs)], " + ") + rhs = isempty(terms) ? "$(l.offset)" : + iszero(l.offset) ? terms : "$terms + $(l.offset)" + print(io, "x[$(l.dof)] = $rhs") +end + +# --------------------------------------------------------------------------- +# ConstraintHandler +# --------------------------------------------------------------------------- + +# DoF numbering: Like FESpaceWithLinearConstraints, we use unsigned DoF indices +# in ConstraintHandler. The function `add_constraint!` accepts signed DoF indices for +# user convenience, which get mapped into DOFs (unsigned). +# Then: +# - `n_dofs`: Total number of DOFs in the space (free + dirichlet) +# - `n_fdofs`: Number of free DOFs (n_dofs >= n_fdofs) +mutable struct ConstraintHandler{T<:Number} + n_dofs::Int + n_fdofs::Int + constraints::Vector{ConstraintLine{T}} + dof_to_constraint::Vector{Int} # 0 if unconstrained + closed::Bool +end + +is_constrained(ch::ConstraintHandler, dof::Int) = + !iszero(ch.dof_to_constraint[_dof_to_DOF(dof, ch.n_fdofs)]) + +num_constrained_dofs(ch::ConstraintHandler) = length(ch.constraints) +num_free_dofs(ch::ConstraintHandler) = ch.n_dofs - num_constrained_dofs(ch) + +function Base.show(io::IO, ch::ConstraintHandler{T}) where T + status = ch.closed ? "closed" : "open" + fdof_str = ch.n_fdofs < ch.n_dofs ? ", $(ch.n_fdofs) free" : "" + print(io, "ConstraintHandler{$T}($(ch.n_dofs) dofs$fdof_str, " * + "$(num_constrained_dofs(ch)) constrained, $status)") +end + +""" + ConstraintHandler(n_dofs, T=Float64) + +Create an empty constraint handler for a space with `n_dofs` degrees of freedom. +All DOFs are treated as free (unsigned indices, no sign conversion). +`T` is the coefficient type. +""" +ConstraintHandler(n_dofs::Int, ::Type{T}=Float64) where T = + ConstraintHandler{T}(n_dofs, n_dofs, ConstraintLine{T}[], zeros(Int, n_dofs), false) + +""" + ConstraintHandler(V::SingleFieldFESpace, T=Float64) + +Create an empty constraint handler for the FE space `V`. DOF indices passed to +`add_constraint!` may use Gridap's signed convention (positive = free, negative = +Dirichlet); they are converted to unsigned indices internally. +`T` is the coefficient type. +""" +function ConstraintHandler(V::SingleFieldFESpace, ::Type{T}=Float64) where T + n_fdofs = num_free_dofs(V) + n_dofs = n_fdofs + num_dirichlet_dofs(V) + ConstraintHandler{T}(n_dofs, n_fdofs, ConstraintLine{T}[], zeros(Int, n_dofs), false) +end + +""" + ConstraintHandler(sDOF_to_dof, sDOF_to_dofs, sDOF_to_coeffs, V::FESpace, T=Float64) + +Reconstruct an open `ConstraintHandler` from a Gridap Constructor-2 constraint +table triple defined on the `SingleFieldFESpace` `V`. This is the inverse of +`constraint_tables`. + +The returned handler is open; call `close!` to chain-resolve and finalise it. +""" +function ConstraintHandler( + sDOF_to_dof, sDOF_to_dofs::Table, sDOF_to_coeffs::Table, + V::SingleFieldFESpace, ::Type{T}=Float64 +) where T + ch = ConstraintHandler(V, T) + for (i, dof) in enumerate(sDOF_to_dof) + add_constraint!(ch, dof, dataview(sDOF_to_dofs, i), dataview(sDOF_to_coeffs, i)) + end + return ch +end + +# --------------------------------------------------------------------------- +# Constraint addition +# --------------------------------------------------------------------------- + +""" + add_constraint!(ch, line; on_conflict=nothing) + +Register a pre-built `ConstraintLine`. The other `add_constraint!` methods +delegate to this one. + +When `on_conflict` is `nothing` (default), constraining an already-constrained +DoF raises an error. Otherwise `on_conflict(existing, incoming) → ConstraintLine` +is called and its return value replaces the stored line. +""" +function add_constraint!( + ch::ConstraintHandler{T}, line::ConstraintLine{T}; + on_conflict = nothing +) where T + @check !ch.closed "Cannot add constraints after close!()" + @check 1 <= line.dof <= ch.n_dofs "DoF $(line.dof) out of range [1, $(ch.n_dofs)]" + @check length(line.masters) == length(line.coeffs) "masters and coeffs must have the same length" + for master in line.masters + @check 1 <= master <= ch.n_dofs "Master DoF $master out of range [1, $(ch.n_dofs)]" + @check master != line.dof "Self-referential constraint on DoF $(line.dof)" + end + + dof = line.dof + if !is_constrained(ch, dof) + push!(ch.constraints, line) + ch.dof_to_constraint[dof] = length(ch.constraints) + elseif !isnothing(on_conflict) + idx = ch.dof_to_constraint[dof] + ch.constraints[idx] = on_conflict(ch.constraints[idx], line) + else + @notimplemented """ + Constraint conflict on DoF $dof: + - existing: $(ch.constraints[ch.dof_to_constraint[dof]]) + - incoming: $line + To resolve conflicts, provide an `on_conflict` function to `add_constraint!` + """ + end + return ch +end + +""" + add_constraint!(ch, dof, masters, coeffs, offset=zero(T); on_conflict=nothing) + +Register the linear constraint x[dof] = Σ coeffs[i] * x[masters[i]] + offset. + +Both `masters` and `coeffs` are copied internally, so the caller can safely reuse +or modify them after the call. +""" +function add_constraint!( + ch::ConstraintHandler{T}, dof::Int, + masters::AbstractVector{<:Integer}, + coeffs::AbstractVector{<:Number}, + offset::Number = zero(T); + kwargs... +) where T + @check !ch.closed "Cannot add constraints to a closed handler" + @check length(masters) == length(coeffs) + DOF = _dof_to_DOF(dof, ch.n_fdofs) + _masters = _dof_to_DOF!(collect(Int, masters), ch.n_fdofs) + line = ConstraintLine{T}(DOF, _masters, collect(T, coeffs), T(offset)) + add_constraint!(ch, line; kwargs...) + return ch +end + +function add_constraint!( + ch::ConstraintHandler{T}, + dofs::AbstractVector{<:Integer}, + masters::AbstractVector{<:Integer}, + coeffs::AbstractMatrix{<:Number}, + offsets::AbstractVector{<:Number} = zeros(T, size(coeffs, 1)); + kwargs... +) where T + n_s, n_m = length(dofs), length(masters) + @check size(coeffs, 1) == n_s "coeffs has $(size(coeffs,1)) rows but dofs has $n_s entries" + @check size(coeffs, 2) == n_m "coeffs has $(size(coeffs,2)) cols but masters has $(length(masters)) entries" + @check length(offsets) == n_s "offsets has $(length(offsets)) entries but dofs has $n_s entries" + @inbounds for j in 1:n_s + add_constraint!( + ch, dofs[j], masters, view(coeffs, j, :), offsets[j]; kwargs... + ) + end + return ch +end + +""" + add_constraint!(ch, dof, offset) + +Register a Dirichlet constraint x[dof] = offset (no master DoFs). +""" +function add_constraint!(ch::ConstraintHandler{T}, dof::Int, offset::Number; kwargs...) where T + DOF = _dof_to_DOF(dof, ch.n_fdofs) + add_constraint!(ch, ConstraintLine{T}(DOF, Int[], T[], T(offset)); kwargs...) +end + +""" + set_offset!(ch, dof, value) + +Update the offset of an already-registered constraint on `dof`. +""" +function set_offset!(ch::ConstraintHandler{T}, dof::Int, value::Number) where T + @check !ch.closed "Cannot modify constraints after close!()" + @check is_constrained(ch, dof) "DoF $dof has no constraint to update" + idx = ch.dof_to_constraint[_dof_to_DOF(dof, ch.n_fdofs)] + line = ch.constraints[idx] + ch.constraints[idx] = ConstraintLine{T}(line.dof, line.masters, line.coeffs, T(value)) + return ch +end + +""" + merge!(a, b; on_conflict=nothing) + +Absorb all constraints from `b` into `a`. Both handlers must be open. +Conflict resolution is optional and controlled by the user-defined function `on_conflict`, +which is called when the same DoF is constrained in both handlers. +Its signature should be + + `on_conflict(existing::ConstraintLine, incoming::ConstraintLine) → ConstraintLine` + +Example resolvers: + - keep existing: `(e, _) -> e` + - keep incoming: `(_, i) -> i` + - average: `(e, i) -> ConstraintLine(e.dof, vcat(e.masters, i.masters), vcat(e.coeffs, i.coeffs) .* 0.5, 0.0)` + +Chain resolution and cycle detection are deferred to `close!(a)`. +""" +function Base.merge!(a::ConstraintHandler{T}, b::ConstraintHandler{T}; on_conflict=nothing) where T + @check !a.closed && !b.closed "Both handlers must be open to merge" + @check a.n_dofs == b.n_dofs "Handlers have incompatible sizes: $(a.n_dofs) vs $(b.n_dofs)" + for line in b.constraints + add_constraint!(a, line; on_conflict) + end + return a +end + +""" + close!(ch; tol=1e-13) + +Finalize the constraint handler: +1. Sort constraints by DoF index; rebuild dof_to_constraint. +2. Detect dependency cycles via topological sort (Kahn's algorithm). +3. Resolve constraint chains by iterative substitution. +4. Merge duplicate master entries, strip near-zero coefficients. +""" +function close!(ch::ConstraintHandler{T}; tol::Float64=1e-13) where T + ch.closed && return ch + + # Step 1: sort, rebuild dof_to_constraint. + sort!(ch.constraints, by=l->l.dof) + fill!(ch.dof_to_constraint, 0) + for (i, line) in enumerate(ch.constraints) + ch.dof_to_constraint[line.dof] = i + end + + # Step 2: cycle detection. + _detect_cycles(ch) + + # Step 3: chain resolution (multi-pass substitution). + # After cycle detection we know the dependency graph is a DAG, so this + # terminates in at most n_slaves passes. + masters, coeffs = Int[], T[] # reusable buffers + constrained = Base.Fix1(is_constrained, ch) + for _ in 1:length(ch.constraints) + changed = false + for i in eachindex(ch.constraints) + any(constrained, ch.constraints[i].masters) || continue + + line = ch.constraints[i] + empty!(masters); empty!(coeffs) + new_offset = line.offset + for (master, coeff) in zip(line.masters, line.coeffs) + if constrained(master) + sub = ch.constraints[ch.dof_to_constraint[master]] + append!(masters, sub.masters) + append!(coeffs, coeff .* sub.coeffs) + new_offset += coeff * sub.offset + else + push!(masters, master); push!(coeffs, coeff) + end + end + resize!(line.masters, length(masters)); copyto!(line.masters, masters) + resize!(line.coeffs, length(coeffs)); copyto!(line.coeffs, coeffs) + ch.constraints[i] = ConstraintLine{T}(line.dof, line.masters, line.coeffs, new_offset) + changed = true + end + !changed && break + end + + # Step 4: merge duplicate masters, strip near-zeros. + # Masters and coeffs are modified in-place — no new ConstraintLine needed. + z = zero(T) + m_to_c = OrderedDict{Int,T}() + for line in ch.constraints + empty!(m_to_c) + for (master, coeff) in zip(line.masters, line.coeffs) + abs(coeff) < tol && continue + m_to_c[master] = get(m_to_c, master, z) + coeff + end + n_new = length(m_to_c) + resize!(line.masters, n_new) + resize!(line.coeffs, n_new) + for (i, (master, coeff)) in enumerate(m_to_c) + line.masters[i] = master + line.coeffs[i] = coeff + end + end + + ch.closed = true + return ch +end + +# Cycle detection via Kahn's topological sort on the slave→slave subgraph. +# Edge direction: t → s means "t must be resolved before s" (s depends on t). +function _detect_cycles(ch::ConstraintHandler) + in_degree = zeros(Int, ch.n_dofs) + adj = Dict{Int,Vector{Int}}() # only populated for slave→slave edges + + for line in ch.constraints + for master in line.masters + if is_constrained(ch, master) + push!(get!(adj, master, Int[]), line.dof) # master must precede line.dof + in_degree[line.dof] += 1 + end + end + end + + q = Queue{Int}() + for line in ch.constraints + iszero(in_degree[line.dof]) && enqueue!(q, line.dof) + end + + processed = 0 + while !isempty(q) + dof = dequeue!(q) + processed += 1 + for t in get(adj, dof, Int[]) + in_degree[t] -= 1 + iszero(in_degree[t]) && enqueue!(q, t) + end + end + + if processed < num_constrained_dofs(ch) + cycle_dofs = sort!([line.dof for line in ch.constraints if in_degree[line.dof] > 0]) + error("Constraint cycle detected among DoFs: $cycle_dofs\n" * + "Cycles are not resolvable by substitution. Check your constraint sources.") + end +end + +""" + get_matrix(ch) → (C, g) + +Return the sparse prolongation matrix `C` (n × n_free) and offset vector `g` +(length n) such that the full solution vector satisfies `x = C * u + g`, +where `u` is the reduced (free-dof) solution. + +Columns of `C` correspond to free DoFs in natural index order. +Free-DoF rows are identity rows; slave-DoF rows carry the constraint coefficients. +""" +function Algebra.get_matrix(ch::ConstraintHandler{T}) where T + @check ch.closed "The ConstraintHandler must be closed before building the constraint matrix." + n = ch.n_dofs + n_free = num_free_dofs(ch) + + # Map each free DoF to its column index in C. + col = 0 + n_entries = 0 + dof_to_col = zeros(Int, n) + for i in 1:n + if !is_constrained(ch, i) + col += 1 + dof_to_col[i] = col + n_entries += 1 + else + line = ch.constraints[ch.dof_to_constraint[i]] + n_entries += length(line.masters) + end + end + + I = Vector{Int}(undef,n_entries) + J = Vector{Int}(undef,n_entries) + V = Vector{T}(undef,n_entries) + g = zeros(T, n) + + k = 1 + for i in 1:n + if !is_constrained(ch, i) + I[k] = i + J[k] = dof_to_col[i] + V[k] = one(T) + k += 1 + else + line = ch.constraints[ch.dof_to_constraint[i]] + g[i] = line.offset + for (master, coeff) in zip(line.masters, line.coeffs) + I[k] = i + J[k] = dof_to_col[master] + V[k] = coeff + k += 1 + end + end + end + + C = sparse(I, J, V, n, n_free) + return C, g +end + +""" + apply_constraints(A, b, ch) → (Ã, b̃) + +Return the reduced system `Ã = C'AC` and `b̃ = C'(b - Ag)` obtained by +eliminating all constrained DoFs via the prolongation matrix `C` and offset +vector `g` from `constraint_matrix(ch)`. + +Works with any matrix type; sparsity is preserved when `A` is sparse. +After solving `Ã * u = b̃`, recover the full solution with `x = C * u + g`. +""" +function apply_constraints(A, b::AbstractVector, ch::ConstraintHandler) + @check ch.closed "Must call close!(ch) before condense" + C, g = get_matrix(ch) + Ã = C' * A * C + b̃ = C' * (b - A * g) + return Ã, b̃ +end + +""" + apply_constraints!(x, ch) + +Overwrite `x[dof]` for every constrained DoF using the stored constraints: + + x[dof] = offset + Σ_j coefficient_j * x[master_j] + +After `close!()`, all masters are free, so the iteration order is immaterial. +""" +function apply_constraints!(x::AbstractVector, ch::ConstraintHandler) + @check ch.closed "Must call close!() before apply_constraints!()" + @check length(x) == ch.n_dofs "Vector length $(length(x)) does not match n_dofs $(ch.n_dofs)" + for line in ch.constraints + val = line.offset + for (master, coeff) in zip(line.masters, line.coeffs) + val += coeff * x[master] + end + x[line.dof] = val + end + return x +end + +# --------------------------------------------------------------------------- +# FESpaceWithLinearConstraints from a ConstraintHandler +# --------------------------------------------------------------------------- + +function constraint_tables(V::SingleFieldFESpace, ch::ConstraintHandler{T}) where T + @check ch.closed "Must call close!(ch) before building constraint tables" + + n_fdofs = num_free_dofs(V) + n_s = length(ch.constraints) + + ptrs = Vector{Int}(undef, n_s + 1) + ptrs[1] = 1 + for (i, line) in enumerate(ch.constraints) + ptrs[i+1] = ptrs[i] + length(line.masters) + end + nnz = ptrs[end] - 1 + master_dofs = Vector{Int}(undef, nnz) + coefficients = Vector{T}(undef, nnz) + for (i, line) in enumerate(ch.constraints) + for j in eachindex(line.masters) + master_dofs[ ptrs[i] + j - 1] = _DOF_to_dof(line.masters[j], n_fdofs) + coefficients[ptrs[i] + j - 1] = line.coeffs[j] + end + end + + sDOF_to_dof = [_DOF_to_dof(line.dof, n_fdofs) for line in ch.constraints] + sDOF_to_dofs = Table(master_dofs, ptrs) + sDOF_to_coeffs = Table(coefficients, copy(ptrs)) + sDOF_to_offsets = T[line.offset for line in ch.constraints] + + return sDOF_to_dof, sDOF_to_dofs, sDOF_to_coeffs, sDOF_to_offsets +end + +""" + FESpaceWithLinearConstraints(V::SingleFieldFESpace, ch::ConstraintHandler) + +Build a `FESpaceWithLinearConstraints` from `V` and a closed `ConstraintHandler`. +""" +function FESpaceWithLinearConstraints(V::SingleFieldFESpace, ch::ConstraintHandler) + sDOF_to_dof, sDOF_to_dofs, sDOF_to_coeffs, sDOF_to_offsets = constraint_tables(V, ch) + return FESpaceWithLinearConstraints(sDOF_to_dof, sDOF_to_dofs, sDOF_to_coeffs, V; sDOF_to_offsets) +end diff --git a/src/FESpaces/FESpaces.jl b/src/FESpaces/FESpaces.jl index 7c44abb06..8a80dfb59 100644 --- a/src/FESpaces/FESpaces.jl +++ b/src/FESpaces/FESpaces.jl @@ -11,6 +11,7 @@ using SparseArrays using LinearAlgebra using StaticArrays using ForwardDiff +using DataStructures using Gridap.Helpers using Gridap.Algebra @@ -209,6 +210,18 @@ export CLagrangianFESpace export DirichletFESpace export FESpaceWithLinearConstraints +export ConstraintLine +export ConstraintHandler +export add_constraint! +export set_offset! +export is_constrained +export num_constrained_dofs +export close! +export merge! +export apply_constraints! +export apply_constraints +export constraint_tables + export FiniteElements export DiscreteModelWithFEMap @@ -266,6 +279,8 @@ include("DirichletFESpaces.jl") include("FESpacesWithLinearConstraints.jl") +include("ConstraintHandlers.jl") + include("DiscreteModelWithFEMaps.jl") include("ConstantFESpaces.jl") diff --git a/src/FESpaces/FESpacesWithLinearConstraints.jl b/src/FESpaces/FESpacesWithLinearConstraints.jl index 1a6be1e79..bc0b09d50 100644 --- a/src/FESpaces/FESpacesWithLinearConstraints.jl +++ b/src/FESpaces/FESpacesWithLinearConstraints.jl @@ -52,519 +52,849 @@ # # - coeff: A coefficient in a linear constraint -# TODO we can optimize memory by only storing info about slave DOFs -# The fields of this struct are private """ struct FESpaceWithLinearConstraints{S<:SingleFieldFESpace} <: SingleFieldFESpace """ struct FESpaceWithLinearConstraints{S<:SingleFieldFESpace} <: SingleFieldFESpace space::S - n_fdofs::Int - n_fmdofs::Int - mDOF_to_DOF::Vector - DOF_to_mDOFs::Table - DOF_to_coeffs::Table - cell_to_lmdof_to_mdof::Table - cell_to_ldof_to_dof::Table + mDOF_to_dof::Vector + sDOF_to_dof::Vector + sDOF_to_mdofs::Table + sDOF_to_coeffs::Table + cell_to_mdofs::Table + dmdof_to_offsets::Vector + n_free::Int end # Constructors -# This is the public constructor +function FESpaceWithLinearConstraints( + space::SingleFieldFESpace, + mDOF_to_dof::AbstractVector{<:Integer}, + sDOF_to_dof::AbstractVector{<:Integer}, + sDOF_to_mdofs::Table, + sDOF_to_coeffs::Table, + n_fmdofs::Int = _count_free_mdofs(mDOF_to_dof,sDOF_to_mdofs), + dmdof_to_offsets::AbstractVector = zeros(Float64, length(mDOF_to_dof) - n_fmdofs) +) + cell_to_mdofs = _generate_cell_to_mdofs( + space, mDOF_to_dof, sDOF_to_dof, sDOF_to_mdofs, n_fmdofs + ) + @check _check_constraints(space, mDOF_to_dof, sDOF_to_dof, sDOF_to_mdofs, n_fmdofs) + return FESpaceWithLinearConstraints( + space, mDOF_to_dof, sDOF_to_dof, sDOF_to_mdofs, + sDOF_to_coeffs, cell_to_mdofs, dmdof_to_offsets, n_fmdofs + ) +end + """ FESpaceWithLinearConstraints( sDOF_to_dof::AbstractVector{<:Integer}, sDOF_to_dofs::Table, sDOF_to_coeffs::Table, - space::SingleFieldFESpace) + space::SingleFieldFESpace + ) """ function FESpaceWithLinearConstraints( sDOF_to_dof::AbstractVector{<:Integer}, sDOF_to_dofs::Table, sDOF_to_coeffs::Table, - space::SingleFieldFESpace) + space::SingleFieldFESpace; + sDOF_to_offsets=nothing +) + mDOF_to_dof, sDOF_to_mdofs, n_fmdofs, n_dmdofs = + _find_master_dofs(sDOF_to_dof, sDOF_to_dofs, space) + if !isnothing(sDOF_to_offsets) + mDOF_to_dof, sDOF_to_mdofs, sDOF_to_coeffs, dmdof_to_offsets = + _attach_offsets(mDOF_to_dof, sDOF_to_mdofs, sDOF_to_coeffs, sDOF_to_offsets, n_dmdofs) + else + dmdof_to_offsets = eltype(sDOF_to_coeffs.data)[] + end + return FESpaceWithLinearConstraints( + space, mDOF_to_dof, sDOF_to_dof, sDOF_to_mdofs, sDOF_to_coeffs, n_fmdofs, dmdof_to_offsets + ) +end - n_fdofs = num_free_dofs(space) - n_ddofs = num_dirichlet_dofs(space) - n_DOFs = n_fdofs + n_ddofs +function FESpaceWithLinearConstraints( + DOF_to_dofs::Table, DOF_to_coeffs::Table, space::SingleFieldFESpace, + DOF_to_offsets=nothing +) + sDOF_to_dof, sDOF_to_dofs, sDOF_to_coeffs, sDOF_to_offsets = _filter_constraints( + DOF_to_dofs, DOF_to_coeffs, DOF_to_offsets, space + ) + return FESpaceWithLinearConstraints( + sDOF_to_dof, sDOF_to_dofs, sDOF_to_coeffs, space; sDOF_to_offsets + ) +end - DOF_to_DOFs, DOF_to_coeffs = _prepare_DOF_to_DOFs( - sDOF_to_dof, sDOF_to_dofs, sDOF_to_coeffs, n_fdofs, n_DOFs) +function FESpaceWithLinearConstraints( + fdof_to_dofs::Table, + fdof_to_coeffs::Table, + ddof_to_dofs::Table, + ddof_to_coeffs::Table, + space::SingleFieldFESpace +) + DOF_to_dofs, DOF_to_coeffs = _merge_free_and_diri_constraints( + fdof_to_dofs, fdof_to_coeffs, ddof_to_dofs, ddof_to_coeffs + ) + return FESpaceWithLinearConstraints(DOF_to_dofs, DOF_to_coeffs, space) +end - FESpaceWithLinearConstraints!(DOF_to_DOFs,DOF_to_coeffs,space) +_dof_to_DOF(dof, n_fdofs) = ifelse(iszero(dof), 0, ifelse(dof > 0, dof, n_fdofs - dof)) +_DOF_to_dof(DOF, n_fdofs) = ifelse(DOF > n_fdofs, -(DOF - n_fdofs), DOF) +function _dof_to_DOF!(dofs::Vector{Int}, n_fdofs::Int) + @inbounds for i in eachindex(dofs) + dofs[i] = _dof_to_DOF(dofs[i], n_fdofs) + end + return dofs end -function _prepare_DOF_to_DOFs( - sDOF_to_dof, sDOF_to_dofs, sDOF_to_coeffs, n_fdofs, n_DOFs) - - Tp = eltype(sDOF_to_dofs.ptrs) - Td = eltype(sDOF_to_dofs.data) - Tc = eltype(sDOF_to_coeffs.data) +# Implementation of FESpace interface - DOF_to_DOFs_ptrs = ones(Tp,n_DOFs+1) +get_triangulation(f::FESpaceWithLinearConstraints) = get_triangulation(f.space) - n_sDOFs = length(sDOF_to_dof) +get_free_dof_ids(f::FESpaceWithLinearConstraints) = Base.OneTo(f.n_free) +get_vector_type(f::FESpaceWithLinearConstraints) = get_vector_type(f.space) - for sDOF in 1:n_sDOFs - a = sDOF_to_dofs.ptrs[sDOF] - b = sDOF_to_dofs.ptrs[sDOF+1] - dof = sDOF_to_dof[sDOF] - DOF = _dof_to_DOF(dof,n_fdofs) - DOF_to_DOFs_ptrs[DOF+1] = b-a - end +get_fe_basis(f::FESpaceWithLinearConstraints) = get_fe_basis(f.space) +get_trial_fe_basis(f::FESpaceWithLinearConstraints) = get_trial_fe_basis(f.space) - length_to_ptrs!(DOF_to_DOFs_ptrs) - ndata = DOF_to_DOFs_ptrs[end]-1 - DOF_to_DOFs_data = zeros(Td,ndata) - DOF_to_coeffs_data = ones(Tc,ndata) +CellField(f::FESpaceWithLinearConstraints,cellvals) = CellField(f.space,cellvals) - for DOF in 1:n_DOFs - q = DOF_to_DOFs_ptrs[DOF] - DOF_to_DOFs_data[q] = DOF - end +# Implementation of the SingleFieldFESpace interface - for sDOF in 1:n_sDOFs - dof = sDOF_to_dof[sDOF] - DOF = _dof_to_DOF(dof,n_fdofs) - q = DOF_to_DOFs_ptrs[DOF]-1 - pini = sDOF_to_dofs.ptrs[sDOF] - pend = sDOF_to_dofs.ptrs[sDOF+1]-1 - for (i,p) in enumerate(pini:pend) - _dof = sDOF_to_dofs.data[p] - _DOF = _dof_to_DOF(_dof,n_fdofs) - coeff = sDOF_to_coeffs.data[p] - DOF_to_DOFs_data[q+i] = _DOF - DOF_to_coeffs_data[q+i] = coeff - end - end +get_cell_dof_ids(f::FESpaceWithLinearConstraints) = f.cell_to_mdofs +get_fe_dof_basis(f::FESpaceWithLinearConstraints) = get_fe_dof_basis(f.space) +get_dof_value_type(f::FESpaceWithLinearConstraints) = get_dof_value_type(f.space) - DOF_to_DOFs = Table(DOF_to_DOFs_data,DOF_to_DOFs_ptrs) - DOF_to_coeffs = Table(DOF_to_coeffs_data,DOF_to_DOFs_ptrs) +get_dirichlet_dof_ids(f::FESpaceWithLinearConstraints) = Base.OneTo(length(f.mDOF_to_dof) - f.n_free) +num_dirichlet_tags(f::FESpaceWithLinearConstraints) = num_dirichlet_tags(f.space) - DOF_to_DOFs, DOF_to_coeffs +function get_dirichlet_dof_tag(f::FESpaceWithLinearConstraints) + ddof_to_tag = get_dirichlet_dof_tag(f.space) + dmdof_to_tag = zeros(eltype(ddof_to_tag),num_dirichlet_dofs(f)) + _ddof_to_dmdof_vals!( + dmdof_to_tag,ddof_to_tag,f.mDOF_to_dof + ) + return dmdof_to_tag end -# Private constructors +function get_dirichlet_dof_values(f::FESpaceWithLinearConstraints) + ddof_to_val = get_dirichlet_dof_values(f.space) + dmdof_to_val = zeros(eltype(ddof_to_val), num_dirichlet_dofs(f)) + _ddof_to_dmdof_vals!(dmdof_to_val, ddof_to_val, f.mDOF_to_dof) + n_offsets = length(f.dmdof_to_offsets) + dmdof_to_val[end-n_offsets+1:end] .= f.dmdof_to_offsets + return dmdof_to_val +end -function FESpaceWithLinearConstraints( - fdof_to_dofs::Table, - fdof_to_coeffs::Table, - ddof_to_dofs::Table, - ddof_to_coeffs::Table, - space::SingleFieldFESpace) +function scatter_free_and_dirichlet_values(f::FESpaceWithLinearConstraints,fmdof_to_val,dmdof_to_val) + fdof_to_val = zero_free_values(f.space) + ddof_to_val = zero_dirichlet_values(f.space) + _mdof_to_dof_vals!( + fdof_to_val, ddof_to_val, fmdof_to_val, dmdof_to_val, + f.mDOF_to_dof, f.sDOF_to_dof, f.sDOF_to_mdofs, f.sDOF_to_coeffs + ) + scatter_free_and_dirichlet_values(f.space,fdof_to_val,ddof_to_val) +end - DOF_to_DOFs, DOF_to_coeffs = _merge_free_and_diri_constraints( - fdof_to_dofs, fdof_to_coeffs, ddof_to_dofs, ddof_to_coeffs) - FESpaceWithLinearConstraints!(DOF_to_DOFs,DOF_to_coeffs,space) +function gather_free_and_dirichlet_values(f::FESpaceWithLinearConstraints,cell_to_ldof_to_val) + fdof_to_val, ddof_to_val = gather_free_and_dirichlet_values(f.space,cell_to_ldof_to_val) + fmdof_to_val = zero_free_values(f) + dmdof_to_val = zero_dirichlet_values(f) + _dof_to_mdof_vals!( + fmdof_to_val, dmdof_to_val, fdof_to_val, ddof_to_val, f.mDOF_to_dof + ) + return fmdof_to_val, dmdof_to_val end -function _merge_free_and_diri_constraints(fdof_to_dofs, fdof_to_coeffs, ddof_to_dofs, ddof_to_coeffs) - n_fdofs = length(fdof_to_dofs) - DOF_to_DOFs = append_tables_globally(fdof_to_dofs,ddof_to_dofs) - for i in 1:length(DOF_to_DOFs.data) - dof = DOF_to_DOFs.data[i] - DOF = _dof_to_DOF(dof,n_fdofs) - DOF_to_DOFs.data[i] = DOF - end - DOF_to_coeffs = Table(vcat(fdof_to_coeffs.data,ddof_to_coeffs.data),DOF_to_DOFs.ptrs) - DOF_to_DOFs, DOF_to_coeffs +function gather_free_and_dirichlet_values!(fmdof_to_val,dmdof_to_val,f::FESpaceWithLinearConstraints,cell_to_ldof_to_val) + fdof_to_val, ddof_to_val = gather_free_and_dirichlet_values(f.space,cell_to_ldof_to_val) + _dof_to_mdof_vals!( + fmdof_to_val, dmdof_to_val, fdof_to_val, ddof_to_val, f.mDOF_to_dof + ) + return fmdof_to_val, dmdof_to_val end +function _mdof_to_dof_vals!( + fdof_to_val,ddof_to_val,fmdof_to_val,dmdof_to_val, + mDOF_to_dof, sDOF_to_dof, sDOF_to_mdofs, sDOF_to_coeffs +) + T = eltype(fdof_to_val) + n_mdofs = length(mDOF_to_dof) + n_fmdofs = length(fmdof_to_val) -function FESpaceWithLinearConstraints(DOF_to_DOFs::Table,DOF_to_coeffs::Table,space::SingleFieldFESpace) - FESpaceWithLinearConstraints!(copy(DOF_to_DOFs),copy(DOF_to_coeffs),space::SingleFieldFESpace) -end + # Map free master dofs + for (mdof,mDOF) in enumerate(1:n_fmdofs) + dof = mDOF_to_dof[mDOF] + if dof > 0 + fdof_to_val[dof] = fmdof_to_val[mdof] + elseif dof < 0 + ddof_to_val[-dof] = fmdof_to_val[mdof] + end + end -function FESpaceWithLinearConstraints!(DOF_to_DOFs::Table,DOF_to_coeffs::Table,space::SingleFieldFESpace) + # Map dirichlet master dofs + for (mdof,mDOF) in enumerate((n_fmdofs+1):n_mdofs) + dof = mDOF_to_dof[mDOF] + if dof < 0 + ddof_to_val[-dof] = dmdof_to_val[mdof] + end + end - n_fdofs = num_free_dofs(space) - mDOF_to_DOF, n_fmdofs = _find_master_dofs(DOF_to_DOFs,n_fdofs) - DOF_to_mDOFs = _renumber_constraints!(DOF_to_DOFs,mDOF_to_DOF) - cell_to_ldof_to_dof = Table(get_cell_dof_ids(space)) - cell_to_lmdof_to_mdof = _setup_cell_to_lmdof_to_mdof(cell_to_ldof_to_dof,DOF_to_mDOFs,n_fdofs,n_fmdofs) - - FESpaceWithLinearConstraints( - space, - n_fdofs, - n_fmdofs, - mDOF_to_DOF, - DOF_to_mDOFs, - DOF_to_coeffs, - cell_to_lmdof_to_mdof, - cell_to_ldof_to_dof) + # Map slave dofs + ptrs = sDOF_to_mdofs.ptrs + for (sDOF,dof) in enumerate(sDOF_to_dof) + val = zero(T) + for k in ptrs[sDOF]:(ptrs[sDOF+1]-1) + mdof = sDOF_to_mdofs.data[k] + coeff = sDOF_to_coeffs.data[k] + if mdof > 0 + val += fmdof_to_val[mdof]*coeff + else + val += dmdof_to_val[-mdof]*coeff + end + end + if dof > 0 + fdof_to_val[dof] = val + else + ddof_to_val[-dof] = val + end + end + return fdof_to_val, ddof_to_val end -function _find_master_dofs(DOF_to_DOFs,n_fdofs) - n_DOFs = length(DOF_to_DOFs) - DOF_to_ismaster = fill(false,n_DOFs) - for DOF in 1:n_DOFs - pini = DOF_to_DOFs.ptrs[DOF] - pend = DOF_to_DOFs.ptrs[DOF+1]-1 - for p in pini:pend - _DOF = DOF_to_DOFs.data[p] - @assert (DOF_to_DOFs.ptrs[_DOF+1]-DOF_to_DOFs.ptrs[_DOF]) == 1 "Rcursive constraints not allowed now" - @assert DOF_to_DOFs.data[DOF_to_DOFs.ptrs[_DOF]] == _DOF "Rcursive constraints not allowed now" - DOF_to_ismaster[_DOF] = true +function _dof_to_mdof_vals!( + fmdof_to_val, dmdof_to_val, fdof_to_val, ddof_to_val, mDOF_to_dof +) + n_mdofs = length(mDOF_to_dof) + n_fmdofs = n_mdofs - length(dmdof_to_val) + + # Map free master dofs + for (mdof,mDOF) in enumerate(1:n_fmdofs) + dof = mDOF_to_dof[mDOF] + if dof > 0 + fmdof_to_val[mdof] = fdof_to_val[dof] + elseif dof < 0 + fmdof_to_val[mdof] = ddof_to_val[-dof] end end - n_fmdofs = 0 - for DOF in 1:n_fdofs - if DOF_to_ismaster[DOF] - n_fmdofs += 1 + + # Map dirichlet master dofs + for (mdof,mDOF) in enumerate((n_fmdofs+1):n_mdofs) + dof = mDOF_to_dof[mDOF] + if dof < 0 + dmdof_to_val[mdof] = ddof_to_val[-dof] end end - mDOF_to_DOF = findall(DOF_to_ismaster) - mDOF_to_DOF, n_fmdofs + + return fmdof_to_val, dmdof_to_val end -function _renumber_constraints!(DOF_to_DOFs,mDOF_to_DOF) - DOF_to_mDOF = zeros(eltype(DOF_to_DOFs.data),length(DOF_to_DOFs)) - DOF_to_mDOF[mDOF_to_DOF] .= 1:length(mDOF_to_DOF) - for i in 1:length(DOF_to_DOFs.data) - DOF = DOF_to_DOFs.data[i] - mDOF = DOF_to_mDOF[DOF] - DOF_to_DOFs.data[i] = mDOF +function _ddof_to_dmdof_vals!( + dmdof_to_val, ddof_to_val, mDOF_to_dof +) + n_mdofs = length(mDOF_to_dof) + n_fmdofs = n_mdofs - length(dmdof_to_val) + for (mdof, mDOF) in enumerate((n_fmdofs+1):n_mdofs) + dof = mDOF_to_dof[mDOF] + if dof < 0 + dmdof_to_val[mdof] = ddof_to_val[-dof] + end end - DOF_to_DOFs + return dmdof_to_val end -function _setup_cell_to_lmdof_to_mdof(cell_to_ldof_to_dof,DOF_to_mDOFs,n_fdofs,n_fmdofs) +# Constraints - n_cells = length(cell_to_ldof_to_dof) - cell_to_lmdof_to_mdof_ptrs = zeros(eltype(cell_to_ldof_to_dof.ptrs),n_cells+1) +ConstraintStyle(::Type{<:FESpaceWithLinearConstraints}) = Constrained() +function get_cell_isconstrained(f::FESpaceWithLinearConstraints) + cell_dofs = get_cell_dof_ids(f.space) + cell_mdofs = get_cell_dof_ids(f) + mDOF_to_dof = f.mDOF_to_dof + + n_cells = length(cell_dofs) + n_mfdofs = num_free_dofs(f) + cell_isconstrained = Vector{Bool}(undef,n_cells) for cell in 1:n_cells - mdofs = Set{Int}() - pini = cell_to_ldof_to_dof.ptrs[cell] - pend = cell_to_ldof_to_dof.ptrs[cell+1]-1 - for p in pini:pend - dof = cell_to_ldof_to_dof.data[p] - DOF = _dof_to_DOF(dof,n_fdofs) - qini = DOF_to_mDOFs.ptrs[DOF] - qend = DOF_to_mDOFs.ptrs[DOF+1]-1 - for q in qini:qend - mDOF = DOF_to_mDOFs.data[q] - mdof = _DOF_to_dof(mDOF,n_fmdofs) - push!(mdofs,mdof) - end + dofs = view(cell_dofs,cell) + mdofs = view(cell_mdofs,cell) + + i = 1 + mask = (length(dofs) != length(mdofs)) + while !mask && (i < length(mdofs)) + mdof = mdofs[i] + mDOF = _dof_to_DOF(mdof,n_mfdofs) + dof = mDOF_to_dof[mDOF] + mask = (dof != dofs[i]) + i += 1 end - cell_to_lmdof_to_mdof_ptrs[cell+1] = length(mdofs) + + cell_isconstrained[cell] = mask end - length_to_ptrs!(cell_to_lmdof_to_mdof_ptrs) - ndata = cell_to_lmdof_to_mdof_ptrs[end]-1 - cell_to_lmdof_to_mdof_data = zeros(eltype(cell_to_ldof_to_dof.data),ndata) + return cell_isconstrained +end - for cell in 1:n_cells - mdofs = Set{Int}() - pini = cell_to_ldof_to_dof.ptrs[cell] - pend = cell_to_ldof_to_dof.ptrs[cell+1]-1 - for p in pini:pend - dof = cell_to_ldof_to_dof.data[p] - DOF = _dof_to_DOF(dof,n_fdofs) - qini = DOF_to_mDOFs.ptrs[DOF] - qend = DOF_to_mDOFs.ptrs[DOF+1]-1 - for q in qini:qend - mDOF = DOF_to_mDOFs.data[q] - mdof = _DOF_to_dof(mDOF,n_fmdofs) - push!(mdofs,mdof) - end - end - o = cell_to_lmdof_to_mdof_ptrs[cell]-1 - for (lmdof, mdof) in enumerate(mdofs) - cell_to_lmdof_to_mdof_data[o+lmdof] = mdof - end - end +function get_cell_constraints(f::FESpaceWithLinearConstraints) + DOF_to_msDOF = generate_DOF_to_msDOF_map(f.space,f.mDOF_to_dof,f.sDOF_to_dof) + k = LinearConstraintsMap( + DOF_to_msDOF, f.sDOF_to_mdofs, f.sDOF_to_coeffs, + length(f.mDOF_to_dof), num_free_dofs(f), num_free_dofs(f.space) + ) + cell_to_mat = get_cell_constraints(f.space) + return lazy_map(k,get_cell_dof_ids(f),get_cell_dof_ids(f.space),cell_to_mat) +end - Table(cell_to_lmdof_to_mdof_data,cell_to_lmdof_to_mdof_ptrs) +struct LinearConstraintsMap{A,B} <: Map + DOF_to_msDOF::Vector{Int} + sDOF_to_mdofs::A + sDOF_to_coeffs::B + n_mdofs::Int + n_fmdofs::Int + n_fdofs::Int end -function _dof_to_DOF(dof,n_fdofs) - if dof > 0 - DOF = dof - else - DOF = n_fdofs - dof - end +function LinearConstraintsMap( + sDOF_to_dof::AbstractVector{<:Integer}, + sDOF_to_dofs::Table, + sDOF_to_coeffs::Table, + space::SingleFieldFESpace +) + mDOF_to_dof, sDOF_to_mdofs, n_fmdofs, _ = + _find_master_dofs(sDOF_to_dof, sDOF_to_dofs, space) + DOF_to_msDOF = generate_DOF_to_msDOF_map(space,mDOF_to_dof,sDOF_to_dof) + return LinearConstraintsMap( + DOF_to_msDOF, sDOF_to_mdofs, sDOF_to_coeffs, + length(mDOF_to_dof), n_fmdofs, num_free_dofs(space) + ) +end + +function return_cache(k::LinearConstraintsMap,mdofs,dofs,mat) + T = typeof(zero(eltype(mat))*zero(eltype(k.sDOF_to_coeffs.data))) + n_lmdofs = length(mdofs) + n_ldofs = length(dofs) + @assert n_ldofs == size(mat,1) + m1 = CachedArray(zeros(T,(n_lmdofs,n_ldofs))) + m2 = CachedArray(zeros(T,(n_lmdofs,size(mat,2)))) + mdof_to_lmdof = Dict{Int,Int}() + return m1, m2, mdof_to_lmdof end -function _DOF_to_dof(DOF,n_fdofs) - if DOF > n_fdofs - dof = -(DOF-n_fdofs) - else - dof = DOF +function evaluate!(cache,k::LinearConstraintsMap,mdofs,dofs,mat) + m1, m2, mdof_to_lmdof = cache + + n_lmdofs = length(mdofs) + n_ldofs = length(dofs) + setsize!(m1,(n_lmdofs,n_ldofs)) + setsize!(m2,(n_lmdofs,size(mat,2))) + a1 = get_array(m1) + a2 = get_array(m2) + fill!(a1,zero(eltype(a1))) + + # Precompute mdof to lmdof map + empty!(mdof_to_lmdof) + for (lmdof,mdof) in enumerate(mdofs) + mdof_to_lmdof[mdof] = lmdof end -end -# Implementation of the SingleFieldFESpace interface + # Precompute DOF to msDOF map + o = one(eltype(a1)) + ptrs = k.sDOF_to_mdofs.ptrs + for (ldof,dof) in enumerate(dofs) + DOF = _dof_to_DOF(dof,k.n_fdofs) + msDOF = k.DOF_to_msDOF[DOF] + if msDOF <= k.n_mdofs # master dof + mDOF = msDOF + mdof = _DOF_to_dof(mDOF,k.n_fmdofs) + lmdof = mdof_to_lmdof[mdof] + a1[lmdof,ldof] = o + else # slave dof + sDOF = msDOF - k.n_mdofs + for i in ptrs[sDOF]:(ptrs[sDOF+1]-1) + mdof = k.sDOF_to_mdofs.data[i] + coeff = k.sDOF_to_coeffs.data[i] + lmdof = mdof_to_lmdof[mdof] + a1[lmdof,ldof] = coeff + end + end + end -function get_cell_dof_ids(f::FESpaceWithLinearConstraints) - f.cell_to_lmdof_to_mdof + #TODO this is not always needed + mul!(a2,a1,mat) + a2 end -function get_fe_dof_basis(f::FESpaceWithLinearConstraints) - get_fe_dof_basis(f.space) -end +# Private methods -get_dof_value_type(f::FESpaceWithLinearConstraints) = get_dof_value_type(f.space) +function _check_constraints( + space, mDOF_to_dof, sDOF_to_dof, sDOF_to_mdofs, n_fmdofs +) + n_fdofs = num_free_dofs(space) + n_dofs = n_fdofs + num_dirichlet_dofs(space) + DOF_is_master = fill(false,n_dofs) + for dof in mDOF_to_dof + iszero(dof) && continue + DOF = _dof_to_DOF(dof,n_fdofs) + DOF_is_master[DOF] = true + end -get_dirichlet_dof_ids(f::FESpaceWithLinearConstraints) = Base.OneTo(length(f.mDOF_to_DOF) - f.n_fmdofs) + mDOF_is_master = fill(true,length(mDOF_to_dof)) + for (mDOF,dof) in enumerate(mDOF_to_dof) + DOF = _dof_to_DOF(dof,n_fdofs) + mDOF_is_master[mDOF] = iszero(DOF) || DOF_is_master[DOF] + end -num_dirichlet_tags(f::FESpaceWithLinearConstraints) = num_dirichlet_tags(f.space) + c = array_cache(sDOF_to_mdofs) + for (sDOF,dof) in enumerate(sDOF_to_dof) + DOF = _dof_to_DOF(dof,n_fdofs) + @check !DOF_is_master[DOF] "A dof cannot be both master and slave." + mdofs = getindex!(c,sDOF_to_mdofs,sDOF) + for mdof in mdofs + mDOF = _dof_to_DOF(mdof,n_fmdofs) + @check mDOF_is_master[mDOF] "Constraints cannot be recursive." + @check (dof > 0) || (mdof < 0) "Dirichlet dofs can only be constrained by Dirichlet dofs." + end + end -function get_dirichlet_dof_tag(f::FESpaceWithLinearConstraints) - ddof_to_tag = get_dirichlet_dof_tag(f.space) - dmdof_to_tag = zeros(eltype(ddof_to_tag),num_dirichlet_dofs(f)) - _setup_ddof_to_tag!( - dmdof_to_tag, - ddof_to_tag, - f.mDOF_to_DOF, - f.n_fdofs, - f.n_fmdofs) - dmdof_to_tag + return true end -function _setup_ddof_to_tag!( - dmdof_to_tag, - ddof_to_tag, - mDOF_to_DOF, - n_fdofs, - n_fmdofs) - - for mDOF in (n_fmdofs+1):length(mDOF_to_DOF) - mdof = _DOF_to_dof(mDOF,n_fmdofs) - @assert mdof < 0 "Dirichlet dofs can only depend on Dirichlet dofs" - dmdof = -mdof - DOF = mDOF_to_DOF[mDOF] - dof = _DOF_to_dof(DOF,n_fdofs) - @assert dof < 0 "Dirichlet dofs can only depend on Dirichlet dofs" - ddof = -dof - dmdof_to_tag[dmdof] = ddof_to_tag[ddof] - end +function _filter_constraints(DOF_to_dofs, DOF_to_coeffs, DOF_to_offsets, space) + n_fdofs = num_free_dofs(space) + isslave(DOF,dofs) = (dofs != [_DOF_to_dof(DOF,n_fdofs)]) + sDOF_to_DOF = findall(map(isslave,eachindex(DOF_to_dofs),DOF_to_dofs)) + sDOF_to_dof = map(Base.Fix2(_DOF_to_dof,n_fdofs),sDOF_to_DOF) + sDOF_to_dofs = DOF_to_dofs[sDOF_to_DOF] + sDOF_to_coeffs = DOF_to_coeffs[sDOF_to_DOF] + sDOF_to_offsets = isnothing(DOF_to_offsets) ? nothing : DOF_to_offsets[sDOF_to_DOF] + return sDOF_to_dof, sDOF_to_dofs, Table(sDOF_to_coeffs.data,sDOF_to_dofs.ptrs), sDOF_to_offsets end -function get_dirichlet_dof_values(f::FESpaceWithLinearConstraints) - ddof_to_tag = get_dirichlet_dof_values(f.space) - dmdof_to_tag = zeros(eltype(ddof_to_tag),num_dirichlet_dofs(f)) - _setup_ddof_to_tag!( - dmdof_to_tag, - ddof_to_tag, - f.mDOF_to_DOF, - f.n_fdofs, - f.n_fmdofs) - dmdof_to_tag +function _merge_free_and_diri_constraints( + fdof_to_dofs, fdof_to_coeffs, ddof_to_dofs, ddof_to_coeffs +) + DOF_to_dofs = append_tables_globally(fdof_to_dofs,ddof_to_dofs) + DOF_to_coeffs = Table( + vcat(fdof_to_coeffs.data,ddof_to_coeffs.data), DOF_to_dofs.ptrs + ) + return DOF_to_dofs, DOF_to_coeffs end -function scatter_free_and_dirichlet_values(f::FESpaceWithLinearConstraints,fmdof_to_val,dmdof_to_val) - fdof_to_val = zero_free_values(f.space) - ddof_to_val = zero_dirichlet_values(f.space) - _setup_dof_to_val!( - fdof_to_val, - ddof_to_val, - fmdof_to_val, - dmdof_to_val, - f.DOF_to_mDOFs, - f.DOF_to_coeffs, - f.n_fdofs, - f.n_fmdofs) - scatter_free_and_dirichlet_values(f.space,fdof_to_val,ddof_to_val) +function _count_free_mdofs(mDOF_to_dof,sDOF_to_mdofs) + a = count(>(0), mDOF_to_dof; init=0) + b = maximum(sDOF_to_mdofs.data) + return max(a, b) end -function _setup_dof_to_val!( - fdof_to_val, - ddof_to_val, - fmdof_to_val, - dmdof_to_val, - DOF_to_mDOFs, - DOF_to_coeffs, - n_fdofs, - n_fmdofs) +function _find_master_dofs(sDOF_to_dof, sDOF_to_dofs, space) + n_fdofs = num_free_dofs(space) + n_dofs = n_fdofs + num_dirichlet_dofs(space) - T = eltype(fdof_to_val) + DOF_ismaster = fill(true, n_dofs) + for dof in sDOF_to_dof + DOF = _dof_to_DOF(dof, n_fdofs) + DOF_ismaster[DOF] = false + end - for DOF in 1:length(DOF_to_mDOFs) - pini = DOF_to_mDOFs.ptrs[DOF] - pend = DOF_to_mDOFs.ptrs[DOF+1]-1 - val = zero(T) - for p in pini:pend - mDOF = DOF_to_mDOFs.data[p] - coeff = DOF_to_coeffs.data[p] - mdof = _DOF_to_dof(mDOF,n_fmdofs) - if mdof > 0 - fmdof = mdof - val += fmdof_to_val[fmdof]*coeff + n_fmdofs = count(view(DOF_ismaster, 1:n_fdofs)) + n_dmdofs = count(view(DOF_ismaster, (n_fdofs+1):n_dofs)) + + kf, kd = 1, 1 + mDOF_to_dof = Vector{Int32}(undef, n_fmdofs + n_dmdofs) + DOF_to_mdof = Vector{Int32}(undef, n_dofs) + for DOF in 1:n_dofs + if DOF_ismaster[DOF] + dof = _DOF_to_dof(DOF, n_fdofs) + if dof > 0 + mDOF_to_dof[kf] = dof + DOF_to_mdof[DOF] = kf + kf += 1 else - dmdof = -mdof - val += dmdof_to_val[dmdof]*coeff + mDOF_to_dof[n_fmdofs + kd] = dof + DOF_to_mdof[DOF] = -kd + kd += 1 end end - dof = _DOF_to_dof(DOF,n_fdofs) - if dof > 0 - fdof = dof - fdof_to_val[fdof] = val - else - ddof = -dof - ddof_to_val[ddof] = val - end end -end + data = zeros(Int32, length(sDOF_to_dofs.data)) + for (i, dof) in enumerate(sDOF_to_dofs.data) + data[i] = DOF_to_mdof[_dof_to_DOF(dof, n_fdofs)] + end -function gather_free_and_dirichlet_values(f::FESpaceWithLinearConstraints,cell_to_ludof_to_val) - fdof_to_val, ddof_to_val = gather_free_and_dirichlet_values(f.space,cell_to_ludof_to_val) - fmdof_to_val = zero_free_values(f) - dmdof_to_val = zero_dirichlet_values(f) - _setup_mdof_to_val!( - fmdof_to_val, - dmdof_to_val, - fdof_to_val, - ddof_to_val, - f.mDOF_to_DOF, - f.n_fdofs, - f.n_fmdofs) - fmdof_to_val, dmdof_to_val -end + return mDOF_to_dof, Table(data, sDOF_to_dofs.ptrs), n_fmdofs, n_dmdofs +end + +# Extends mDOF_to_dof, sDOF_to_mdofs, and sDOF_to_coeffs with fictitious Dirichlet +# masters encoding affine inhomogeneities. One fictitious master per nonzero offset, +# appended with coefficient 1. Returns (mDOF_to_dof, sDOF_to_mdofs, sDOF_to_coeffs, dmdof_to_offsets). +function _attach_offsets(mDOF_to_dof, sDOF_to_mdofs, sDOF_to_coeffs, sDOF_to_offsets, n_dmdofs) + n_slaves = length(sDOF_to_mdofs) + n_offsets = count(!iszero, sDOF_to_offsets) + OT = eltype(sDOF_to_offsets) + CT = eltype(sDOF_to_coeffs.data) + iszero(n_offsets) && return mDOF_to_dof, sDOF_to_mdofs, sDOF_to_coeffs, OT[] + + dmdof_to_offsets = Vector{OT}(undef, n_offsets) + mdofs_data = Vector{Int32}(undef, n_offsets) + coeffs_data = Vector{CT}(undef, n_offsets) + ptrs = zeros(eltype(sDOF_to_mdofs.ptrs), n_slaves + 1) + ptrs[1] = 1 + kf = 0 + for s in 1:n_slaves + b = sDOF_to_offsets[s] + ptrs[s+1] = ptrs[s] + if !iszero(b) + kf += 1 + dmdof_to_offsets[kf] = OT(b) + mdofs_data[kf] = Int32(-(n_dmdofs + kf)) + coeffs_data[kf] = one(CT) + ptrs[s+1] += 1 + end + end -function gather_free_and_dirichlet_values!(fmdof_to_val,dmdof_to_val,f::FESpaceWithLinearConstraints,cell_to_ludof_to_val) - fdof_to_val, ddof_to_val = gather_free_and_dirichlet_values(f.space,cell_to_ludof_to_val) - _setup_mdof_to_val!( - fmdof_to_val, - dmdof_to_val, - fdof_to_val, - ddof_to_val, - f.mDOF_to_DOF, - f.n_fdofs, - f.n_fmdofs) - fmdof_to_val, dmdof_to_val -end + return ( + vcat(mDOF_to_dof, zeros(Int32, n_offsets)), + append_tables_locally(sDOF_to_mdofs, Table(mdofs_data, ptrs)), + append_tables_locally(sDOF_to_coeffs, Table(coeffs_data, ptrs)), + dmdof_to_offsets + ) +end + +function _generate_cell_to_mdofs( + space, mDOF_to_dof, sDOF_to_dof, sDOF_to_mdofs, n_fmdofs +) + DOF_to_msDOF = generate_DOF_to_msDOF_map( + space,mDOF_to_dof,sDOF_to_dof + ) + cell_dofs = get_cell_dof_ids(space) + c1 = array_cache(cell_dofs) + c2 = array_cache(sDOF_to_mdofs) + + n_cells = length(cell_dofs) + n_fdofs = num_free_dofs(space) + n_mDOFs = length(mDOF_to_dof) -function _setup_mdof_to_val!( - fmdof_to_val, - dmdof_to_val, - fdof_to_val, - ddof_to_val, - mDOF_to_DOF, - n_fdofs, - n_fmdofs) - - for mDOF in 1:length(mDOF_to_DOF) - DOF = mDOF_to_DOF[mDOF] - dof = _DOF_to_dof(DOF,n_fdofs) - if dof > 0 - fdof = dof - val = fdof_to_val[fdof] - else - ddof = -dof - val = ddof_to_val[ddof] - end - mdof = _DOF_to_dof(mDOF,n_fmdofs) - if mdof > 0 - fmdof = mdof - fmdof_to_val[fmdof] = val - else - dmdof = -mdof - dmdof_to_val[dmdof] = val + acc = OrderedSet{Int32}() + ptrs = zeros(Int32,n_cells+1) + for cell in 1:n_cells + dofs = getindex!(c1,cell_dofs,cell) + for dof in dofs + DOF = _dof_to_DOF(dof,n_fdofs) + msDOF = DOF_to_msDOF[DOF] + if msDOF <= n_mDOFs # master dof + mdof = _DOF_to_dof(msDOF,n_fmdofs) + push!(acc, mdof) + else # slave dof + sDOF = msDOF - n_mDOFs + mdofs = getindex!(c2,sDOF_to_mdofs,sDOF) + push!(acc, mdofs...) + end end + ptrs[cell+1] += length(acc) + empty!(acc) end + Arrays.length_to_ptrs!(ptrs) -end - -# Implementation of FESpace interface + data = zeros(Int32,ptrs[end]-1) + for cell in 1:n_cells + dofs = getindex!(c1,cell_dofs,cell) + for dof in dofs + DOF = _dof_to_DOF(dof,n_fdofs) + msDOF = DOF_to_msDOF[DOF] + if msDOF <= n_mDOFs # master dof + mdof = _DOF_to_dof(msDOF,n_fmdofs) + push!(acc, mdof) + else # slave dof + sDOF = msDOF - n_mDOFs + mdofs = getindex!(c2,sDOF_to_mdofs,sDOF) + push!(acc, mdofs...) + end + end + data[ptrs[cell]:(ptrs[cell+1]-1)] .= collect(acc) + empty!(acc) + end -function get_triangulation(f::FESpaceWithLinearConstraints) - get_triangulation(f.space) + return Table(data,ptrs) end -get_free_dof_ids(f::FESpaceWithLinearConstraints) = Base.OneTo(f.n_fmdofs) - -function get_vector_type(f::FESpaceWithLinearConstraints) - get_vector_type(f.space) -end +function generate_DOF_to_msDOF_map(space, mDOF_to_dof, sDOF_to_dof) + n_mdofs = length(mDOF_to_dof) + n_fdofs = num_free_dofs(space) + n_dofs = n_fdofs + num_dirichlet_dofs(space) + DOF_to_msDOF = Vector{Int}(undef,n_dofs) -function get_fe_basis(f::FESpaceWithLinearConstraints) - get_fe_basis(f.space) -end + for (mDOF,dof) in enumerate(mDOF_to_dof) + iszero(dof) && continue + DOF = _dof_to_DOF(dof,n_fdofs) + DOF_to_msDOF[DOF] = mDOF + end -function get_trial_fe_basis(f::FESpaceWithLinearConstraints) - get_trial_fe_basis(f.space) -end + for (sDOF,dof) in enumerate(sDOF_to_dof) + DOF = _dof_to_DOF(dof,n_fdofs) + DOF_to_msDOF[DOF] = sDOF + n_mdofs + end -function CellField(f::FESpaceWithLinearConstraints,cellvals) - CellField(f.space,cellvals) + return DOF_to_msDOF end -ConstraintStyle(::Type{<:FESpaceWithLinearConstraints}) = Constrained() +# --------------------------------------------------------------------------- +# Table-level merge and chain-resolution for constraint tables +# --------------------------------------------------------------------------- -function get_cell_isconstrained(f::FESpaceWithLinearConstraints) - #TODO this can be heavily optimized - n = length(get_cell_dof_ids(f)) - Fill(true,n) -end +""" + merge_slave_constraint_tables( + space, s1_dof, s1_dofs, s1_coeffs, s2_dof, s2_dofs, s2_coeffs, + s1_offsets, s2_offsets; on_conflict) + → (sDOF_to_dof, sDOF_to_dofs, sDOF_to_coeffs, sDOF_to_offsets) + +Merge two Constructor-2 constraint table triples defined on the same `space`. +Each triple lists only slave DOFs: `s_dof[i]` is the signed dof of slave `i`, +`s_dofs[i]` its master signed dofs, `s_coeffs[i]` the corresponding coefficients. +`s1_offsets` / `s2_offsets` are optional per-slave inhomogeneity vectors (default zero). + +A conflict arises when the same slave dof appears in both sets: + - `on_conflict=nothing` (default): raise an error. + - otherwise: `on_conflict(dof, dofs1, coeffs1, offset1, dofs2, coeffs2, offset2) → (dofs, coeffs, offset)` + is called and its return value replaces the stored row. +""" +function merge_slave_constraint_tables( + space::SingleFieldFESpace, + s1_dof, s1_dofs::Table, s1_coeffs::Table{T}, + s2_dof, s2_dofs::Table, s2_coeffs::Table{T}, + s1_offsets::AbstractVector = zeros(T, length(s1_dof)), + s2_offsets::AbstractVector = zeros(T, length(s2_dof)); + on_conflict=nothing +) where T + n_fdofs = num_free_dofs(space) + n_dofs = n_fdofs + num_dirichlet_dofs(space) + n_s1 = length(s1_dof) + n_s2 = length(s2_dof) + + s1_DOF = [_dof_to_DOF(dof, n_fdofs) for dof in s1_dof] + DOF_to_s1 = find_inverse_index_map(s1_DOF, n_dofs) + + # conflicted_s1[i1] = i2 (s2 index of the conflicting entry) or 0 (no conflict). + conflicted_s1 = zeros(Int, n_s1) + new_s2_indices = Int[] + sizehint!(new_s2_indices, n_s2) + + for i2 in 1:n_s2 + i1 = DOF_to_s1[_dof_to_DOF(s2_dof[i2], n_fdofs)] + if iszero(i1) + push!(new_s2_indices, i2) + elseif isnothing(on_conflict) + error("Constraint conflict on dof $(s2_dof[i2]):\n" * + " existing: $(collect(dataview(s1_dofs, i1)))\n" * + " incoming: $(collect(dataview(s2_dofs, i2)))") + else + conflicted_s1[i1] = i2 + end + end -function get_cell_constraints(f::FESpaceWithLinearConstraints) + # Fast path: no conflicts. + if all(iszero, conflicted_s1) + return vcat(s1_dof, s2_dof), + append_tables_globally(s1_dofs, s2_dofs), + append_tables_globally(s1_coeffs, s2_coeffs), + vcat(s1_offsets, s2_offsets) + end - k = LinearConstraintsMap( - f.DOF_to_mDOFs, - f.DOF_to_coeffs, - length(f.mDOF_to_DOF), - f.n_fmdofs, - f.n_fdofs) + n_new_s2 = length(new_s2_indices) + n_out = n_s1 + n_new_s2 - cell_to_mat = get_cell_constraints(f.space) - lazy_map(k,f.cell_to_lmdof_to_mdof,f.cell_to_ldof_to_dof,cell_to_mat) + # Pessimistic data size: conflicts get len(s1_row) + len(s2_row) as upper bound. + n_data = 0 + for i in 1:n_s1 + i2 = conflicted_s1[i] + n_data += s1_dofs.ptrs[i+1] - s1_dofs.ptrs[i] + !iszero(i2) && (n_data += s2_dofs.ptrs[i2+1] - s2_dofs.ptrs[i2]) + end + for i2 in new_s2_indices + n_data += s2_dofs.ptrs[i2+1] - s2_dofs.ptrs[i2] + end -end + out_dof = Vector{Int}(undef, n_out) + out_offsets = Vector{T}(undef, n_out) + out_data_d = Vector{Int}(undef, n_data) + out_data_c = Vector{T}(undef, n_data) + ptrs = zeros(Int32, n_out + 1) + copyto!(out_dof, 1, s1_dof, 1, n_s1) + + k = 0 + for i1 in 1:n_s1 + i2 = conflicted_s1[i1] + if iszero(i2) + rin = datarange(s1_dofs, i1) + nl = length(rin) + out_data_d[k+1:k+nl] = view(s1_dofs.data, rin) + out_data_c[k+1:k+nl] = view(s1_coeffs.data, rin) + out_offsets[i1] = s1_offsets[i1] + else + r1 = datarange(s1_dofs, i1); r2 = datarange(s2_dofs, i2) + od, oc, oo = on_conflict(s1_dof[i1], + view(s1_dofs.data, r1), view(s1_coeffs.data, r1), s1_offsets[i1], + view(s2_dofs.data, r2), view(s2_coeffs.data, r2), s2_offsets[i2]) + nl = length(od) + out_data_d[k+1:k+nl] = od + out_data_c[k+1:k+nl] = oc + out_offsets[i1] = oo + end + ptrs[i1+1] = nl + k += nl + end + for (j, i2) in enumerate(new_s2_indices) + ii2 = n_s1 + j + rin = datarange(s2_dofs, i2) + nl = length(rin) + out_dof[ii2] = s2_dof[i2] + out_offsets[ii2] = s2_offsets[i2] + out_data_d[k+1:k+nl] = view(s2_dofs.data, rin) + out_data_c[k+1:k+nl] = view(s2_coeffs.data, rin) + ptrs[ii2+1] = nl + k += nl + end -struct LinearConstraintsMap{A,B} <: Map - DOF_to_mDOFs::A - DOF_to_coeffs::B - n_mDOFs::Int - n_fmdofs::Int - n_fdofs::Int + Arrays.length_to_ptrs!(ptrs) + resize!(out_data_d, k) + resize!(out_data_c, k) + return out_dof, Table(out_data_d, ptrs), Table(out_data_c, ptrs), out_offsets end -function return_cache(k::LinearConstraintsMap,lmdof_to_mdof,ldof_to_dof,mat) - n_lmdofs = length(lmdof_to_mdof) - n_ldofs = length(ldof_to_dof) - n_ludofs = size(mat,2) - @assert n_ldofs == size(mat,1) - m1 = CachedArray(zeros(n_lmdofs,n_ldofs)) - m2 = CachedArray(zeros(n_lmdofs,n_ludofs)) - mDOF_to_lmdof = zeros(Int16,k.n_mDOFs) - m1, m2, mDOF_to_lmdof -end +""" + close_slave_constraint_tables( + space, sDOF_to_dof, sDOF_to_dofs, sDOF_to_coeffs, sDOF_to_offsets; tol) + → (sDOF_to_dof, sDOF_to_dofs, sDOF_to_coeffs, sDOF_to_offsets) -function evaluate!(cache,k::LinearConstraintsMap,lmdof_to_mdof,ldof_to_dof,mat) - m1, m2, mDOF_to_lmdof = cache - n_lmdofs = length(lmdof_to_mdof) - n_ldofs = length(ldof_to_dof) - n_ludofs = size(mat,2) +Chain-resolve a Constructor-2 constraint table triple: substitute any slave master +with its own constraint row until all masters are free. Errors on cycles. +Merges duplicate masters and strips near-zero coefficients (threshold `tol`). +The slave set is invariant under resolution; output rows are in topological order. - setsize!(m1,(n_lmdofs,n_ldofs)) - setsize!(m2,(n_lmdofs,n_ludofs)) - a1 = m1.array - a2 = m2.array - fill!(a1,zero(eltype(a1))) +`sDOF_to_offsets` (optional, default zero) carries inhomogeneities: if slave i +depends on slave j the offset accumulates as `offset_i += coeff_ij * offset_j`. +""" +function close_slave_constraint_tables( + space::SingleFieldFESpace, + sDOF_to_dof, sDOF_to_dofs::Table, sDOF_to_coeffs::Table{T}, + sDOF_to_offsets::AbstractVector = zeros(T, length(sDOF_to_dof)); + tol::T = 1000*eps(T) +) where T + n_fdofs = num_free_dofs(space) + n_dofs = n_fdofs + num_dirichlet_dofs(space) + n_slaves = length(sDOF_to_dof) + + s_DOF = [_dof_to_DOF(dof, n_fdofs) for dof in sDOF_to_dof] + dof_to_sDOF = find_inverse_index_map(s_DOF, n_dofs) + + # Build dependency graph (CSR): edge t → s means slave s depends on slave t. + dep_ptrs = zeros(Int32, n_slaves + 1) + in_degree = zeros(Int, n_slaves) + for s in 1:n_slaves + for d in dataview(sDOF_to_dofs, s) + t = dof_to_sDOF[_dof_to_DOF(d, n_fdofs)] + if t != 0 + dep_ptrs[t+1] += Int32(1) + in_degree[s] += 1 + end + end + end + all(iszero, in_degree) && return sDOF_to_dof, sDOF_to_dofs, sDOF_to_coeffs, copy(sDOF_to_offsets) + + Arrays.length_to_ptrs!(dep_ptrs) + dep_data = Vector{Int}(undef, dep_ptrs[end] - 1) + fill_pos = copy(dep_ptrs) + for s in 1:n_slaves + for d in dataview(sDOF_to_dofs, s) + t = dof_to_sDOF[_dof_to_DOF(d, n_fdofs)] + if t != 0 + dep_data[fill_pos[t]] = s + fill_pos[t] += 1 + end + end + end + dependents = Table(dep_data, dep_ptrs) + + # Kahn's topological sort; cycle detection implicit (processed < n_slaves ⟹ cycle). + head = 0 + queue = findall(iszero, in_degree) + order = Vector{Int}(undef, n_slaves) + order_rev = Vector{Int}(undef, n_slaves) + while !isempty(queue) + t = pop!(queue) + head += 1 + order[head] = t + order_rev[t] = head + for s in dataview(dependents, t) + in_degree[s] -= 1 + iszero(in_degree[s]) && push!(queue, s) + end + end + if head < n_slaves + cycle_dofs = sort!([sDOF_to_dof[s] for s in 1:n_slaves if in_degree[s] > 0]) + error("Constraint cycle detected among dofs: $cycle_dofs.\n" * + "Cycles are not resolvable by substitution.") + end - for (lmdof,mdof) in enumerate(lmdof_to_mdof) - mDOF = _dof_to_DOF(mdof,k.n_fmdofs) - mDOF_to_lmdof[mDOF] = lmdof + # Exact pre-dedup/strip output size: single forward pass in topo order. + n_data = 0 + ub = zeros(Int, n_slaves) + for i in 1:n_slaves + s = order[i] + for d in dataview(sDOF_to_dofs, s) + t = dof_to_sDOF[_dof_to_DOF(d, n_fdofs)] + ub[i] += (t != 0) ? ub[order_rev[t]] : 1 + end + n_data += ub[i] end - for (ldof,dof) in enumerate(ldof_to_dof) - DOF = _dof_to_DOF(dof,k.n_fdofs) - qini = k.DOF_to_mDOFs.ptrs[DOF] - qend = k.DOF_to_mDOFs.ptrs[DOF+1]-1 - for q in qini:qend - mDOF = k.DOF_to_mDOFs.data[q] - coeff = k.DOF_to_coeffs.data[q] - lmdof = mDOF_to_lmdof[mDOF] - a1[lmdof,ldof] = coeff + out_offsets = copy(sDOF_to_offsets) + out_data_d = Vector{Int}(undef, n_data) + out_data_c = Vector{T}(undef, n_data) + ptrs = zeros(Int32, n_slaves + 1) + ptrs[1] = Int32(1) + acc = Dict{Int, T}() + k = 0 + + for i in 1:n_slaves + s = order[i] + rd = dataview(sDOF_to_dofs, s) + rc = dataview(sDOF_to_coeffs, s) + + any_sub = any(dof_to_sDOF[_dof_to_DOF(d, n_fdofs)] != 0 for d in rd) + + if any_sub + empty!(acc) + for kk in eachindex(rd) + d = rd[kk]; c = rc[kk] + t = dof_to_sDOF[_dof_to_DOF(d, n_fdofs)] + if t != 0 + j = order_rev[t] + for m in ptrs[j]:ptrs[j+1]-1 + dd = out_data_d[m] + acc[dd] = get(acc, dd, zero(T)) + c * out_data_c[m] + end + out_offsets[s] += c * out_offsets[t] + else + acc[d] = get(acc, d, zero(T)) + c + end + end + end + + for (d, c) in (any_sub ? acc : zip(rd, rc)) + abs(c) > tol || continue + k += 1; out_data_d[k] = d; out_data_c[k] = c end + ptrs[i+1] = Int32(k + 1) end - #TODO this is not always needed - mul!(a2,a1,mat) - a2 + resize!(out_data_d, k) + resize!(out_data_c, k) + out_dof = [sDOF_to_dof[order[i]] for i in 1:n_slaves] + out_off = [out_offsets[order[i]] for i in 1:n_slaves] + return out_dof, Table(out_data_d, ptrs), Table(out_data_c, ptrs), out_off end diff --git a/test/FESpacesTests/ConstraintHandlerTests.jl b/test/FESpacesTests/ConstraintHandlerTests.jl new file mode 100644 index 000000000..ebf623929 --- /dev/null +++ b/test/FESpacesTests/ConstraintHandlerTests.jl @@ -0,0 +1,487 @@ +module ConstraintHandlerTests + +using Test +using LinearAlgebra +using SparseArrays +using Gridap +using Gridap.FESpaces +using Gridap.Arrays: Table, datarange + +# ----------------------------------------------------------------------- +# Construction and basic queries +# ----------------------------------------------------------------------- + +ch = ConstraintHandler(10) +@test num_constrained_dofs(ch) == 0 +@test num_free_dofs(ch) == 10 +@test !is_constrained(ch, 3) + +add_constraint!(ch, 3, [1, 2], [0.5, 0.5]) +@test is_constrained(ch, 3) +@test num_constrained_dofs(ch) == 1 +@test num_free_dofs(ch) == 9 + +# ----------------------------------------------------------------------- +# Dirichlet (pure offset) +# ----------------------------------------------------------------------- + +ch = ConstraintHandler(5) +add_constraint!(ch, 2, 3.0) +close!(ch) +x = zeros(5) +apply_constraints!(x, ch) +@test x[2] ≈ 3.0 + +# ----------------------------------------------------------------------- +# Homogeneous linear constraint +# ----------------------------------------------------------------------- + +ch = ConstraintHandler(5) +add_constraint!(ch, 3, [1, 2], [0.5, 0.5]) +close!(ch) +x = [2.0, 4.0, 0.0, 0.0, 0.0] +apply_constraints!(x, ch) +@test x[3] ≈ 3.0 + +# ----------------------------------------------------------------------- +# Affine constraint (linear + offset) +# ----------------------------------------------------------------------- + +ch = ConstraintHandler(4) +add_constraint!(ch, 4, [1], [2.0], 1.0) +close!(ch) +x = [3.0, 0.0, 0.0, 0.0] +apply_constraints!(x, ch) +@test x[4] ≈ 7.0 + +# ----------------------------------------------------------------------- +# set_offset! +# ----------------------------------------------------------------------- + +ch = ConstraintHandler(3) +add_constraint!(ch, 2, 1.0) +set_offset!(ch, 2, 5.0) +close!(ch) +x = zeros(3) +apply_constraints!(x, ch) +@test x[2] ≈ 5.0 + +# ----------------------------------------------------------------------- +# Chain resolution: one level +# ----------------------------------------------------------------------- + +ch = ConstraintHandler(5) +add_constraint!(ch, 2, [1], [1.0]) +add_constraint!(ch, 3, [2], [1.0]) # x3 = x2 = x1 +close!(ch) + +line3 = ch.constraints[findfirst(l -> l.dof == 3, ch.constraints)] +@test length(line3.masters) == 1 +@test line3.masters[1] == 1 +@test line3.coeffs[1] ≈ 1.0 + +x = [7.0, 0.0, 0.0, 0.0, 0.0] +apply_constraints!(x, ch) +@test x[2] ≈ 7.0 +@test x[3] ≈ 7.0 + +# ----------------------------------------------------------------------- +# Chain resolution: weighted chain +# ----------------------------------------------------------------------- + +ch = ConstraintHandler(4) +add_constraint!(ch, 2, [1], [2.0]) +add_constraint!(ch, 3, [2], [3.0]) # x3 = 6*x1 +close!(ch) + +line3 = ch.constraints[findfirst(l -> l.dof == 3, ch.constraints)] +@test length(line3.masters) == 1 +@test line3.masters[1] == 1 +@test line3.coeffs[1] ≈ 6.0 + +x = [1.0, 0.0, 0.0, 0.0] +apply_constraints!(x, ch) +@test x[2] ≈ 2.0 +@test x[3] ≈ 6.0 + +# ----------------------------------------------------------------------- +# Chain resolution: offset propagation +# ----------------------------------------------------------------------- + +# x2 = x1 + 1, x3 = x2 + 2 → x3 = x1 + 3 +ch = ConstraintHandler(5) +add_constraint!(ch, 2, [1], [1.0], 1.0) +add_constraint!(ch, 3, [2], [1.0], 2.0) +close!(ch) + +x = [10.0, 0.0, 0.0, 0.0, 0.0] +apply_constraints!(x, ch) +@test x[2] ≈ 11.0 +@test x[3] ≈ 13.0 + +# ----------------------------------------------------------------------- +# Duplicate master entry merging +# ----------------------------------------------------------------------- + +ch = ConstraintHandler(3) +add_constraint!(ch, 3, [1, 1, 2], [0.3, 0.2, 0.5]) +close!(ch) + +masters = Dict(zip(ch.constraints[1].masters, ch.constraints[1].coeffs)) +@test masters[1] ≈ 0.5 +@test masters[2] ≈ 0.5 + +# ----------------------------------------------------------------------- +# Near-zero coefficient stripping +# ----------------------------------------------------------------------- + +ch = ConstraintHandler(3) +add_constraint!(ch, 3, [1, 2], [1.0, 1e-15]) +close!(ch) +@test length(ch.constraints[1].masters) == 1 +@test ch.constraints[1].masters[1] == 1 + +# ----------------------------------------------------------------------- +# Idempotent close! +# ----------------------------------------------------------------------- + +ch = ConstraintHandler(3) +add_constraint!(ch, 2, [1], [1.0]) +close!(ch) +close!(ch) +@test ch.closed +@test num_constrained_dofs(ch) == 1 + +# ----------------------------------------------------------------------- +# Cycle detection +# ----------------------------------------------------------------------- + +ch = ConstraintHandler(5) +add_constraint!(ch, 2, [3], [1.0]) +add_constraint!(ch, 3, [2], [1.0]) # cycle: 2 ↔ 3 +@test_throws ErrorException close!(ch) +@test occursin("cycle", sprint(showerror, try close!(ch) catch e e end)) + +# ----------------------------------------------------------------------- +# merge!: no conflict +# ----------------------------------------------------------------------- + +ch1 = ConstraintHandler(6) +add_constraint!(ch1, 1, [3], [1.0]) +ch2 = ConstraintHandler(6) +add_constraint!(ch2, 2, [4, 5], [0.5, 0.5]) +merge!(ch1, ch2) +close!(ch1) +@test is_constrained(ch1, 1) +@test is_constrained(ch1, 2) +x = [0.0, 0.0, 7.0, 2.0, 4.0, 0.0] +apply_constraints!(x, ch1) +@test x[1] ≈ 7.0 +@test x[2] ≈ 3.0 + +# ----------------------------------------------------------------------- +# merge!: chain across two handlers +# ----------------------------------------------------------------------- + +ch1 = ConstraintHandler(4) +add_constraint!(ch1, 3, [2], [1.0]) # x3 = x2 +ch2 = ConstraintHandler(4) +add_constraint!(ch2, 2, [1], [2.0]) # x2 = 2*x1 +merge!(ch1, ch2) +close!(ch1) # should resolve x3 = x2 = 2*x1 + +line3 = ch1.constraints[findfirst(l -> l.dof == 3, ch1.constraints)] +@test length(line3.masters) == 1 +@test line3.masters[1] == 1 +@test line3.coeffs[1] ≈ 2.0 + +# ----------------------------------------------------------------------- +# merge!: conflict raises error by default +# ----------------------------------------------------------------------- + +ch1 = ConstraintHandler(4) +add_constraint!(ch1, 2, [1], [1.0]) +ch2 = ConstraintHandler(4) +add_constraint!(ch2, 2, [3], [1.0]) # conflict on dof 2 +@test_throws ErrorException merge!(ch1, ch2) + +# ----------------------------------------------------------------------- +# merge!: user-supplied conflict resolver +# ----------------------------------------------------------------------- + +ch1 = ConstraintHandler(4) +add_constraint!(ch1, 2, [1], [1.0]) +ch2 = ConstraintHandler(4) +add_constraint!(ch2, 2, [3], [1.0]) + +# Average the two constraints: x2 = 0.5*x1 + 0.5*x3 +resolve = (existing, incoming) -> + ConstraintLine(existing.dof, + vcat(existing.masters, incoming.masters), + vcat(existing.coeffs, incoming.coeffs) .* 0.5, + 0.0) + +merge!(ch1, ch2; on_conflict=resolve) +close!(ch1) +x = [4.0, 0.0, 6.0, 0.0] +apply_constraints!(x, ch1) +@test x[2] ≈ 5.0 + +# ----------------------------------------------------------------------- +# apply_constraints (system reduction): pure Dirichlet +# ----------------------------------------------------------------------- + +# 2 dofs: x1 free, x2 = 3.0 (Dirichlet) +# Ã = [2], b̃ = [5 - 1*3] = [2] → u = [1] → x = [1, 3] +A = sparse([2.0 1.0; 1.0 2.0]) +b = [5.0, 6.0] +ch = ConstraintHandler(2) +add_constraint!(ch, 2, 3.0) +close!(ch) + +Ã, b̃ = apply_constraints(A, b, ch) +C, g = get_matrix(ch) + +@test size(Ã) == (1, 1) +u = Ã \ b̃ +x = C * u + g +@test x[1] ≈ 1.0 +@test x[2] ≈ 3.0 + +# ----------------------------------------------------------------------- +# apply_constraints (system reduction): linear constraint +# ----------------------------------------------------------------------- + +# 3 dofs: x3 = 0.5*x1 + 0.5*x2, A = I +# b = A * [1, 2, 1.5] = [1, 2, 1.5] (consistent with constraint) +n = 3 +A = sparse(1.0I, n, n) +x_target = [1.0, 2.0, 1.5] +b = A * x_target + +ch = ConstraintHandler(n) +add_constraint!(ch, 3, [1, 2], [0.5, 0.5]) +close!(ch) + +Ã, b̃ = apply_constraints(A, b, ch) +C, g = get_matrix(ch) + +@test size(Ã) == (2, 2) +u = Ã \ b̃ +x = C * u + g +@test x[1] ≈ 1.0 +@test x[2] ≈ 2.0 +@test x[3] ≈ 1.5 + +# ----------------------------------------------------------------------- +# apply_constraints (system reduction): full system round-trip +# ----------------------------------------------------------------------- + +# Poisson-like system with a known analytical solution. +# 4 dofs: symmetric tridiagonal stiffness with constraints: +# x1 = 0 (Dirichlet left) +# x4 = 1 (Dirichlet right) +n = 4 +A = sparse([2.0 -1.0 0.0 0.0; + -1.0 2.0 -1.0 0.0; + 0.0 -1.0 2.0 -1.0; + 0.0 0.0 -1.0 2.0]) +b = zeros(n) + +ch = ConstraintHandler(n) +add_constraint!(ch, 1, 0.0) +add_constraint!(ch, 4, 1.0) +close!(ch) + +Ã, b̃ = apply_constraints(A, b, ch) +C, g = get_matrix(ch) + +@test size(Ã) == (2, 2) +u = Ã \ b̃ +x = C * u + g + +# Linear interpolation between 0 and 1 is exact for this 1D Poisson. +@test x[1] ≈ 0.0 atol=1e-12 +@test x[2] ≈ 1/3 atol=1e-12 +@test x[3] ≈ 2/3 atol=1e-12 +@test x[4] ≈ 1.0 atol=1e-12 + +# ----------------------------------------------------------------------- +# Error cases +# ----------------------------------------------------------------------- + +ch = ConstraintHandler(5) +close!(ch) +@test_throws AssertionError add_constraint!(ch, 1, 0.0) + +ch = ConstraintHandler(5) +add_constraint!(ch, 2, [1], [1.0]) +@test_throws ErrorException add_constraint!(ch, 2, [3], [1.0]) + +ch = ConstraintHandler(5) +@test_throws AssertionError add_constraint!(ch, 0, 0.0) +@test_throws AssertionError add_constraint!(ch, 6, 0.0) + +ch = ConstraintHandler(5) +@test_throws AssertionError add_constraint!(ch, 3, [3], [1.0]) + +ch = ConstraintHandler(3) +add_constraint!(ch, 2, [1], [1.0]) +@test_throws AssertionError apply_constraints!(zeros(3), ch) + +ch1 = ConstraintHandler(3) +close!(ch1) +ch2 = ConstraintHandler(3) +@test_throws AssertionError merge!(ch1, ch2) + +ch1 = ConstraintHandler(4) +ch2 = ConstraintHandler(5) +@test_throws AssertionError merge!(ch1, ch2) + +# ----------------------------------------------------------------------- +# constraint_tables — requires a real FESpace +# ----------------------------------------------------------------------- + +# Space: 6 free + 3 Dirichlet dofs (2×2 Cartesian mesh, Lagrangian order 1, +# Dirichlet on left edge x=0). +# ConstraintHandler DoF convention: 1..6 free, 7..9 Dirichlet. +# Gridap signed-dof convention: 1..6 free, -1..-3 Dirichlet. + +_model = CartesianDiscreteModel((0,1,0,1), (2,2)) +_reffe = ReferenceFE(lagrangian, Float64, 1) +let labels = get_face_labeling(_model) + add_tag_from_tags!(labels, "left", ["tag_7", "tag_1", "tag_3"]) +end +_space = FESpace(_model, _reffe; dirichlet_tags="left") + +# Free master +ch = ConstraintHandler(num_free_dofs(_space) + num_dirichlet_dofs(_space)) +add_constraint!(ch, 3, [1, 2], [0.5, 0.5]) +close!(ch) +sDOF_to_dof, sDOF_to_dofs, sDOF_to_coeffs = constraint_tables(_space, ch) +@test sDOF_to_dof == [3] +@test collect(sDOF_to_dofs[1]) == [1, 2] +@test collect(sDOF_to_coeffs[1]) ≈ [0.5, 0.5] + +# Dirichlet master (ConstraintHandler DoF 7 → signed dof -1) +ch = ConstraintHandler(num_free_dofs(_space) + num_dirichlet_dofs(_space)) +add_constraint!(ch, 3, [7], [1.0]) +close!(ch) +sDOF_to_dof, sDOF_to_dofs, sDOF_to_coeffs = constraint_tables(_space, ch) +@test sDOF_to_dof == [3] +@test collect(sDOF_to_dofs[1]) == [-1] + +# Slave is a Dirichlet dof (ConstraintHandler DoF 7 → signed dof -1) +ch = ConstraintHandler(num_free_dofs(_space) + num_dirichlet_dofs(_space)) +add_constraint!(ch, 7, [1], [1.0]) +close!(ch) +sDOF_to_dof, sDOF_to_dofs, sDOF_to_coeffs = constraint_tables(_space, ch) +@test sDOF_to_dof == [-1] +@test collect(sDOF_to_dofs[1]) == [1] + +# Chain resolved before output: x3 = x2 = x1 +ch = ConstraintHandler(num_free_dofs(_space) + num_dirichlet_dofs(_space)) +add_constraint!(ch, 2, [1], [1.0]) +add_constraint!(ch, 3, [2], [1.0]) +close!(ch) +sDOF_to_dof, sDOF_to_dofs, sDOF_to_coeffs = constraint_tables(_space, ch) +i3 = findfirst(==(3), sDOF_to_dof) +@test collect(sDOF_to_dofs[i3]) == [1] +@test collect(sDOF_to_coeffs[i3]) ≈ [1.0] + +# Round-trip: constraint_tables → ConstraintHandler +ch_orig = ConstraintHandler(num_free_dofs(_space) + num_dirichlet_dofs(_space)) +add_constraint!(ch_orig, 3, [1, 2], [0.5, 0.5]) +add_constraint!(ch_orig, 5, [4], [1.0]) +close!(ch_orig) +sDOF_to_dof, sDOF_to_dofs, sDOF_to_coeffs = constraint_tables(_space, ch_orig) +ch_rt = ConstraintHandler(sDOF_to_dof, sDOF_to_dofs, sDOF_to_coeffs, _space) +close!(ch_rt) +@test num_constrained_dofs(ch_rt) == num_constrained_dofs(ch_orig) +@test is_constrained(ch_rt, 3) +@test is_constrained(ch_rt, 5) + + +# ----------------------------------------------------------------------- +# merge_slave_constraint_tables / close_slave_constraint_tables +# ----------------------------------------------------------------------- + +using Gridap.FESpaces: merge_slave_constraint_tables, close_slave_constraint_tables + +# Helper: slave_dof → Dict(master_dof → coeff), order-independent comparison. +_cdict(dofs, dtable, ctable) = Dict( + dofs[i] => Dict(zip(collect(dtable[i]), collect(ctable[i]))) + for i in eachindex(dofs)) + +_n_dofs = num_free_dofs(_space) + num_dirichlet_dofs(_space) + +# Merge: disjoint sets → fast path, all slaves present in output. +let + ch1 = ConstraintHandler(_n_dofs); add_constraint!(ch1, 3, [1,2], [0.5,0.5]); close!(ch1) + ch2 = ConstraintHandler(_n_dofs); add_constraint!(ch2, 5, [4], [1.0]); close!(ch2) + d1,t1,c1,o1 = constraint_tables(_space, ch1) + d2,t2,c2,o2 = constraint_tables(_space, ch2) + out_dof, out_dofs, out_coeffs, _ = merge_slave_constraint_tables( + _space, d1,t1,c1, d2,t2,c2, o1,o2) + @test Set(out_dof) == Set([3, 5]) + @test _cdict(out_dof, out_dofs, out_coeffs)[3] == Dict(1 => 0.5, 2 => 0.5) + @test _cdict(out_dof, out_dofs, out_coeffs)[5] == Dict(4 => 1.0) +end + +# Merge: conflict → on_conflict (first-wins keeps s1 row). +let + ch1 = ConstraintHandler(_n_dofs); add_constraint!(ch1, 3, [1], [1.0]); close!(ch1) + ch2 = ConstraintHandler(_n_dofs); add_constraint!(ch2, 3, [2], [1.0]); close!(ch2) + d1,t1,c1,o1 = constraint_tables(_space, ch1) + d2,t2,c2,o2 = constraint_tables(_space, ch2) + first_wins = (dof,md1,mc1,mo1,md2,mc2,mo2) -> (collect(md1), collect(mc1), mo1) + out_dof, out_dofs, out_coeffs, _ = merge_slave_constraint_tables( + _space, d1,t1,c1, d2,t2,c2, o1,o2; on_conflict=first_wins) + @test _cdict(out_dof, out_dofs, out_coeffs)[3] == Dict(1 => 1.0) +end + +# Hard test: path A (merge ConstraintHandlers → close → tables) must equal +# path B (tables per handler → merge_slave → close_slave). +# +# ch1: dof3=dof2, dof5=0.5*dof1+0.5*dof4 +# ch2: dof2=dof1, dof5=dof4 (conflict on dof5; ch1/s1 wins) +# +# After merge+close: dof2=dof1, dof3=dof1 (chain dof3→dof2→dof1), dof5=0.5*dof1+0.5*dof4. +# Path B exercises close_slave_constraint_tables on the cross-handler chain dof3→dof2. +let + first_wins_ch = (existing, _) -> existing + first_wins_tbl = (dof,md1,mc1,mo1,md2,mc2,mo2) -> (collect(md1), collect(mc1), mo1) + + # Path A + ch_a1 = ConstraintHandler(_n_dofs) + add_constraint!(ch_a1, 3, [2], [1.0]) + add_constraint!(ch_a1, 5, [1,4], [0.5, 0.5]) + ch_a2 = ConstraintHandler(_n_dofs) + add_constraint!(ch_a2, 2, [1], [1.0]) + add_constraint!(ch_a2, 5, [4], [1.0]) + merge!(ch_a1, ch_a2; on_conflict=first_wins_ch) + close!(ch_a1) + dof_A, dofs_A, coeffs_A, _ = constraint_tables(_space, ch_a1) + + # Path B + ch_b1 = ConstraintHandler(_n_dofs) + add_constraint!(ch_b1, 3, [2], [1.0]) + add_constraint!(ch_b1, 5, [1,4], [0.5, 0.5]) + close!(ch_b1) + d1,t1,c1,o1 = constraint_tables(_space, ch_b1) + + ch_b2 = ConstraintHandler(_n_dofs) + add_constraint!(ch_b2, 2, [1], [1.0]) + add_constraint!(ch_b2, 5, [4], [1.0]) + close!(ch_b2) + d2,t2,c2,o2 = constraint_tables(_space, ch_b2) + + dof_B, dofs_B, coeffs_B, offs_B = merge_slave_constraint_tables( + _space, d1,t1,c1, d2,t2,c2, o1,o2; on_conflict=first_wins_tbl) + dof_B, dofs_B, coeffs_B, _ = close_slave_constraint_tables( + _space, dof_B, dofs_B, coeffs_B, offs_B) + + @test _cdict(dof_A, dofs_A, coeffs_A) == _cdict(dof_B, dofs_B, coeffs_B) +end + +end # module \ No newline at end of file diff --git a/test/FESpacesTests/FESpacesWithLinearConstraintsTests.jl b/test/FESpacesTests/FESpacesWithLinearConstraintsTests.jl index fc5f6222d..b7008601a 100644 --- a/test/FESpacesTests/FESpacesWithLinearConstraintsTests.jl +++ b/test/FESpacesTests/FESpacesWithLinearConstraintsTests.jl @@ -1,18 +1,18 @@ module FESpacesWithLinearConstraintsTests +using Test +using LinearAlgebra +using Gridap + using Gridap.Algebra using Gridap.Arrays using Gridap.Fields using Gridap.Geometry using Gridap.FESpaces -using Test -using LinearAlgebra using Gridap.CellData using Gridap.ReferenceFEs -domain = (0,1,0,1) -partition = (2,2) -model = CartesianDiscreteModel(domain,partition) +model = CartesianDiscreteModel((0,1,0,1),(2,2)) labels = get_face_labeling(model) add_tag_from_tags!(labels,"dirichlet",[1,2,5]) @@ -26,7 +26,9 @@ dΓ = Measure(Γ,2) dΛ = Measure(Λ,2) V = FESpace( - model,ReferenceFE(lagrangian,Float64,1), conformity=:H1, dirichlet_tags="dirichlet") + model,ReferenceFE(lagrangian,Float64,1); + conformity=:H1, dirichlet_tags="dirichlet" +) test_single_field_fe_space(V) fdof_to_val = collect(Float64,1:num_free_dofs(V)) @@ -38,18 +40,14 @@ sDOF_to_dofs = Table([[-1,4],[4,6],[-1,-3]]) sDOF_to_coeffs = Table([[0.5,0.5],[0.5,0.5],[0.5,0.5]]) Vc = FESpaceWithLinearConstraints( - sDOF_to_dof, - sDOF_to_dofs, - sDOF_to_coeffs, - V) + sDOF_to_dof, sDOF_to_dofs, sDOF_to_coeffs, V +) test_single_field_fe_space(Vc) @test has_constraints(Vc) @test isa(get_cell_constraints(Vc,Λ)[1],ArrayBlock) - -@test Vc.n_fdofs == 6 -@test Vc.n_fmdofs == 4 +@test num_free_dofs(Vc) == 4 fmdof_to_val = collect(Float64,1:num_free_dofs(Vc)) dmdof_to_val = -collect(Float64,1:num_dirichlet_dofs(Vc)) @@ -60,7 +58,6 @@ r = [[-1.0, -1.5, 1.0, 1.0], [-1.5, -2.0, 1.0, 2.0], [1.0, 1.0, 3.0, 3.5], [1.0, v(x) = sin(4*x[1]+0.4)*cos(5*x[2]+0.7) vch = interpolate(v,Vc) - #using Gridap.Visualization #writevtk(Ω,"trian",nsubcells=10,cellfields=["vh"=>vh,"vch"=>vch]) @@ -71,7 +68,6 @@ Uc = TrialFESpace(Vc,u) uch = interpolate(u,Uc) n_Γ = get_normal_vector(Γ) - a(u,v) = ∫( ∇(v)⋅∇(u) )*dΩ + ∫( jump(u)*jump(v) )*dΛ l(v) = ∫( v*f )*dΩ + ∫( v*(n_Γ⋅∇(u)) )*dΓ @@ -82,7 +78,6 @@ uch = solve(op) #writevtk(trian,"trian",nsubcells=10,cellfields=["uch"=>uch]) e = u - uch - e_l2 = sqrt(sum(∫( e*e )*dΩ)) e_h1 = sqrt(sum(∫( e*e + ∇(e)⋅∇(e) )*dΩ)) @@ -90,15 +85,91 @@ tol = 1.e-9 @test e_l2 < tol @test e_h1 < tol +# Test with complex values + V2 = FESpace( - model,ReferenceFE(lagrangian,Float64,1), conformity=:H1, dirichlet_tags="dirichlet", vector_type=ComplexF64) + model,ReferenceFE(lagrangian,Float64,1); + conformity=:H1, dirichlet_tags="dirichlet", vector_type=ComplexF64 +) Vc2 = FESpaceWithLinearConstraints( sDOF_to_dof, sDOF_to_dofs, sDOF_to_coeffs, - V2) + V2 +) @test get_dof_value_type(Vc2) <: ComplexF64 +# Alternative constructor + +fdof_to_dofs = Table([[-1,4],[2],[3],[4],[4,6],[6]]) +fdof_to_coeffs = Table([[0.5,0.5],[1.0],[1.0],[1.0],[0.5,0.5],[1.0]]) +ddof_to_dofs = Table([[-1],[-1,-3],[-3]]) +ddof_to_coeffs = Table([[1.0],[0.5,0.5],[1.0]]) + +Vc3 = FESpaceWithLinearConstraints( + fdof_to_dofs, fdof_to_coeffs, ddof_to_dofs, ddof_to_coeffs, V +) + +test_single_field_fe_space(Vc3) + +@test Vc3.mDOF_to_dof == Vc.mDOF_to_dof +@test Vc3.sDOF_to_dof == Vc.sDOF_to_dof +@test Vc3.sDOF_to_mdofs == Vc.sDOF_to_mdofs +@test Vc3.sDOF_to_coeffs == Vc.sDOF_to_coeffs + +################################################################################### +# Affine constraints — nonzero offsets + +# All-zero offsets: no fictitious masters should be added, result identical to Vc +Vca0 = FESpaceWithLinearConstraints( + sDOF_to_dof, sDOF_to_dofs, sDOF_to_coeffs, V; sDOF_to_offsets=zeros(3) +) +@test num_free_dofs(Vca0) == num_free_dofs(Vc) +@test num_dirichlet_dofs(Vca0) == num_dirichlet_dofs(Vc) +@test isempty(Vca0.dmdof_to_offsets) + +# Two nonzero offsets: slaves 1 and 3 get shifts of +2 and -1 respectively +sDOF_to_offsets_affine = [2.0, 0.0, -1.0] +Vca = FESpaceWithLinearConstraints( + sDOF_to_dof, sDOF_to_dofs, sDOF_to_coeffs, V; sDOF_to_offsets=sDOF_to_offsets_affine +) + +test_single_field_fe_space(Vca) +@test has_constraints(Vca) +@test num_free_dofs(Vca) == 4 +@test num_dirichlet_dofs(Vca) == 4 # 2 real Dirichlet masters + 2 fictitious +@test Vca.dmdof_to_offsets ≈ [2.0, -1.0] # nonzero offsets, in slave order + +dvals_ca = get_dirichlet_dof_values(Vca) # get_dirichlet_dof_values injects offsets into the last entries +@test length(dvals_ca) == 4 +@test dvals_ca[end-1:end] ≈ [2.0, -1.0] + +# Scatter test: same master values as the homogeneous case, plus offset values for the fictitious masters. +# dmdof layout: [Dir master 1, Dir master 2, fictitious for slave-1 offset, fictitious for slave-3 offset] +fmdof_to_val_a = collect(Float64, 1:num_free_dofs(Vca)) # [1, 2, 3, 4] +dmdof_to_val_a = vcat(-collect(Float64, 1:2), [2.0, -1.0]) # [-1, -2, 2.0, -1.0] +vca = FEFunction(Vca, fmdof_to_val_a, dmdof_to_val_a) + +# Expected slave values: +# free DOF 1: 0.5*(-1.0) + 0.5*3.0 + 1.0*2.0 = 3.0 (homogeneous was 1.0) +# free DOF 5: 0.5*3.0 + 0.5*4.0 = 3.5 (offset 0 → unchanged) +# Dir DOF 2: 0.5*(-1.0) + 0.5*(-2.0) + 1.0*(-1.0) = -2.5 (homogeneous was -1.5) +r_affine = [[-1.0, -2.5, 3.0, 1.0], [-2.5, -2.0, 1.0, 2.0], + [3.0, 1.0, 3.0, 3.5], [1.0, 2.0, 3.5, 4.0]] +@test get_cell_dof_values(vca) ≈ r_affine + +################################################################################### +# ConstraintHandler constructor +n_dofs = num_free_dofs(V) + num_dirichlet_dofs(V) +ch = ConstraintHandler(n_dofs) +add_constraint!(ch, 1, [4], [0.5], 0.0) # x1 = 0.5*x4 (free slave, free master) +add_constraint!(ch, 5, [4, 6], [0.5, 0.5]) # x5 = 0.5*(x4+x6) +close!(ch) + +Vc4 = FESpaceWithLinearConstraints(V, ch) +test_single_field_fe_space(Vc4) +@test num_free_dofs(Vc4) == num_free_dofs(V) - 2 + end # module diff --git a/test/FESpacesTests/runtests_2.jl b/test/FESpacesTests/runtests_2.jl index eefab398b..0678eac4f 100644 --- a/test/FESpacesTests/runtests_2.jl +++ b/test/FESpacesTests/runtests_2.jl @@ -12,6 +12,8 @@ using Test @testset "ExtendedFESpaces" begin include("ExtendedFESpacesTests.jl") end +@testset "ConstraintHandler" begin include("ConstraintHandlerTests.jl") end + @testset "FESpacesWithLinearConstraints" begin include("FESpacesWithLinearConstraintsTests.jl") end @testset "FEAutodiff" begin include("FEAutodiffTests.jl") end