# OCC Periodic Mesh of geometry with multiple shapes on the boundary faces

Hi, Im using OCC, NGsolve and right now I’m struggling with meshing problem of the periodic geometry of patterns. I need to create periodic mesh of my geometry. The utility I’m using can prepare periodic mesh under the condition that there is only one closed shape on the boundary surface. If there are more shapes, it will make a nonperiodic mesh. Could someone please advice me, if there is any way to create such a periodic mesh in OCC?
Below is my code and an example geometry, where the only one direction (Y) has a periodic faces.
Thanks

``````def Cell50Block(nX, nY, nZ, L, r, mesh_order=2, minh=0.2, maxh=0.4):
for i in trange(nX, desc='cells generating...'):
for j in range(nY):
for k in range(nZ):
lp1x = np.array([-L/2. + L*i, L*j, L*k])
lp2x = np.array([ L/2. + L*i, L*j, L*k])
linex = line(lp1x, lp2x, r)

lp1y = np.array([L*i, -L/2 + L*j, L*k])
lp2y = np.array([L*i,  L/2 + L*j, L*k])
liney = line(lp1y, lp2y, r)

lp1z = np.array([L*i, L*j, -L/2 + L*k])
lp2z = np.array([L*i, L*j,  L/2 + L*k])
linez = line(lp1z, lp2z, r)

if i==0 and j==0 and k==0:
RVE = linex + liney + linez
else:
RVE += linex + liney + linez

RVE = RVE.Move((-L*(nX-2)+1,-L*(nY-2),-L*(nZ-2)))

RVE.faces.Min(occ.X).Identify(RVE.faces.Max(occ.X), "periodicX",
IdentificationType.PERIODIC)
RVE.faces.Min(occ.Y).Identify(RVE.faces.Max(occ.Y), "periodicY",
IdentificationType.PERIODIC)
RVE.faces.Min(occ.Z).Identify(RVE.faces.Max(occ.Z), "periodicZ",
IdentificationType.PERIODIC)

posX = occ.X >= L/2.
negX = occ.X <= -L/2.

posY = occ.Y >= L/2.
negY = occ.Y <= -L/2.

posZ = occ.Z >= L/2.
negZ = occ.Z <= -L/2.

RVE.faces[posX].name = 'pos_x'
RVE.faces[posX].col = (0, 0, 1)
RVE.faces[negX].name = 'neg_x'
RVE.faces[negX].col = (0, 0, 1)

RVE.faces[posY].name = 'pos_y'
RVE.faces[posY].col = (0, 0, 1)
RVE.faces[negY].name = 'neg_y'
RVE.faces[negY].col = (0, 0, 1)

RVE.faces[posZ].name = 'pos_z'
RVE.faces[posZ].col = (0, 0, 1)
RVE.faces[negZ].name = 'neg_z'
RVE.faces[negZ].col = (0, 0, 1)

# DrawGeo(RVE);

geo = occ.OCCGeometry(RVE)
mesh = Mesh(geo.GenerateMesh(minh=minh, maxh=maxh))
mesh.Curve(mesh_order)
return mesh
``````

1 Like

Can you provide a minimal example code that is not working as expected?

Surely, below is functional code and the area of interest is the function “Cell50Block()”. The mesh should be periodic in every direction, but it appears to be periodic only in the Y axis, bcs on the X and Z boundaries are multiple closed shapes (see picture abowe). Below is screenshot of periodicity check.

Again, thanks a lot for any advice

``````import dolfin as df
import meshio
import netgen.occ as occ
import ngsolve as ng
import numpy as np
from dolfin import near
from mpi4py import MPI
from netgen.meshing import IdentificationType
from ngsolve import Mesh
from tqdm import trange

# The indices of the full stiffness matrix of (orthorhombic) interest
Voigt_notation = [(0, 0), (1, 1), (2, 2), (1, 2), (0, 2), (0, 1)]

