Solid mesh generation

Hello, I’m using Netgen for generating the cube mesh (Matrix + ellipsoids).

Now I can generate the mesh based on examples.

However, when importing the mesh into COMSOL, there are not “domains” in the mesh.

I thought “GenerateVolumeMesh()” gave me the mesh, but solid mesh.

I want the solid mesh. In other words, solid ellipsoids and solid matrix are what I wanted, not hollowed ones.

Here is the code that I use.

If you don’t mind, please tell me which part should be modified to generate the solid mesh?

# Generate brick and mesh it
geo1 = CSGeometry()
cube = OrthoBrick( Pnt(0,0,0), Pnt(64,64,64) )
geo1.Add(cube)
m1 = geo1.GenerateMesh (maxh=5)
# m1.Refine()

# Read Coordinates_XYZ and AxisLength_abc from *.txt files
Coordinates_XYZ = np.loadtxt("outputs/PackingAlgorithm/Centroid_XYZ_periodic.txt")
AxisLength_abc  = np.loadtxt("outputs/PackingAlgorithm/AxisLength_Manual_abc_31p_periodic.txt")

# Make Ellipsoids based on given information
t0_start = time.time()
# Generate sphere and mesh it
geo2 = CSGeometry()
ellips = [0] * len(Coordinates_XYZ[:,0])
for i in range(0,len(Coordinates_XYZ[:,0])):
    ellips[i] = Ellipsoid(Pnt(Coordinates_XYZ[i,0],Coordinates_XYZ[i,1],Coordinates_XYZ[i,2]),
                          Vec(AxisLength_abc[i,0],0,0),
                          Vec(0,AxisLength_abc[i,1],0),
                          Vec(0,0,AxisLength_abc[i,2]))
    # geo2.Add(ellips[i])
    print(i)

# Store all ellipsoids.
all_ellips = ellips[0]
for i in range(1,len(Coordinates_XYZ[:,0])):
    all_ellips = all_ellips + ellips[i]
# Intersection between all ellipsoids and cube
all_ellips = all_ellips * cube

geo2.Add(all_ellips,maxh=5)

m2 = geo2.GenerateMesh (maxh=5)
t0_end = time.time()
print("Making ellipsoids took --- %s seconds ---" % (t0_end - t0_start))

# create an empty mesh
ngmesh = Mesh()

# a face-descriptor stores properties associated with a set of surface elements
# bc .. boundary condition marker,
# domin/domout .. domain-number in front/back of surface elements (0 = void),
# surfnr .. number of the surface described by the face-descriptor
fd_outside = ngmesh.Add (FaceDescriptor(bc=1,domin=1,surfnr=1))
fd_inside = ngmesh.Add (FaceDescriptor(bc=2,domin=2,domout=1,surfnr=2))

# copy all boundary points from first mesh to new mesh.
# pmap1 maps point-numbers from old to new mesh
pmap1 = { }
for e in m1.Elements2D():
    for v in e.vertices:
        if (v not in pmap1):
            pmap1[v] = ngmesh.Add (m1[v])


# copy surface elements from first mesh to new mesh
# we have to map point-numbers:
for e in m1.Elements2D():
    ngmesh.Add (Element2D (fd_outside, [pmap1[v] for v in e.vertices]))

# same for the second mesh:
pmap2 = { }
for e in m2.Elements2D():
    for v in e.vertices:
        if (v not in pmap2):
            pmap2[v] = ngmesh.Add (m2[v])

for e in m2.Elements2D():
    ngmesh.Add (Element2D (fd_inside, [pmap2[v] for v in e.vertices]))

t1_start = time.time()
ngmesh.GenerateVolumeMesh()
print("Generating Volume mesh done.")
# for visualization we need a NGSolve mesh
import ngsolve
mesh = ngsolve.Mesh(ngmesh)
t1_end = time.time()
print("GenerateMesh took --- %s seconds ---" % (t1_end - t1_start))

Hi,

you do not need to generate two different meshes and then merge them. Adding the ellipsoids into the main geometry is possible after you take the difference.

I attached a jupyter-notebook generating the mesh with OCC and assigned domain names.

With COMSOL I have no experience.

cube_particle.ipynb (2.5 KB)

Best,
Michael