diff --git a/py_tutorials/merge.py b/py_tutorials/merge.py index 931d095a..abacd38e 100644 --- a/py_tutorials/merge.py +++ b/py_tutorials/merge.py @@ -3,28 +3,35 @@ 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() -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 ("** start merging **") -print ("********************") +print ("***************************") +print ("** merging suface meshes **") +print ("***************************") +# create an empty mesh mesh = Mesh() -fd_outside = mesh.Add (FaceDescriptor(surfnr=1,domin=1)) -fd_inside = mesh.Add (FaceDescriptor(surfnr=2,domin=2,domout=1)) +# 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(): @@ -32,11 +39,17 @@ for e in m1.Elements2D(): 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: