This is indeed a good point. I’ve checked the surf_dofs for the different faces, and they are in fact not disjoint (which they should be, as far as I understood). What would be the way to go? I fear that I need to tag somehow all the edges and then assign the dofs each egde to one face uniquely…
Below is a minimal working example. Putting all of them in one directory, generating the msh-file with GMSH, running conv_mesh.py and solve.py should work within the dolfinx_mpc-0.4.1-docker.
Thanks again.
mymesh.geo
// Gmsh project created on Mon Jun 20 19:24:08 2022
SetFactory("OpenCASCADE");
//+
Box(1) = {-0.5, -0.5, -0.5, 1, 1, 1};
//+
Sphere(2) = {0, 0, 0, 0.3, -Pi/2, Pi/2, 2*Pi};
//+
Physical Surface("face22") = {4};
//+
Physical Surface("face12") = {6};
//+
Physical Surface("face21") = {2};
//+
Physical Surface("face11") = {1};
//+
Physical Surface("face31") = {5};
//+
Physical Surface("face32") = {3};
//+
Physical Surface("active_mat") = {7};
//+
BooleanDifference(99) = { Volume{1}; Delete; }{ Volume{2}; Delete; };
//+
Physical Volume("vol") = {99};
solve.py
import numpy as np
import meshio
import dolfinx
from dolfinx.fem import (Constant, dirichletbc, FunctionSpace, locate_dofs_topological)
from dolfinx.io import XDMFFile
from ufl import (TestFunction, TrialFunction, dx, grad, inner)
import dolfinx_mpc
from dolfinx_mpc import LinearProblem
from mpi4py import MPI
from petsc4py.PETSc import ScalarType
# Read meshes
with XDMFFile(MPI.COMM_WORLD, "mymesh_tetra.xdmf", "r") as xdmf:
mesh = xdmf.read_mesh(name="Grid")
ct = xdmf.read_meshtags(mesh, name="Grid")
mesh.topology.create_connectivity(mesh.topology.dim-1, 0)
with XDMFFile(MPI.COMM_WORLD, "mymesh_triangle.xdmf", "r") as xdmf:
ft = xdmf.read_meshtags(mesh, name="Grid")
meshiomesh = meshio.read("mymesh.msh")
# Topology stuff / Dirichlet BC
V = FunctionSpace(mesh, ("CG", 1))
surf_groups = {key: meshiomesh.field_data[key][0] for key in meshiomesh.field_data if meshiomesh.field_data[key][1] == 2}
surf_faces = {}
surf_dofs = {}
bcs = {}
for key in surf_groups:
surf_faces[key] = ft.indices[ft.values==surf_groups[key]]
surf_dofs[key] = locate_dofs_topological(V, mesh.topology.dim-1, surf_faces[key])
bcs["active_mat"] = dirichletbc(ScalarType(2.0), surf_dofs["active_mat"], V)
bcs_list = [bcs[key] for key in bcs]
# Periodic BC
def periodic_relation_face1i(x):
out_x = np.zeros(x.shape)
out_x[0] = x[0] - 1
out_x[1] = x[1]
out_x[2] = x[2]
return out_x
def periodic_relation_face2i(x):
out_x = np.zeros(x.shape)
out_x[0] = x[0]
out_x[1] = x[1] - 1
out_x[2] = x[2]
return out_x
def periodic_relation_face3i(x):
out_x = np.zeros(x.shape)
out_x[0] = x[0]
out_x[1] = x[1]
out_x[2] = x[2] - 1
return out_x
mpc = dolfinx_mpc.MultiPointConstraint(V)
mpc.create_periodic_constraint_topological(V, ft, surf_groups["face12"], periodic_relation_face1i, bcs_list, 1)
mpc.create_periodic_constraint_topological(V, ft, surf_groups["face22"], periodic_relation_face2i, bcs_list, 1)
mpc.create_periodic_constraint_topological(V, ft, surf_groups["face32"], periodic_relation_face3i, bcs_list, 1)
mpc.finalize()
# Solve
u, v = TrialFunction(V), TestFunction(V)
a = inner(grad(u), grad(v)) * dx
rhs = Constant(mesh, ScalarType(1)) * v * dx
problem = LinearProblem(a, rhs, mpc, bcs=bcs_list, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh = problem.solve()
# Save
with dolfinx.io.XDMFFile(MPI.COMM_WORLD, "output.xdmf", "w") as xdmf:
xdmf.write_mesh(mesh)
xdmf.write_function(uh)
xdmf.close()
conv_mesh.py
import meshio
def create_mesh(mesh, cell_type, prune_z=False):
cells = mesh.get_cells_type(cell_type)
cell_data = mesh.get_cell_data("gmsh:physical", cell_type)
points = mesh.points[:,:2] if prune_z else mesh.points
out_mesh = meshio.Mesh(points=points, cells={cell_type: cells}, cell_data={"name_to_read":[cell_data]})
return out_mesh
mesh = meshio.read("mymesh.msh")
newmesh1 = create_mesh(mesh, "triangle")
meshio.write("mymesh_triangle.xdmf", newmesh1)
newmesh2 = create_mesh(mesh, "tetra")
meshio.write("mymesh_tetra.xdmf", newmesh2)
I’m trying to implement periodic BCs (Poisson.-eq., of course 
)?
