I think you are misunderstanding what should be modified in your code.
Here you are only changing the index map, but not the dofmap.
You need to change both in a consistent way, which is quite some work.
Here is a code that creates a periodic BC for a 1D problem (with a single dof being replaced in the dofmap)
# A DOFMAP reduction code for periodic boundary conditions in 1D
# SPDX-License-Identifier: MIT
# Author: Jørgen S. Dokken
from petsc4py import PETSc
from mpi4py import MPI
import ufl
from dolfinx import fem, mesh, cpp, io
from dolfinx.fem.petsc import LinearProblem
from ufl import dx
import basix
import numpy as np
msh = mesh.create_unit_interval(
MPI.COMM_WORLD, 10, ghost_mode=mesh.GhostMode.none)
V = fem.functionspace(msh, ("Lagrange", 1))
def left(x):
return np.isclose(x[0], 0.0)
def right(x):
return np.isclose(x[0], 1.0)
puppet = fem.locate_dofs_geometrical(V, left)
master = fem.locate_dofs_geometrical(V, right)
assert len(puppet) <= 2
assert len(master) <= 2
# Currently assuming single puppet
V_index_map = V.dofmap.index_map
has_master = np.full(msh.comm.size, -1, dtype=np.int64)
has_puppet = np.zeros(msh.comm.size, dtype=np.int32)
owned_puppet = False
if len(puppet) > 0 and puppet[0] < V_index_map.size_local:
has_puppet[:] = 1
owned_puppet = True
# Send master to do slave process
owned_master = False
if len(master) > 0 and master[0] < V_index_map.size_local:
has_master[:] = V_index_map.local_range[0] + master[0]
owned_master = True
puppet_to_master = np.zeros_like(has_puppet)
master_to_puppet = np.zeros_like(has_master)
msh.comm.Alltoall(has_puppet, puppet_to_master)
msh.comm.Alltoall(has_master, master_to_puppet)
puppet_ranks = np.flatnonzero(puppet_to_master)
master_indicator = (master_to_puppet > -1).astype(np.int32)
master_ranks = np.flatnonzero(master_indicator)
# If has puppet, this process should communicate with the master
V_index_map = V.dofmap.index_map
i_to_dest = V_index_map.index_to_dest_ranks()
src = np.unique(i_to_dest.array[:i_to_dest.offsets[V_index_map.size_local]])
dest = np.unique(V_index_map.owners)
local_range = V_index_map.local_range
local_size = V_index_map.size_local
# Map ghosts to a lower index to account for the removed dof
old_ghosts = V.dofmap.index_map.ghosts
new_ghosts = []
new_owners = []
for owner, ghost in zip(V_index_map.owners, old_ghosts):
new_ghosts.append(ghost - sum(puppet_to_master[:owner]))
new_owners.append(owner)
if not owned_master:
assert not owned_master
if (gd := (master_to_puppet[master_ranks[0]] - sum(puppet_to_master[:master_ranks[0]]))) not in new_ghosts:
ghost_pos = len(new_ghosts)
new_ghosts.append(gd)
new_owners.append(master_ranks[0])
else:
ghost_pos = np.where(np.asarray(new_ghosts) == gd)[0][0]
if owned_puppet:
local_size -= 1
new_index_map = cpp.common.IndexMap(
MPI.COMM_WORLD, local_size, np.asarray(new_ghosts, dtype=np.int64), np.asarray(new_owners, dtype=np.int32))
old_dofs = V.dofmap.list.copy()
shape = old_dofs.shape
new_dofs = old_dofs.reshape(-1).copy()
pupp = -1
if len(puppet) > 0:
pupp = puppet[0]
for i, dof in enumerate(old_dofs.reshape(-1)):
if dof > pupp:
new_dofs[i] = dof - int(pupp > -1)
elif dof == pupp:
if owned_master:
new_dofs[i] = master[0] - int(master[0] > pupp)*int(pupp > -1)
else:
new_dofs[i] = local_size + ghost_pos
dofs = cpp.graph.AdjacencyList_int32(new_dofs.reshape(shape))
new_dofmap = cpp.fem.DofMap(V.dofmap.dof_layout,
new_index_map, V.dofmap.index_map_bs,
dofs, V.dofmap.bs)
# Creat the functionspace
cppV = cpp.fem.FunctionSpace_float64(
msh._cpp_object, V.element, new_dofmap, V.value_shape)
e = fem.ElementMetaData("Lagrange", 1)
ufl_e = basix.ufl.element(e.family, msh.basix_cell(), e.degree, shape=e.shape,
symmetry=e.symmetry)
V_pbc = fem.FunctionSpace(msh, ufl_e, cppV)
def solve_problem(V):
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
x = ufl.SpatialCoordinate(msh)
f = x[0]*ufl.sin(x[0])
a = ufl.dot(ufl.grad(u), ufl.grad(v)) * dx
L = f * v * dx
problem = LinearProblem(a, L, bcs=[], petsc_options={
"ksp_type": "preonly", "pc_type": "lu",
"pc_factor_mat_solver_type": "mumps",
"mat_mumps_icntl_24": 1, # Option to support solving a singular matrix
"mat_mumps_icntl_25": 0 # Option to support solving a singular matrix
})
nullspace = PETSc.NullSpace().create(constant=True) # type: ignore
PETSc.Mat.setNearNullSpace(problem.A, nullspace)
uh = problem.solve()
print(problem.A[:, :])
print(problem.b[:])
with io.XDMFFile(msh.comm, "test.xdmf", "w") as xdmf:
xdmf.write_mesh(msh)
xdmf.write_function(uh, 0.0)
solve_problem(V)
solve_problem(V_pbc)
yielding (non-periodic)
[[ 10. -10. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
[-10. 20. -10. 0. 0. 0. 0. 0. 0. 0. 0.]
[ 0. -10. 20. -10. 0. 0. 0. 0. 0. 0. 0.]
[ 0. 0. -10. 20. -10. 0. 0. 0. 0. 0. 0.]
[ 0. 0. 0. -10. 20. -10. 0. 0. 0. 0. 0.]
[ 0. 0. 0. 0. -10. 20. -10. 0. 0. 0. 0.]
[ 0. 0. 0. 0. 0. -10. 20. -10. 0. 0. 0.]
[ 0. 0. 0. 0. 0. 0. -10. 20. -10. 0. 0.]
[ 0. 0. 0. 0. 0. 0. 0. -10. 20. -10. 0.]
[ 0. 0. 0. 0. 0. 0. 0. 0. -10. 20. -10.]
[ 0. 0. 0. 0. 0. 0. 0. 0. 0. -10. 10.]]
[8.32777934e-05 1.16322600e-03 4.13331213e-03 9.01733729e-03
1.57171653e-02 2.40974738e-02 3.39877899e-02 4.51850600e-02
5.74567199e-02 7.05442245e-02 3.97830923e-02]
Periodic
[[ 20. -10. 0. 0. 0. 0. 0. 0. 0. -10.]
[-10. 20. -10. 0. 0. 0. 0. 0. 0. 0.]
[ 0. -10. 20. -10. 0. 0. 0. 0. 0. 0.]
[ 0. 0. -10. 20. -10. 0. 0. 0. 0. 0.]
[ 0. 0. 0. -10. 20. -10. 0. 0. 0. 0.]
[ 0. 0. 0. 0. -10. 20. -10. 0. 0. 0.]
[ 0. 0. 0. 0. 0. -10. 20. -10. 0. 0.]
[ 0. 0. 0. 0. 0. 0. -10. 20. -10. 0.]
[ 0. 0. 0. 0. 0. 0. 0. -10. 20. -10.]
[-10. 0. 0. 0. 0. 0. 0. 0. -10. 20.]]
[0.00116323 0.00413331 0.00901734 0.01571717 0.02409747 0.03398779
0.04518506 0.05745672 0.07054422 0.03986637]