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))