#NETGEN GEOMETRY DATA (PYTHON) GENERATED BY PROGRAM ARIGEN/Ss ON 02/12/26
#NumSpheres= 10     VolFract= 4.0001E-01
#                   CellSize= 1.0000E+00     RefElSize= 1.0000E-01
#Periodic arrangement of spheres, basic geometry
#
# part 1: MATRIX          (xi=0.6000)
# part 2: INHOMOGENEITIES (xi=0.4000)
#
#
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

#
# import coordinates and radii of spheres
#
numSph = 10
spCc   = Matrix(numSph,3)
spCc[0,0] =  0.003104
spCc[0,1] =  0.329353
spCc[0,2] =  0.916689
spCc[1,0] =  0.724646
spCc[1,1] =  0.017073
spCc[1,2] =  0.264546
spCc[2,0] =  0.435834
spCc[2,1] =  0.429475
spCc[2,2] =  0.274924
spCc[3,0] =  0.802864
spCc[3,1] =  0.617388
spCc[3,2] =  0.470076
spCc[4,0] =  0.484647
spCc[4,1] =  0.931435
spCc[4,2] =  0.851611
spCc[5,0] =  0.556863
spCc[5,1] =  0.381654
spCc[5,2] =  0.837692
spCc[6,0] =  0.075071
spCc[6,1] =  0.227175
spCc[6,2] =  0.467438
spCc[7,0] =  0.931254
spCc[7,1] =  0.848801
spCc[7,2] =  0.894271
spCc[8,0] =  0.289118
spCc[8,1] =  0.881403
spCc[8,2] =  0.260720
spCc[9,0] =  0.229195
spCc[9,1] =  0.623627
spCc[9,2] =  0.632293
spRdP  = Vector(numSph)
spRdP[0] =  0.212157
spRdP[1] =  0.212157
spRdP[2] =  0.212157
spRdP[3] =  0.212157
spRdP[4] =  0.212157
spRdP[5] =  0.212157
spRdP[6] =  0.212157
spRdP[7] =  0.212157
spRdP[8] =  0.212157
spRdP[9] =  0.212157

#
# 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 spheres
sphrsP = None
for iP in range(numSph):
    cp  = Pnt(spCc[iP,0],spCc[iP,1],spCc[iP,2])
    spP = DuplicateSphere (cp,spRdP[iP])
    sphrsP = spP + sphrsP if sphrsP else spP

# set up geometry
geo.Add (VolEl-sphrsP)
geo.Add (sphrsP*VolEl)
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.1000)
mesh.SecondOrder()

mesh.Save ("TEST4008_1.vol")