7. Create CSG 3D-geometry

This demo is implemented in a single Python file, demo_csg-3D.py, and demonstrates usage of 3D geometries in DOLFIN.

This demo illustrates how to:

• Define Constructive Solid Geometry (CSG)
• Generate meshes from geometries

The domain looks as follows:

and the mesh:

7.1. Problem definition

We do not have any equation in this demo. The demo focuses purely on the geometry and how to generate a mesh from it.

7.2. Implementation

This description goes through how to make 3-dimensional geometries and meshes in DOLFIN, as implemented in demo_csg-3D.py.

First, the dolfin module is imported:

from dolfin import *

Then we check if CGAL is installed, as it is needed to run this demo:

if not has_cgal():
        print "DOLFIN must be compiled with CGAL to run this demo."
        exit(0)

Now we define 3D geometries. We start with defining a box by sending the coordinates of two opposite corners as arguments to the class Box.

box = Box(0, 0, 0, 1, 1, 1)

We create a sphere by sending the center and radius to Sphere. The center \((x_0, y_0, z_0)\) is represented by an instance of Point.

sphere = Sphere(Point(0, 0, 0), 0.3)

We define a Cone by four arguments, the center at one end Point (\(x_1,y_1,z_1\)), the center at the other end Point (\(x_2,y_2,z_2\)), and the radii at these points.

cone = Cone(Point(0, 0, -1), Point(0, 0, 1), 1., .5)

Now we have some geometries that we can play with:

g3d = box + cone - sphere

This simple line makes a geometry of the union of the box and the cone, from which the sphere is subtracted.

To get information about our new geometry we use the function info. This function takes a string or a DOLFIN object as argument, and optionally we can give a second argument to indicate whether verbose object data should be printed. If the second argument is False (which is default), a one-line summary is printed. If True, verbose and sometimes very exhaustive object data are printed.

# Test printing
info("\nCompact output of 3D geometry:")
info(g3d)
info("\nVerbose output of 3D geometry:")
info(g3d, True)

To visualize our geometry we plot it:

plot(g3d, "3D geometry (surface)")

The second argument is optional and it specifies title of the plot.

Finally, we generate a mesh using Mesh and plot it.

mesh3d = Mesh(g3d, 32)
info(mesh3d)
plot(mesh3d, "3D mesh")

Note that when we create a mesh from a CSG geometry, the resolution must be specified as a second argument to the Mesh constructor.

7.3. Complete code

from dolfin import *

if not has_cgal():
    print "DOLFIN must be compiled with CGAL to run this demo."
    exit(0)


# Define 3D geometry
box = Box(0, 0, 0, 1, 1, 1)
sphere = Sphere(Point(0, 0, 0), 0.3)
cone = Cone(Point(0, 0, -1), Point(0, 0, 1), 1., .5)

g3d = box + cone - sphere

# Test printing
info("\nCompact output of 3D geometry:")
info(g3d)
info("\nVerbose output of 3D geometry:")
info(g3d, True)

# Plot geometry
plot(g3d, "3D geometry (surface)")

# Generate and plot mesh
mesh3d = Mesh(g3d, 32)
info(mesh3d)
plot(mesh3d, "3D mesh")

interactive()