Stop meshing since boundary mesh is overlapping

Main Forum Old-Forum - do not post here
1

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.

2

could you attach the calling code as well?
Best
Christopher

3

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.

800×600 106 KB

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