# 6.1.1 Generating surface meshes¶

[1]:

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.

[2]:

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

[2]:

BaseWebGuiScene


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

[3]:

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)

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

[3]:

BaseWebGuiScene


Curving the resulting meshes is done as for volume meshes by

[4]:

mesh.Curve(2)
Draw(mesh)

[4]:

BaseWebGuiScene


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

[5]:

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

[6]:

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

[7]:

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

[8]:

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

[9]:

geo = CSGeometry()
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)

[10]:

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.

[11]:

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

mesh = Mesh(geo.GenerateMesh(maxh=0.5))
Draw(mesh)

[11]:

BaseWebGuiScene


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

[12]:

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


## 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$
[13]:

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.

[14]:

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

[14]:

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.

[15]:

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

[15]:

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.

[16]:

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.

[17]:

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

('force',)

[ ]: