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 :slight_smile:

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

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 :slight_smile:

Screenshot from 2024-01-12 13-09-41

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:
        infile.read(mesh)
    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 :slight_smile:

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.