class ExternalCell2Dolfin(object):
def __init__(self, mesh):
self.mesh = mesh
self.centre = (self.mesh.coordinates().max(
0) + self.mesh.coordinates().min(0)) / 2
# centering
self.mesh.translate(df.Point(-self.centre))
self.xmin, self.ymin, self.zmin = self.mesh.coordinates().min(0)
self.xmax, self.ymax, self.zmax = self.mesh.coordinates().max(0)

self.lx = self.xmax - self.xmin
self.ly = self.ymax - self.ymin
self.lz = self.zmax - self.zmin

self.EPS = 1e-5
self.get_markers()

def get_left_corner(self):
EPS = self.EPS
xmin, ymin, zmin = self.xmin, self.ymin, self.zmin

class get_corner(df.SubDomain):
def inside(self, x, on_boundary):
return df.near(x[0], xmin, EPS) \
and df.near(x[1], ymin, EPS) \
and df.near(x[2], zmin, EPS) and on_boundary
return get_corner()

def get_markers(self):
EPS = self.EPS
xmin, ymin, zmin = self.xmin, self.ymin, self.zmin
xmax, ymax, zmax = self.xmax, self.ymax, self.zmax

class neg_x(df.SubDomain):
def inside(self, x, on_boundary):
return df.near(x[0], xmin, EPS) and on_boundary

class neg_y(df.SubDomain):
def inside(self, x, on_boundary):
return df.near(x[1], ymin, EPS) and on_boundary

class pos_x(df.SubDomain):
def inside(self, x, on_boundary):
return df.near(x[0], xmax, EPS) and on_boundary

class pos_y(df.SubDomain):
def inside(self, x, on_boundary):
return df.near(x[1], ymax, EPS) and on_boundary

class neg_z(df.SubDomain):
def inside(self, x, on_boundary):
return df.near(x[2], zmin, EPS) and on_boundary

class pos_z(df.SubDomain):
def inside(self, x, on_boundary):
return df.near(x[2], zmax, EPS) and on_boundary

self.boundaries = {'neg_x': neg_x(), 'pos_x': pos_x(),
'neg_y': neg_y(), 'pos_y': pos_y(),
'neg_z': neg_z(), 'pos_z': pos_z()}

def save_mesh(self, fname="reentrant_cell.pvd"):
file = df.File(fname)
file << self.mesh

def check_periodic(self):
points = self.mesh.coordinates()
pcs = [[self.xmin, self.xmax],
[self.ymin, self.ymax],
[self.zmin, self.zmax]]
nms = ['x', 'y', 'z']
for i in range(len(pcs)):
pc_min, pc_max = pcs[i]
print(f'Checking {nms[i]}-periodicity...', end=" ", flush=True)
points_min = points[(points[:, i] >= (pc_min-self.EPS))
& (points[:, i] <= (pc_min + self.EPS))]
points_max = points[(points[:, i] >= (pc_max-self.EPS))
& (points[:, i] <= (pc_max + self.EPS))]
if len(points_min) == len(points_max):
print('PASS')
else:
print('FAILED')

def line(point1, point2, r, cup=False):

lp1 = occ.Sphere(occ.Pnt(point1), r=r)
lp2 = occ.Sphere(occ.Pnt(point2), r=r)
h = np.linalg.norm(point2-point1)
direction = (point2-point1)/h
line = occ.Cylinder(occ.Pnt(point1),
occ.gp_Vec(tuple(direction)),
r=r, h=h)
if cup:
return line + lp1 + lp2
else:
return line

def Cell50Block(nX, nY, nZ, L, r, mesh_order=2, minh=0.2, maxh=0.4):
for i in trange(nX, desc='cells generating...'):
for j in range(nY):
for k in range(nZ):
lp1x = np.array([-L/2. + L*i, L*j, L*k])
lp2x = np.array([ L/2. + L*i, L*j, L*k])
linex = line(lp1x, lp2x, r)

lp1y = np.array([L*i, -L/2 + L*j, L*k])
lp2y = np.array([L*i,  L/2 + L*j, L*k])
liney = line(lp1y, lp2y, r)

lp1z = np.array([L*i, L*j, -L/2 + L*k])
lp2z = np.array([L*i, L*j,  L/2 + L*k])
linez = line(lp1z, lp2z, r)

if i==0 and j==0 and k==0:
RVE = linex + liney + linez
else:
RVE += linex + liney + linez

RVE = RVE.Move((-L*(nX-2)+1,-L*(nY-2),-L*(nZ-2)))

RVE.faces.Min(occ.X).Identify(RVE.faces.Max(occ.X), "periodicX",
IdentificationType.PERIODIC)
RVE.faces.Min(occ.Y).Identify(RVE.faces.Max(occ.Y), "periodicY",
IdentificationType.PERIODIC)
RVE.faces.Min(occ.Z).Identify(RVE.faces.Max(occ.Z), "periodicZ",
IdentificationType.PERIODIC)

posX = occ.X >= L/2.
negX = occ.X <= -L/2.
posY = occ.Y >= L/2.
negY = occ.Y <= -L/2.
posZ = occ.Z >= L/2.
negZ = occ.Z <= -L/2.

RVE.faces[posX].name = 'pos_x'
RVE.faces[posX].col = (0, 0, 1)
RVE.faces[negX].name = 'neg_x'
RVE.faces[negX].col = (0, 0, 1)

RVE.faces[posY].name = 'pos_y'
RVE.faces[posY].col = (0, 0, 1)
RVE.faces[negY].name = 'neg_y'
RVE.faces[negY].col = (0, 0, 1)

RVE.faces[posZ].name = 'pos_z'
RVE.faces[posZ].col = (0, 0, 1)
RVE.faces[negZ].name = 'neg_z'
RVE.faces[negZ].col = (0, 0, 1)

geo = occ.OCCGeometry(RVE)
mesh = Mesh(geo.GenerateMesh(minh=minh, maxh=maxh))
mesh.Curve(mesh_order)
return mesh

def netgen2dolfin(netgen_mesh):
MESH_FILE_NAME = 'mesh.xdmf'
# only linear elements
vertices = np.array([list(v.point) for v in netgen_mesh.vertices])

connectivity = []
mat_ids = []
materials = netgen_mesh.GetMaterials()
for i, el in enumerate(netgen_mesh.Elements(ng.VOL)):
connectivity.append([v.nr for v in el.vertices])
mat_ids.append(0)
for j, material in enumerate(materials):
if el.mat == material:
mat_ids[i] = j

mesh = meshio.Mesh(vertices, {"tetra": connectivity},
cell_data={'materials': [mat_ids]})
mesh.write(MESH_FILE_NAME)

mesh = df.Mesh()
with df.XDMFFile(MESH_FILE_NAME) as infile:
return mesh

if __name__ == '__main__':
mesh = Cell50Block(1, 4, 1, 1.0, 0.15, mesh_order=1, minh=0.004, maxh=0.06)
mesh = netgen2dolfin(mesh)
domain = ExternalCell2Dolfin(mesh)
domain.save_mesh('mesh.pvd')
domain.check_periodic()
``````

for multiple faces you need to also give all the faces and provide the trafo between the faces. Or you identify all pair separately.

faces.Min(X) only gives one minimal face (depending on numeric rounding errros which).
you can do

``````    minX = RVE.faces.Min(occ.X) + 1e-8
maxX = RVE.faces.Max(occ.X) - 1e-8
trf = gp_Transf.Translation(distance * occ.X)
RVE.faces.RVE.faces[occ.X < minX].Identify(RVE.faces[occ.X > maxX], "periodicX",
IdentificationType.PERIODIC, trf)
``````

to identify all X faces.

and what should be the purpose of the function gp_Transf? Thanks

