From 464b5bf2c9e1a7550c2b5fe4ab9332d85adce0bc Mon Sep 17 00:00:00 2001 From: Jordi Manyer Date: Sat, 22 Nov 2025 19:37:32 +1100 Subject: [PATCH 1/2] Refactor of FESpaceWithLinearConstraints --- src/FESpaces/FESpaces.jl | 1 + src/FESpaces/FESpacesWithLinearConstraints.jl | 794 +++++++++--------- .../FESpacesWithLinearConstraintsTests.jl | 54 +- 3 files changed, 445 insertions(+), 404 deletions(-) diff --git a/src/FESpaces/FESpaces.jl b/src/FESpaces/FESpaces.jl index c8ee6f871..40a5a97fd 100644 --- a/src/FESpaces/FESpaces.jl +++ b/src/FESpaces/FESpaces.jl @@ -12,6 +12,7 @@ using SparseArrays using LinearAlgebra using StaticArrays using ForwardDiff +using DataStructures using Gridap.Helpers using Gridap.Algebra diff --git a/src/FESpaces/FESpacesWithLinearConstraints.jl b/src/FESpaces/FESpacesWithLinearConstraints.jl index a1a5771fd..db9435d2e 100644 --- a/src/FESpaces/FESpacesWithLinearConstraints.jl +++ b/src/FESpaces/FESpacesWithLinearConstraints.jl @@ -52,90 +52,60 @@ # # - 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 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 + 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_dofs::Table, + sDOF_to_mdofs::Table, sDOF_to_coeffs::Table, - space::SingleFieldFESpace) - - n_fdofs = num_free_dofs(space) - n_ddofs = num_dirichlet_dofs(space) - n_DOFs = n_fdofs + n_ddofs - - DOF_to_DOFs, DOF_to_coeffs = _prepare_DOF_to_DOFs( - sDOF_to_dof, sDOF_to_dofs, sDOF_to_coeffs, n_fdofs, n_DOFs) - - FESpaceWithLinearConstraints!(DOF_to_DOFs,DOF_to_coeffs,space) - + n_fmdofs::Int = _count_free_mdofs(mDOF_to_dof,sDOF_to_mdofs) +) + 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, n_fmdofs + ) 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) - - DOF_to_DOFs_ptrs = ones(Tp,n_DOFs+1) - - n_sDOFs = length(sDOF_to_dof) - - 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 - - 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) - - for DOF in 1:n_DOFs - q = DOF_to_DOFs_ptrs[DOF] - DOF_to_DOFs_data[q] = DOF - end - - 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 - - DOF_to_DOFs = Table(DOF_to_DOFs_data,DOF_to_DOFs_ptrs) - DOF_to_coeffs = Table(DOF_to_coeffs_data,DOF_to_DOFs_ptrs) - - DOF_to_DOFs, DOF_to_coeffs +function FESpaceWithLinearConstraints( + sDOF_to_dof::AbstractVector{<:Integer}, + sDOF_to_dofs::Table, + sDOF_to_coeffs::Table, + space::SingleFieldFESpace +) + mDOF_to_dof, sDOF_to_mdofs, nfmdofs = _find_master_dofs( + sDOF_to_dof, sDOF_to_dofs, space + ) + return FESpaceWithLinearConstraints( + space, mDOF_to_dof, sDOF_to_dof, sDOF_to_mdofs, sDOF_to_coeffs, nfmdofs + ) end -# Private constructors +function FESpaceWithLinearConstraints( + DOF_to_dofs::Table, DOF_to_coeffs::Table, space::SingleFieldFESpace +) + sDOF_to_dof, sDOF_to_dofs, sDOF_to_coeffs = _filter_constraints( + DOF_to_dofs, DOF_to_coeffs, space + ) + return FESpaceWithLinearConstraints( + sDOF_to_dof, sDOF_to_dofs, sDOF_to_coeffs, space + ) +end function FESpaceWithLinearConstraints( fdof_to_dofs::Table, @@ -144,131 +114,10 @@ function FESpaceWithLinearConstraints( 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) - FESpaceWithLinearConstraints!(DOF_to_DOFs,DOF_to_coeffs,space) -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 -end - - -function FESpaceWithLinearConstraints(DOF_to_DOFs::Table,DOF_to_coeffs::Table,space::SingleFieldFESpace) - FESpaceWithLinearConstraints!(copy(DOF_to_DOFs),copy(DOF_to_coeffs),space::SingleFieldFESpace) -end - -function FESpaceWithLinearConstraints!(DOF_to_DOFs::Table,DOF_to_coeffs::Table,space::SingleFieldFESpace) - - 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) - -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 - end - end - n_fmdofs = 0 - for DOF in 1:n_fdofs - if DOF_to_ismaster[DOF] - n_fmdofs += 1 - end - end - mDOF_to_DOF = findall(DOF_to_ismaster) - mDOF_to_DOF, n_fmdofs -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 - end - DOF_to_DOFs -end - -function _setup_cell_to_lmdof_to_mdof(cell_to_ldof_to_dof,DOF_to_mDOFs,n_fdofs,n_fmdofs) - - 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) - - 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 - cell_to_lmdof_to_mdof_ptrs[cell+1] = length(mdofs) - 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) - - 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 - - Table(cell_to_lmdof_to_mdof_data,cell_to_lmdof_to_mdof_ptrs) + 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 function _dof_to_DOF(dof,n_fdofs) @@ -287,270 +136,261 @@ function _DOF_to_dof(DOF,n_fdofs) end end -# Implementation of the SingleFieldFESpace interface +# Implementation of FESpace interface -function get_cell_dof_ids(f::FESpaceWithLinearConstraints) - f.cell_to_lmdof_to_mdof -end +get_triangulation(f::FESpaceWithLinearConstraints) = get_triangulation(f.space) -function get_fe_dof_basis(f::FESpaceWithLinearConstraints) - get_fe_dof_basis(f.space) -end +get_free_dof_ids(f::FESpaceWithLinearConstraints) = Base.OneTo(f.n_free) +get_vector_type(f::FESpaceWithLinearConstraints) = get_vector_type(f.space) -get_dof_value_type(f::FESpaceWithLinearConstraints) = get_dof_value_type(f.space) +get_fe_basis(f::FESpaceWithLinearConstraints) = get_fe_basis(f.space) +get_trial_fe_basis(f::FESpaceWithLinearConstraints) = get_trial_fe_basis(f.space) + +CellField(f::FESpaceWithLinearConstraints,cellvals) = CellField(f.space,cellvals) + +# Implementation of the SingleFieldFESpace interface -get_dirichlet_dof_ids(f::FESpaceWithLinearConstraints) = Base.OneTo(length(f.mDOF_to_DOF) - f.n_fmdofs) +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) +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) 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 -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 + _ddof_to_dmdof_vals!( + dmdof_to_tag,ddof_to_tag,f.mDOF_to_dof + ) + return dmdof_to_tag 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 + _ddof_to_dmdof_vals!( + dmdof_to_tag,ddof_to_tag,f.mDOF_to_dof + ) + return dmdof_to_tag 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) + _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 -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 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 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) + + # 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 + + # 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 - for DOF in 1:length(DOF_to_mDOFs) - pini = DOF_to_mDOFs.ptrs[DOF] - pend = DOF_to_mDOFs.ptrs[DOF+1]-1 + # Map slave dofs + ptrs = sDOF_to_mdofs.ptrs + for (sDOF,dof) in enumerate(sDOF_to_dof) 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) + for k in ptrs[sDOF]:(ptrs[sDOF+1]-1) + mdof = sDOF_to_mdofs.data[k] + coeff = sDOF_to_coeffs.data[k] if mdof > 0 - fmdof = mdof - val += fmdof_to_val[fmdof]*coeff + val += fmdof_to_val[mdof]*coeff else - dmdof = -mdof - val += dmdof_to_val[dmdof]*coeff + val += dmdof_to_val[-mdof]*coeff end end - dof = _DOF_to_dof(DOF,n_fdofs) if dof > 0 - fdof = dof - fdof_to_val[fdof] = val + fdof_to_val[dof] = val else - ddof = -dof - ddof_to_val[ddof] = val + ddof_to_val[-dof] = val end end + return fdof_to_val, ddof_to_val 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 +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) -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 - -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) + # Map free master dofs + for (mdof,mDOF) in enumerate(1:n_fmdofs) + dof = mDOF_to_dof[mDOF] if dof > 0 - fdof = dof - val = fdof_to_val[fdof] - else - ddof = -dof - val = ddof_to_val[ddof] + fmdof_to_val[mdof] = fdof_to_val[dof] + elseif dof < 0 + fmdof_to_val[mdof] = ddof_to_val[-dof] 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 + end + + # 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 + return fmdof_to_val, dmdof_to_val end -# Implementation of FESpace interface - -function get_triangulation(f::FESpaceWithLinearConstraints) - get_triangulation(f.space) +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 + return dmdof_to_val end -get_free_dof_ids(f::FESpaceWithLinearConstraints) = Base.OneTo(f.n_fmdofs) - -function get_vector_type(f::FESpaceWithLinearConstraints) - get_vector_type(f.space) -end +# Constraints -function get_fe_basis(f::FESpaceWithLinearConstraints) - get_fe_basis(f.space) -end +ConstraintStyle(::Type{<:FESpaceWithLinearConstraints}) = Constrained() -function get_trial_fe_basis(f::FESpaceWithLinearConstraints) - get_trial_fe_basis(f.space) -end +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 -function CellField(f::FESpaceWithLinearConstraints,cellvals) - CellField(f.space,cellvals) -end + n_cells = length(cell_dofs) + n_mfdofs = num_free_dofs(f) + cell_isconstrained = Vector{Bool}(undef,n_cells) + for cell in 1:n_cells + 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 -ConstraintStyle(::Type{<:FESpaceWithLinearConstraints}) = Constrained() + cell_isconstrained[cell] = mask + end -function get_cell_isconstrained(f::FESpaceWithLinearConstraints) - #TODO this can be heavily optimized - n = length(get_cell_dof_ids(f)) - Fill(true,n) + return cell_isconstrained 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( - f.DOF_to_mDOFs, - f.DOF_to_coeffs, - length(f.mDOF_to_DOF), - f.n_fmdofs, - f.n_fdofs) - + 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) - lazy_map(k,f.cell_to_lmdof_to_mdof,f.cell_to_ldof_to_dof,cell_to_mat) - + return lazy_map(k,get_cell_dof_ids(f),get_cell_dof_ids(f.space),cell_to_mat) end struct LinearConstraintsMap{A,B} <: Map - DOF_to_mDOFs::A - DOF_to_coeffs::B - n_mDOFs::Int + DOF_to_msDOF::Vector{Int} + sDOF_to_mdofs::A + sDOF_to_coeffs::B + n_mdofs::Int n_fmdofs::Int n_fdofs::Int 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) +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(n_lmdofs,n_ldofs)) - m2 = CachedArray(zeros(n_lmdofs,n_ludofs)) - mDOF_to_lmdof = zeros(Int16,k.n_mDOFs) - m1, m2, mDOF_to_lmdof + 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 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) - +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,n_ludofs)) - a1 = m1.array - a2 = m2.array + setsize!(m2,(n_lmdofs,size(mat,2))) + a1 = get_array(m1) + a2 = get_array(m2) fill!(a1,zero(eltype(a1))) - for (lmdof,mdof) in enumerate(lmdof_to_mdof) - mDOF = _dof_to_DOF(mdof,k.n_fmdofs) - mDOF_to_lmdof[mDOF] = lmdof + # Precompute mdof to lmdof map + empty!(mdof_to_lmdof) + for (lmdof,mdof) in enumerate(mdofs) + mdof_to_lmdof[mdof] = lmdof end - for (ldof,dof) in enumerate(ldof_to_dof) + # 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) - 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 + 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 @@ -558,3 +398,185 @@ function evaluate!(cache,k::LinearConstraintsMap,lmdof_to_mdof,ldof_to_dof,mat) mul!(a2,a1,mat) a2 end + +# Private methods + +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 + DOF = _dof_to_DOF(dof,n_fdofs) + DOF_is_master[DOF] = true + end + + 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] = DOF_is_master[DOF] + end + + 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 + + return true +end + +function _filter_constraints(DOF_to_dofs, DOF_to_coeffs, 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] + return sDOF_to_dof, sDOF_to_dofs, Table(sDOF_to_coeffs.data,sDOF_to_dofs.ptrs) +end + +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 _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 _find_master_dofs( + sDOF_to_dof, sDOF_to_dofs, space +) + n_fdofs = num_free_dofs(space) + n_dofs = n_fdofs + num_dirichlet_dofs(space) + + # Flag master dofs + DOF_ismaster = fill(true,n_dofs) + for dof in sDOF_to_dof + DOF = _dof_to_DOF(dof,n_fdofs) + DOF_ismaster[DOF] = false + end + + # Counting mdofs + n_fmdofs = count(view(DOF_ismaster,1:n_fdofs)) + n_dmdofs = count(view(DOF_ismaster,(n_fdofs+1):n_dofs)) + n_mdofs = n_fmdofs + n_dmdofs + + # DOF to mdof mapping + kf, kd = 1, 1 + mDOF_to_dof = Vector{Int32}(undef,n_mdofs) + 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 + mDOF_to_dof[n_fmdofs + kd] = dof + DOF_to_mdof[DOF] = -kd + kd += 1 + end + end + end + + # Map constraints + data = zeros(Int32,length(sDOF_to_dofs.data)) + for (i,dof) in enumerate(sDOF_to_dofs.data) + DOF = _dof_to_DOF(dof,n_fdofs) + data[i] = DOF_to_mdof[DOF] + end + sDOF_to_mdofs = Table(data,sDOF_to_dofs.ptrs) + + return mDOF_to_dof, sDOF_to_mdofs, n_fmdofs +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) + + 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) + + 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 + + return Table(data,ptrs) +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) + + for (mDOF,dof) in enumerate(mDOF_to_dof) + DOF = _dof_to_DOF(dof,n_fdofs) + DOF_to_msDOF[DOF] = mDOF + end + + for (sDOF,dof) in enumerate(sDOF_to_dof) + DOF = _dof_to_DOF(dof,n_fdofs) + DOF_to_msDOF[DOF] = sDOF + n_mDOFs + end + + return DOF_to_msDOF +end diff --git a/test/FESpacesTests/FESpacesWithLinearConstraintsTests.jl b/test/FESpacesTests/FESpacesWithLinearConstraintsTests.jl index fc5f6222d..309f498c5 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,38 @@ 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 +# Test other constructors + +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 + end # module From 137e76b0f69daef8981bd4dd48fadb9085dbddc9 Mon Sep 17 00:00:00 2001 From: Jordi Manyer Date: Mon, 24 Nov 2025 15:26:49 +1100 Subject: [PATCH 2/2] Minor changes --- docs/src/dev-notes/constraints.md | 31 +++++++++++++++++++ src/FESpaces/FESpacesWithLinearConstraints.jl | 29 +++++++---------- .../FESpacesWithLinearConstraintsTests.jl | 2 +- 3 files changed, 43 insertions(+), 19 deletions(-) create mode 100644 docs/src/dev-notes/constraints.md 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/FESpacesWithLinearConstraints.jl b/src/FESpaces/FESpacesWithLinearConstraints.jl index db9435d2e..c604f52ae 100644 --- a/src/FESpaces/FESpacesWithLinearConstraints.jl +++ b/src/FESpaces/FESpacesWithLinearConstraints.jl @@ -121,19 +121,11 @@ function FESpaceWithLinearConstraints( end function _dof_to_DOF(dof,n_fdofs) - if dof > 0 - DOF = dof - else - DOF = n_fdofs - dof - end + ifelse(iszero(dof), 0, ifelse(dof > 0, dof, n_fdofs - dof)) end function _DOF_to_dof(DOF,n_fdofs) - if DOF > n_fdofs - dof = -(DOF-n_fdofs) - else - dof = DOF - end + ifelse(DOF > n_fdofs, -(DOF - n_fdofs), DOF) end # Implementation of FESpace interface @@ -167,12 +159,12 @@ function get_dirichlet_dof_tag(f::FESpaceWithLinearConstraints) 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)) + 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_tag,ddof_to_tag,f.mDOF_to_dof + dmdof_to_val,ddof_to_val,f.mDOF_to_dof ) - return dmdof_to_tag + return dmdof_to_val end function scatter_free_and_dirichlet_values(f::FESpaceWithLinearConstraints,fmdof_to_val,dmdof_to_val) @@ -334,7 +326,6 @@ function get_cell_constraints(f::FESpaceWithLinearConstraints) 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 - struct LinearConstraintsMap{A,B} <: Map DOF_to_msDOF::Vector{Int} sDOF_to_mdofs::A @@ -408,6 +399,7 @@ function _check_constraints( 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 @@ -415,7 +407,7 @@ function _check_constraints( 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] = DOF_is_master[DOF] + mDOF_is_master[mDOF] = iszero(DOF) || DOF_is_master[DOF] end c = array_cache(sDOF_to_mdofs) @@ -563,19 +555,20 @@ function _generate_cell_to_mdofs( end function generate_DOF_to_msDOF_map(space, mDOF_to_dof, sDOF_to_dof) - n_mDOFs = length(mDOF_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) for (mDOF,dof) in enumerate(mDOF_to_dof) + iszero(dof) && continue DOF = _dof_to_DOF(dof,n_fdofs) DOF_to_msDOF[DOF] = mDOF end for (sDOF,dof) in enumerate(sDOF_to_dof) DOF = _dof_to_DOF(dof,n_fdofs) - DOF_to_msDOF[DOF] = sDOF + n_mDOFs + DOF_to_msDOF[DOF] = sDOF + n_mdofs end return DOF_to_msDOF diff --git a/test/FESpacesTests/FESpacesWithLinearConstraintsTests.jl b/test/FESpacesTests/FESpacesWithLinearConstraintsTests.jl index 309f498c5..d25232cb4 100644 --- a/test/FESpacesTests/FESpacesWithLinearConstraintsTests.jl +++ b/test/FESpacesTests/FESpacesWithLinearConstraintsTests.jl @@ -101,7 +101,7 @@ Vc2 = FESpaceWithLinearConstraints( @test get_dof_value_type(Vc2) <: ComplexF64 -# Test other constructors +# 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]])