]> https://gitweb.dealii.org/ - dealii.git/commitdiff
step-54 is now completed (it only needs comments). three different projectors have...
authorAndrea Mola <mola.andrea@gmail.com>
Tue, 14 Oct 2014 15:45:03 +0000 (17:45 +0200)
committerLuca Heltai <luca.heltai@sissa.it>
Fri, 31 Oct 2014 10:07:09 +0000 (11:07 +0100)
Fixed indentation.

fixed some problems with residual calls to Pnt() (now Point()) and axis_ intersection() (now line_intersection()) functions in boundary_lib. files coarse_sphere and coarse_circle removed from step-54 directory. also parameters.prm has been removed, and n_cycles value is asigned in an hard coded way

all input and output files have been converted to .vtk

added comments to step-54.cc

fixed all the notes by luca exept for that on the closest point function

Fixed indentation.

examples/step-54/coarse_circle.inp [deleted file]
examples/step-54/coarse_sphere.inp [deleted file]
examples/step-54/initial_mesh_3d.inp [deleted file]
examples/step-54/initial_mesh_3d.vtk [new file with mode: 0644]
examples/step-54/parameters.prm [deleted file]
examples/step-54/step-54.cc
include/deal.II/opencascade/boundary_lib.h
include/deal.II/opencascade/utilities.h
source/opencascade/boundary_lib.cc
source/opencascade/boundary_lib.inst.in
source/opencascade/utilities.cc

diff --git a/examples/step-54/coarse_circle.inp b/examples/step-54/coarse_circle.inp
deleted file mode 100644 (file)
index dda97cd..0000000
+++ /dev/null
@@ -1,21 +0,0 @@
-10 10 0 0 0
-1      0 1.0 0
-2      0.587785252292 0.809016994375 0
-3      0.951056516295 0.309016994375 0
-4      0.951056516295 -0.309016994375 0
-5      0.587785252292 -0.809016994375 0
-6      0 -1.0 0
-7      -0.587785252292 -0.809016994375 0
-8      -0.951056516295 -0.309016994375 0
-9      -0.951056516295 0.309016994375 0
-10     -0.587785252292 0.809016994375 0
-1 1 line 1 2
-2 1 line 2 3
-3 1 line 3 4
-4 1 line 4 5
-5 1 line 5 6
-6 1 line 6 7
-7 1 line 7 8
-8 1 line 8 9
-9 1 line 9 10
-10 1 line 10 1
diff --git a/examples/step-54/coarse_sphere.inp b/examples/step-54/coarse_sphere.inp
deleted file mode 100644 (file)
index 5ea249c..0000000
+++ /dev/null
@@ -1,15 +0,0 @@
-8 6 0 0 0
-1      -0.577350269 -0.577350269 -0.577350269
-2       0.577350269 -0.577350269 -0.577350269
-3      -0.577350269  0.577350269 -0.577350269
-4       0.577350269  0.577350269 -0.577350269
-5      -0.577350269 -0.577350269  0.577350269
-6       0.577350269 -0.577350269  0.577350269
-7      -0.577350269  0.577350269  0.577350269
-8       0.577350269  0.577350269  0.577350269
-1 1 quad 3 4 2 1
-2 1 quad 5 6 8 7
-3 1 quad 1 2 6 5 
-4 1 quad 3 7 8 4
-5 1 quad 5 7 3 1
-6 1 quad 2 4 8 6
diff --git a/examples/step-54/initial_mesh_3d.inp b/examples/step-54/initial_mesh_3d.inp
deleted file mode 100644 (file)
index 8080f67..0000000
+++ /dev/null
@@ -1,11 +0,0 @@
-4 1 0 0 0
-1      -3.0595482 0 0.4043209
-2      -2.63758308 0 -0.34670038
-3      -2.010663 0 -0.249154631
-4      -2.010663 0.314450974 0.331320309
-1 1 quad 1 2 3 4
-2 2 line 1 2
-3 3 line 2 3
-4 4 line 3 4
-5 5 line 4 1
-
diff --git a/examples/step-54/initial_mesh_3d.vtk b/examples/step-54/initial_mesh_3d.vtk
new file mode 100644 (file)
index 0000000..ddf377e
--- /dev/null
@@ -0,0 +1,17 @@
+# vtk DataFile Version 3.0
+#This file was generated by the deal.II library on 2014/10/15 at 14:32:28
+ASCII
+DATASET UNSTRUCTURED_GRID
+
+POINTS 4 double
+-3.05955 0 0.404321
+-2.63758 0 -0.3467
+-2.01066 0.314451 0.33132
+-2.01066 0 -0.249155
+
+CELLS 1 5
+4      0       1       3       2
+
+CELL_TYPES 1
+ 9
+POINT_DATA 4
diff --git a/examples/step-54/parameters.prm b/examples/step-54/parameters.prm
deleted file mode 100644 (file)
index 534c102..0000000
+++ /dev/null
@@ -1,12 +0,0 @@
-# Listing of Parameters
-# ---------------------
-set Number of cycles                = 7
-
-
-subsection Quadrature rules
-  set Quadrature order          = 4
-  set Quadrature type           = gauss
-end
-
-
-
index 7fecec4a9e3a52178f055dcc960609eae74179c4..2c030087ea39a6e81f7c2244c08aee8634a25e8c 100644 (file)
@@ -43,6 +43,7 @@
 #include <deal.II/grid/grid_out.h>
 #include <deal.II/grid/tria_boundary_lib.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>
 
@@ -73,29 +74,28 @@ namespace Step54
 
   // @sect3{The TriangulationOnCAD class}
 
-  // The structure of a boundary element method code is very similar to the
-  // structure of a finite element code, and so the member functions of this
-  // class are like those of most of the other tutorial programs. In
-  // particular, by now you should be familiar with reading parameters from an
-  // external file, and with the splitting of the different tasks into
-  // different modules. The same applies to boundary element methods, and we
-  // won't comment too much on them, except on the differences.
+  // 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
+  // other tutorial programs deal with the setup of the grid for the
+  // simulations.
 
   class TriangulationOnCAD
   {
   public:
-    TriangulationOnCAD(const unsigned int fe_degree = 1,
-               const unsigned int mapping_degree = 1);
+    TriangulationOnCAD(const std::string &initial_mesh_filename,
+                       const std::string &output_filename,
+                       const unsigned int &surface_projection_kind = 0);
+
 
-    
     ~TriangulationOnCAD();
 
     void run();
 
   private:
 
-    void read_parameters (const std::string &filename);
-
     void read_domain();
 
     void refine_and_resize();
@@ -103,146 +103,115 @@ namespace Step54
     void output_results(const unsigned int cycle);
 
     Triangulation<2, 3>   tria;
-    FE_Q<2,3>             fe;
-    DoFHandler<2,3>       dh;
-    MappingQ<2,3>      mapping;
-
-    std_cxx11::shared_ptr<Quadrature<2> > quadrature;
-
-    SolverControl solver_control;
 
-    unsigned int n_cycles;
-
-    OpenCASCADE::ArclengthProjectionLineManifold<2,3> *line_projector;
-    OpenCASCADE::NormalProjectionBoundary<2,3> *normal_projector;
-    OpenCASCADE::DirectionalProjectionBoundary<2,3> *directional_projector;
+    const std::string &initial_mesh_filename;
+    const std::string &output_filename;
+    const unsigned int &surface_projection_kind;
 
   };
 
 
-  // @sect4{TriangulationOnCAD::TriangulationOnCAD and TriangulationOnCAD::read_parameters}
+  // @sect4{TriangulationOnCAD::TriangulationOnCAD }
 
-  // The constructor initializes the various object in much the same way as
-  // done in the finite element programs such as step-4 or step-6. The only
-  // new ingredient here is the ParsedFunction object, which needs, at
-  // construction time, the specification of the number of components.
-  //
-  // For the exact solution the number of vector components is one, and no
-  // action is required since one is the default value for a ParsedFunction
-  // object. The wind, however, requires dim components to be
-  // specified. Notice that when declaring entries in a parameter file for the
-  // expression of the Functions::ParsedFunction, we need to specify the
-  // number of components explicitly, since the function
-  // Functions::ParsedFunction::declare_parameters is static, and has no
-  // knowledge of the number of components.
-
-  TriangulationOnCAD::TriangulationOnCAD(const unsigned int fe_degree,
-                              const unsigned int mapping_degree)
-    :
-    fe(fe_degree),
-    dh(tria),
-    mapping(mapping_degree, true)
-  {}
+  // The constructor of the TriangulationOnCAD class is very simple.
+  // The input arguments are strings for the input and output file
+  // names, and an unsigned int flag which can only assume values
+  // 0,1,2 and determines which kind of surface projector is used in
+  // the mesh refinement cycles (see below for details).
 
-  TriangulationOnCAD::~TriangulationOnCAD()
+  TriangulationOnCAD::TriangulationOnCAD(const std::string &initial_mesh_filename,
+                                         const std::string &output_filename,
+                                         const unsigned int &surface_projection_kind)
+    :
+    initial_mesh_filename(initial_mesh_filename),
+    output_filename(output_filename),
+    surface_projection_kind(surface_projection_kind)
   {
-  tria.set_manifold(1);
-  tria.set_manifold(2);
-  delete line_projector;
-  delete normal_projector;
-  delete directional_projector;
   }
 
-  void TriangulationOnCAD::read_parameters (const std::string &filename)
+  TriangulationOnCAD::~TriangulationOnCAD()
   {
-    deallog << std::endl << "Parsing parameter file " << filename << std::endl
-            << "for a three dimensional geometry. " << std::endl;
-
-    ParameterHandler prm;
-
-    prm.declare_entry("Number of cycles", "4",
-                      Patterns::Integer());
-
-    prm.enter_subsection("Quadrature rules");
-    {
-      prm.declare_entry("Quadrature type", "gauss",
-                        Patterns::Selection(QuadratureSelector<(2)>::get_quadrature_names()));
-      prm.declare_entry("Quadrature order", "4", Patterns::Integer());
-    }
-    prm.leave_subsection();
-
-    // After declaring all these parameters to the ParameterHandler object,
-    // let's read an input file that will give the parameters their values. We
-    // then proceed to extract these values from the ParameterHandler object:
-    prm.read_input(filename);
-
-    n_cycles = prm.get_integer("Number of cycles");
-
-    prm.enter_subsection("Quadrature rules");
-    {
-      quadrature =
-        std_cxx11::shared_ptr<Quadrature<2> >
-        (new QuadratureSelector<2> (prm.get("Quadrature type"),
-                                        prm.get_integer("Quadrature order")));
-    }
-    prm.leave_subsection();
-
   }
 
 
+
   // @sect4{TriangulationOnCAD::read_domain}
 
 
-  // Some of the mesh formats supported in deal.II use by default three
-  // dimensional points to describe meshes. These are the formats which are
-  // compatible with the boundary element method capabilities of deal.II. In
-  // particular we can use either UCD or GMSH formats. In both cases, we have
-  // to be particularly careful with the orientation of the mesh, because,
-  // unlike in the standard finite element case, no reordering or
-  // compatibility check is performed here.  All meshes are considered as
-  // oriented, because they are embedded in a higher dimensional space. (See
-  // the documentation of the GridIn and of the Triangulation for further
-  // details on orientation of cells in a triangulation.) In our case, the
-  // normals to the mesh are external to both the circle in 2d or the sphere
-  // in 3d.
+  // The following function represents the core of the present tutorial 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.
+  // The presence of several convex and concave high curvature regions makes this
+  // 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
+  // refinement.
   //
-  // The other detail that is required for appropriate refinement of the
-  // boundary element mesh, is an accurate description of the manifold that
-  // the mesh is approximating. We already saw this several times for the
-  // boundary of standard finite element meshes (for example in step-5 and
-  // step-6), and here the principle and usage is the same, except that the
-  // HyperBallBoundary class takes an additional template parameter that
-  // specifies the embedding space dimension. The function object still has to
-  // be static to live at least as long as the triangulation object to which
-  // it is attached.
+  // 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.
+  //
+  // 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
+  // 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
+  // implemented in a different class, which can be assigned to the set_manifold method
+  // of a deal.ii 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
+    // 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
+    // class, namely a TopoDS_Shape.
 
-    std::ifstream in;
-
-    in.open ("initial_mesh_3d.inp");
+    TopoDS_Shape bow_surface = OpenCASCADE::read_IGES("DTMB-5415_bulbous_bow.iges",1e-3);
 
-    GridIn<2,3> gi;
-    gi.attach_triangulation (tria);
-    gi.read_ucd (in);
+    // Each CAD geometrical object is defined along with a tolerance, which indicates
+    // possible inaccuracy of its placement. For instance, the tolerance tol of a vertex
+    // indicates that it can be located in any point contained in a sphere centered
+    // in the nominal position and having radius tol. While projecting a point onto a
+    // surface (which will in turn have its tolerance) we must keep in mind that the
+    // precision of the projection will be limited by the tolerance with which the
+    // surface is built.
 
+    // Thus, we use a method that extracts the tolerance of a desired shape
 
-    Triangulation<2,3>::active_cell_iterator cell = tria.begin_active();
-    cell->set_manifold_id(1);
+    double tolerance = OpenCASCADE::get_shape_tolerance(bow_surface);
 
-    for (unsigned int f=0; f<GeometryInfo<2>::faces_per_cell; ++f)
-        cell->face(f)->set_manifold_id(2);
 
-    TopoDS_Shape bow_surface = OpenCASCADE::read_IGES("DTMB-5415_bulbous_bow.iges",1e-3);
+    // To stay out of trouble, we make this tolerance a bit bigger
+    tolerance*=5.0;
 
+    // We now want to extract from the generic shape, a set of composite sub-shapes (we are
+    // in particular interested in the single wire contained in the CAD file, which will
+    // allow us to define a line projector).
+    // To extract all these sub-shapes, we resort to a method of the OpenCASCADE namespace.
+    // The input of extract_compound_shapes is a shape and a set of empty std::vectors
+    // of subshapes.
     std::vector<TopoDS_Compound> compounds;
     std::vector<TopoDS_CompSolid> compsolids;
     std::vector<TopoDS_Solid> solids;
     std::vector<TopoDS_Shell> shells;
-    std::vector<TopoDS_Wire> wires;    
+    std::vector<TopoDS_Wire> wires;
 
     OpenCASCADE::extract_compound_shapes(bow_surface,
                                          compounds,
@@ -251,30 +220,91 @@ namespace Step54
                                          shells,
                                          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.
+    std::ifstream in;
+
+    in.open(initial_mesh_filename.c_str());
+
+    GridIn<2,3> gi;
+    gi.attach_triangulation(tria);
+    gi.read_vtk(in);
+
+    // 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.
+    // and assgin it the manifold_id 1
+    Triangulation<2,3>::active_cell_iterator cell = tria.begin_active();
+    cell->set_manifold_id(1);
+
+    // We also get an iterator to its faces, and assign each of them to manifold_id 2.
 
-    line_projector = new OpenCASCADE::ArclengthProjectionLineManifold<2,3>(wires[0]);
-    normal_projector = new OpenCASCADE::NormalProjectionBoundary<2,3>(bow_surface);
-    directional_projector = new OpenCASCADE::DirectionalProjectionBoundary<2,3>(bow_surface, Point<3>(0.0,1.0,0.0));
+    for (unsigned int f=0; f<GeometryInfo<2>::faces_per_cell; ++f)
+      cell->face(f)->set_manifold_id(2);
+
+    // Once both the CAD geometry and the initial mesh have been imported and digested, we
+    // use the CAD surfaces and curves to define the projectors and assign them to the
+    // manifold ids just specified.
+
+    // A first projector is defined using the single wire contained in our CAD file.
+    // The ArclengthProjectionLineManifold will make sure that every mesh edge located
+    // on the wire is refined with a point that lies on the wire and splits in two equal arcs
+    // the wire portion lying between the edge vertices.
+
+    static OpenCASCADE::ArclengthProjectionLineManifold<2,3> line_projector(wires[0], tolerance);
+
+    // Once the projector is created, we assign it to all the edges with manifold_id = 2
+    tria.set_manifold(2, line_projector);
+
+    // The surface projector is created according to what specified with the surface_projection_kind
+    // option of the constructor.
+    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
+        // 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.
+        tria.set_manifold(1,normal_projector);
+        break;
+      case 1:
+        // If the 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.
+        static OpenCASCADE::DirectionalProjectionBoundary<2,3> directional_projector(bow_surface, Point<3>(0.0,1.0,0.0), tolerance);
+        tria.set_manifold(1,directional_projector);
+        break;
+      case 2:
+        // If the 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)
+        // and a tolerance.
+        static OpenCASCADE::NormalToMeshProjectionBoundary<2,3> normal_to_mesh_projector(bow_surface, tolerance);
+        tria.set_manifold(1,normal_to_mesh_projector);
+        break;
+      default:
+        AssertThrow(false, ExcMessage("No valid projector selected: surface_projection_kind must be 0,1 or 2."));
+        break;
+      }
 
-    tria.set_manifold(2, *line_projector);
-    tria.set_manifold(1,*directional_projector);
   }
 
 
   // @sect4{TriangulationOnCAD::refine_and_resize}
 
-  // This function globally refines the mesh, distributes degrees of freedom,
-  // and resizes matrices and vectors.
+  // 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.
 
 
   void TriangulationOnCAD::refine_and_resize()
   {
     tria.refine_global(1);
-
-    dh.distribute_dofs(fe);
-
-    const unsigned int n_dofs =  dh.n_dofs();
-
   }
 
 
@@ -288,13 +318,12 @@ namespace Step54
   void TriangulationOnCAD::output_results(const unsigned int cycle)
   {
 
-    std::string filename = ( Utilities::int_to_string(3) +
-                             "d_meshhh_" +
+    std::string filename = ( output_filename + "_" +
                              Utilities::int_to_string(cycle) +
-                             ".inp" );
-  std::ofstream logfile(filename.c_str());
-  GridOut grid_out;
-  grid_out.write_ucd(tria, logfile);
+                             ".vtk" );
+    std::ofstream logfile(filename.c_str());
+    GridOut grid_out;
+    grid_out.write_vtk(tria, logfile);
 
 
   }
@@ -308,14 +337,13 @@ namespace Step54
   void TriangulationOnCAD::run()
   {
 
-    read_parameters("parameters.prm");
 
     read_domain();
-
+    unsigned int n_cycles = 5;
     for (unsigned int cycle=0; cycle<n_cycles; ++cycle)
       {
         refine_and_resize();
-        output_results(cycle);
+        output_results(cycle+1);
       }
 
   }
@@ -333,12 +361,42 @@ int main ()
       using namespace dealii;
       using namespace Step54;
 
-      const unsigned int degree = 1;
-      const unsigned int mapping_degree = 1;
 
       deallog.depth_console (3);
-      TriangulationOnCAD triangulation_on_cad(degree, mapping_degree);
-      triangulation_on_cad.run();
+
+      std::string in_mesh_filename = "initial_mesh_3d.vtk";
+
+      cout<<"----------------------------------------------------------"<<endl;
+      cout<<"Testing projection in direction normal to CAD surface"<<endl;
+      cout<<"----------------------------------------------------------"<<endl;
+      std::string out_mesh_filename = ( "3d_mesh_normal_projection" );
+      TriangulationOnCAD tria_on_cad_norm(in_mesh_filename,out_mesh_filename,0);
+      tria_on_cad_norm.run();
+      cout<<"----------------------------------------------------------"<<endl;
+      cout<<endl;
+      cout<<endl;
+
+      cout<<"----------------------------------------------------------"<<endl;
+      cout<<"Testing projection in y-axis direction"<<endl;
+      cout<<"----------------------------------------------------------"<<endl;
+      out_mesh_filename = ( "3d_mesh_directional_projection" );
+      TriangulationOnCAD tria_on_cad_dir(in_mesh_filename,out_mesh_filename,1);
+      tria_on_cad_dir.run();
+      cout<<"----------------------------------------------------------"<<endl;
+      cout<<endl;
+      cout<<endl;
+
+      cout<<"----------------------------------------------------------"<<endl;
+      cout<<"Testing projection in direction normal to mesh elements"<<endl;
+      cout<<"----------------------------------------------------------"<<endl;
+      out_mesh_filename = ( "3d_mesh_normal_to_mesh_projection" );
+      TriangulationOnCAD tria_on_cad_norm_to_mesh(in_mesh_filename,out_mesh_filename,2);
+      tria_on_cad_norm_to_mesh.run();
+      cout<<"----------------------------------------------------------"<<endl;
+      cout<<endl;
+      cout<<endl;
+
+
 
     }
   catch (std::exception &exc)
index f7ba94d8ebed95fa6e90d7fa36ec79346c5dd859..4a0553d91f81ca8f7740839b002b54c4c0b88845 100644 (file)
@@ -40,7 +40,7 @@ namespace OpenCASCADE
   /**
    * A Boundary object based on OpenCASCADE TopoDS_Shape where where
    * new points are first computed by averaging the surrounding points
-   * in the same way as FlatManifold does, and then projecting in the
+   * in the same way as FlatManifold does, and are then projected in the
    * normal direction using OpenCASCADE utilities.
    *
    * This class makes no assumptions on the shape you pass to it, and
@@ -107,7 +107,7 @@ namespace OpenCASCADE
   /**
    * A Boundary object based on OpenCASCADE TopoDS_Shape where new
    * points are first computed by averaging the surrounding points in
-   * the same way as FlatManifold does, and then projecting in onto
+   * the same way as FlatManifold does, and then projecting them onto
    * the manifold along the direction specified at construction time
    * using OpenCASCADE utilities.
    *
@@ -173,6 +173,99 @@ namespace OpenCASCADE
     const double tolerance;
   };
 
+
+  /**
+   * A Boundary object based on OpenCASCADE TopoDS_Shape where new
+   * points are first computed by averaging the surrounding points in
+   * the same way as FlatManifold does, and then projecting it using
+   * OpenCASCADE utilities onto the manifold along a direction which
+   * is an estimation of the surrounding points (hence mesh cell) normal.
+   * The direction normal to the mesh is particularly useful because
+   * it is the direction in which the mesh is missing nodes. For
+   * instance, during the refinement of a cell a new node is initially
+   * created around the baricenter of the cell. This location somehow
+   * ensures a uniform distance from the nodes of the old cell.
+   * Projecting such cell baricenter onto the CAD surface in the direction
+   * normal to the original cell will then retain uniform distance from
+   * the points of the original cell. Of course, at the stage of mesh
+   * generation, no dof handler nor finite element are defined, and
+   * such direction has to be estimated.
+   * For the case in which 8 surrounding points are present, 4 different
+   * triangles are identified with the points assigned, and the normals
+   * of such triangles are averaged to obtain the approximation of
+   * the normal to the cell.
+   * The case in which 2 surrounding points are present (i.e.:a cell
+   * edge is being refined) is of course more tricky. The average of
+   * the CAD surface normals at the 2 surrounding points is first
+   * computed, and then projected onto the plane normal to the
+   * segment linking the surrounding points. This again is an attempt
+   * to have the new point with equal distance with respect to the
+   * surrounding points
+   *
+   * This class only operates with CAD faces and makes the
+   * assumption that the shape you pass to it
+   * contains at least a face. If that is not the case, an Exeption
+   * is thrown. In debug mode there is a sanity check to
+   * make sure that the surrounding points (the ones used in
+   * project_to_manifold()) actually live on the Manifold, i.e.,
+   * calling OpenCASCADE::closest_point() on those points leaves them
+   * untouched. If this is not the case, an ExcPointNotOnManifold is
+   * thrown.
+   *
+   *
+   * Notice that this type of Boundary descriptor may fail to give
+   * results if the triangulation to be refined is close to the
+   * boundary of the given TopoDS_Shape, or when the normal direction
+   * estimated from the surrounding points does not intersect the shape.
+   * An exception
+   * is thrown when this appens.
+   *
+   * @author Luca Heltai, Andrea Mola, 2011--2014.
+   */
+  template <int dim, int spacedim>
+  class NormalToMeshProjectionBoundary : public Boundary<dim,spacedim>
+  {
+  public:
+    /**
+     * Construct a Boundary object which will project points on the
+     * TopoDS_Shape @p sh, along the given @p direction.
+     */
+    NormalToMeshProjectionBoundary(const TopoDS_Shape &sh,
+                                   const double tolerance=1e-7);
+
+    /**
+     * Perform the actual projection onto the manifold. This function,
+     * in debug mode, checks that each of the @p surrounding_points is
+     * within tolerance from the given TopoDS_Shape. If this is not
+     * the case, an exception is thrown.
+     *
+     * The projected point is computed using OpenCASCADE directional
+     * projection algorithms.
+     */
+    virtual Point<spacedim>
+    project_to_manifold (const std::vector<Point<spacedim> > &surrounding_points,
+                         const Point<spacedim> &candidate) const;
+
+  private:
+    /**
+     * The topological shape which is used internally to project
+     * points. You can construct such a shape by calling the
+     * OpenCASCADE::read_IGES() function, which will create a
+     * TopoDS_Shape with the geometry contained in the IGES file.
+     */
+    const TopoDS_Shape sh;
+
+    /**
+     * Direction used to project new points on the shape.
+     */
+    const Point<3> direction;
+
+    /**
+     * Relative tolerance used by this class to compute distances.
+     */
+    const double tolerance;
+  };
+
   /**
    * A Boundary object based on OpenCASCADE TopoDS_Shape objects which
    * have topological dimension equal to one (TopoDS_Edge or
index b3d033fe4e5bbdb9eaf17227051120c92030fe5d..8b6c27c3d0923a838ef33cc7db8f0b4006470c40 100644 (file)
@@ -121,7 +121,20 @@ namespace OpenCASCADE
   void write_IGES(const TopoDS_Shape &shape,
                   const std::string &filename);
 
-
+  /**
+    * This function returns the tolerance associated with the shape.
+    * Each CAD geometrical object is defined along with a tolerance, which indicates
+    * possible inaccuracy of its placement. For instance, the tolerance tol of a vertex
+    * indicates that it can be located in any point contained in a sphere centered
+    * in the nominal position and having radius tol. While carrying out an operation
+    * such as projecting a point onto a
+    * surface (which will in turn have its tolerance) we must keep in mind that the
+    * precision of the projection will be limited by the tolerance with which the
+    * surface is built.
+    * The tolerance is computed taking the maximum tolerance among the subshapes
+    * composing the shape.
+    */
+  double get_shape_tolerance(const TopoDS_Shape &shape);
 
   /**
    * Perform the intersection of the given topological shape with the
index f8110886be32e94757753751ed6b804be7154540..5a1864753092bcba29d84a51d804486d5b968017 100644 (file)
@@ -79,6 +79,7 @@ namespace OpenCASCADE
     return closest_point(sh, candidate, out_shape, u, v);
   }
 
+
   /*============================== DirectionalProjectionBoundary ==============================*/
   template <int dim, int spacedim>
   DirectionalProjectionBoundary<dim,spacedim>::DirectionalProjectionBoundary(const TopoDS_Shape &sh,
@@ -94,6 +95,46 @@ namespace OpenCASCADE
 
   template <int dim, int spacedim>
   Point<spacedim>  DirectionalProjectionBoundary<dim,spacedim>::
+  project_to_manifold (const std::vector<Point<spacedim> > &surrounding_points,
+                       const Point<spacedim> &candidate) const
+  {
+    TopoDS_Shape out_shape;
+    double u=0, v=0;
+    for (unsigned int i=0; i<surrounding_points.size(); ++i)
+      Assert(closest_point(sh, surrounding_points[i], out_shape, u, v)
+             .distance(surrounding_points[i]) <
+             std::max(tolerance*surrounding_points[i].norm(), tolerance),
+             ExcPointNotOnManifold(surrounding_points[i]));
+
+    return line_intersection(sh, candidate, direction, tolerance);
+
+  }
+
+
+
+  /*============================== NormalToMeshProjectionBoundary ==============================*/
+  template <int dim, int spacedim>
+  NormalToMeshProjectionBoundary<dim,spacedim>::NormalToMeshProjectionBoundary(const TopoDS_Shape &sh,
+      const double tolerance) :
+    sh(sh),
+    tolerance(tolerance)
+  {
+    Assert(spacedim == 3, ExcNotImplemented());
+
+    unsigned int n_faces;
+    unsigned int n_edges;
+    unsigned int n_vertices;
+    count_elements(sh,
+                   n_faces,
+                   n_edges,
+                   n_vertices);
+
+    Assert(n_faces > 0, ExcMessage("NormalToMeshProjectionBoundary needs a shape containing faces to operate."));
+  }
+
+
+  template <int dim, int spacedim>
+  Point<spacedim>  NormalToMeshProjectionBoundary<dim,spacedim>::
   project_to_manifold (const std::vector<Point<spacedim> > &surrounding_points,
                        const Point<spacedim> &candidate) const
   {
@@ -101,57 +142,76 @@ namespace OpenCASCADE
     double u=0, v=0;
     Point<3> average_normal(0.0,0.0,0.0);
     for (unsigned int i=0; i<surrounding_points.size(); ++i)
-        {
-        Point<3> surface_normal;
-        double mean_curvature;
-        Assert(closest_point_and_differential_forms(sh, surrounding_points[i], surface_normal, mean_curvature)
+      {
+        Assert(closest_point(sh, surrounding_points[i], out_shape, u, v)
                .distance(surrounding_points[i]) <
-             1e3*std::max(tolerance*surrounding_points[i].norm(), tolerance),
-             ExcPointNotOnManifold(surrounding_points[i]));
-        average_normal += surface_normal;
-        }
-
-    average_normal/=surrounding_points.size();
-
-    if (surrounding_points.size() == 2)
-       {
-       Point<3> P = (surrounding_points[0]+surrounding_points[1])/2;
-       Point<3> N = surrounding_points[0]-surrounding_points[1];
-       N = N/sqrt(N.square());
-       average_normal = average_normal-(average_normal*N)*N;
-       average_normal = average_normal/sqrt(average_normal.square());
-       }
-    else if (surrounding_points.size() == 8)
-       {
-       //cout<<"Ps = ["<<endl;
-       //for (unsigned int i=0; i<surrounding_points.size(); ++i)
-       //    cout<<surrounding_points[i]<<endl;
-       //cout<<"]"<<endl;
-       Point<3> u = surrounding_points[1]-surrounding_points[0];
-       Point<3> v = surrounding_points[2]-surrounding_points[0];
-       Point<3> n1(u(1)*v(2)-u(2)*v(1),u(2)*v(0)-u(0)*v(2),u(1)*v(1)-u(1)*v(0));
-       n1 = n1/n1.norm();
-       u = surrounding_points[2]-surrounding_points[3];
-       v = surrounding_points[1]-surrounding_points[3];
-       Point<3> n2(u(1)*v(2)-u(2)*v(1),u(2)*v(0)-u(0)*v(2),u(1)*v(1)-u(1)*v(0));
-       n2 = n2/n2.norm();
-
-       average_normal = (n1+n2)/2.0;
-       average_normal = average_normal/average_normal.norm();
-       }
-
-  
-    average_normal = average_normal/sqrt(average_normal.square());
-    // if for any reason the normals have zero average, just use the direction
-    // specified at the construction of the projector. Otherwise use "local" normal estimate    
-    if (average_normal.norm() < 0.9)
-       return axis_intersection(sh, candidate, direction, tolerance);
-    else
-       return axis_intersection(sh, candidate, average_normal, tolerance);
+               std::max(tolerance*surrounding_points[i].norm(), tolerance),
+               ExcPointNotOnManifold(surrounding_points[i]));
+      }
+
+    switch (surrounding_points.size())
+      {
+      case 2:
+      {
+        for (unsigned int i=0; i<surrounding_points.size(); ++i)
+          {
+            Point<3> surface_normal;
+            double mean_curvature;
+            closest_point_and_differential_forms(sh, surrounding_points[i], surface_normal, mean_curvature);
+            average_normal += surface_normal;
+          }
+
+        average_normal/=2.0;
+
+        Assert(average_normal.norm() > 1e-4,
+               ExcMessage("Failed to refine cell: the average of the surface normals at the surrounding edge turns out to be a null vector, making the projection direction undetermined."));
+
+        Point<3> P = (surrounding_points[0]+surrounding_points[1])/2;
+        Point<3> N = surrounding_points[0]-surrounding_points[1];
+        N = N/sqrt(N.square());
+        average_normal = average_normal-(average_normal*N)*N;
+        average_normal = average_normal/average_normal.norm();
+        break;
+      }
+      case 8:
+      {
+        Point<3> u = surrounding_points[1]-surrounding_points[0];
+        Point<3> v = surrounding_points[2]-surrounding_points[0];
+        Point<3> n1(u(1)*v(2)-u(2)*v(1),u(2)*v(0)-u(0)*v(2),u(1)*v(1)-u(1)*v(0));
+        n1 = n1/n1.norm();
+        u = surrounding_points[2]-surrounding_points[3];
+        v = surrounding_points[1]-surrounding_points[3];
+        Point<3> n2(u(1)*v(2)-u(2)*v(1),u(2)*v(0)-u(0)*v(2),u(1)*v(1)-u(1)*v(0));
+        n2 = n2/n2.norm();
+        u = surrounding_points[4]-surrounding_points[7];
+        v = surrounding_points[6]-surrounding_points[7];
+        Point<3> n3(u(1)*v(2)-u(2)*v(1),u(2)*v(0)-u(0)*v(2),u(1)*v(1)-u(1)*v(0));
+        n3 = n3/n3.norm();
+        u = surrounding_points[6]-surrounding_points[7];
+        v = surrounding_points[5]-surrounding_points[7];
+        Point<3> n4(u(1)*v(2)-u(2)*v(1),u(2)*v(0)-u(0)*v(2),u(1)*v(1)-u(1)*v(0));
+        n4 = n4/n4.norm();
+
+        average_normal = (n1+n2+n3+n4)/4.0;
+
+        Assert(average_normal.norm() > 1e-4,
+               ExcMessage("Failed to refine cell: the average of the surface normals at the surrounding edge turns out to be a null vector, making the projection direction undetermined."));
+
+        average_normal = average_normal/average_normal.norm();
+        break;
+      }
+      default:
+      {
+        AssertThrow(false, ExcNotImplemented());
+        break;
+      }
+      }
+
+    return line_intersection(sh, candidate, average_normal, tolerance);
   }
 
-  /*============================== ArclengthProjectionLineManifold ==============================*/
 
+  /*============================== ArclengthProjectionLineManifold ==============================*/
   template <int dim, int spacedim>
   ArclengthProjectionLineManifold<dim,spacedim>::ArclengthProjectionLineManifold(const TopoDS_Shape &sh,
       const double tolerance):
index 8194c9a707f35a089a0227646e8d8c4fbf510dda..0b8bb6e8221251585f9c77e840b08926b02e955d 100644 (file)
@@ -20,6 +20,7 @@ for (deal_II_dimension : DIMENSIONS)
   {
        template class NormalProjectionBoundary<deal_II_dimension, 3>;
        template class DirectionalProjectionBoundary<deal_II_dimension, 3>;
+        template class NormalToMeshProjectionBoundary<deal_II_dimension, 3>;
        template class ArclengthProjectionLineManifold<deal_II_dimension, 3>; 
   }
 
index 348e4ee0c99b25d599fe581f6a9501d23c61f8bc..c44b63ade471976c9aa9b5243869ce722a8fd791 100644 (file)
@@ -218,6 +218,31 @@ namespace OpenCASCADE
     AssertThrow(OK, ExcMessage("Failed to write IGES file."));
   }
 
+  double get_shape_tolerance(const TopoDS_Shape &shape)
+  {
+    double tolerance = 0.0;
+
+    std::vector<TopoDS_Face> faces;
+    std::vector<TopoDS_Edge> edges;
+    std::vector<TopoDS_Vertex> vertices;
+
+    extract_geometrical_shapes(shape,
+                               faces,
+                               edges,
+                               vertices);
+
+    for (unsigned int i=0; i<vertices.size(); ++i)
+      tolerance = fmax(tolerance,BRep_Tool::Tolerance(vertices[i]));
+
+    for (unsigned int i=0; i<edges.size(); ++i)
+      tolerance = fmax(tolerance,BRep_Tool::Tolerance(edges[i]));
+
+    for (unsigned int i=0; i<faces.size(); ++i)
+      tolerance = fmax(tolerance,BRep_Tool::Tolerance(faces[i]));
+
+
+    return tolerance;
+  }
 
   TopoDS_Shape intersect_plane(const TopoDS_Shape &in_shape,
                                const double c_x,
@@ -321,20 +346,20 @@ namespace OpenCASCADE
     Assert(Inters.IsDone(), ExcMessage("Could not project point."));
 
     double minDistance = 1e7;
-    double distance; 
+    double distance;
     int lowest_dist_int = 0;
     Point<3> result;
     for (int i=0; i<Inters.NbPnt(); ++i)
-        {
-        distance = Pnt(origin).Distance(Inters.Pnt(i+1));
-        //cout<<"Point "<<i<<": "<<Pnt(Inters.Pnt(i+1))<<"  distance: "<<distance<<endl;
+      {
+        distance = point(origin).Distance(Inters.Pnt(i+1));
+        //cout<<"Point "<<i<<": "<<point(Inters.Pnt(i+1))<<"  distance: "<<distance<<endl;
         if (distance < minDistance)
-           {
-           minDistance = distance;
-           result = Pnt(Inters.Pnt(i+1));
-           lowest_dist_int = i+1;
-           }
-        }
+          {
+            minDistance = distance;
+            result = point(Inters.Pnt(i+1));
+            lowest_dist_int = i+1;
+          }
+      }
 
     return result;
   }
@@ -420,12 +445,12 @@ namespace OpenCASCADE
     // to loop on edges. Even if the closest point lies on the boundary of a parametric surface,
     // we need in fact to retain the face and both u and v, if we want to use this method to
     // retrieve the surface normal
-    if (face_counter==0)    
-       for (exp.Init(in_shape, TopAbs_EDGE); exp.More(); exp.Next())
-         {
-           TopoDS_Edge edge = TopoDS::Edge(exp.Current());
-           if (!BRep_Tool::Degenerated(edge))
-              {
+    if (face_counter==0)
+      for (exp.Init(in_shape, TopAbs_EDGE); exp.More(); exp.Next())
+        {
+          TopoDS_Edge edge = TopoDS::Edge(exp.Current());
+          if (!BRep_Tool::Degenerated(edge))
+            {
               TopLoc_Location L;
               Standard_Real First;
               Standard_Real Last;
@@ -434,18 +459,18 @@ namespace OpenCASCADE
               // curve upon which the edge is defined
               Handle(Geom_Curve) CurveToProj = BRep_Tool::Curve(edge,L,First,Last);
 
-              GeomAPI_ProjectPointOnCurve Proj(Pnt(origin),CurveToProj);
+              GeomAPI_ProjectPointOnCurve Proj(point(origin),CurveToProj);
               unsigned int num_proj_points = Proj.NbPoints();
               if ((num_proj_points > 0) && (Proj.LowerDistance() < minDistance))
-                 {
-                 minDistance = Proj.LowerDistance();
-                 Pproj = Proj.NearestPoint();
-                 out_shape = edge;
-                 u=Proj.LowerDistanceParameter();
-                 ++counter;
-                 }
-              }
-         }
+                {
+                  minDistance = Proj.LowerDistance();
+                  Pproj = Proj.NearestPoint();
+                  out_shape = edge;
+                  u=Proj.LowerDistanceParameter();
+                  ++counter;
+                }
+            }
+        }
 
     Assert(counter > 0, ExcMessage("Could not find projection points."));
     return point(Pproj);
@@ -472,15 +497,15 @@ namespace OpenCASCADE
 
     // just a check here: the number of faces in out_shape must be 1, otherwise
     // something is wrong
-    unsigned int n_faces, n_edges, n_vertices;   
+    unsigned int n_faces, n_edges, n_vertices;
     count_elements(out_shape,
                    n_faces,
                    n_edges,
                    n_vertices);
-    
+
     Assert(n_faces > 0, ExcMessage("Could not find normal: the shape containing the closest point has 0 faces."));
     Assert(n_faces < 2, ExcMessage("Could not find normal: the shape containing the closest point has more than 1 face."));
-    
+
 
     TopExp_Explorer exp;
     exp.Init(in_shape, TopAbs_FACE);
@@ -488,7 +513,7 @@ namespace OpenCASCADE
 
     Handle(Geom_Surface) SurfToProj = BRep_Tool::Surface(face);
     GeomLProp_SLProps props(SurfToProj, u, v, 1, tolerance);
-    gp_Dir Normal = props.Normal(); 
+    gp_Dir Normal = props.Normal();
     Standard_Real Mean_Curvature = props.MeanCurvature();
 
     surface_normal = Point<3>(Normal.X(),Normal.Y(),Normal.Z());

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.