This page was generated from unit-6.1.1-surfacemeshes/surface_meshes.ipynb.

6.1.1 Generating surface meshes

from ngsolve import *
from netgen.csg import *
from netgen.meshing import MeshingStep
from ngsolve.webgui import Draw

Closed surfaces of volume meshes

When creating a volume mesh NETGEN starts meshing the surface and afterwards the volume. By setting the perstepsend argument in the CreateMesh method we can tell NETGEN to stop after the surface is meshed resulting in a valid surface mesh.

geo = CSGeometry()
geo.Add(Sphere(Pnt(0,0,0),1))
mesh = Mesh(geo.GenerateMesh(maxh=0.5, perfstepsend=MeshingStep.MESHSURFACE))
Draw(mesh)

 Start Findpoints
 main-solids: 1
 Found points 0
 Analyze spec points
 Find edges
 0 edges found
 Check intersecting edges
 CalcLocalH: 2 Points 0 Elements 0 Surface Elements
 Start Findpoints
 main-solids: 1
 Found points 0
 Analyze spec points
 Find edges
 0 edges found
 Check intersecting edges
 Start Findpoints
 main-solids: 1
 Found points 0
 Analyze spec points
 Find edges
 0 edges found
 Check intersecting edges
 Surface 1 / 1
 load internal triangle rules
 Surface meshing done
 0 illegal triangles
 Optimize Surface
 Edgeswapping, topological
 Smoothing
 Combine improve
 Smoothing
 Edgeswapping, metric
 Smoothing
 Combine improve
 Smoothing
 Edgeswapping, metric
 Smoothing
 Combine improve
 Smoothing
 126 elements, 65 points
 SplitSeparateFaces
 Update mesh topology
 Update clusters

BaseWebGuiScene

With this arbitrary complex closed surface meshes based on CSG volume meshes can be generated.

geo = CSGeometry()
sphere_bot = Sphere(Pnt(0,0,0),1)
sphere_middle = Sphere(Pnt(0,1,0),0.7)
sphere_top = Sphere(Pnt(0,1.9,0),0.4)
geo.Add(sphere_bot+sphere_middle+sphere_top)

mesh = Mesh(geo.GenerateMesh(maxh=0.5, perfstepsend=MeshingStep.MESHSURFACE))
Draw(mesh)

 Start Findpoints
 main-solids: 1
 Found points 4
 Analyze spec points
 Find edges
 2 edges found
 Check intersecting edges
 CalcLocalH: 21 Points 0 Elements 0 Surface Elements
 Start Findpoints
 Analyze spec points
 Find edges
 2 edges found
 Check intersecting edges
 Start Findpoints
 Analyze spec points
 Find edges
 2 edges found
 Check intersecting edges
 Surface 1 / 3
 Surface meshing done
 0 illegal triangles
 Optimize Surface
 Edgeswapping, topological
 Smoothing
 Combine improve
 Smoothing
 Edgeswapping, metric
 Smoothing
 Combine improve
 Smoothing
 Edgeswapping, metric
 Smoothing
 Combine improve
 Smoothing
 120 elements, 76 points
 Surface 2 / 3
 Surface meshing done
 0 illegal triangles
 Optimize Surface
 Edgeswapping, topological
 Smoothing
 Combine improve
 Smoothing
 Edgeswapping, metric
 Smoothing
 Combine improve
 Smoothing
 Edgeswapping, metric
 Smoothing
 Combine improve
 Smoothing
 105 elements, 118 points
 Surface 3 / 3
 Surface meshing done
 0 illegal triangles
 Optimize Surface
 Edgeswapping, topological
 Smoothing
 Combine improve
 Smoothing
 Edgeswapping, metric
 Smoothing
 Combine improve
 Smoothing
 Edgeswapping, metric
 Smoothing
 Combine improve
 Smoothing
 79 elements, 154 points
 SplitSeparateFaces
 Update mesh topology
 Update clusters

BaseWebGuiScene

Curving the resulting meshes is done as for volume meshes by

mesh.Curve(2)
Draw(mesh)

 Curve elements, order = 2
 Update clusters
 Curving edges
 Curving faces

BaseWebGuiScene

Also boundary names are specified in the same manner as for volume meshes

