Here is the minimal working example:
import gmsh
import os
import math
import meshio
import dolfin
###################################################################################################################################################################
fname = 'cube'
width = 1
r = width/2.5
shift_x = width/2
shift_y = width/2
shift_z = width/2
center = width/2
xmin = 0.
ymin = 0.
zmin = 0.
xmax = width
ymax = width
zmax = width
x0 = center
y0 = center
z0 = center
r0 = r
l = width/10
e = 1e-6
def setPeriodic(coord):
# From https://gitlab.onelab.info/gmsh/gmsh/-/issues/744
smin = gmsh.model.getEntitiesInBoundingBox(xmin - e, ymin - e, zmin - e,
(xmin + e) if (coord == 0) else (xmax + e),
(ymin + e) if (coord == 1) else (ymax + e),
(zmin + e) if (coord == 2) else (zmax + e),
2)
dx = (xmax - xmin) if (coord == 0) else 0
dy = (ymax - ymin) if (coord == 1) else 0
dz = (zmax - zmin) if (coord == 2) else 0
for i in smin:
bb = gmsh.model.getBoundingBox(i[0], i[1])
bbe = [bb[0] - e + dx, bb[1] - e + dy, bb[2] - e + dz,
bb[3] + e + dx, bb[4] + e + dy, bb[5] + e + dz]
smax = gmsh.model.getEntitiesInBoundingBox(bbe[0], bbe[1], bbe[2],
bbe[3], bbe[4], bbe[5])
for j in smax:
bb2 = list(gmsh.model.getBoundingBox(j[0], j[1]))
bb2[0] -= dx; bb2[1] -= dy; bb2[2] -= dz;
bb2[3] -= dx; bb2[4] -= dy; bb2[5] -= dz;
if ((abs(bb2[0] - bb[0]) < e) and (abs(bb2[1] - bb[1]) < e) and
(abs(bb2[2] - bb[2]) < e) and (abs(bb2[3] - bb[3]) < e) and
(abs(bb2[4] - bb[4]) < e) and (abs(bb2[5] - bb[5]) < e)):
gmsh.model.mesh.setPeriodic(2, [j[1]], [i[1]], [1, 0, 0, dx,\
0, 1, 0, dy,\
0, 0, 1, dz,\
0, 0, 0, 1 ])
gmsh.initialize()
box_tag = 1
hole_tag1 = 2
rve_tag = 3
gmsh.model.occ.addBox(x=xmin, y=ymin, z=zmin, dx=xmax-xmin, dy=ymax-ymin, dz=zmax-zmin, tag=box_tag)
gmsh.model.occ.addSphere(xc=xmin+shift_x, yc=ymin+shift_y, zc=zmin+shift_z, radius=r0, tag=hole_tag1)
gmsh.model.occ.cut(objectDimTags=[(3, box_tag)], toolDimTags=[(3, hole_tag1)], tag=rve_tag)
gmsh.model.occ.synchronize()
gmsh.model.addPhysicalGroup(dim=3, tags=[rve_tag])
setPeriodic(0)
setPeriodic(1)
setPeriodic(2)
gmsh.model.mesh.setSize(dimTags=gmsh.model.getEntities(0), size=l)
gmsh.model.mesh.generate(dim=3)
gmsh.write(fname+"-mesh.vtk")
gmsh.finalize()
mesh = meshio.read(fname+"-mesh.vtk")
mesh.points = mesh.points[:, :3]
meshio.write(fname+"-mesh.xdmf", mesh)
mesh = dolfin.Mesh()
dolfin.XDMFFile(fname+"-mesh.xdmf").read(mesh)
vol = width**3
phi = (4/3*math.pi*r**3)/vol
print("porosity:" +str(phi))
#####################################
tol = 1E-8
xmin_sd = dolfin.CompiledSubDomain("near(x[0], x0, tol) && on_boundary", x0=xmin, tol=tol)
xmax_sd = dolfin.CompiledSubDomain("near(x[0], x0, tol) && on_boundary", x0=xmax, tol=tol)
ymin_sd = dolfin.CompiledSubDomain("near(x[1], x0, tol) && on_boundary", x0=ymin, tol=tol)
ymax_sd = dolfin.CompiledSubDomain("near(x[1], x0, tol) && on_boundary", x0=ymax, tol=tol)
zmin_sd = dolfin.CompiledSubDomain("near(x[2], x0, tol) && on_boundary", x0=zmin, tol=tol)
zmax_sd = dolfin.CompiledSubDomain("near(x[2], x0, tol) && on_boundary", x0=zmax, tol=tol)
xmin_id = 1
xmax_id = 2
ymin_id = 3
ymax_id = 4
zmin_id = 5
zmax_id = 6
boundaries_mf = dolfin.MeshFunction("size_t", mesh, mesh.topology().dim()-1)
boundaries_mf.set_all(0)
xmin_sd.mark(boundaries_mf, xmin_id)
xmax_sd.mark(boundaries_mf, xmax_id)
ymin_sd.mark(boundaries_mf, ymin_id)
ymax_sd.mark(boundaries_mf, ymax_id)
zmin_sd.mark(boundaries_mf, zmin_id)
zmax_sd.mark(boundaries_mf, zmax_id)
xdmf_file_boundaries = dolfin.XDMFFile(fname+"-boundaries.xdmf")
xdmf_file_boundaries.write(boundaries_mf)
xdmf_file_boundaries.close()
It works with finer mesh, but I don’t know why this problem exists with coarser mesh.
Thanks!
