Open Cascade Technology Geometry (OCC Geometry)

2.4. Open Cascade Technology Geometry (OCC Geometry)#

Netgen provides a Python wrapper around the Open Cascade Technology (OCCT) geometry kernel. It allows to model complex geometric objects. It also allows to import models via STEP format, explore, and modify the geometry.

In case you are building Netgen/NGsolve from scratch you need the cmake-flag -DUSE_OCC=ON.

In the wrapper we aim in bringing most of the structure of OCCT to Python. If you are familiar with the C++ interface of Open Cascade you will recognize many classes.

from ngsolve import *
from ngsolve.webgui import Draw
from netgen.occ import *

2.4.1. Construction of 2D objects (Workplanes)#

Workplanes are the basic building blocks of 2D geometry, they can be understood as a plane in 3D space.

There are several methods that are already implemented to define standard shapes (Rectangle, Circle, Ellipse…).

After “sketching” on the workplane you need to extact the object of interest.

wp = WorkPlane()
wp.Rectangle(2, 5) # you still cannot draw a rectangle with the same syntax as in the GUI

# Draw(topo_rect) # try ... it fails

<netgen.libngpy._NgOCC.WorkPlane at 0x7f7c704824b0>

Before drawing the object you need to extract the face from the workplane. This is done by the Face method.

rect = wp.Face()
Draw(rect);

2.4.1.1. Orientation of the workplane#

The concept of Workplane is quite powerful, it allows to draw on a 2D plane and automatically map it to the 3D space.

The WorkPlane is initialized with:

class WorkPlane(
    axes: Axes = ...,
    pos: gp_Ax2d = ...
)

where Axes is a set a point and two vectors that define the orientation of the workplane.

Axes(
    p: gp_Pnt = (0, 0, 0), # origin 
    n: gp_Dir = (0, 0, 1), # normal vector
    h: gp_Dir = (1, 0, 0)  # horizontal vector
    ) -> None

The gp_Ax2d is a set of a point in 2D and the main direction.

gp_Ax2d(
    p: gp_Pnt2d = (0, 0), # origin
    v: gp_Dir2d = (1, 0)  # direction
    ) -> None

axes = Axes( (1,1,0) ,(0,1,0) , (-1,1,0) ) 

# axes = Axes((1, 1, 0), Y, Y-X) # also works

pos = gp_Ax2d()

wp2 = WorkPlane( axes = axes, pos = pos )
rect2 = wp2.Rectangle(2, 5).Face()
Draw(rect2)

BaseWebGuiScene

union = rect + rect2
Draw(union);

Let’s create the a workplane, the default position is the center of the plane and the default direction is the x-axis.

To create a figure think of it like giving directions to a driver. For example:

1. Go straight for 1 meter
2. immediately turn (left) 60 degrees

then repeat the previous steps for 6 times in total.

wp = WorkPlane()
for i in range(6):
    wp.Line(1).Rotate(60)
Draw(wp.Face());

Suppose to change slightly the above instructions to:

1. Go straight for 1 meter
2. Turn (right) 60 degrees while still driving so that you make an arc of 0.5 meters

then repeat the previous steps for 6 times in total.

for i in range(6):
    wp.Line(1).Arc(0.5, 60)
face = wp.Face()
Draw(face);

The drawing needs to be positively oriented, otherwise the face will not be created.

If you have drawn the object in the wrong direction you can use the Reverse method to change the orientation.

A negative oriented Wire is the same of subtracting to the figure a positive oriented Wire.

wp = WorkPlane().RectangleC(2, 1) \
    .Circle(0.5, 0, 0.2).Reverse() \
    .Circle(-0.5, 0, 0.2).Reverse()
Draw(wp.Face());

There are two main ways to obtain a 3D object from a 2D object:

Extrusion in a certain direction
Revolution around an axis

face = WorkPlane().MoveTo(10, 0).Line(18).Arc(2, 90).Line(6).Arc(2, 90) \
    .Line(18).Rotate(90).Line(2).Rotate(90).Line(10) \
    .Arc(3, -180).Line(10, "nice line").Close().Face()

Draw(face);

ext = face.Extrude(-40).Move((10, 0, 0))

rev = face.Revolve(Axis((0, 0, 0), Y+Z+X), 280).Move((-10, 0, 0))

Draw(rev+ext );

2.4.1.2. Modify 3D objects#

fillet = ext.MakeFillet(ext.edges, 0.3)

Draw(fillet);

chamfer = ext.MakeChamfer(ext.edges, 0.3)

Draw(chamfer);

There are several other ready-to-use 2D shapes. Look into the module to find more.

2.4.1.3. Construction of 3D objects#

We define a box aligned with the Cartesian coordinates given by two points with minimal and maximal x/y/z coordinates. A cylinder is given by a point on the axis, a direction vector, the radius \(r\), and the height \(h\). The symbols X, Y, and Z are predefined basis vectors for the Cartesian coordinates.

box = Box(Pnt(0,0,0), Pnt(1,1,1))
cyl = Cylinder(Pnt(1,0.5,0.5), X, r=0.3, h=0.5)

The Boolean operations fuse, common, and cut provided by OCCT are made available by the operators +, * and -. More can be found in the OCCT-documentation. Note that the fuse generates one new solid, there is no interface face where the cylinder is touching the box.

# name the faces 
cyl.faces.Min(X).name = "cyl_minx"
box.faces.Min(X).name = "box_minx"
box.faces.Max(X).name = "box_maxx"

# create the union
fused = box+cyl

# impose beforehand the mesh size
fused.faces["box_.*"].maxh = 0.05

# che check the names
Draw(fused);

The generated objects are Py-wrapped OCCT objects derived from TopoDS_Shape.

from ngsolve import Mesh
from ngsolve.webgui import Draw

geo = OCCGeometry(fused)
mesh = Mesh(geo.GenerateMesh(maxh=0.2, grading=0.9))
mesh.Curve(4)
Draw (mesh, clipping={"z":-1});

Instead of fusing, we can glue together shapes. Then, the resulting composite solid contains the interface face between the solids. This is important when dealing with separate material regions:

geo = Glue( [ cyl, box]) 

mesh = Mesh(OCCGeometry(geo).GenerateMesh(maxh=0.2, grading=0.9))
mesh.Curve(4)

Draw (mesh, clipping={"z": -1});

be aware of the order
change the order and look at the names of the common face

A third option is to form a compound of shapes, then the component shapes are meshed independently:

geo = Compound( [box, cyl])

mesh = Mesh(OCCGeometry(geo).GenerateMesh(maxh=0.2)).Curve(3)
Draw (mesh, clipping=True);

2.4.2. Transformation of shapes#

We can translate, rotate and mirror shapes. The transformation preserves the original shape, and returns a transformed copy. The translation Move takes a vector as argument, the rotation Rotation needs an axis given by a point and a direction, and an angle. The sign of the angle reflects the right hand screw rule.

solid = Box((0,0,0), (5,3,1)) + Sphere((0,0,0), 0.3)
solid2 = solid.Move((5,0,2))
solid3 = solid.Move((0,0,4)).Rotate( Axis((0,0,4), X), 45)
Draw (solid + solid2 + solid3);

2.4.3. Selection of shapes#

Sometimes we want to set properties of shapes, boundary conditions or mesh-size. For that we can use shape selectors.

The Max and Min selectors finds the sub-shapes where the center of gravity has maximal or minimal coordinates in a given direction

More of them will come soon …

box = Box((0,0,0), (1,1,1))
box.faces.Max(Z).col = (1,0,0)
box.faces.Min(Y).col = (0,0,1)
Draw(box);

Faces where center of gravity has \(x\) coordinate less than 0.8:

Draw(Compound(box.faces[X<0.8]));

2.4.3.1. Last remark on the workflow OCC-Netgen-NGSolve#

The workflow on how to create a domain is as follows:

Create the geometry using OCC
Pass the geometry to Netgen with the meshing parameters
Mesh the geometry using Netgen
Pass the netgen mesh to NGSolve