geo = CSGeometry()
sphere_bot = Sphere(Pnt(0,0,0),1).bc("bottom")
sphere_middle = Sphere(Pnt(0,1,0),0.7).bc("middle")
sphere_top = Sphere(Pnt(0,1.9,0),0.4).bc("top")
geo.Add(sphere_bot+sphere_middle+sphere_top)
mesh = Mesh(geo.GenerateMesh(maxh=0.5, perfstepsend=MeshingStep.MESHSURFACE))
print("Bnd = ", mesh.GetBoundaries())

 Start Findpoints
 main-solids: 1
 Found points 4
 Analyze spec points
 Find edges
 2 edges found
 Check intersecting edges
 CalcLocalH: 21 Points 0 Elements 0 Surface Elements
 Start Findpoints
 Analyze spec points
 Find edges
 2 edges found
 Check intersecting edges
 Start Findpoints
 Analyze spec points
 Find edges
 2 edges found
 Check intersecting edges
 Surface 1 / 3
 Surface meshing done
 0 illegal triangles
 Optimize Surface
 Edgeswapping, topological
 Smoothing
 Combine improve
 Smoothing
 Edgeswapping, metric
 Smoothing
 Combine improve
 Smoothing
 Edgeswapping, metric
 Smoothing
 Combine improve
 Smoothing
 120 elements, 76 points
 Surface 2 / 3
 Surface meshing done
 0 illegal triangles
 Optimize Surface
 Edgeswapping, topological
 Smoothing
 Combine improve
 Smoothing
 Edgeswapping, metric
 Smoothing
 Combine improve
 Smoothing
 Edgeswapping, metric
 Smoothing
 Combine improve
 Smoothing
 105 elements, 118 points
 Surface 3 / 3
 Surface meshing done
 0 illegal triangles
 Optimize Surface
 Edgeswapping, topological
 Smoothing
 Combine improve
 Smoothing
 Edgeswapping, metric
 Smoothing
 Combine improve
 Smoothing
 Edgeswapping, metric
 Smoothing
 Combine improve
 Smoothing
 79 elements, 154 points
 SplitSeparateFaces
Bnd =  Update mesh topology
 ('bottom', 'middle', 'top')
 Update clusters

or by calling after the mesh generation the method SetBCName

mesh.ngmesh.SetBCName(0, "other_name")
print("Bnd = ", mesh.GetBoundaries())

Bnd =  ('other_name', 'middle', 'top')

where the index is 0-based. To label edges, so-called BBoundaries, if available, one can use the NameEdge method

geo = CSGeometry()
sphere_bot = Sphere(Pnt(0,0,0),1).bc("bottom")
sphere_middle = Sphere(Pnt(0,1,0),0.7).bc("middle")
sphere_top = Sphere(Pnt(0,1.9,0),0.4).bc("top")
geo.Add(sphere_bot+sphere_middle+sphere_top)

geo.NameEdge(sphere_bot,sphere_middle, "lower")
geo.NameEdge(sphere_middle,sphere_top, "upper")

mesh = Mesh(geo.GenerateMesh(maxh=0.5, perfstepsend=MeshingStep.MESHSURFACE))
print("BBnd = ", mesh.GetBBoundaries())

 Start Findpoints
 main-solids: 1
 Found points 4
 Analyze spec points
 Find edges
 2 edges found
 Check intersecting edges
 CalcLocalH: 21 Points 0 Elements 0 Surface Elements
 Start Findpoints
 Analyze spec points
 Find edges
 2 edges found
 Check intersecting edges
 Start Findpoints
 Analyze spec points
 Find edges
 2 edges found
 Check intersecting edges
 Surface 1 / 3
 Surface meshing done
 0 illegal triangles
 Optimize Surface
 Edgeswapping, topological
 Smoothing
 Combine improve
 Smoothing
 Edgeswapping, metric
 Smoothing
 Combine improve
 Smoothing
 Edgeswapping, metric
 Smoothing
 Combine improve
 Smoothing
 120 elements, 76 points
 Surface 2 / 3
 Surface meshing done
 0 illegal triangles
 Optimize Surface
 Edgeswapping, topological
 Smoothing
 Combine improve
 Smoothing
 Edgeswapping, metric
 Smoothing
 Combine improve
 Smoothing
 Edgeswapping, metric
 Smoothing
 Combine improve
 Smoothing
 105 elements, 118 points
 Surface 3 / 3
 Surface meshing done
 0 illegal triangles
 Optimize Surface
 Edgeswapping, topological
 Smoothing
 Combine improve
 Smoothing
 Edgeswapping, metric
 Smoothing
 Combine improve
 Smoothing
 Edgeswapping, metric
 Smoothing
 Combine improve
 Smoothing
 79 elements, 154 points
 SplitSeparateFaces
BBnd =  Update mesh topology
 ('lower', 'upper')
 Update clusters

