Stop meshing since boundary mesh is overlapping

Hi,
I am trying to mesh a unit cube containing several randomly distributed spherical inclusions. Since the cube is triply periodic (along x, y and z) I create matching spheres on the faces/edges to be able to generate a mesh with matching pairs of vertices.

I am currently running two cases:

  1. Not including the spheres (basically meshing cube - all_spheres object). This runs fine and produces a matching mesh as desired.

  2. Including the spheres. In this case I get the error, " Stop meshing since boundary mesh is overlapping ". I am assuming that it could arise from the spheres on two faces not being exactly mirrored. But then why would case 1. run perfectly fine? Any suggestions?

I am attaching the portion of the code that generates the mesh

def generateGmshMesh(xc, pc, rc):
    """
    docstring
    """    
    left = Plane(Pnt(0,0,0),Vec(-1,0,0))
    bottom = Plane(Pnt(0,0,0),Vec(0,-1,0))
    back = Plane(Pnt(0,0,0),Vec(0,0,-1))
    right = Plane(Pnt(1,0,0),Vec(1,0,0))
    top = Plane(Pnt(0,1,0),Vec(0,1,0))
    front = Plane(Pnt(0,0,1),Vec(0,0,1))

    cube = left * right * bottom * top * back * front
    all_spheres = [Sphere(Pnt(sph[0], sph[1], sph[2]), rc[i]) for i, sph in enumerate(xc)]
    all_particles = [Sphere(Pnt(sph[0], sph[1], sph[2]), rc[i]) for i, sph in enumerate(pc)]

    for i, sph in enumerate(all_spheres):
        if i==0:
            sphAll = all_spheres[i]
        else:
            sphAll += all_spheres[i]

    for i, sph in enumerate(all_particles):
        if i==0:
            partAll = all_particles[i]
        else:
            partAll += all_particles[i]
     
    # reduce(operator.mul, all_spheres)
    # partAll = reduce(operator.mul, all_particles, 1)

    matrix = cube - sphAll
    partsinCube = cube * partAll
    matrix.maxh(0.01)
    partsinCube.maxh(0.005)
    matrix.transp()
    partsinCube.col([1,0,0])

    matrix.mat("mat")
    partsinCube.mat("parts")
    geo = CSGeometry()
    geo.Add(matrix)
    geo.Add(partsinCube)
    geo.PeriodicSurfaces(left,right)
    geo.PeriodicSurfaces(bottom,top)
    geo.PeriodicSurfaces(back,front)
    geo.Draw()
    Redraw()
    msh = geo.GenerateMesh()
    msh.Export("trialMesh.msh", "Gmsh2 Format")
    msh.Export("trialMesh.inp", "Abaqus Format")

I tried playing with the maxh paramter but to no avail.

could you attach the calling code as well?
Best
Christopher

Hi Christopher,
Thanks for your reply. Here is the calling code:

rad = np.loadtxt("radii_10.txt")
sph = np.loadtxt("spheres_10.txt")
generateGmshMesh(sph, sph, rad) 

Note that if I do

rad = np.loadtxt("radii_10.txt")
sph = np.loadtxt("spheres_10.txt")
generateGmshMesh(sph,[], rad) 

and comment out the relevant lines in the function, it generates the mesh as desired. The microstructure looks something like this. The spheres are close but non-intersecting, and the unit cell (cube with spheres) is periodic by construction.

https://ngsolve.org/media/kunena/attachments/1242/spheres_10.txt

https://ngsolve.org/media/kunena/attachments/1242/radii_10.txt

Attachment: spheres_10_2020-03-02.txt

Attachment: radii_10_2020-03-02.txt