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:
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()
