This program was contributed by Andrea Mola and Luca Heltai. @note This program elaborates on concepts of industrial geometry, using tools that interface with the OpenCASCADE library (http://www.opencascade.org) that allow the specification of arbitrary IGES files to describe the boundaries for your geometries.

Installing and Running

If you have a recent version of deal.II installed, you can configure, compile, and run using the following commands: $ cmake -D DEAL_II_DIR=/path/to/dealii . $ make $ make run

Introduction

In some of the previous tutorial programs (step-1, step-3, step-5, step-6 and step-49 among others) we have learned how to use the mesh refinement methods provided in deal.II. These tutorials have shown how to employ such tools to produce a fine grid for a single simulation, as done in step-3; or to start from a coarse grid and carry out a series of simulations on adaptively refined grids, as is the case of step-6. Regardless of which approach is taken, the mesh refinement requires a suitable geometrical description of the computational domain boundary in order to place, at each refinement, the new mesh nodes onto the boundary surface. For instance, step-5 shows how to assign a circular shape to the boundary of the computational domain, so that the faces lying on the boundary are refined onto the circle. step-53 shows how to do this with a boundary defined by experimentally obtained data. But, at least as far as elementary boundary shapes are concerned, deal.II really only provides circles, spheres, boxes and various combinations. In this tutorial, we will show how to use a set of classes developed to import arbitrary CAD geometries, assign them to the desired boundary of the computational domain, and refine a computational grid on such complex shapes.

CAD surfaces

In the most common industrial practice, the geometrical models of arbitrarily shaped objects are realized by means of Computer Aided Design (CAD) tools. The use of CAD modelers has spread in the last decades, as they allow for the generation of a full virtual model of each designed object, which through a computer can be visualized, inspected, and analyzed in its finest details well before it is physically crafted. From a mathematical perspective, the engine lying under the hood of CAD modelers is represented by analytical geometry, and in particular by parametric curves and surfaces such as B-splines and NURBS that are rich enough that they can represent most surfaces of practical interest. Once a virtual model is ready, all the geometrical features of the desired object are stored in files which materially contain the coefficients of the parametric surfaces and curves composing the object. Depending on the specific CAD tool used to define the geometrical model, there are of course several different file formats in which the information of a CAD model can be organized. To provide a common ground to exchange data across CAD tools, the U.S. National Bureau of Standards published in 1980 the Initial Graphics Exchange Representation (IGES) neutral file format, which is used in this example.

The boundary projector classes

To import and interrogate CAD models, the deal.II library contains a series of wrapper functions for the OpenCASCADE open source library for CAD modeling. These functions allow to import IGES files into OpenCASCADE native objects, and wrap them inside a series of Manifold classes. Once imported from an IGES file, the model is stored in a TopoDS_Shape, which is the generic topological entity defined in the OpenCASCADE framework. From a TopoDS_Shape, it is then possible to access all the sub-shapes (such as vertices, edges and faces) composing it, along with their geometrical description. In the deal.II framework, the topological entities composing a shape are used to create objects of the Manifold or Boundary classes. In Step-6 we saw how to build a HyperBallBoundary and assign it to a set of faces (or cells, for co-dimension 1) of a Triangulation, to have cells and faces refined on a sphere or circle. The functions of the CAD modeling interface have been designed to retain the same structure, allowing the user to build a projector object using the imported CAD shapes, maintaining the very same procedure we use with HyperBallBoundary, i.e., assigning such projector object to cells, faces or edges of a coarse mesh. At each refinement cycle, the new mesh nodes will be then automatically generated by projecting a midpoint of an existing object onto the specified geometry. Differently from a spherical or circular boundary, a boundary with a complex geometry poses problems as to where it is best to place the new nodes created upon refinement on the prescribed shape. HyperBallBoundary first creates the new nodes on the face or edge to be refined by averaging the surrounding points in the same way as FlatManifold does. Then, it goes on to project such nodes on the circle or sphere along the radial direction. On such a geometry, the radial direction ensures that the newly generated nodes remain evenly spaced when remaining on a given refinement level. In the case of an arbitrary and complex shape though, the best direction of the projection cannot be identified that easily. The OpenCASCADE wrappers in deal.II provide several projector classes that employ different projection strategies. A first projector, implemented in the OpenCASCADE::ArclengthProjectionLineManifold class, is to be used only for edge refinement. It is built assigning it a topological shape of dimension one, either a TopoDS_Edge or a TopoDS_Wire (which is a compound shape, made of several connected TopoDS_Edges) and refines a mesh edge finding the new vertex as the point splitting in two even parts the curvilinear length of the CAD curve portion that lies between the vertices of the original edge.

A different projection strategy has been implemented in the OpenCASCADE::NormalProjectionBoundary class. The TopoDS_Shape assigned at construction time can be arbitrary (a collection of shapes, faces, edges or a single face or edge will all work). The new cell nodes are first computed by averaging the surrounding points in the same way as FlatManifold does. In a second step, all the new nodes will be projected onto the TopoDS_Shape along the direction normal to the shape. If no normal projection is available, the point which is closest to the shape---typically lying on the shape boundary---is selected. If the shape is composed of several sub-shapes, the projection is carried out onto every single sub-shape, and the closest projection point point is selected.

As we are about to experience, for some shapes, setting the projection direction as that normal to the CAD surface will not lead to surface mesh elements of suitable quality. This is because the direction normal to the CAD surface has in principle nothing to do with the direction along which the mesh needs the new nodes to be located. The OpenCASCADE::DirectionalProjectionBoundary class, in this case, can help. This class is constructed assigning a TopoDS_Shape (containing at least a face) and a direction along which all the projections will be carried out. New points will be computed by first averaging the surrounding points (as in the FlatManifold case), and then taking the closest intersection between the topological shape and the line passing through the resulting point, along the direction used at construction time. In this way, the user will have a higher control on the projection direction to be enforced to ensure good mesh quality.

Of course the latter approach is effective only when the orientation of the surface is rather uniform, so that a single projection direction can be identified. In cases in which the surface direction is approaching the projection direction, it is even possible that the directional projection is not found. To overcome these problems, the OpenCASCADE::NormalToMeshProjectionBoundary class implements a third projection algorithm. The OpenCASCADE::NormalToMeshProjectionBoundary class is built assigning a TopoDS_Shape (containing at least one face) to the constructor, and works exactly like a OpenCASCADE::DirectionalProjection. But, as the name of the class suggests, OpenCASCADE::NormalToMeshProjectionBoundary tries to come up with a suitable estimate of the direction normal to the mesh elements to be refined, and uses it for the projection of the new nodes onto the CAD surface. If we consider a mesh edge in a 2D space, the direction of its axis is a direction along which to split it in order to give rise to two new cells of the same length. We here extended this concept in 3D, and project all new nodes in a direction that approximates the cell normal. In the next figure, which is inspired by the geometry considered in this tutorial, we make an attempt to compare the behavior of the three projectors considered. As can be seen on the left, given the original cell (in blue), the new point found with the normal projection is in a position which does not allow for the generation of evenly spaced new elements (in red). The situation will get worse in further refinement steps. Since the geometry we considered is somehow perpendicular to the horizontal direction, the directional projection (central image) defined with horizontal direction as the projection direction, does a rather good job in getting the new mesh point. Yet, since the surface is almost horizontal at the bottom of the picture, we can expect problems in those regions when further refinement steps are carried out. Finally, the picture on the right shows that a node located on the cell axis will result in two new cells having the same length. Of course the situation in 3D gets a little more complicated than that described in this simple 2D case. Nevertheless, the results of this test confirm that the normal to the mesh direction is the best approach among the three tested, when arbitrarily shaped surfaces are considered, and unless you have a geometry for which a more specific approach is known to be appropriate.

The testcase

In this program, we will consider creating a surface mesh for a real geometry describing the bow of a ship (this geometry is frequently used in CAD and mesh generation comparisons and is freely available). The surface mesh we get from this could then be used to solve a boundary element equation to simulate the flow of water around the ship (in a way similar to step-34) but we will not try to do this here. To already give you an idea of the geometry we consider, here is a picture:

In the program, we read both the geometry and a coarse mesh from files, and then employ several of the options discussed above to place new vertices for a sequence of mesh refinement steps.

Results

The program execution produces a series of mesh files 3d_mesh_*.vtk that we can visualize with any of the usual visualization programs that can read the VTK file format. The following table illustrates the results obtained employing the normal projection strategy. The first two rows of the table show side views of the grids obtained for progressive levels of refinement, overlain on a very fine rendering of the exact geometry. The dark and light red areas simply indicate whether the current mesh or the fine geometry is closer to the observer; the distinction does not carry any particularly deep meaning. The last row of pictures depict front views (mirrored to both sides of the geometry) of the same grids shown in the second row.

As can be seen in the pictures---and as we anticipated---the normal refinement strategy is unable to produce nicely shaped elements when applied to surfaces with significant curvature changes. This is particularly apparent at the bulb of the hull where all new points have been placed in the upper part of the bulb and the lower part remains completely unresolved. The following table, which is arranged as the previous one, illustrates the results obtained adopting the directional projection approach, in which the projection direction selected was the y-axis (which is indicated with a small yellow arrow at the bottom left of each image).

The images confirm that the quality of the mesh obtained with a directional projection is sensibly higher than that obtained projecting along the surface normal. Yet, a number of elements elongated in the y-direction are observed around the bottom of the bulb, where the surface is almost parallel to the direction chosen for the projection. The final test shows results using instead the projection normal to the faces:

The pictures confirm that the normal to mesh projection approach leads to grids that remain evenly spaced throughtout the refinement steps. At the same time, these meshes represent rather well the original geometry even in the bottom region of the bulb, which is not well recovered employing the directional projector or the normal projector.