OCC shape	Netgen shape	Netgen mesh	NGSolve mesh
TopoDS_Shape	OCCGeometry(…)	.GenerateMesh(…)	Mesh(…)

2.4.4. Metamaterials:#

using the OCC module it is possible to create complex geometries for metamaterials. Let’s createa simle geometry with a certain chirality.

in this case a crown with some additional arms that pop out of the figure.

def BaseShape(r_in, r_out, l):
    c_in = Circle(Pnt(0, 0), r_in).Face()
    c_out = Circle(Pnt(0, 0), r_out).Face()
    h = r_out - r_in

    r_top = MoveTo(0, r_in).Rectangle(l/2, h).Face()
    t_right = r_top.Rotate(Axis((0, 0, 0), (0, 0, 1)), 90)
    r_bot = r_top.Rotate(Axis((0, 0, 0), (0, 0, 1)), 180)
    t_left = r_top.Rotate(Axis((0, 0, 0), (0, 0, 1)), 270)
    shape = c_out - c_in + r_top + r_bot + t_right + t_left

    return shape


l = 6
r_in = 1
r_out = 2
base = BaseShape(r_in, r_out, l)

Draw(base);

figure = base
for i in range(10):
    for j in range(5):
        if (i+j) % 2 == 0:
            figure += base.Move((i*6, j*6, 0))
        else:
            figure += base.Mirror(Axis((0, 0, 0), (1, 0, 0))
                                  ).Move((i*6, j*6, 0))


Draw(figure)

BaseWebGuiScene

figure = figure.Extrude(2)

# add a left handle to the figure 1 x 5*l x 2
lbar = Box((0, 0, 0), (1, 5*l, 2)).Move((-l/2-1, -l/2, 0))
figure += lbar

# add a right handle to the figure 1 x 5*l x 2
rbar = Box((0, 0, 0), (1, 5*l, 2)).Move((9*l + l/2, -l/2, 0))
figure += rbar

figure.faces.Min(X).name = "left"
figure.faces.Max(X).name = "right"


# Draw(base)
Draw(figure)

BaseWebGuiScene

Example of elasticity problem with metamaterials

from ngsolve.solvers import Newton


E, nu = 210, 0.2
mu = E / 2 / (1+nu)
lam = E * nu / ((1+nu)*(1-2*nu))


def C(u):
    F = Id(3) + Grad(u)
    return F.trans * F


def NeoHooke(C):
    return 0.5*mu*(Trace(C-Id(3)) + 2*mu/lam*Det(C)**(-lam/2/mu)-1)


factor = Parameter(0)
force = CoefficientFunction((factor, 0, 0))


mesh = Mesh(OCCGeometry(figure).GenerateMesh(maxh=2)).Curve(4)

fes = H1(mesh, order=2, dirichlet="left", dim=mesh.dim)
u = fes.TrialFunction()

a = BilinearForm(fes)
a += Variation(NeoHooke(C(u)).Compile()*dx)
a += Variation((-InnerProduct(force, u)).Compile()*ds("right"))

gfu = GridFunction(fes)
gfu.vec[:] = 0


scene = Draw(gfu, mesh, deformation=True, scale=2)


numsteps = 2
for step in range(numsteps):
    print(step)
    factor.Set(0.25*(step+1)/numsteps)
    with TaskManager():
        Newton(a, gfu, printing=True, dampfactor=0.9)
    scene.Redraw()

0
Newton iteration  0

err =  3.7671381034659084
Newton iteration  1

err =  1.3623778861407454
Newton iteration  2

err =  0.02762639706179149
Newton iteration  3

err =  0.0009191364582180783
Newton iteration  4

err =  2.120839539324405e-06
Newton iteration  5

err =  1.2240630960776989e-11
Newton iteration  6

err =  1.0479399806498996e-13

1
Newton iteration  0

err =  3.4505115334315652
Newton iteration  1

err =  0.9868738993964846
Newton iteration  2

err =  0.017656445692363357
Newton iteration  3

err =  0.00023902478664546687
Newton iteration  4

err =  1.3738727817032852e-07
Newton iteration  5

err =  1.840911034197887e-13