It needs to be the transformation from the first set of faces to the second one. The faces are not necessarily ordered the same way, so we need to find the matching pairs by the transformation.

If there is only one face mapped to the other we automatically generate the transformation as a translation between the centers of the faces.

I was trying to implement your code and I’m having trouble importing gp_Transf function. Netgen doesn’t know this function, or did I miss something? Below is my code. Thank you again !

``````import meshio
import netgen.occ as occ
from netgen import gp_Transf
import ngsolve as ng
import dolfin as df
import numpy as np
from dolfin import near
from mpi4py import MPI
from netgen.meshing import IdentificationType
from ngsolve import Mesh
from tqdm import trange

# The indices of the full stiffness matrix of (orthorhombic) interest
Voigt_notation = [(0, 0), (1, 1), (2, 2), (1, 2), (0, 2), (0, 1)]

class ExternalCell2Dolfin(object):
def __init__(self, mesh):
self.mesh = mesh
self.centre = (self.mesh.coordinates().max(
0) + self.mesh.coordinates().min(0)) / 2
# centering
self.mesh.translate(df.Point(-self.centre))
self.xmin, self.ymin, self.zmin = self.mesh.coordinates().min(0)
self.xmax, self.ymax, self.zmax = self.mesh.coordinates().max(0)

self.lx = self.xmax - self.xmin
self.ly = self.ymax - self.ymin
self.lz = self.zmax - self.zmin

self.EPS = 1e-5
self.get_markers()

def get_left_corner(self):
EPS = self.EPS
xmin, ymin, zmin = self.xmin, self.ymin, self.zmin

class get_corner(df.SubDomain):
def inside(self, x, on_boundary):
return df.near(x[0], xmin, EPS) \
and df.near(x[1], ymin, EPS) \
and df.near(x[2], zmin, EPS) and on_boundary
return get_corner()

def get_markers(self):
EPS = self.EPS
xmin, ymin, zmin = self.xmin, self.ymin, self.zmin
xmax, ymax, zmax = self.xmax, self.ymax, self.zmax

class neg_x(df.SubDomain):
def inside(self, x, on_boundary):
return df.near(x[0], xmin, EPS) and on_boundary

class neg_y(df.SubDomain):
def inside(self, x, on_boundary):
return df.near(x[1], ymin, EPS) and on_boundary

class pos_x(df.SubDomain):
def inside(self, x, on_boundary):
return df.near(x[0], xmax, EPS) and on_boundary

class pos_y(df.SubDomain):
def inside(self, x, on_boundary):
return df.near(x[1], ymax, EPS) and on_boundary

class neg_z(df.SubDomain):
def inside(self, x, on_boundary):
return df.near(x[2], zmin, EPS) and on_boundary

class pos_z(df.SubDomain):
def inside(self, x, on_boundary):
return df.near(x[2], zmax, EPS) and on_boundary

self.boundaries = {'neg_x': neg_x(), 'pos_x': pos_x(),
'neg_y': neg_y(), 'pos_y': pos_y(),
'neg_z': neg_z(), 'pos_z': pos_z()}

def save_mesh(self, fname="reentrant_cell.pvd"):
file = df.File(fname)
file << self.mesh

def check_periodic(self):
points = self.mesh.coordinates()
pcs = [[self.xmin, self.xmax],
[self.ymin, self.ymax],
[self.zmin, self.zmax]]
nms = ['x', 'y', 'z']
for i in range(len(pcs)):
pc_min, pc_max = pcs[i]
print(f'Checking {nms[i]}-periodicity...', end=" ", flush=True)
points_min = points[(points[:, i] >= (pc_min-self.EPS))
& (points[:, i] <= (pc_min + self.EPS))]
points_max = points[(points[:, i] >= (pc_max-self.EPS))
& (points[:, i] <= (pc_max + self.EPS))]
if len(points_min) == len(points_max):
print('PASS')
else:
print('FAILED')

