Periodic BC in dolfinx

Dear all,

it’s me again :wink: I’m trying to implement periodic BCs (Poisson.-eq., of course :wink: ) in dolfinx 0.4.1.

I have a rather complicated mesh, which is generated by GMSH and imported from an xdmf-file. We can assume the mesh to be a cube with some complicated holes. I have acces to the facets and dofs of each face of the cube, see my last question.

In prior versions, it seems that SubDomain has been the way to go, but apparently it has been removed? Besides, I import my subdomain from an external mesh anyway, and further I can imagine it’s not so easy to map one side to the other, since the mesh is rather unstructered and wild (the nodes may not coincide/be not parallel)

I’ve found dolfinx_mpc, but tbh it seems a little bit overpowered for my purposes (I just have a scalar solution, no mixed space, and I really just need “vanilla” periodic BCs) and I would rather avoid the troubles of its compilation/installation on computers where I have no root access.

So, is there an easy way?

Thanks in advance.

1 Like

The issue with periodic BCs is that you have to alter the sparsity pattern and index maps underlying your FunctionSpaces, no matter how simple or complicated those may be. A general framework for enforcing periodic boundary conditions for scalable solution via MPI is not trivial.

dolfinx_mpc is very likely the best option for you; its compilation and installation should be straightforward without root access. The syntax here for the periodic Poisson problem is quite similar to old DOLFIN.

If your domain is not a simple shape, and the periodic boundaries are best defined topologically then dolfinx_mpc is absolutely the way to go as topological searches are available.

And finally: if your periodic boundary’s mesh does not conform on both “sides”, dolfinx_mpc takes care of this for you with linear combinations of the appropriate finite element bases over non-matching elements. Provided your “support” degrees of freedom are as fine or coarser than your “main” degrees of freedom, you should converge at optimal rates.

1 Like

Thank you very much. Then I will try to stick to dolfinx_mpc. Unfortuantely, the installation is not straightforward on my machine even with root access…
The docker image can’t for some reason be pulled on my machine, and the Readme says the docker image only exists for dolfinx 0.4.0. My docker-fenics has 0.4.2, however, and my conda-fenics (which I usually use) 0.4.1. And compiling on my own fails because of hdf5-libraries, which apparently can’t be found.

But this forum is probably not the right place for my problem. I’ll play around a little bit more and open an issue on github as a last resort. Or maybe @dokken can recommend the “official” way to install dolfinx_mpc (sorry for tagging :slight_smile: )?

I’ve now pushed a v0.4.1 docker image, containing both DOLFINx and DOLFINx_mpc v0.4.1.

I’ve also added a release tag: Release Release compatible with DOLFINx v0.4.1 · jorgensd/dolfinx_mpc · GitHub

Without an error message, it is hard to give any guidance on such errors. DOLFINx_mpc does not directly depend on hdf5, it only uses the python wrapper (h5py) for outputting data.

Thank you very much. The docker version seems to work!

That’s really nice, before that I went went through the git history and just checked out the last commit on May, 2nd :sweat_smile:

I was able to solve the first problem, but there is another one now. I’ll open an issue on github.

Thanks again everyone.

Ok, so I got dolfinx_mpc running using the 0.4.1-docker-container. I’ve implemented my problem, but the solutions is not really periodic… It is possible to impose mutliple periodic BCs, isn’t it?

I have a cube with a hole in it. On the surface of the hole (calles ‘active_mat’) I’ve imposed DIrichlet BCs, and I want all surfaces of the cube to be periodic. The faces face11 and face12 are parallel, so are face21 and face22, etc.

Here’s the code. If needed, I would provide all files/meshes so that it’s possible to run the code on your system. I would be very thankful for any tips.

Load 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")

Load/create tags

V = FunctionSpace(mesh, ("CG", 1))

# Basically a dict, mapping the string names to  the resp. tags
# surf_groups = {'face22': 1,  'face12': 2, 'face21': 3,
# 'face11': 4,  'face31': 5,  'face32': 6, 'active_mat': 7}
surf_groups = {key: meshiomesh.field_data[key][0] for key in meshiomesh.field_data if meshiomesh.field_data[key][1] == 2}
vol_groups = {key: meshiomesh.field_data[key][0] for key in meshiomesh.field_data if meshiomesh.field_data[key][1] == 3}

# faces and dofs
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])

# Dirichlet on the hole
bcs["active_mat"] = dirichletbc(ScalarType(2.0), surf_dofs["active_mat"], V)
bcs_list = [bcs[key] for key in bcs]

Periodic BCs

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)
# Use ft for tag ids and get the correct number from surf_groups-dict.
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()

You need to take some care if you want all sides of the problem to be periodic, as you need to avoid imposing the same constraint on the same DOF twice (which would be at all edges of your cube)…

Without having the mesh (or the equivalent gmsh file) it is hard to give you any further guidance.

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)

Hi welahi,
did you solve your problem? If yes, I would be interested in how to implement multiple periodic boundary conditions on a cube. Thanks!

Dear tuderic,

unfortunately no. I’ve implemented a version, which takes care of the problem that dokken pointed out (imposing contraints on DOFs twice), though. I basically assign each edge to only one face (but I think I did not take care of the corner points of the cube). It still doesn’t really work.

After that I didn’t have much time to work on this problem, but it’s still on my agenda. If you and others are also interested, I could create a github repo and push what I’ve coded till now and we could try to figure it out or so :slight_smile: I think that would be better than posting code snippets here.

Best,
welahi

Thank you for your fast answer. Of course, you could set up a repository and anyone who is interested can work on this problem.
In the old dolfin version 2019.1.0 it was relatively easy to construct periodic boundary conditions using SubDomain class, but - as you said - this doesn’t seem to be longer available.

Probably, I will be working on this problem in a few weeks. If I find a solution, I’ll let you know.

Oh, really? You don’t happen to have a working example of periodic cube with SubDomain-classes and fenics 2019.1.0, do you? :smiley:

I mean, dolfinx_mpc does technically work, but my solution is just not periodic. Maybe I’ve just messed up somewhere else.

See for instance my recent comment at: Nodes belonging to multiple periodic boundary conditions are improperly handled · Issue #21 · jorgensd/dolfinx_mpc · GitHub

1 Like