allow the specification of arbitrary IGES files as boundaries for your
geometries.
-
<a name="Intro"></a>
<h1>Introduction</h1>
-In some of the previous tutorial programs (step-1 and step-5 among others)
-we have seen that the deal.II mesh refinement is able to
--mesh refinement -> rule to refine boundary cells
--set_boundary
--so far only elementary shapes were considered
+In some of the previous tutorial programs (step-1, step-3, step-5 and step-6 among others)
+we have learned how to use the mesh refinement methods provided in deal.II. The 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.
+So far, only elementary boundary shapes such as circles, spheres and boxes were available in the
+library. In this tutorial, we will show how to use a new 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.
+
+
+<h3> CAD surfaces </h3>
+
+In the most common industrial practice, the geometrical models of arbitrary 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.
+By 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.
+Thus, 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.
+
+<h3> The boundary projector classes </h3>
+
+To import and interrogate CAD models, in the current deal.II release we developed a series of
+classes that wrap functions of the OpenCASCADE open source library for CAD modeling. Once imported
+from an IGES file, the model is stored in a TopoDS_Shape, which is the generic topological entity
+defined into the OpenCASCADE framework. From a TopoDS_Shape, it is then possible to access to 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 an object
+of the Boundary class. 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 new functions of the CAD modeling interface have been designed to retain the same structure. Hence, we build a
+projector object using the imported CAD shapes, and using the
+very same procedure we were previously using with HyperBallBoundary we assign 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
+on the specified geometry. But differently from a spherical or circular boundary, a boundary with a complex geometry
+poses proplems as to where exactly place on the prescribed shape the new nodes created upon refinement.
+HyperBallBoundary first creates the new nodes on the cell to be refined by averaging the
+surrounding points in the same way as FlatManifold does. Successively, it goes on to project such nodes on the circle or
+sphere along the radial direction. On such geometry, the radial direction ensures that the newly generated nodes remain evenly
+spaced across the refinement levels required.
+In the case of an arbitrary and complex shape though, the correct direction of projection cannot be identified that easily.
+That is why we have implemented several different projector classes that employ different projection strategies. A first
+projector, the ArclengthProjectionLineManifold, is to be used only for edge refinement. It is built assigning it a TopoDS_Edge,
+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.
+
+<img src="http://www.dealii.org/images/steps/developer/step-54.CurveSplit.png" alt="" width="500">
+
+
+A different projection stategy has been implemented in the NormalProjectionBoundary. The TopoDS_Shape assigned in construction can be arbitrary
+(a collection of shapes, faces, edges or a single face or edge willall work). The new cell nodes are first computed by averaging the
+surrounding points in the same way as FlatManifold does. In a second moment 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 closest point ---typically lying on the shape boundary---
+to the original one is selected. In cases instead where more than one normal projection is available, the closest of them to the original
+point is again selected. If the shape is composed by several sub-shapes, the projection is carried out onto every single
+sub-shape, and the closest projection point point is selected.
+
+<img src="http://www.dealii.org/images/steps/developer/step-54.NormalProjectionEdge.png" alt="" width="500">
+<img src="http://www.dealii.org/images/steps/developer/step-54.NormalProjection.png" alt="" width="500">
+
+As we are about to experience, in presence of arbitrary shapes setting the projection direction as that normal to the CAD surface will not
+lead to surface mesh elements of suitable quality. This 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. For this reason we implemented the DirectionalProjectionBoundary, which
+is constructed assigning a TopoDS_Shape (containing at least a face) and a direction along which all the projections will be carried out.
+In this way, the user will have a higher control on the projection direction to be enforced to ensure good mesh quality.
+
+<img src="http://www.dealii.org/images/steps/developer/step-54.DirectionalProjection.png" alt="" width="500">
+
+
+Of course the latter approach results 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, in the NormalToMeshProjectionBoundary class we also implemented a third
+projection algorithm. The NormalToMeshProjectionBoundary is built assigning a TopoDS_Shape (containing at least a face) to the contructor,
+and works exactly like a DirectionalProjection. But, as the name of the class suggests, 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 inpired to 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 found
+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 here considered is somehow perpendicular to the horizontal direction, we notice that 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 such region with a directional projection, if 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.
+
--this tutorial shows how to refine a mesh on top of a complex surface
-generated with a CAD modeler and saved in an IGES file.
+<img src="http://www.dealii.org/images/steps/developer/step-54.ProjectionComparisons.png" alt="" width="700">
-<h2>Where geometry and meshes intersect</h2>
<h1>Results</h1>
-Running the program produces a series of mesh files <code>mesh.vtu</code>
-that we can
-visualize with any of the usual visualization programs that can read the VTK
+The program execution produces a series of mesh files <code>3d_mesh_*.vtk</code>
+that we can visualize with any of the usual visualization programs that can read the VTK
file format.
-<img src="http://www.dealii.org/images/steps/developer/step-54.mesh.png" alt="">
-<img src="http://www.dealii.org/images/steps/developer/step-54.mesh.png" alt="">
-<img src="http://www.dealii.org/images/steps/developer/step-54.mesh.png" alt="">
-<img src="http://www.dealii.org/images/steps/developer/step-54.mesh.png" alt="">
\ No newline at end of file
+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. The last row
+of pictures depict front views of the same grids shown in the second row.
+
+
+<table style="width:90%">
+ <tr>
+ <td><img src="http://www.dealii.org/images/steps/developer/step-54.common_0.png" alt="" width="500"></td>
+ <td><img src="http://www.dealii.org/images/steps/developer/step-54.normal_1.png" alt="" width="500"></td>
+ <td><img src="http://www.dealii.org/images/steps/developer/step-54.normal_2.png" alt="" width="500"></td>
+ </tr>
+ <tr>
+ <td><img src="http://www.dealii.org/images/steps/developer/step-54.normal_3.png" alt="" width="500"></td>
+ <td><img src="http://www.dealii.org/images/steps/developer/step-54.normal_4.png" alt="" width="500"></td>
+ <td><img src="http://www.dealii.org/images/steps/developer/step-54.normal_5.png" alt="" width="500"></td>
+ </tr>
+ <tr>
+ <td><img src="http://www.dealii.org/images/steps/developer/step-54.normal_front_3.png" alt="" width="500"></td>
+ <td><img src="http://www.dealii.org/images/steps/developer/step-54.normal_front_4.png" alt="" width="500"></td>
+ <td><img src="http://www.dealii.org/images/steps/developer/step-54.normal_front_5.png" alt="" width="500"></td>
+ </tr>
+</table>
+
+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.
+
+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).
+
+
+<table style="width:90%">
+ <tr>
+ <td><img src="http://www.dealii.org/images/steps/developer/step-54.common_0.png" alt="" width="500"></td>
+ <td><img src="http://www.dealii.org/images/steps/developer/step-54.directional_1.png" alt="" width="500"></td>
+ <td><img src="http://www.dealii.org/images/steps/developer/step-54.directional_2.png" alt="" width="500"></td>
+ </tr>
+ <tr>
+ <td><img src="http://www.dealii.org/images/steps/developer/step-54.directional_3.png" alt="" width="500"></td>
+ <td><img src="http://www.dealii.org/images/steps/developer/step-54.directional_4.png" alt="" width="500"></td>
+ <td><img src="http://www.dealii.org/images/steps/developer/step-54.directional_5.png" alt="" width="500"></td>
+ </tr>
+ <tr>
+ <td><img src="http://www.dealii.org/images/steps/developer/step-54.directional_front_3.png" alt="" width="500"></td>
+ <td><img src="http://www.dealii.org/images/steps/developer/step-54.directional_front_4.png" alt="" width="500"></td>
+ <td><img src="http://www.dealii.org/images/steps/developer/step-54.directional_front_5.png" alt="" width="500"></td>
+ </tr>
+</table>
+
+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 carried out results instead in the meshes illustrated in the following table.
+
+<table style="width:90%">
+ <tr>
+ <td><img src="http://www.dealii.org/images/steps/developer/step-54.common_0.png" alt="" width="500"></td>
+ <td><img src="http://www.dealii.org/images/steps/developer/step-54.normal_to_mesh_1.png" alt="" width="500"></td>
+ <td><img src="http://www.dealii.org/images/steps/developer/step-54.normal_to_mesh_2.png" alt="" width="500"></td>
+ </tr>
+ <tr>
+ <td><img src="http://www.dealii.org/images/steps/developer/step-54.normal_to_mesh_3.png" alt="" width="500"></td>
+ <td><img src="http://www.dealii.org/images/steps/developer/step-54.normal_to_mesh_4.png" alt="" width="500"></td>
+ <td><img src="http://www.dealii.org/images/steps/developer/step-54.normal_to_mesh_5.png" alt="" width="500"></td>
+ </tr>
+ <tr>
+ <td><img src="http://www.dealii.org/images/steps/developer/step-54.normal_to_mesh_front_3.png" alt="" width="500"></td>
+ <td><img src="http://www.dealii.org/images/steps/developer/step-54.normal_to_mesh_front_4.png" alt="" width="500"></td>
+ <td><img src="http://www.dealii.org/images/steps/developer/step-54.normal_to_mesh_front_5.png" alt="" width="500"></td>
+ </tr>
+</table>
+
+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.
+
+
// @sect3{Include files}
-// The program starts with including a bunch of include files that we will use
+// We start with including a bunch of files that we will use
// in the various parts of the program. Most of them have been discussed in
// previous tutorials already:
-#include <deal.II/base/smartpointer.h>
-#include <deal.II/base/convergence_table.h>
-#include <deal.II/base/quadrature_lib.h>
-#include <deal.II/base/quadrature_selector.h>
-#include <deal.II/base/parsed_function.h>
-#include <deal.II/base/utilities.h>
-
-#include <deal.II/lac/full_matrix.h>
-#include <deal.II/lac/vector.h>
-#include <deal.II/lac/solver_control.h>
-#include <deal.II/lac/solver_gmres.h>
-#include <deal.II/lac/precondition.h>
-
#include <deal.II/grid/tria.h>
#include <deal.II/grid/tria_iterator.h>
#include <deal.II/grid/tria_accessor.h>
#include <deal.II/grid/grid_in.h>
#include <deal.II/grid/grid_out.h>
#include <deal.II/grid/tria_boundary_lib.h>
+#include <deal.II/numerics/data_out.h>
+#include <deal.II/numerics/vector_tools.h>
// And here are the headers of the opencascade support classes and functions:
#include <deal.II/opencascade/boundary_lib.h>
#include <deal.II/opencascade/utilities.h>
-#include <deal.II/dofs/dof_handler.h>
-#include <deal.II/dofs/dof_accessor.h>
-#include <deal.II/dofs/dof_tools.h>
-
-#include <deal.II/fe/fe_q.h>
-#include <deal.II/fe/fe_values.h>
-#include <deal.II/fe/mapping_q.h>
-#include <deal.II/numerics/data_out.h>
-#include <deal.II/numerics/vector_tools.h>
-
-// And here are a few C++ standard header files that we will need:
+// Finally, a few C++ standard header files that we will need:
#include <cmath>
#include <iostream>
#include <fstream>
// The structure of this class is very small. Since we only want
// to show how a triangulation can be refined onto a CAD surface, the
// arguments of this class are basically just the input and output file
- // names and the triangulation we want to play with.
- // The member functions of this class are like those that in most of the
+ // names, and a flag indicating the surface projection strategy we want to
+ // test. Along with the input argument, the only other member of the class
+ // is the triangulation we want to play with.
+ // The member functions of this class are similar to those that in most of the
// other tutorial programs deal with the setup of the grid for the
// simulations.
// @sect4{TriangulationOnCAD::read_domain}
- // The following function represents the core of the present tutorial program.
+ // The following function represents the core of thhis program.
// In this function we in fact import the CAD shape upon which we want to generate
// and refine our triangulation. Such CAD surface is contained in the IGES
// file "DTMB-5415_bulbous_bow.iges", and represents the bulbous bow of a ship.
// geometry a particularly meaningful example.
//
// So, after importing the hull bow surface, we extract some of the curves and surfaces
- // comosing it, and use them to generate a set of projectors. Such projectors substantially
- // define the rules deal.ii has to follow to position each new node during the cell
+ // composing it, and use them to generate a set of projectors. Such projectors substantially
+ // define the rules the Triangulation has to follow to position each new node during the cell
// refinement.
//
- // As for the triangulation, as done in previous tutorial programs, we import a
- // pre-existing grid saved in .vtk format. The imported mesh is composed of a single
- // quadrilateral cell the vertices of which have been located on the CAD shape. In
- // this tutorial, we chose to import our mesh in vtk format.
+ // To initialize the Triangulation, as done in previous tutorial programs, we import a
+ // pre-existing grid saved in VTK format. The imported mesh is composed of a single
+ // quadrilateral cell the vertices of which have been placed on the CAD shape.
//
- // So, after importing both the initial mesh, we assign the projectors
- // previously generated to each of the edges and cells which will have to be
+ // So, after importing both the IGES geometry and the initial mesh, we assign the projectors
+ // previously discussed to each of the edges and cells which will have to be
// refined on the CAD surface.
//
- // In this tutorial, we will test three different ways to project new mesh nodes onto
- // the CAD surface, and will analyze the results obtained with each surface projection
- // strategy. A first approach consists in projecting each node in the direction which
- // is normal to the surface. A second possibility is represented by chosing a single
- // direction along which to project all the nodes on the surface. The third strategy
- // consists in projecting the new nodes on the surface along a direction which represents
- // an estimate of the mesh cell normal. Each of such projection strategies has been
+ // In this tutorial, we will test the three different CAD surface projectors described in
+ // the introduction, and will analyze the results obtained with each of them.
+ // As mentioned, each of such projection strategiy has been
// implemented in a different class, which can be assigned to the set_manifold method
- // of a deal.ii triangulation class.
+ // of a Triangulation class.
//
void TriangulationOnCAD::read_domain()
{
- // this function allows for the CAD file of interest (in IGES format) to be imported.
- // The function input parameters are a string containing the desired file name, and
+ // The following function allows for the CAD file of interest (in IGES format) to be imported.
+ // The function arguments are a string containing the desired file name, and
// a scale factor. In this example, such scale factor is set to 1e-3, as the original
// geometry is written in millimeters, while we prefer to work in meters.
- // The output of the function is an object of opencascade generic topological shape
+ // The output of the function is an object of OpenCASCADE generic topological shape
// class, namely a TopoDS_Shape.
TopoDS_Shape bow_surface = OpenCASCADE::read_IGES("DTMB-5415_bulbous_bow.iges",1e-3);
wires);
// The next few steps are more familiar, and allow us to import an existing
- // mesh from an external vtk file, and convert it to a deal triangulation.
+ // mesh from an external VTK file, and convert it to a deal triangulation.
std::ifstream in;
in.open(initial_mesh_filename.c_str());
// We output this initial mesh saving it as the refinement step 0.
output_results(0);
- // The mesh imported has a single cell. so, we get an iterator to that cell.
+ // The mesh imported has a single cell. So, we get an iterator to that cell.
// and assgin it the manifold_id 1
Triangulation<2,3>::active_cell_iterator cell = tria.begin_active();
cell->set_manifold_id(1);
switch (surface_projection_kind)
{
case 0:
- // If the value is 0, we select the NormalProjectionBoundary. The new mesh points will initially
- // generated at the baricenter of the cell/edge considere, and then projected
+ // If surface_projection_kind value is 0, we select the NormalProjectionBoundary. The new mesh points will initially
+ // generated at the baricenter of the cell/edge considered, and then projected
// on the CAD surface along its normal direction.
// The NormalProjectionBoundary constructor only needs a shape and a tolerance.
static OpenCASCADE::NormalProjectionBoundary<2,3> normal_projector(bow_surface, tolerance);
- // The normal projector is assigned to the manifold having id 1.
+ // Once created, the normal projector is assigned to the manifold having id 1.
tria.set_manifold(1,normal_projector);
break;
case 1:
- // If the value is 1, we select the DirectionalProjectionBoundary. The new mesh points will initially
+ // If surface_projection_kind value is 1, we select the DirectionalProjectionBoundary. The new mesh points will initially
// generated at the baricenter of the cell/edge considere, and then projected
// on the CAD surface along a direction that is specified to the DirectionalProjectionBoundary
// constructor. In this case, the projection is done along the y-axis.
tria.set_manifold(1,directional_projector);
break;
case 2:
- // If the value is 2, we select the NormaToMeshlProjectionBoundary. The new mesh points will initially
+ // If surface_projection_kind value is 2, we select the NormaToMeshlProjectionBoundary. The new mesh points will initially
// generated at the baricenter of the cell/edge considere, and then projected
// on the CAD surface along a direction that is an estimate of the mesh normal direction.
// The NormalToMeshProjectionBoundary constructor only requires a shape (containing at least a face)
// This function globally refines the mesh. In other tutorials, it tipically also distributes degrees
// of freedom, and resizes matrices and vectors. These tasks are not carried out
- // here, since we are not running any simulation on the triangolation produced.
+ // here, since we are not running any simulation on the Triangulation produced.
void TriangulationOnCAD::refine_and_resize()
* The standard constructor takes a generic TopoDS_Shape @p sh,
* and a tolerance used to compute distances internally.
*
- * The TopoDS_Shape can be of arbitrary, i.e., a collection of
+ * The TopoDS_Shape can be arbitrary, i.e., a collection of
* shapes, faces, edges or a single face or edge.
*/
NormalProjectionBoundary(const TopoDS_Shape &sh,
projects / dealii.git / commitdiff