def line(point1, point2, r, cup=False):

lp1 = occ.Sphere(occ.Pnt(point1), r=r)
lp2 = occ.Sphere(occ.Pnt(point2), r=r)
h = np.linalg.norm(point2-point1)
direction = (point2-point1)/h
line = occ.Cylinder(occ.Pnt(point1),
occ.gp_Vec(tuple(direction)),
r=r, h=h)
if cup:
return line + lp1 + lp2
else:
return line

def Cell50Block(nX, nY, nZ, L, r, mesh_order=2, minh=0.2, maxh=0.4):
for i in trange(nX, desc='cells generating...'):
for j in range(nY):
for k in range(nZ):
lp1x = np.array([-L/2. + L*i, L*j, L*k])
lp2x = np.array([ L/2. + L*i, L*j, L*k])
linex = line(lp1x, lp2x, r)

lp1y = np.array([L*i, -L/2 + L*j, L*k])
lp2y = np.array([L*i,  L/2 + L*j, L*k])
liney = line(lp1y, lp2y, r)

lp1z = np.array([L*i, L*j, -L/2 + L*k])
lp2z = np.array([L*i, L*j,  L/2 + L*k])
linez = line(lp1z, lp2z, r)

if i==0 and j==0 and k==0:
RVE = linex + liney + linez
else:
RVE += linex + liney + linez

RVE = RVE.Move((-L*(nX-2)+1,-L*(nY-2),-L*(nZ-2)))
distance = L
epsx = occ.gp_Vec((1e-8, 0, 0))
epsy = occ.gp_Vec((0, 1e-8, 0))
epsz = occ.gp_Vec((0, 0, 1e-8))

minX = RVE.faces.Min(occ.X + epsx)
trf = occ.gp_Transf.Translation(-distance/2 * occ.X)
RVE.faces.RVE.faces[occ.X < minX].Identify(RVE.faces[occ.X > maxX], "periodicX",
IdentificationType.PERIODIC, trf)

minY = RVE.faces.Min(occ.Y + epsy)
maxY = RVE.faces.Max(occ.Y - epsy)
trf = occ.gp_Transf.Translation(-distance/2 * occ.Y)
RVE.faces.RVE.faces[occ.Y < minY].Identify(RVE.faces[occ.Y > maxY], "periodicY",
IdentificationType.PERIODIC, trf)

minZ = RVE.faces.Min(occ.Z + epsz)
maxZ = RVE.faces.Max(occ.Z - epsz)
trf = occ.gp_Transf.Translation(-distance/2 * occ.Z)
RVE.faces.RVE.faces[occ.Z < minZ].Identify(RVE.faces[occ.Z > maxZ], "periodicZ",
IdentificationType.PERIODIC, trf)

# RVE.faces.Min(occ.X).Identify(RVE.faces.Max(occ.X), "periodicX",
#                               IdentificationType.PERIODIC)
# RVE.faces.Min(occ.Y).Identify(RVE.faces.Max(occ.Y), "periodicY",
#                               IdentificationType.PERIODIC)
# RVE.faces.Min(occ.Z).Identify(RVE.faces.Max(occ.Z), "periodicZ",
#                               IdentificationType.PERIODIC)

posX = occ.X >= L/2.
negX = occ.X <= -L/2.
posY = occ.Y >= L/2.
negY = occ.Y <= -L/2.
posZ = occ.Z >= L/2.
negZ = occ.Z <= -L/2.

RVE.faces[posX].name = 'pos_x'
RVE.faces[posX].col = (0, 0, 1)
RVE.faces[negX].name = 'neg_x'
RVE.faces[negX].col = (0, 0, 1)

RVE.faces[posY].name = 'pos_y'
RVE.faces[posY].col = (0, 0, 1)
RVE.faces[negY].name = 'neg_y'
RVE.faces[negY].col = (0, 0, 1)

