4.2 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

[1]:

# import netgen.gui
from ngsolve 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.

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

[3]:

geo = CSGeometry()

mesh = geo.GenerateMesh(maxh=0.25)
Redraw()
# mesh.Save("cube.vol")

[3]:

True

[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()
mesh = geo.GenerateMesh(maxh=0.1)
Redraw()

[4]:

True


Setting properties of solids¶

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

[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

[6]:

sphere.bc("sphere").maxh(0.25)

[6]:

<netgen.libngpy._csg.Solid at 0x7fdf35c92230>


and define a material for the volume

[7]:

sphere.mat("iron")

[7]:

<netgen.libngpy._csg.Solid at 0x7fdf35c92230>


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

[8]:

sphere.col([1,0,0])#.transp()

[8]:

<netgen.libngpy._csg.Solid at 0x7fdf35c92230>

[9]:

geo = CSGeometry()
geo.Draw()

[10]:

ngmesh = geo.GenerateMesh()
print(type(ngmesh))
Redraw()

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

[10]:

True


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.

[11]:

from ngsolve.comp import Mesh
mesh = Mesh(ngmesh)
print(type(mesh))
Redraw()

<class 'ngsolve.comp.Mesh'>

[11]:

False

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

[13]:

unit_cube.GenerateMesh(maxh=0.4)

[13]:

<netgen.libngpy._meshing.Mesh at 0x7fdf35c89fb0>

[14]:

from netgen.meshing import MeshingParameters
mp = MeshingParameters(maxh=0.4)
unit_cube.GenerateMesh(mp = mp)

[14]:

<netgen.libngpy._meshing.Mesh at 0x7fdf35c9dcf0>


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

[15]:

geo = CSGeometry()

brick = OrthoBrick(Pnt(-2,-2,-2),Pnt(2,2,2))
sphere = Sphere(Pnt(0,0,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.

[16]:

geo = CSGeometry()

brick = OrthoBrick(Pnt(-2,-2,-2),Pnt(2,2,2))
sphere = Sphere(Pnt(0,0,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.

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

ngmesh = geo.GenerateMesh(mp = mp)


Anisotropic meshes¶

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

[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

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)

surface ids1 = 0: 21

surface ids2 = 0: 13



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.

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

[20]:

halfsphere = sphere * Plane(Pnt(0,0,0),Vec(1,0,0)).bc('plane')
box = brick-sphere
geo = CSGeometry()
geo.Draw()

[21]:

ngmesh = geo.GenerateMesh()
mesh = Mesh(ngmesh)
mesh.GetBoundaries()

[21]:

('outer',
'outer',
'outer',
'outer',
'outer',
'outer',
'sphere',
'halfsphere',
'plane')

[ ]: