]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Added OCC::NormalProjectionBoundary
authorLuca Heltai <luca.heltai@sissa.it>
Tue, 12 Aug 2014 13:54:12 +0000 (15:54 +0200)
committerLuca Heltai <luca.heltai@sissa.it>
Fri, 31 Oct 2014 10:06:28 +0000 (11:06 +0100)
Added interpolation curve boundary test.

Updated documentation.

include/deal.II/grid/occ_boundary_lib.h [new file with mode: 0644]
include/deal.II/grid/occ_utilities.h
source/grid/CMakeLists.txt
source/grid/occ_boundary_lib.cc [new file with mode: 0644]
source/grid/occ_boundary_lib.inst.in [new file with mode: 0644]
source/grid/occ_utilities.cc
tests/opencascade/circle_normal_projection.cc [new file with mode: 0644]
tests/opencascade/circle_normal_projection.output [new file with mode: 0644]
tests/opencascade/interpolation_curve_boundary.cc [new file with mode: 0644]
tests/opencascade/interpolation_curve_boundary.output [new file with mode: 0644]

diff --git a/include/deal.II/grid/occ_boundary_lib.h b/include/deal.II/grid/occ_boundary_lib.h
new file mode 100644 (file)
index 0000000..e46be21
--- /dev/null
@@ -0,0 +1,110 @@
+// ---------------------------------------------------------------------
+// $Id$
+//
+// Copyright (C) 2014 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+#ifndef __deal2__occ_boundary_lib_h
+#define __deal2__occ_boundary_lib_h
+
+#include <deal.II/base/config.h>
+
+#ifdef DEAL_II_WITH_OPENCASCADE
+
+#include <deal.II/grid/occ_utilities.h>
+#include <deal.II/grid/tria_boundary.h>
+
+DEAL_II_NAMESPACE_OPEN
+
+namespace OpenCASCADE 
+{
+  /**
+   * @addtogroup OpenCASCADE
+   * @{
+   * 
+   * A Boundary object based on OpenCASCADE TopoDS_Shape where new
+   * points are first computed using the FlatManifold class, and then
+   * projected in the normal direction using OpenCASCADE utilities.
+   *
+   * This class makes no assumptions on the shape you pass to it, and
+   * the topological dimension of the Manifold is inferred from the
+   * TopoDS_Shape itself. 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.
+   *
+   * This could happen, for example, if you are trying to use a shape
+   * of type TopoDS_Edge when projecting on a face. In this case, the
+   * vertices of the face would be collapsed to the edge, and your
+   * surrounding points would not be on lying on the given shape,
+   * raising an exception.
+   * 
+   * @author Luca Heltai, Andrea Mola, 2011--2014.
+   */
+  template <int dim, int spacedim>
+  class NormalProjectionBoundary : public Boundary<dim,spacedim> {
+    public:
+    NormalProjectionBoundary(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 #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 normal
+     * projection algorithms.
+     */
+    virtual Point<spacedim>
+    project_to_manifold (const std::vector<Point<spacedim> > &surrounding_points,
+                        const Point<spacedim> &candidate) const;
+
+    /**
+     * Exception thrown when the point specified as argument does not
+     * lie between #tolerance from the given TopoDS_Shape.
+     */
+    DeclException1 (ExcPointNotOnManifold,
+                   Point<spacedim>,
+                   <<"The point [ "<<arg1<<" ] is not on the manifold.");
+
+  private:
+    /**
+     * The topological shape which is used internally to project
+     * points. You can construct one 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;
+
+    /**
+     * Relative tolerance used by this class to compute distances.
+     */
+    const double tolerance;
+  };
+
+}
+
+/*@}*/
+
+DEAL_II_NAMESPACE_CLOSE
+
+
+#endif // DEAL_II_WITH_OPENCASCADE
+
+/*------------------------------ occ_boundary_lib.h ------------------------------*/
+#endif
+/*------------------------------ occ_boundary_lib.h ------------------------------*/
index e6f9b5c15bb91925b99451bc6b02fa85256f24d1..45b1d52b0f2ba289d176b5cf183b465b55762b34 100644 (file)
@@ -18,6 +18,8 @@
 #ifndef __deal2__occ_utilities_h
 #define __deal2__occ_utilities_h
 
+#include <deal.II/base/config.h>
+
 #ifdef DEAL_II_WITH_OPENCASCADE
 
 #include <string>
@@ -27,7 +29,7 @@
 #include <Geom_Curve.hxx>
 #include <gp_Pnt.hxx>
 
-#include <base/point.h>
+#include <deal.II/base/point.h>
 
 DEAL_II_NAMESPACE_OPEN
 
@@ -123,18 +125,28 @@ namespace OpenCASCADE
   
   /**
    * Creates a 3D smooth BSpline curve passing through the points in
-   * the assigned vector. The points are reordered internally
-   * according to their scalar product with the direction, if
-   * direction is different from zero, otherwise they are used as
+   * the assigned vector, and store it in the returned TopoDS_Shape
+   * (which is of type TopoDS_Edge). The points are reordered
+   * internally according to their scalar product with the direction,
+   * if direction is different from zero, otherwise they are used as
    * passed. Notice that this function changes the input points if
    * required by the algorithm.
    *
    * This class is used to interpolate a BsplineCurve passing through
-   * an array of points, with a C2 Continuity.
+   * an array of points, with a C2 Continuity. If the optional
+   * parameter #closed is set to true, then the curve will be C2 at
+   * all points execpt the first (where only C1 continuity will be
+   * given), and it will be a closed curve.
+   *
+   * The curve is garanteed to be at distance #tolerance from the
+   * input points. If the algorithm fails in generating such a curve,
+   * an exception is thrown.
    */
-  TopoDS_Shape interpolation_curve(std::vector<dealii::Point<3> >  &curve_points,
-                                  const dealii::Point<3> direction=dealii::Point<3>(), 
-                                  const double tolerance=1e-7);
+  TopoDS_Edge interpolation_curve(std::vector<dealii::Point<3> >  &curve_points,
+                                 const dealii::Point<3> direction=dealii::Point<3>(), 
+                                 const bool closed=false,
+                                 const double tolerance=1e-7);
+
 
   /**
    * Get the closest point to the given topological shape. If the
@@ -155,19 +167,20 @@ namespace OpenCASCADE
   /**
    * Convert OpenCASCADE point into a Point<3>.
    */
-  Point<3> Pnt(const gp_Pnt &p);
+  inline Point<3> Pnt(const gp_Pnt &p);
 
 
   /**
    * Convert Point<3> into OpenCASCADE point.
    */
-  gp_Pnt Pnt(const Point<3> &p);
+  inline gp_Pnt Pnt(const Point<3> &p);
 
   
   /**
    * Sort two points according to their scalar product with
    * direction. If the norm of the direction is zero, then use
-   * lexycographical ordering.
+   * lexycographical ordering. The optional parameter is used as a
+   * relative tolerance when comparing objects.
    */
   inline bool point_compare(const dealii::Point<3> &p1, const dealii::Point<3> &p2,
                            const dealii::Point<3> direction=Point<3>(),
index 738fc8c43b6e09973505f94a3bb4c98d5a900a26..8fa6c75d027b27fdb4284ae8fe0105aed5475f81 100644 (file)
@@ -26,6 +26,7 @@ SET(_src
   manifold.cc
   manifold_lib.cc
   occ_utilities.cc
+  occ_boundary_lib.cc
   persistent_tria.cc
   tria_accessor.cc
   tria_boundary.cc
@@ -45,6 +46,7 @@ SET(_inst
   intergrid_map.inst.in
   manifold.inst.in
   manifold_lib.inst.in
+  occ_boundary_lib.inst.in
   tria_accessor.inst.in
   tria_boundary.inst.in
   tria_boundary_lib.inst.in
diff --git a/source/grid/occ_boundary_lib.cc b/source/grid/occ_boundary_lib.cc
new file mode 100644 (file)
index 0000000..2ef10fa
--- /dev/null
@@ -0,0 +1,42 @@
+#include <deal.II/grid/occ_boundary_lib.h>
+
+#ifdef DEAL_II_WITH_OPENCASCADE
+
+DEAL_II_NAMESPACE_OPEN
+
+namespace OpenCASCADE
+{
+  template <int dim, int spacedim>
+  NormalProjectionBoundary<dim,spacedim>::NormalProjectionBoundary(const TopoDS_Shape sh, 
+                                                                  const double tolerance) :
+    sh(sh),
+    tolerance(tolerance) 
+  {
+    Assert(spacedim == 3, ExcImpossibleInDim(spacedim));
+  }
+  
+  
+  template <int dim, int spacedim>
+  Point<spacedim>  NormalProjectionBoundary<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]) < (surrounding_points[i].norm()>0 ? 
+                                                tolerance*surrounding_points[i].norm() :
+                                                tolerance),
+            ExcPointNotOnManifold(surrounding_points[i]));
+    
+    return closest_point(sh, candidate, out_shape, u, v);
+  }
+  
+  // Explicit instantiations
+#include "occ_boundary_lib.inst"  
+  
+} // end namespace OpenCASCADE
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif
diff --git a/source/grid/occ_boundary_lib.inst.in b/source/grid/occ_boundary_lib.inst.in
new file mode 100644 (file)
index 0000000..6e97bc7
--- /dev/null
@@ -0,0 +1,25 @@
+// ---------------------------------------------------------------------
+// $Id$
+//
+// Copyright (C) 1998 - 2013 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+
+for (deal_II_dimension : DIMENSIONS)
+  {
+       template class NormalProjectionBoundary<deal_II_dimension, 3>;
+  }
+
+
+
index 61c62445c9d16fb7d7d42fe67a9c5ab9c24afd7c..73b04d3ed129c9c0a2e8b79b3748ccb14cff4997 100644 (file)
@@ -1,5 +1,4 @@
 #include <deal.II/grid/occ_utilities.h>
-#include <deal.II/base/config.h>
 
 #ifdef DEAL_II_WITH_OPENCASCADE
 
@@ -106,10 +105,10 @@ namespace OpenCASCADE
   }
   
   inline bool point_compare(const dealii::Point<3> &p1, const dealii::Point<3> &p2,
-                           const dealii::Point<3> direction=Point<3>(),
-                           const double tolerance=1e-10)
+                           const dealii::Point<3> direction,
+                           const double tolerance)
   {
-    const double rel_tol=std::max(p1.norm(), p2.norm())*tolerance;
+    const double rel_tol=std::max(tolerance, std::max(p1.norm(), p2.norm())*tolerance);
     if(direction.norm() > 0.0)
       return (p1*direction < p2*direction-rel_tol);
     else 
@@ -245,8 +244,9 @@ namespace OpenCASCADE
 
 
 
-  TopoDS_Shape interpolation_curve(std::vector<Point<3> > &curve_points,
+  TopoDS_Edge interpolation_curve(std::vector<Point<3> > &curve_points,
                                   const Point<3> direction,
+                                  const bool closed,
                                   const double tolerance)
   {
 
@@ -266,12 +266,12 @@ namespace OpenCASCADE
       }
 
 
-    GeomAPI_Interpolate bspline_generator(vertices, false, tolerance);
+    GeomAPI_Interpolate bspline_generator(vertices, closed, tolerance);
     bspline_generator.Perform();
     Assert( (bspline_generator.IsDone()), ExcMessage("Interpolated bspline generation failed"));
     
     Handle(Geom_BSplineCurve) bspline = bspline_generator.Curve();
-    TopoDS_Shape out_shape = BRepBuilderAPI_MakeEdge(bspline);
+    TopoDS_Edge out_shape = BRepBuilderAPI_MakeEdge(bspline);
     return out_shape;
   }
 
diff --git a/tests/opencascade/circle_normal_projection.cc b/tests/opencascade/circle_normal_projection.cc
new file mode 100644 (file)
index 0000000..f40c3ed
--- /dev/null
@@ -0,0 +1,102 @@
+
+//----------------------------  occ_circle_projection.cc  ---------------------------
+//    $Id: testsuite.html 13373 2006-07-13 13:12:08Z kanschat $
+//    Version: $Name$ 
+//
+//    Copyright (C) 2005 by the deal.II authors 
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//----------------------------  occ_circle_projection.cc  ---------------------------
+
+
+// Create a circle, a Triangulation, and try to project normally on
+// it.
+
+#include "../tests.h"
+
+#include <deal.II/grid/occ_boundary_lib.h>
+
+#include <fstream>
+#include <base/logstream.h>
+#include <grid/tria.h>
+#include <grid/tria_accessor.h>
+#include <grid/grid_out.h>
+#include <grid/tria_iterator.h>
+#include <grid/grid_generator.h>
+#include <grid/tria_boundary_lib.h>
+
+#include <gp_Pnt.hxx>
+#include <gp_Dir.hxx>
+#include <gp_Ax2.hxx>
+#include <GC_MakeCircle.hxx>
+#include <BRepBuilderAPI_MakeEdge.hxx>
+#include <BRepBuilderAPI_MakeFace.hxx>
+#include <BRepBuilderAPI_MakeWire.hxx>
+#include <TopoDS_Wire.hxx>
+#include <TopoDS_Edge.hxx>
+#include <TopoDS_Face.hxx>
+
+using namespace OpenCASCADE;
+
+int main () 
+{
+  std::ofstream logfile("output");
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+
+                                  // The circle passing through the
+                                  // vertices of the unit square
+  gp_Dir z_axis(0.,0.,1.);
+  gp_Pnt center(.5,.5,0.);
+  gp_Ax2 axis(center, z_axis);
+  Standard_Real radius(std::sqrt(2.)/2.);
+  
+  Handle(Geom_Curve) circle = GC_MakeCircle(axis, radius);
+  TopoDS_Edge edge = BRepBuilderAPI_MakeEdge(circle);
+  
+  
+                                  // Create a boundary projector.
+  NormalProjectionBoundary<2,3> boundary_line(edge);
+  
+                                  // This one is for checking: This
+                                  // is what deal.II would do for a
+                                  // circle.
+  HyperBallBoundary<2,3> boundary_line_deal (Point<3>(.5,.5,0),
+                                            std::sqrt(2.)/2.);
+  
+  
+  
+                                  // The unit square.
+  Triangulation<2,3> tria;
+  GridGenerator::hyper_cube(tria);
+
+                                  // Set the exterior boundary
+  tria.set_boundary(0, boundary_line);
+                                  // Uncomment the following line to
+                                  // see the difference.
+// tria.set_boundary(0, boundary_line_deal);
+  
+
+                                  // This is here to ignore the
+                                  // points created in the interior
+                                  // of the face.
+  tria.begin()->set_material_id(1);
+
+                                  // We refine twice, and expect the
+                                  // outer points to end up on the
+                                  // circle.
+  tria.refine_global(2);
+
+
+                                  // You can open the generated file
+                                  // with paraview.
+  GridOut gridout;
+  gridout.write_ucd (tria, logfile);
+  
+  return 0;
+}
+                  
diff --git a/tests/opencascade/circle_normal_projection.output b/tests/opencascade/circle_normal_projection.output
new file mode 100644 (file)
index 0000000..3291307
--- /dev/null
@@ -0,0 +1,43 @@
+
+25 16 0 0 0
+1  0.00000 0.00000 0.00000
+2  1.00000 0.00000 0.00000
+3  0.00000 1.00000 0.00000
+4  1.00000 1.00000 0.00000
+5  0.500000 -0.207107 0.00000
+6  -0.207107 0.500000 0.00000
+7  1.20711 0.500000 0.00000
+8  0.500000 1.20711 0.00000
+9  0.500000 0.500000 0.00000
+10  0.229402 -0.153281 0.00000
+11  0.770598 -0.153281 0.00000
+12  -0.153281 0.229402 0.00000
+13  -0.153281 0.770598 0.00000
+14  1.15328 0.229402 0.00000
+15  1.15328 0.770598 0.00000
+16  0.229402 1.15328 0.00000
+17  0.770598 1.15328 0.00000
+18  0.500000 0.146447 0.00000
+19  0.500000 0.853553 0.00000
+20  0.146447 0.500000 0.00000
+21  0.853553 0.500000 0.00000
+22  0.189433 0.189433 0.00000
+23  0.810567 0.189433 0.00000
+24  0.189433 0.810567 0.00000
+25  0.810567 0.810567 0.00000
+1 1 quad    1 10 22 12 
+2 1 quad    10 5 18 22 
+3 1 quad    12 22 20 6 
+4 1 quad    22 18 9 20 
+5 1 quad    5 11 23 18 
+6 1 quad    11 2 14 23 
+7 1 quad    18 23 21 9 
+8 1 quad    23 14 7 21 
+9 1 quad    6 20 24 13 
+10 1 quad    20 9 19 24 
+11 1 quad    13 24 16 3 
+12 1 quad    24 19 8 16 
+13 1 quad    9 21 25 19 
+14 1 quad    21 7 15 25 
+15 1 quad    19 25 17 8 
+16 1 quad    25 15 4 17 
diff --git a/tests/opencascade/interpolation_curve_boundary.cc b/tests/opencascade/interpolation_curve_boundary.cc
new file mode 100644 (file)
index 0000000..9d43fa3
--- /dev/null
@@ -0,0 +1,87 @@
+
+//----------------------------  occ_circle_projection.cc  ---------------------------
+//    $Id: testsuite.html 13373 2006-07-13 13:12:08Z kanschat $
+//    Version: $Name$ 
+//
+//    Copyright (C) 2005 by the deal.II authors 
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//----------------------------  occ_circle_projection.cc  ---------------------------
+
+
+// Create a Triangulation, interpolate its boundary points to a close
+// smooth BSpline, and use that as a Boundary Descriptor.
+
+#include "../tests.h"
+
+#include <deal.II/grid/occ_boundary_lib.h>
+
+#include <fstream>
+#include <base/logstream.h>
+#include <grid/tria.h>
+#include <grid/tria_accessor.h>
+#include <grid/grid_out.h>
+#include <grid/tria_iterator.h>
+#include <grid/grid_generator.h>
+#include <grid/tria_boundary_lib.h>
+
+#include <gp_Pnt.hxx>
+#include <gp_Dir.hxx>
+#include <gp_Ax2.hxx>
+#include <GC_MakeCircle.hxx>
+#include <BRepBuilderAPI_MakeEdge.hxx>
+#include <BRepBuilderAPI_MakeFace.hxx>
+#include <BRepBuilderAPI_MakeWire.hxx>
+#include <TopoDS_Wire.hxx>
+#include <TopoDS_Edge.hxx>
+#include <TopoDS_Face.hxx>
+
+using namespace OpenCASCADE;
+
+int main () 
+{
+  std::ofstream logfile("output");
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+
+  // The curve passing through the vertices of the unit square.
+  std::vector<Point<3> > vertices;
+  vertices.push_back(Point<3>(0,0,0));
+  vertices.push_back(Point<3>(1,0,0));
+  vertices.push_back(Point<3>(1,1,0));
+  vertices.push_back(Point<3>(0,1,0));
+  TopoDS_Shape shape = interpolation_curve(vertices, Point<3>(), true);
+  
+  // Create a boundary projector.
+  NormalProjectionBoundary<2,3> boundary_line(shape);
+
+  // The unit square.
+  Triangulation<2,3> tria;
+  GridGenerator::hyper_cube(tria);
+
+  // Set the exterior boundary
+  tria.set_boundary(0, boundary_line);
+
+  // This is here to ignore the points created in the interior of the
+  // face.
+  tria.begin()->set_material_id(1);
+
+  // We refine twice, and expect the outer points to end up on a
+  // smooth curve interpolating the square vertices.
+  tria.refine_global(2);
+
+
+  // You can open the generated file with paraview.
+  GridOut gridout;
+  gridout.write_msh (tria, logfile);
+
+  std::ofstream of("/tmp/grid.msh");
+  gridout.write_msh(tria, of);
+  
+  return 0;
+}
+                  
diff --git a/tests/opencascade/interpolation_curve_boundary.output b/tests/opencascade/interpolation_curve_boundary.output
new file mode 100644 (file)
index 0000000..894bc26
--- /dev/null
@@ -0,0 +1,48 @@
+
+$NOD
+25
+1  0.00000 0.00000 0.00000
+2  1.00000 0.00000 0.00000
+3  0.00000 1.00000 0.00000
+4  1.00000 1.00000 0.00000
+5  0.439917 -0.236739 0.00000
+6  -0.310022 0.473832 0.00000
+7  1.14617 0.492348 0.00000
+8  0.488859 1.16307 0.00000
+9  0.470616 0.486564 0.00000
+10  0.207735 -0.140071 0.00000
+11  0.764326 -0.218873 0.00000
+12  -0.202060 0.202837 0.00000
+13  -0.224611 0.783020 0.00000
+14  1.10302 0.237503 0.00000
+15  1.12328 0.758622 0.00000
+16  0.234185 1.11236 0.00000
+17  0.760208 1.13858 0.00000
+18  0.455267 0.124913 0.00000
+19  0.479738 0.824816 0.00000
+20  0.0802967 0.480198 0.00000
+21  0.808395 0.489456 0.00000
+22  0.142719 0.173942 0.00000
+23  0.773464 0.171896 0.00000
+24  0.152383 0.790482 0.00000
+25  0.784659 0.794181 0.00000
+$ENDNOD
+$ELM
+16
+1 3 1 0 4 1 10 22 12 
+2 3 1 0 4 10 5 18 22 
+3 3 1 0 4 12 22 20 6 
+4 3 1 0 4 22 18 9 20 
+5 3 1 0 4 5 11 23 18 
+6 3 1 0 4 11 2 14 23 
+7 3 1 0 4 18 23 21 9 
+8 3 1 0 4 23 14 7 21 
+9 3 1 0 4 6 20 24 13 
+10 3 1 0 4 20 9 19 24 
+11 3 1 0 4 13 24 16 3 
+12 3 1 0 4 24 19 8 16 
+13 3 1 0 4 9 21 25 19 
+14 3 1 0 4 21 7 15 25 
+15 3 1 0 4 19 25 17 8 
+16 3 1 0 4 25 15 4 17 
+$ENDELM

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.