Constructive Solid Geometry (CSG)

These geometries are bases on primitives (e.g. sphere, cylinder, plane) which are used to build solids by performing boolean operations. Netgen offers the following primitives

primitive csg syntax meaning
half-space Plane(Pnt a,Vec n) point p in plane, normal vector
sphere Sphere(Pnt c,float r) sphere with center c and radius r
cylinder Cylinder(Pnt a, Pnt b, float r) points a and b define the axes of a infinite cylinder of radius r
brick OrthoBrick (Pnt a, Pnt b) axes parallel brick with minimal coordinates a and maximal coordinates b

and the boolean operators

operator set operation
* intersection
+ union
- intersection with complement
In [1]:
import netgen.gui
%gui tk
from ngsolve.solve import Draw, Redraw # just for visualization

Using these primitives and operations, we can easily construct a cube. First we import the netgen.csg module, create 6 plane and intersect them to get the solid cube.

In [2]:
from netgen.csg import *

left  = Plane (Pnt(0,0,0), Vec(-1,0,0) )
right = Plane (Pnt(1,1,1), Vec( 1,0,0) )
front = Plane (Pnt(0,0,0), Vec(0,-1,0) )
back  = Plane (Pnt(1,1,1), Vec(0, 1,0) )
bot   = Plane (Pnt(0,0,0), Vec(0,0,-1) )
top   = Plane (Pnt(1,1,1), Vec(0,0, 1) )

cube = left * right * front * back * bot * top

Then we create a CSGeometry object and add the solid.

In [3]:
geo = CSGeometry()
geo.Add (cube)

mesh = geo.GenerateMesh(maxh=0.25)
Redraw()
# mesh.Save("cube.vol")
In [4]:
from netgen.csg import *

cube = OrthoBrick( Pnt(0,0,0), Pnt(1,1,1) )
hole = Cylinder ( Pnt(0.5, 0.5, 0), Pnt(0.5, 0.5, 1), 0.2)

geo = CSGeometry()
geo.Add (cube-hole.maxh(0.05))
mesh = geo.GenerateMesh(maxh=0.1)
Redraw()

Setting properties of solids

A solid has members which we can set to define the desired properties.

In [5]:
sphere = Sphere(Pnt(0,0,0),1)

Now we can set a boundary name and a maximal mesh size on the surface of this sphere

In [6]:
sphere.bc("sphere").maxh(0.25)
Out[6]:
<netgen.libngpy._csg.Solid at 0x7fbb64148bc8>

and define a material for the volume

In [7]:
sphere.mat("iron")
Out[7]:
<netgen.libngpy._csg.Solid at 0x7fbb64148bc8>

In case we want to visualize the geometry we can define the color (using rgb values) and transparency of the solid.

In [8]:
sphere.col([1,0,0])#.transp()
Out[8]:
<netgen.libngpy._csg.Solid at 0x7fbb64148bc8>
In [9]:
geo = CSGeometry()
geo.Add(sphere)
geo.Draw()
Redraw()
In [10]:
ngmesh = geo.GenerateMesh()
print(type(ngmesh))
Redraw()

<class 'netgen.libngpy._meshing.Mesh'>

To improve the approximation of curved geometries it is possible to use curved elements. This can be done within NGSolve. Thus we have to convert the Netgen mesh to a NGSolve mesh before curving it.

In [11]:
from ngsolve.comp import Mesh
mesh = Mesh(ngmesh)
print(type(mesh))
Redraw()

<class 'ngsolve.comp.Mesh'>
In [12]:
mesh.Curve(3)
Draw(mesh)

Setting the mesh size

There are the following options to set the mesh size:

  • globally as argument maxh of GenerateMesh
  • to the surface of one solid (maxh property as above mentioned)
  • for the volume of a solid as optional argument when adding it to the geometry Add(...,bc)
  • restrict the mesh size for one point using RestrictH
  • use CloseSurfaces to generate anisotropic meshes

Global mesh size