RVE.faces[posZ].name = 'pos_z'
RVE.faces[posZ].col = (0, 0, 1)
RVE.faces[negZ].name = 'neg_z'
RVE.faces[negZ].col = (0, 0, 1)

geo = occ.OCCGeometry(RVE)
mesh = Mesh(geo.GenerateMesh(minh=minh, maxh=maxh))
mesh.Curve(mesh_order)
return mesh

def netgen2dolfin(netgen_mesh):
MESH_FILE_NAME = 'mesh.xdmf'
# only linear elements
vertices = np.array([list(v.point) for v in netgen_mesh.vertices])

connectivity = []
mat_ids = []
materials = netgen_mesh.GetMaterials()
for i, el in enumerate(netgen_mesh.Elements(ng.VOL)):
connectivity.append([v.nr for v in el.vertices])
mat_ids.append(0)
for j, material in enumerate(materials):
if el.mat == material:
mat_ids[i] = j

mesh = meshio.Mesh(vertices, {"tetra": connectivity},
cell_data={'materials': [mat_ids]})
mesh.write(MESH_FILE_NAME)

mesh = df.Mesh()
with df.XDMFFile(MESH_FILE_NAME) as infile:
return mesh

if __name__ == '__main__':
mesh = Cell50Block(1, 4, 1, 1.0, 0.15, mesh_order=1, minh=0.004, maxh=0.06)
mesh = netgen2dolfin(mesh)
domain = ExternalCell2Dolfin(mesh)
domain.save_mesh('mesh.pvd')
domain.check_periodic()
``````

the function should be in `netgen.occ`

Best Christopher

It seems that the function is gp_Trsf form netgen.occ. I adapted your code a little bit, as follows

``````    epsx = occ.gp_Vec((1e-8, 0, 0))
epsy = occ.gp_Vec((0, 1e-8, 0))
epsz = occ.gp_Vec((0, 0, 1e-8))

minX = RVE.faces.Min(occ.X + epsx)
maxX = RVE.faces.Max(occ.X - epsx)
trf = occ.gp_Trsf.Translation(-distance/2 * occ.X)
RVE.faces.RVE.faces[occ.X < minX].Identify(RVE.faces[occ.X > maxX], "periodicX",
IdentificationType.PERIODIC, trf)
``````

The netget.occ can’t work with defining faceses like in this code `RVE.faces.RVE.faces[occ.x < minX]`… Is there any way to fix this line? Or it can be fixed by for-loop for example? Below is shown error message. Thank you

``````--> 145 RVE.faces.RVE.faces[occ.X < minX].Identify(RVE.faces[occ.X > maxX], "periodicX",
146                               IdentificationType.PERIODIC, trf)

AttributeError: 'netgen.libngpy._NgOCC.ListOfShapes' object has no attribute 'RVE'
``````

ah yes sorry typo. but you having RVE.faces twice is also a typo or?

oh, i see.

``````    RVE.faces[occ.X < minX].Identify(RVE.faces[occ.X > maxX], "periodicX",
IdentificationType.PERIODIC, trf)
``````

But there is still problem with arguments as follows:

``````--> 145 RVE.faces[occ.X < minX].Identify(RVE.faces[occ.X > maxX], "periodicX",
146                               IdentificationType.PERIODIC, trf)

TypeError: __lt__(): incompatible function arguments. The following argument types are supported:
1. (self: netgen.libngpy._NgOCC.gp_Vec, arg0: float) -> netgen::DirectionalInterval

Invoked with: (1, 0, 0), <netgen.libngpy._NgOCC.Face object at 0x7377c2e32370>
``````

in `occ.X<minX` the `minX` should be a number, but in your case it is a shape

Thank you for your prompt and significant support. I really appreciate it! I finally figure it out! Below is my code, which works for me.

``````    minX = -distance/2 + 1e-3
maxX = distance/2 - 1e-3
trf = occ.gp_Trsf.Translation(distance * occ.X)
RVE.faces[occ.X < minX].Identify(RVE.faces[occ.X > maxX], "periodicX",
IdentificationType.PERIODIC, trf)
"--> likewise for other directions"
``````

It connects different geometries in the right way, but still has a problem with the following geometry:

The code above can create mesh with periodic faces only for 2 directions. In this case it’s not able to create X direction periodic faces. Below are the face mashes of interest in the overlay for X (not periodic) and Y (periodic) direction.

Do you have any idea what’s wrong? Again, thank you very much for your help.

To reference joachim:

I think this can be a bug, but if you post a full new minimal example that we can run, we have a change to check it

Yes, sorry. Here is the new executable code:

``````import meshio
import netgen.occ as occ
from netgen.occ import gp_Trsf
import ngsolve as ng
import dolfin as df
import numpy as np
from dolfin import near
from mpi4py import MPI
from netgen.meshing import IdentificationType
from ngsolve import Mesh
from tqdm import trange

# The indices of the full stiffness matrix of (orthorhombic) interest
Voigt_notation = [(0, 0), (1, 1), (2, 2), (1, 2), (0, 2), (0, 1)]

class ExternalCell2Dolfin(object):
def __init__(self, mesh):
self.mesh = mesh
self.centre = (self.mesh.coordinates().max(
0) + self.mesh.coordinates().min(0)) / 2
# centering
self.mesh.translate(df.Point(-self.centre))
self.xmin, self.ymin, self.zmin = self.mesh.coordinates().min(0)
self.xmax, self.ymax, self.zmax = self.mesh.coordinates().max(0)

self.lx = self.xmax - self.xmin
self.ly = self.ymax - self.ymin
self.lz = self.zmax - self.zmin

self.EPS = 1e-5
self.get_markers()

def get_left_corner(self):
EPS = self.EPS
xmin, ymin, zmin = self.xmin, self.ymin, self.zmin

class get_corner(df.SubDomain):
def inside(self, x, on_boundary):
return df.near(x[0], xmin, EPS) \
and df.near(x[1], ymin, EPS) \
and df.near(x[2], zmin, EPS) and on_boundary
return get_corner()

def get_markers(self):
EPS = self.EPS
xmin, ymin, zmin = self.xmin, self.ymin, self.zmin
xmax, ymax, zmax = self.xmax, self.ymax, self.zmax

class neg_x(df.SubDomain):
def inside(self, x, on_boundary):
return df.near(x[0], xmin, EPS) and on_boundary

class neg_y(df.SubDomain):
def inside(self, x, on_boundary):
return df.near(x[1], ymin, EPS) and on_boundary

class pos_x(df.SubDomain):
def inside(self, x, on_boundary):
return df.near(x[0], xmax, EPS) and on_boundary

class pos_y(df.SubDomain):
def inside(self, x, on_boundary):
return df.near(x[1], ymax, EPS) and on_boundary

class neg_z(df.SubDomain):
def inside(self, x, on_boundary):
return df.near(x[2], zmin, EPS) and on_boundary

class pos_z(df.SubDomain):
def inside(self, x, on_boundary):
return df.near(x[2], zmax, EPS) and on_boundary

self.boundaries = {'neg_x': neg_x(), 'pos_x': pos_x(),
'neg_y': neg_y(), 'pos_y': pos_y(),
'neg_z': neg_z(), 'pos_z': pos_z()}

def save_mesh(self, fname="reentrant_cell.pvd"):
file = df.File(fname)
file << self.mesh

def check_periodic(self):
points = self.mesh.coordinates()
pcs = [[self.xmin, self.xmax],
[self.ymin, self.ymax],
[self.zmin, self.zmax]]
nms = ['x', 'y', 'z']
for i in range(len(pcs)):
pc_min, pc_max = pcs[i]
print(f'Checking {nms[i]}-periodicity...', end=" ", flush=True)
points_min = points[(points[:, i] >= (pc_min-self.EPS))
& (points[:, i] <= (pc_min + self.EPS))]
points_max = points[(points[:, i] >= (pc_max-self.EPS))
& (points[:, i] <= (pc_max + self.EPS))]
if len(points_min) == len(points_max):
print('PASS')
else:
print('FAILED')

