3D Periodic Boundary Conditions with Gmsh

1

I am trying to make a 3D simulation that has periodic boundaries on the x and y direction. The geometry on one side is a polygon, so I have been using gmsh to create the meshes. The issue I am running into is that fenics is not correctly applying periodic boundary conditions to the mesh. Some of the points are not correctly mapped from the slave side to the master side.

A simple working example uses a square mesh using pygmsh, a python interface to gmsh.

import dolfin as d
import pygmsh as pg
import meshio
import matplotlib.pyplot as plt
plt.rcParams['image.cmap'] = 'jet'

# Define boundary conditions
class PeriodicBoundary(d.SubDomain):

    def inside(self, x, on_boundary):
        # return True if on left or bottom boundary AND NOT on one of the two slave edges
        if on_boundary:
            if d.near(x[0], 0,eps= 1E-8) and not d.near(x[1], 1, eps= 1E-8):
                return True
            if d.near(x[1], 0,eps= 1E-8) and not d.near(x[0],1, eps= 1E-8):
                return True
            return False

    def map(self, x, y):
        if d.near(x[0],1,eps= 1E-8):
            y[0] = x[0] - 1
        else:
            y[0] = x[0]
            
        if d.near(x[1], 1,eps= 1E-8):
            y[1] = x[1] - 1
        else:
            y[1] = x[1]
        
        y[2] = x[2]




# Create a mesh using pygmsh
geom = pg.built_in.Geometry()

# Make points for a cube
pt1 = geom.add_point([0,0,0])
pt2 = geom.add_point([1,0,0])
pt3 = geom.add_point([1,1,0])
pt4 = geom.add_point([0,1,0])
pt5 = geom.add_point([0,0,1])
pt6 = geom.add_point([1,0,1])
pt7 = geom.add_point([1,1,1])
pt8 = geom.add_point([0,1,1])
# Create lines
line1 = geom.add_line(pt1,pt2)
line2 = geom.add_line(pt2,pt3)
line3 = geom.add_line(pt3,pt4)
line4 = geom.add_line(pt4,pt1)
line5 = geom.add_line(pt5,pt6)
line6 = geom.add_line(pt6,pt7)
line7 = geom.add_line(pt7,pt8)
line8 = geom.add_line(pt8,pt5)
line9 = geom.add_line(pt1,pt5)
line10 = geom.add_line(pt2,pt6)
line11 = geom.add_line(pt3,pt7)
line12 = geom.add_line(pt4,pt8)
# Create line loops
LineLoop1 = geom.add_line_loop([line1,line2,line3,line4])
LineLoop2 = geom.add_line_loop([line5,line6,line7,line8])
LineLoop3 = geom.add_line_loop([line1,line10,-line5,-line9])
LineLoop4 = geom.add_line_loop([line2,line11,-line6,-line10])
LineLoop5 = geom.add_line_loop([line3,line12,-line7,-line11])
LineLoop6 = geom.add_line_loop([line4,line9,-line8,-line12])
# Create Surfaces
surf1 = geom.add_surface(LineLoop1)
surf2 = geom.add_surface(LineLoop2)
surf3 = geom.add_surface(LineLoop3)
surf4 = geom.add_surface(LineLoop4)
surf5 = geom.add_surface(LineLoop5)
surf6 = geom.add_surface(LineLoop6)
# Create Surface Loop
surfLoop = geom.add_surface_loop([surf1,surf2,surf3,surf4,surf5,surf6])
# Create Volume
vol = geom.add_volume(surfLoop)
# Make two sides periodic
geom.add_raw_code(f"Periodic Surface {{rs2}} = {{rs4}} Translate{{0, -1, 0}};") #This defines surfXZ0Y and surfXZPosY as periodic boundaries 
geom.add_raw_code(f"Periodic Surface {{rs3}} = {{rs5}} Translate{{1, 0, 0}};") # This defines surfZYNegX and surfZYPosX as periodic boundaries
# Print gmsh code for reference
print(geom.get_code())
# Make Mesh 
mesh = pg.generate_mesh(geom)

# Save mesh in xdmf format using meshio
meshio.write('tmp/msh.xdmf', meshio.Mesh(points=mesh.points, cells={"tetra": mesh.cells["tetra"]}))
# meshio.write("tmp/mf.xdmf", meshio.Mesh(points=mesh.points, cells={"triangle": mesh.cells["triangle"]}, cell_data={"triangle": {"name_to_read": mesh.cell_data["triangle"]["gmsh:physical"]}}))

# Create problem using doflin

# Load mesh
mesh = d.Mesh()
with d.XDMFFile('tmp/msh.xdmf') as meshFile:
    meshFile.read(mesh)


# create function space and problem
pbc = PeriodicBoundary()    
V = d.FunctionSpace(mesh, "Lagrange", 1, constrained_domain=pbc)
f = d.Function(V)
f.vector()[:] = 1

# Apply a consant to the periodic boundary
pc = d.DirichletBC(V,d.Constant(-1),pbc)
pc.apply(f.vector())
d.plot(f)
plt.show()


# Using a build in mesh
mesh = d.BoxMesh(d.Point(0,0,0),d.Point(1,1,1),10,10,10)
# create function space and problem
pbc = PeriodicBoundary()    
V = d.FunctionSpace(mesh, "Lagrange", 1, constrained_domain=pbc)
f = d.Function(V)
f.vector()[:] = 1

# Apply a consant to the periodic boundary
pc = d.DirichletBC(V,d.Constant(-1),pbc)
pc.apply(f.vector())
d.plot(f)
plt.show()

Plotting the marked boundaries for both a gmsh mesh and a dolfin mesh shows that dolfin is not mapping one of the side points correctly for the gmsh mesh. When I use a larger more complicated mesh the number of points that fail to map increase. I have tried adding tolerances to the near functions with no affect.

2

I was able to figure it out. The issue is that the tolerance for mapping the slave side to the master was too low. Looking at the C++ code for the subdomain function I found the following:

/// Constructor
///
/// @param map_tol (double)
///         The tolerance used when identifying mapped points using
///         the function SubDomain::map.
SubDomain(const double map_tol=1.0e-10);

By calling the subDomain function with a larger tolerance I was able to fix the mapping issue.

pbc = PeriodicBoundary(1E-5)

1 Like