The global mesh size can be set with the named argument maxh. The following two versions are equivalent since all arguments of the of the GenerateMesh function are parsed to the MeshingParameters if no named argument mp is given.

In [13]:
unit_cube.GenerateMesh(maxh=0.4)
Out[13]:
<netgen.libngpy._meshing.Mesh at 0x7fbb641606c0>
In [14]:
mp = MeshingParameters(maxh=0.4)
unit_cube.GenerateMesh(mp = mp)
Out[14]:
<netgen.libngpy._meshing.Mesh at 0x7fbb641608f0>

Mesh size for one solid

To set the mesh size for one domain of the mesh we have to add the desired maxh as argument when adding the solid to the geometry

In [15]:
geo = CSGeometry()

brick = OrthoBrick(Pnt(-2,-2,-2),Pnt(2,2,2))
sphere = Sphere(Pnt(0,0,0),1)

geo.Add(brick-sphere)
geo.Add(sphere,maxh=0.1)
ngmesh = geo.GenerateMesh(maxh=0.4)

Mesh size on a surface

If we want to refine just on a surface we define it as property of the solid.

In [16]:
geo = CSGeometry()

brick = OrthoBrick(Pnt(-2,-2,-2),Pnt(2,2,2))
sphere = Sphere(Pnt(0,0,0),1)

geo.Add(brick-sphere)
geo.Add(sphere.maxh(0.1))
ngmesh = geo.GenerateMesh()

Mesh size in points

This can be done with the MeshingParameters. Using RestrictH we can define the mesh size in an arbitrary point.

In [17]:
geo = CSGeometry()

brick = OrthoBrick(Pnt(-2,-2,-2),Pnt(2,2,2))
sphere = Sphere(Pnt(0,0,0),1)

mp = MeshingParameters(maxh=0.4)
mp.RestrictH (x=0, y=0, z=1, h=0.025)
        
geo.Add(brick-sphere)
geo.Add(sphere)
ngmesh = geo.GenerateMesh(mp = mp)

Anisotropic meshes

If the geometry contains thin layers we can use CloseSurfaces to avoid elements with small angles.

In [18]:
from netgen.csg import *

geo = CSGeometry()

box = OrthoBrick(Pnt(0,0,0),Pnt(1,1,1))
top = Plane(Pnt(0,0,0.52),Vec(0,0,1))
bot = Plane(Pnt(0,0,0.48),Vec(0,0,-1))
plate = box * top * bot

geo.Add((box-top).mat("air"))
geo.Add(plate.mat("plate"))
geo.Add((box-bot).mat("air"))

slices = [2**(-i) for i in reversed(range(1,6))]
# define the close surfaces
geo.CloseSurfaces(bot,top)#,slices)
nmesh = geo.GenerateMesh(maxh=0.3)
# refine the mesh between the close surfaces
# ZRefinement(nmesh,geo)

Setting boundary conditions

Boundary condition on the surface of a solid

Setting a boundary condition on the whole surface of solid can be achieved by adding it as property to the solid.

In [19]:
brick = OrthoBrick(Pnt(-2,-2,-2),Pnt(2,2,2)).bc('outer')
sphere = Sphere(Pnt(0,0,0),1).bc('sphere')

Modify boundary between two solids

This can be done by adding the named argument bcmod when adding the solid to the geometry. Here we change the boundary condition on the surface between the halfsphere and the already added box.

In [20]:
halfsphere = sphere * Plane(Pnt(0,0,0),Vec(1,0,0)).bc('plane')
box = brick-sphere
geo = CSGeometry()
geo.Add(box.col([1,0,0]).transp())
geo.Add(halfsphere.col([0,0,1]),bcmod=[(box,"halfsphere")])
geo.Draw()
In [21]:
ngmesh = geo.GenerateMesh()
mesh = Mesh(ngmesh)
mesh.GetBoundaries()
Out[21]:
['outer',
 'outer',
 'outer',
 'outer',
 'outer',
 'outer',
 'sphere',
 'halfsphere',
 'plane']