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)

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)

BaseWebGuiScene

Curving the resulting meshes is done as for volume meshes by

mesh.Curve(2)
Draw(mesh)

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

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

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

BBnd =  ('lower', 'upper')

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

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)

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

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)

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)

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