Periodic BC in dolfinx

mesh
8

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)