def line(point1, point2, r, cup=False):

lp1 = occ.Sphere(occ.Pnt(point1), r=r)
lp2 = occ.Sphere(occ.Pnt(point2), r=r)
h = np.linalg.norm(point2-point1)
direction = (point2-point1)/h
line = occ.Cylinder(occ.Pnt(point1),
occ.gp_Vec(tuple(direction)),
r=r, h=h)
if cup:
return line + lp1 + lp2
else:
return line

def Cell21X(L, r, mesh_order=2, minh=0.2, maxh=0.4):
lp1 = np.array([- L/2., - L/2., - L/2.])
lp2 = np.array([- L/2., + L/2., - L/2.])
lp3 = np.array([+ L/2., + L/2., - L/2.])
lp3out = np.array([+ L/2. + r, + L/2. + r, - L/2. - r])
lp4 = np.array([+ L/2., - L/2., - L/2.])
lp5 = np.array([- L/2., - L/2., + L/2.])
lp5out = np.array([- L/2. - r, - L/2. - r, + L/2. + r])
lp6 = np.array([- L/2., + L/2., + L/2.])
lp7 = np.array([+ L/2., + L/2., + L/2.])
lp8 = np.array([+ L/2., - L/2., + L/2.])
line1 = line(lp1, lp7, r)
line2 = line(lp2, lp8, r)
line3 = line(lp3, lp5, r)
line4 = line(lp4, lp6, r)

RVE = line1 + line2 + line3 + line4
lp_outer = occ.Box(occ.Pnt(lp3out), occ.Pnt(lp5out))
lp_inner = occ.Box(occ.Pnt(lp3), occ.Pnt(lp5))
GG = lp_outer-lp_inner
RVE -= GG

minX = -L/2 + 1e3
maxX = L/2 - 1e-3
trf = occ.gp_Trsf.Translation(L * occ.X)
RVE.faces[occ.X < minX].Identify(RVE.faces[occ.X > maxX], "periodicX",
IdentificationType.PERIODIC, trf)

minY = -L/2 + 1e-3
maxY = L/2 - 1e-3
trf = occ.gp_Trsf.Translation(L * occ.Y)
RVE.faces[occ.Y < minY].Identify(RVE.faces[occ.Y > maxY], "periodicY",
IdentificationType.PERIODIC, trf)

minZ = -L/2 + 1e-3
maxZ = L/2 - 1e-3
trf = occ.gp_Trsf.Translation(L * occ.Z)
RVE.faces[occ.Z < minZ].Identify(RVE.faces[occ.Z > maxZ], "periodicZ",
IdentificationType.PERIODIC, trf)

geo = occ.OCCGeometry(RVE)
mesh = Mesh(geo.GenerateMesh(minh=minh, maxh=maxh))
mesh.Curve(mesh_order)
return mesh

def netgen2dolfin(netgen_mesh):
MESH_FILE_NAME = 'mesh.xdmf'
# only linear elements
vertices = np.array([list(v.point) for v in netgen_mesh.vertices])

connectivity = []
mat_ids = []
materials = netgen_mesh.GetMaterials()
for i, el in enumerate(netgen_mesh.Elements(ng.VOL)):
connectivity.append([v.nr for v in el.vertices])
mat_ids.append(0)
for j, material in enumerate(materials):
if el.mat == material:
mat_ids[i] = j

mesh = meshio.Mesh(vertices, {"tetra": connectivity},
cell_data={'materials': [mat_ids]})
mesh.write(MESH_FILE_NAME)

mesh = df.Mesh()
with df.XDMFFile(MESH_FILE_NAME) as infile: