From 6a2b90cfa8b4ccf41fa1b605f20d35d567585d5a Mon Sep 17 00:00:00 2001 From: Andrea Mola Date: Tue, 21 Oct 2014 16:59:20 +0200 Subject: [PATCH] comments added to step-54 more comments added to step-54. and docs removed some files some improvements to step-54 comments, intro and results some more improvements to step-54 comments, intro and results pictures path in html files of step-54 is now changed to its final location (pictures will need to be uploaded on the deal website) --- examples/step-54/doc/intro.dox | 113 +++++++++++++++++++-- examples/step-54/doc/results.dox | 87 ++++++++++++++-- examples/step-54/step-54.cc | 86 ++++++---------- include/deal.II/opencascade/boundary_lib.h | 2 +- 4 files changed, 216 insertions(+), 72 deletions(-) diff --git a/examples/step-54/doc/intro.dox b/examples/step-54/doc/intro.dox index e8b6366ebe..b18b820c45 100644 --- a/examples/step-54/doc/intro.dox +++ b/examples/step-54/doc/intro.dox @@ -7,19 +7,114 @@ tools that interface with OpenCASCADE (http://www.opencascade.org) that allow the specification of arbitrary IGES files as boundaries for your geometries. -

Introduction

-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. + + +

CAD surfaces

+ +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. + +

The boundary projector classes

+ +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. + + + + +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. + + + + +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. + + + + +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. + -

Where geometry and meshes intersect

diff --git a/examples/step-54/doc/results.dox b/examples/step-54/doc/results.dox index ef0e3651db..1e975d7386 100644 --- a/examples/step-54/doc/results.dox +++ b/examples/step-54/doc/results.dox @@ -1,11 +1,84 @@

Results

-Running the program produces a series of mesh files mesh.vtu -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 3d_mesh_*.vtk +that we can visualize with any of the usual visualization programs that can read the VTK file format. - - - - \ 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. + + + + + + + + + + + + + + + + + + +
+ +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). + + + + + + + + + + + + + + + + + + +
+ +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. + + + + + + + + + + + + + + + + + +
+ +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. + + diff --git a/examples/step-54/step-54.cc b/examples/step-54/step-54.cc index 2c030087ea..24be8b69fb 100644 --- a/examples/step-54/step-54.cc +++ b/examples/step-54/step-54.cc @@ -19,22 +19,9 @@ // @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 -#include -#include -#include -#include -#include - -#include -#include -#include -#include -#include - #include #include #include @@ -42,23 +29,15 @@ #include #include #include +#include +#include // And here are the headers of the opencascade support classes and functions: #include #include -#include -#include -#include - -#include -#include -#include -#include -#include - -// 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 #include #include @@ -77,8 +56,10 @@ namespace Step54 // 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. @@ -138,7 +119,7 @@ namespace Step54 // @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. @@ -146,28 +127,23 @@ namespace Step54 // 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. // @@ -176,11 +152,11 @@ namespace Step54 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); @@ -221,7 +197,7 @@ namespace Step54 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()); @@ -233,7 +209,7 @@ namespace Step54 // 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); @@ -262,16 +238,16 @@ namespace Step54 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. @@ -279,7 +255,7 @@ namespace Step54 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) @@ -299,7 +275,7 @@ namespace Step54 // 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() diff --git a/include/deal.II/opencascade/boundary_lib.h b/include/deal.II/opencascade/boundary_lib.h index cffff0fe79..9d0af6833c 100644 --- a/include/deal.II/opencascade/boundary_lib.h +++ b/include/deal.II/opencascade/boundary_lib.h @@ -69,7 +69,7 @@ namespace OpenCASCADE * 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, -- 2.39.5