1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
|
from netgen.meshing import *
from netgen.csg import *
from ngsolve import ngsglobals
ngsglobals.msg_level = 2
# generate brick and mesh it
geo1 = CSGeometry()
geo1.Add (OrthoBrick( Pnt(0,0,0), Pnt(1,1,1) ))
m1 = geo1.GenerateMesh (maxh=0.1)
m1.Refine()
# generate sphere and mesh it
geo2 = CSGeometry()
geo2.Add (Sphere (Pnt(0.5,0.5,0.5), 0.1))
m2 = geo2.GenerateMesh (maxh=0.05)
m2.Refine()
m2.Refine()
print ("****************************")
print ("** merging surface meshes **")
print ("****************************")
# create an empty mesh
mesh = 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 = mesh.Add (FaceDescriptor(bc=1,domin=1,surfnr=1))
fd_inside = mesh.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] = mesh.Add (m1[v])
# copy surface elements from first mesh to new mesh
# we have to map point-numbers:
for e in m1.Elements2D():
mesh.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] = mesh.Add (m2[v])
for e in m2.Elements2D():
mesh.Add (Element2D (fd_inside, [pmap2[v] for v in e.vertices]))
print ("******************")
print ("** merging done **")
print ("******************")
mesh.GenerateVolumeMesh()
mesh.Save ("newmesh.vol")
|