Here, the edge resulting from intersecting two geometries is labeled. After the mesh was generated the 1-based SetCD2Name (CD = co-dimension) method can be used.

mesh.ngmesh.SetCD2Name(1, "other_bbnd_name")
print("BBnd = ", mesh.GetBBoundaries())

BBnd =  ('other_bbnd_name', 'upper')

It is possible to add and label single points, called BBBoundaries, e.g. for defining a point load, by

geo = CSGeometry()
geo.Add(Sphere(Pnt(0,0,0),1))
geo.AddPoint(Pnt(0,0,1), "pntload")
mesh = Mesh(geo.GenerateMesh(maxh=0.5, perfstepsend=MeshingStep.MESHSURFACE))
Draw(mesh)
print("BBBnd = ", mesh.GetBBBoundaries())

 Start Findpoints
 main-solids: 1
 Found points 1
 Analyze spec points
 Find edges
 0 edges found
 Check intersecting edges
 CalcLocalH: 3 Points 0 Elements 0 Surface Elements
 Start Findpoints
 Analyze spec points
 Find edges
 0 edges found
 Check intersecting edges
 Start Findpoints
 Analyze spec points
 Find edges
 0 edges found
 Check intersecting edges
 Surface 1 / 1
 Surface meshing done
 0 illegal triangles
 Optimize Surface
 Edgeswapping, topological
 Smoothing
 Combine improve
 Smoothing
 Edgeswapping, metric
 Smoothing
 Combine improve
 Smoothing
 Edgeswapping, metric
 Smoothing
 Combine improve
 Smoothing
 126 elements, 66 points
 SplitSeparateFaces
 Update mesh topology
 Update clusters

BBBnd =  ('pntload',)

Again, renaming after mesh creation is done by (1-based index)

mesh.ngmesh.SetCD3Name(1, "other_bbbnd_name")
print("BBBnd = ", mesh.GetBBBoundaries())

BBBnd =  ('other_bbbnd_name',)

Non-closed surfaces of volume mesh

To generate non-closed surface meshes based on a volume geometry one has to use the AddSurface method of CSGeometry. Therefore, a (finite) volume geometry needs to be defined, which acts as base geometry. Next we use AddSurface to add a surface geometry, which is cut automatically with the base geometry. Here, no additional flag needs to be set in the GenerateMesh method.

geo       = CSGeometry()
cyl       = Cylinder(Pnt(0,0,0), Pnt(1,0,0), 1)
bot       = Plane(Pnt(0,0,0), Vec(0,0,-1))
right     = Plane( Pnt(3,0,0), Vec(1,0,0))
left      = Plane(Pnt(0,0,0), Vec(-1,0,0))
finitecyl = cyl * bot * left * right

geo.AddSurface(cyl, finitecyl)
mesh = Mesh(geo.GenerateMesh(maxh=0.5))
Draw(mesh)

 Start Findpoints
 main-solids: 1
 Found points 6
 Analyze spec points
 Find edges
 4 edges found
 Check intersecting edges
 CalcLocalH: 24 Points 0 Elements 0 Surface Elements
 Start Findpoints
 Analyze spec points
 Find edges
 4 edges found
 Check intersecting edges
 Start Findpoints
 Analyze spec points
 Find edges
 4 edges found
 Check intersecting edges
 Surface 1 / 1
 Surface meshing done
 0 illegal triangles
 Optimize Surface
 Edgeswapping, topological
 Smoothing
 Combine improve
 Smoothing
 Edgeswapping, metric
 Smoothing
 Combine improve
 Smoothing
 Edgeswapping, metric
 Smoothing
 Combine improve
 Smoothing
 94 elements, 60 points
 SplitSeparateFaces
 Remove Illegal Elements
 Volume Optimization
 CombineImprove
 ImproveMesh
 SplitImprove
 ImproveMesh
 SwapImprove
 SwapImprove2
 ImproveMesh
 CombineImprove
 ImproveMesh
 SplitImprove
 ImproveMesh
 SwapImprove
 SwapImprove2
 ImproveMesh
 CombineImprove
 ImproveMesh
 SplitImprove
 ImproveMesh
 SwapImprove
 SwapImprove2
 ImproveMesh
 Update mesh topology
 Update clusters

BaseWebGuiScene

Curving, adding points, and labeling boundaries, BBoundaries, and BBBoundaries is done in the same matter as before, e.g.

geo       = CSGeometry()
cyl       = Cylinder(Pnt(0,0,0), Pnt(1,0,0), 1)
bot       = Plane(Pnt(0,0,0), Vec(0,0,-1))
right     = Plane( Pnt(3,0,0), Vec(1,0,0))
left      = Plane(Pnt(0,0,0), Vec(-1,0,0))
finitecyl = cyl * bot * left * right

geo.AddSurface(cyl, finitecyl.bc("surface"))
geo.AddPoint(Pnt(0,0,1), "pntload")

geo.NameEdge(cyl,bot, "sym")
geo.NameEdge(cyl,left, "left")
geo.NameEdge(cyl,right, "right")

mesh = Mesh(geo.GenerateMesh(maxh=0.5))
mesh.Curve(2)
Draw(mesh)
print("Bnd   = ", mesh.GetBoundaries())
print("BBnd  = ", mesh.GetBBoundaries())
print("BBBnd = ", mesh.GetBBBoundaries())

 Start Findpoints
 main-solids: 1
 Found points 6
 Analyze spec points
 Find edges
 5 edges found
 Check intersecting edges
 CalcLocalH: 24 Points 0 Elements 0 Surface Elements
 Start Findpoints
 Analyze spec points
 Find edges
 5 edges found
 Check intersecting edges
 Start Findpoints
 Analyze spec points
 Find edges
 5 edges found
 Check intersecting edges
 Surface 1 / 1
 Surface meshing done
 0 illegal triangles
 Optimize Surface
 Edgeswapping, topological
 Smoothing
 Combine improve
 Smoothing
 Edgeswapping, metric
 Smoothing
 Combine improve
 Smoothing
 Edgeswapping, metric
 Smoothing
 Combine improve
 Smoothing
 94 elements, 60 points
 SplitSeparateFaces
 Remove Illegal Elements
 Volume Optimization
 CombineImprove
 ImproveMesh
 SplitImprove
 ImproveMesh
 SwapImprove
 SwapImprove2
 ImproveMesh
 CombineImprove
 ImproveMesh
 SplitImprove
 ImproveMesh
 SwapImprove
 SwapImprove2
 ImproveMesh
 CombineImprove
 ImproveMesh
 SplitImprove
 ImproveMesh
 SwapImprove
 SwapImprove2
 ImproveMesh
 Update mesh topology
 Update clusters
 Curve elements, order = 2
 Update clusters
 Curving edges
 Curving faces

Bnd   =  ('surface',)
BBnd  =  ('sym', 'sym', 'left', 'right', 'left')
BBBnd =  ('pntload',)

Non-closed surfaces via embedding

One can define a (possibly complicated) surfaces by a mapping

\[\Phi:[0,1]\times[0,1]\times\{0\}\to\mathbb{R}^3\]

from math import pi
L  = 12
W  = 1.1
mapping = lambda x,y,z : (x*L, -W*(y-0.5)*sin(pi/2*x), W*(y-0.5)*cos(pi/2*x) )

Then the MakeStructuredSurfaceMesh function is used, which takes the mapping and generates the corresponding structured surface mesh. If no mapping is given, the unit-square \([0,1]\times[0,1]\times\{0\}\) is generated.

from ngsolve.meshes import MakeStructuredSurfaceMesh
mesh = MakeStructuredSurfaceMesh(quads=True, nx=12, ny=2, mapping=mapping)
Draw(mesh)

 Update mesh topology
 Update clusters

BaseWebGuiScene

The number of elements in \(x\) and \(y\) direction (on the unit-square) can be specified, by nx and ny, as well as if a structured quadrilateral or triangular mesh should be used. The triangle orientation can be changed via the flip_triangles argument.

mesh = MakeStructuredSurfaceMesh(quads=False, nx=12, ny=2, mapping=mapping, flip_triangles=True)
Draw(mesh)

 Update mesh topology
 Update clusters

BaseWebGuiScene

Note, that the boundary names are always set to “left, right, top, bottom” according the unit-square and can be changed via the SetCD2Name method.

print(mesh.GetBBoundaries())
mesh.ngmesh.SetCD2Name(1, "other_name")
print(mesh.GetBBoundaries())

('bottom', 'right', 'top', 'left')
('other_name', 'right', 'top', 'left')

To name points as BBBoundaries one can specify them in the array bbbpts and add the corresponding names in the bbbnames list. Note, that the number of elements together with the mapping have to be specified such that the points in the list are already Gridpoints. Otherwise an exception is thrown.

mesh = MakeStructuredSurfaceMesh(nx=12, ny=2, mapping=mapping, bbbpts=[[L,0,0]], bbbnames=["force"])
print(mesh.GetBBBoundaries())

('force',)
 Update mesh topology
 Update clusters