Meshing coated spheres periodically arranged in unit box

I’m trying to mesh unit cubes containing spherically coated, periodically arranged spheres; the geometries are defined via python. I’ve run into a meshing difficulty I do not understand.

In a simple case the unit cube contains just 2 compound spheres, one in its interior and one crossing a surface of the cube. If I’m content with having three parts (“matrix”, “core” and “interphase”, where core and interphase make up the coated sphere), the following Pyton sequence works fine:

from netgen.csg import *

geo = CSGeometry()
celsz = 1.000000

generate cell

left = Plane( Pnt(0.00,0.00,0.00), Vec(-1.0,0.0,0.0) )
front = Plane( Pnt(0.00,0.00,0.00), Vec(0.0,-1.0,0.0) )
bot = Plane( Pnt(0.00,0.00,0.00), Vec(0.0,0.0,-1.0) )
right = Plane( Pnt(celsz,celsz,celsz), Vec(1.0,0.0,0.0) )
back = Plane( Pnt(celsz,celsz,celsz), Vec(0.0,1.0,0.0) )
top = Plane( Pnt(celsz,celsz,celsz), Vec(0.0,0.0,1.0) )
VolEl = left * right * bot * top * front * back

center coordinates and radii of spheres

numSph = 2

centers of spheres

spCc = Matrix(numSph,3)
spCc[0,0] = 0.440789
spCc[0,1] = 0.466176
spCc[0,2] = 0.373534
spCc[1,0] = 0.423044
spCc[1,1] = 0.034480
spCc[1,2] = 0.948371

radii of outer spheres

spRdP = Vector(numSph)
spRdP[0] = 0.212157
spRdP[1] = 0.212157

radii of inner spheres

spRdL = Vector(numSph)
spRdL[0] = 0.180000
spRdL[1] = 0.180000

function for duplicating spheres for periodicity

def DuplicateSphere (c,r):
spheres = None
for x in [-1,0,1]:
for y in [-1,0,1]:
for z in [-1,0,1]:
sph = Sphere(c+celsz*Vec(x,y,z),r)
if spheres:
spheres = spheres + sph
else:
spheres = sph
return spheres

set up compound spheres corresponding to particles

sphere00

cp = Pnt(spCc[0,0],spCc[0,1],spCc[0,2])
sphP00 = DuplicateSphere (cp,spRdP[0])
sphL00 = DuplicateSphere (cp,spRdL[0])
sphrsP = sphP00
sphrsL = sphL00

sphere01

cp = Pnt(spCc[1,0],spCc[1,1],spCc[1,2])
sphP01 = DuplicateSphere (cp,spRdP[1])
sphL01 = DuplicateSphere (cp,spRdL[1])
sphrsP = sphrsP + sphP01
sphrsL = sphrsL + sphL01
#==============================================================================

set up geometry - parts correspond to phases

geo.Add (VolEl-sphrsP)
geo.Add (sphrsL*VolEl)
geo.Add ((sphrsP-sphrsL)*VolEl, maxh=0.030000)
#==============================================================================

geo.PeriodicSurfaces(left, right)
geo.PeriodicSurfaces(top, bot)
geo.PeriodicSurfaces(front, back)
geo.SetBoundingBox(Pnt(-1,-1,-1), Pnt(2,2,2))

mesh = geo.GenerateMesh(maxh=0.0500)
mesh.SecondOrder()

mesh.Save (“test.vol”)

If, however, I try to obtain separate parts for the core and the coating of each sphere by using

#============================================================================== # set up geometry - core and layer parts correspond to individual spheres geo.Add (VolEl-sphP00-sphP01) geo.Add (sphP00*VolEl) geo.Add ((sphP00-sphL00)*VolEl, maxh=0.030000) geo.Add (sphP01*VolEl) geo.Add ((sphP01-sphL01)*VolEl, maxh=0.030000) #==============================================================================
NetGen returns an error. If I mesh only sphere00, I get the expected mesh, but if I try meshing only sphere01, NetGen throws a
netgen.libngpy._meshing.NgException: Meshing failed!

Can anybody help me understand what I did wrong?