I tried to get the MPC strategy to work with Nedelec elements, but got into similar issues as you did (it is probably related to the facet orientation and dof orientation).
Therefore, I decided to go down another route that I have wanted to try for a long time, which is to make the mesh topology periodic, as one would no longer need specialized assembly routines.
Following is a prototype that works in serial:
# Create a periodic mesh in serial
# SPDX-License-Identifier: MIT
# Author: Jørgen S. Dokken
from mpi4py import MPI
import numpy as np
import dolfinx
import ufl
def create_periodic_mesh(mesh, indicator, mapping_function):
"""
Create a periodic mesh that takes all facets that satisfy the `indicator` function,
and map the vertices of these facets to the vertices that satisfies the mapping function.
Example:
.. code-block:: python
.. highlight:: python
mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 7, 19)
def indicator(x):
return numpy.isclose(x[1], 1)
def map(x):
values = x.copy()
values[1] -= 1
return values
periodic_mesh = create_periodic_mesh(mesh, indicator, map)
"""
geometry = mesh.geometry._cpp_object
topology = mesh.topology
mesh.topology.create_connectivity(mesh.topology.dim, mesh.topology.dim-1)
mesh.topology.create_connectivity(mesh.topology.dim-1, mesh.topology.dim)
# Map left side to right side
left_facets = dolfinx.mesh.locate_entities_boundary(mesh, mesh.topology.dim - 1, indicator)
left_vertices = dolfinx.mesh.compute_incident_entities(mesh.topology, left_facets, mesh.topology.dim -1, 0)
left_midpoints = dolfinx.mesh.compute_midpoints(mesh, 0, left_vertices)
right_midpoints = mapping_function(left_midpoints.T).T
right_vertices = dolfinx.mesh.locate_entities_boundary(mesh, 0, lambda x: x[0] > 1.0 - 1e-10)
# Find closest vertex on right side
bb_tree = dolfinx.geometry.bb_tree(mesh,0, right_vertices)
mid_tree = dolfinx.geometry.create_midpoint_tree(mesh, 0, right_vertices)
closest_vertex = dolfinx.geometry.compute_closest_entity(bb_tree, mid_tree, mesh, right_midpoints)
# Keep only left side vertices
num_vertices_local = mesh.topology.index_map(0).size_local + mesh.topology.index_map(0).num_ghosts
keep_vertices = np.ones(num_vertices_local, dtype=np.bool_)
keep_vertices[right_vertices] = False
# Create submap
new_vertices = np.flatnonzero(keep_vertices)
new_vertex_map, sub_to_parent = dolfinx.cpp.common.create_sub_index_map(mesh.topology.index_map(0), new_vertices, allow_owner_change=False)
# Invert map
num_vertices_local = mesh.topology.index_map(0).size_local + mesh.topology.index_map(0).num_ghosts
parent_to_sub = np.full(num_vertices_local, -1, dtype=np.int32)
parent_to_sub[sub_to_parent] = np.arange(sub_to_parent.size, dtype=np.int32)
# Create map from right to left vertices for replacement
replace_map = np.arange(num_vertices_local, dtype=np.int32)
replace_map[closest_vertex] = left_vertices
# First map vertices from right to left vertices
c_to_v = mesh.topology.connectivity(mesh.topology.dim, 0)
new_c = replace_map[c_to_v.array]
# Then map vertices to new indices
new_c = parent_to_sub[new_c]
new_o = c_to_v.offsets.copy()
new_c_to_v = dolfinx.graph.adjacencylist(new_c, new_o)
new_v_to_v = dolfinx.graph.adjacencylist(np.arange(len(sub_to_parent), dtype=np.int32))
topology = dolfinx.cpp.mesh.Topology(MPI.COMM_WORLD, mesh.topology.cell_type)
topology.set_index_map(0, new_vertex_map)
topology.set_index_map(mesh.topology.dim, mesh.topology.index_map(mesh.topology.dim))
topology.set_connectivity(new_v_to_v, 0,0)
topology.set_connectivity(new_c_to_v, mesh.topology.dim, 0)
c_el = dolfinx.fem.coordinate_element(mesh._ufl_domain.ufl_coordinate_element().basix_element)
geometry = dolfinx.mesh.create_geometry(mesh.geometry.index_map(), mesh.geometry.dofmap, c_el._cpp_object, mesh.geometry.x[:, :mesh.geometry.dim].copy(), mesh.geometry.input_global_indices)
if mesh.geometry.x.dtype == np.float64:
cpp_mesh = dolfinx.cpp.mesh.Mesh_float64(mesh.comm, topology, geometry._cpp_object)
elif mesh.geometry.x.dtype == np.float32:
cpp_mesh = dolfinx.cpp.mesh.Mesh_float32(mesh.comm, topology, geometry._cpp_object)
else:
raise RuntimeError(f"Unsupported dtype for mesh {mesh.geometry.x.dtype}")
new_mesh = dolfinx.mesh.Mesh(cpp_mesh, domain = ufl.Mesh(mesh._ufl_domain.ufl_coordinate_element()))
return new_mesh
mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 19, 13)
def indicator(x):
return np.isclose(x[0], 0.0)
def mapping(x):
values = x.copy()
values[0] += 1
return values
new_mesh = create_periodic_mesh(mesh, indicator, mapping)
with dolfinx.io.XDMFFile(MPI.COMM_WORLD, "periodic_mesh.xdmf", "w") as xdmf:
xdmf.write_mesh(new_mesh)
V = dolfinx.fem.functionspace(new_mesh, ("Lagrange", 2, (new_mesh.geometry.dim, )))
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
a = ufl.inner(u, v) * ufl.dx
x = ufl.SpatialCoordinate(new_mesh)
f = ufl.as_vector([5*x[0]*ufl.sin(3*np.pi * x[0])+2*x[1], x[0]*ufl.sin(3*np.pi * x[1])])
L = ufl.inner(f, v) * ufl.dx
import dolfinx.fem.petsc
problem = dolfinx.fem.petsc.LinearProblem(a, L, bcs=[], petsc_options = {"ksp_type": "preonly", "pc_type": "lu", "pc_factor_mat_solver_type": "mumps"})
uh = problem.solve()
V_out = dolfinx.fem.functionspace(new_mesh, ("DG", 2, (new_mesh.geometry.dim, )))
u_out = dolfinx.fem.Function(V_out)
u_out.interpolate(uh)
with dolfinx.io.VTXWriter(new_mesh.comm, "u_periodic.bp", [u_out]) as writer:
writer.write(0.0)
yielding:
for second order Lagrange, and with N1curl (considering tangential continuity)
