Moving towards 3D

My initial plan was to develop the 2D solver first with sufficient features to make a business case before even attempting to get myself into the difficulties of 3D geometry. Especially because I was slow to realise that GMSH has a full CAD engine inside of it.

I couldn’t shake the the urge to see if 3D would be possible and how it would be so tonight I have been experimenting with GMSH to get myself a bit acquainted anted with it. The biggest reason is that the success of this entire project could fall apart if 3D geometry and the associated difficulties aren’t solved. Why even continue if its not possible to realistically make 3D work?

GMSH’s 3D support

I knew boolean geometry worked in Shapely and always thought 3D geometry was going to be a huge summit to climb. Then I looked a bit better in the GMSH tutorials and discovered that they included the openCASCADE CAD kernal. Great! This actually means there is ample reason to cut out Shapely from EMerge entirely.

The hardest parts of boolean geomtry still apply. For those who are lucky enough to have never tried, there is a reason that making CAD kernals is an absolute nightmare. It is in fact so difficult to do that it is my firm believe that the limited amount of FEM solvers on the market is more due to the nightmar of 3D CAD software than actually implementing the physics of the Finite Element Method.

In theory the idea of geometric boolean operations is very elegant. You can define shapes like boxes, cyllinders and spheres and then subtract them and add them as if they where set-operations.

What makes it so hard is that while these are in fact, set operations, the elemtents of the sets are continuous spaces that have a topology defined on them which means that you can’t consider them als singular objects, ‘sets’ and instead as collections of distinct bounded subsets of R3. What this all boils down to for the programmer is that where usually the union or difference of two sets are other sets, now they can be collections of sets with and without holes. In other words, the return products depend a lot on what you put into it. A union may produce one or more geometries and so may an intersection, or difference. Then domains may have holes in them or sub-volumes etc. That makes it very difficult to generalize operations on geometries.

Lets say I define two blocks in programming terms:

block1 = Box(args, kwargs)

block2 = Box(args, kwargs)

output = subtract(block1, block2)

The difficulty with this desired code implementation is that output can either be one or more objects and its not at all clear how many this would be. If i want to add another geometry to the output of the subtraction, how should I do that? Take the following example:

outputs = subtract(block1, block2)

block3 = Box(args, kwargs)

output2 = [union(subpart, block3) for subpart in outputs]

This cannot work as the union tool object block3 now exists multiple times and on top of that, the objects in outputs unioninzed with block3 may form one new object all together.

The generic solution would be to then make sure that boolean operation operate on collections of domains of which the inputs can just be single element lists. Ok sure, this sounds good in principle but this only solves one problem. Its also possible that intersections of multiple objects return objects of a lower dimension. On top of that it may return no objects at all and the union of two objects that should touch may no longer touch due to floating point rounding errors.

Current implementation

I am currently working on a class structure that differentiates between two different contexts for geometries and domains (even though in GMSH they are not distinct). The first is a set of 3D volumes called GMSHVolume objects that generalize domains in GMSH and mange the tags that refer to them. The available methods will call the openCASCADE methods on these objects to generate new GMSHVolume objects from them (or a list of them).

These objects are then flattened() again to generate a final simulation domain.

Within these domains you have, again, a set of domains, boundaries, edges and points which are the constituent parts of the GMSH model.

A domain here is very similar to the GMSHVolume in terms of what it is but is different in the sense that it is final and disjoint from other domains. A Domain object in that sense is directly associated to something GMSH meshes instead of an object that may be deformed or unionized.

Domain objects will thus also directly relate to the resultant mesh.

I am not entirely convinced this is the best way forward but for know I will keep it this way.

A code example and the resultant GMSH mesh is shown below.

The code generating some arbitrary geometry

The output mesh generated from the EMerge input commands

Concluding remarks

A lot of the difficulties associated with these boolean operations is managed by the openCASCADE kernel. The challenge for EMerge will be to offer an interface that allows users to profit off of the features offered by GMSH without having to know how to make CAD objects with GMSH. Maybe at some point I will just interface with the GMSH geometry coding language. However, an important aspect of the FEM solver is interacting with these domains when defining the physics and materials. Thus having the entire geometry definition in Python as well helps with the link between the two especially when leveraging the auto-complete features of the IDE.