From: Luca Heltai <luca.heltai@sissa.it>
Date: Tue, 12 Aug 2014 15:58:41 +0000 (+0200)
Subject: Added Axis Projection.
X-Git-Tag: v8.2.0-rc1~81^2~10
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=d0442dbfec9bdeaf9068191b81d9d990a83b9b40;p=dealii.git

Added Axis Projection.
---

diff --git a/include/deal.II/grid/occ_boundary_lib.h b/include/deal.II/grid/occ_boundary_lib.h
index e46be218fe..eb53542afa 100644
--- a/include/deal.II/grid/occ_boundary_lib.h
+++ b/include/deal.II/grid/occ_boundary_lib.h
@@ -73,13 +73,58 @@ namespace OpenCASCADE
     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 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;
+
     /**
-     * Exception thrown when the point specified as argument does not
-     * lie between #tolerance from the given TopoDS_Shape.
+     * Relative tolerance used by this class to compute distances.
      */
-    DeclException1 (ExcPointNotOnManifold,
-		    Point<spacedim>,
-		    <<"The point [ "<<arg1<<" ] is not on the manifold.");
+    const double tolerance;
+  };
+
+  /**
+   * A Boundary object based on OpenCASCADE TopoDS_Shape where new
+   * points are first computed using the FlatManifold class, and then
+   * projected along the direction given at construction time, 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.
+   * 
+   * @author Luca Heltai, Andrea Mola, 2011--2014.
+   */
+  template <int dim, int spacedim>
+  class AxisProjectionBoundary : public Boundary<dim,spacedim> {
+    public:
+    AxisProjectionBoundary(const TopoDS_Shape sh, 
+			   const Point<3> direction, 
+			   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;
 
   private:
     /**
@@ -90,6 +135,11 @@ namespace OpenCASCADE
      */ 
     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.
      */
diff --git a/include/deal.II/grid/occ_utilities.h b/include/deal.II/grid/occ_utilities.h
index 45b1d52b0f..d483b58e11 100644
--- a/include/deal.II/grid/occ_utilities.h
+++ b/include/deal.II/grid/occ_utilities.h
@@ -164,6 +164,19 @@ namespace OpenCASCADE
 			 double &v, 
 			 const double tolerance=1e-7);
 
+  /**
+   * Intersect a line passing through the given #origin point along
+   * #direction and the given topological shape. If there is more than
+   * one intersection, it will return the closest one.
+   *
+   * The optional #tolerance parameter is used to compute distances.
+   */
+  Point<3> axis_intersection(const TopoDS_Shape in_shape, 
+			     const Point<3> origin, 
+			     const Point<3> direction, 
+			     const double tolerance=1e-7);
+  
+
   /**
    * Convert OpenCASCADE point into a Point<3>.
    */
@@ -185,8 +198,34 @@ 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);
-				      
-}
+
+
+  /**
+   * Exception thrown when the point specified as argument does not
+   * lie between #tolerance from the given TopoDS_Shape.
+   */
+  DeclException1 (ExcPointNotOnManifold,
+		  Point<3>,
+		  <<"The point [ "<<arg1<<" ] is not on the manifold.");
+
+  /**
+   * Exception thrown when the point specified as argument cannot be
+   * projected to the manifold.
+   */			      
+  DeclException1 (ExcProjectionFailed, 
+		  Point<3>,
+		  <<"Projection of point [ "<< arg1
+		  << " ] failed.");
+
+  /**
+   * Thrown when internal OpenCASCADE utilities fail to return the OK
+   * status.
+   */
+  DeclException1 (ExcOCCError, 
+		  IFSelect_ReturnStatus, 
+		  <<"An OpenCASCADE routine failed with return status "
+		  <<arg1);
+} 
 /*@}*/
 
 DEAL_II_NAMESPACE_CLOSE
diff --git a/source/grid/occ_boundary_lib.cc b/source/grid/occ_boundary_lib.cc
index 2ef10faf2e..64724e9148 100644
--- a/source/grid/occ_boundary_lib.cc
+++ b/source/grid/occ_boundary_lib.cc
@@ -6,6 +6,8 @@ DEAL_II_NAMESPACE_OPEN
 
 namespace OpenCASCADE
 {
+
+  /*============================== NormalProjectionBoundary ==============================*/
   template <int dim, int spacedim>
   NormalProjectionBoundary<dim,spacedim>::NormalProjectionBoundary(const TopoDS_Shape sh, 
 								   const double tolerance) :
@@ -27,10 +29,41 @@ namespace OpenCASCADE
 	     .distance(surrounding_points[i]) < (surrounding_points[i].norm()>0 ? 
 						 tolerance*surrounding_points[i].norm() :
 						 tolerance),
-	     ExcPointNotOnManifold(surrounding_points[i]));
+	     // ExcPointNotOnManifold(surrounding_points[i]));
+	     ExcMessage("Points are not on manifold"));
     
     return closest_point(sh, candidate, out_shape, u, v);
   }
+
+  /*============================== AxisProjectionBoundary ==============================*/
+  template <int dim, int spacedim>
+  AxisProjectionBoundary<dim,spacedim>::AxisProjectionBoundary(const TopoDS_Shape sh, 
+							       const Point<3> direction,
+							       const double tolerance) :
+    sh(sh),
+    direction(direction),
+    tolerance(tolerance) 
+  {
+    Assert(spacedim == 3, ExcImpossibleInDim(spacedim));
+  }
+  
+  
+  template <int dim, int spacedim>
+  Point<spacedim>  AxisProjectionBoundary<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 axis_intersection(sh, candidate, direction, tolerance);
+  }
+
   
   // Explicit instantiations
 #include "occ_boundary_lib.inst"  
diff --git a/source/grid/occ_boundary_lib.inst.in b/source/grid/occ_boundary_lib.inst.in
index 6e97bc7391..3471c0aea8 100644
--- a/source/grid/occ_boundary_lib.inst.in
+++ b/source/grid/occ_boundary_lib.inst.in
@@ -19,6 +19,7 @@
 for (deal_II_dimension : DIMENSIONS)
   {
 	template class NormalProjectionBoundary<deal_II_dimension, 3>;
+	template class AxisProjectionBoundary<deal_II_dimension, 3>;
   }
 
 
diff --git a/source/grid/occ_utilities.cc b/source/grid/occ_utilities.cc
index 73b04d3ed1..1bf4571e1b 100644
--- a/source/grid/occ_utilities.cc
+++ b/source/grid/occ_utilities.cc
@@ -26,10 +26,12 @@
 #include <TColStd_SequenceOfTransient.hxx>
 #include <TColStd_HSequenceOfTransient.hxx>
 #include <TopExp_Explorer.hxx>
-// #include <gp_Pnt.hxx>
-// #include <gp_Vec.hxx>
+#include <gp_Pnt.hxx>
+#include <gp_Lin.hxx>
+#include <gp_Vec.hxx>
 #include <GeomAPI_ProjectPointOnSurf.hxx>
 #include <GeomAPI_ProjectPointOnCurve.hxx>
+#include <IntCurvesFace_ShapeIntersector.hxx>
 // #include <Standard_Real.hxx>
 // #include <Standard_Integer.hxx>
 // #include <BRep_Tool.hxx>
@@ -243,6 +245,29 @@ namespace OpenCASCADE
    }
 
 
+  Point<3> axis_intersection(const TopoDS_Shape in_shape, 
+			     const Point<3> origin, 
+			     const Point<3> direction, 
+			     const double tolerance)
+  {
+    // translating original Point<dim> to gp point
+
+    gp_Pnt P0 = Pnt(origin);
+    gp_Ax1 gpaxis(P0, gp_Dir(direction(0), direction(1), direction(2)));
+    gp_Lin line(gpaxis);
+
+    // destination point
+    gp_Pnt Pproj(0.0,0.0,0.0);
+
+    // we prepare now the surface for the projection we get the whole
+    // shape from the iges model
+    IntCurvesFace_ShapeIntersector Inters;
+    Inters.Load(in_shape,tolerance);
+	
+    Inters.PerformNearest(line,-RealLast(),+RealLast());
+    Assert(Inters.IsDone(), ExcMessage("Could not project point."));
+    return Pnt(Inters.Pnt(1));
+  }
 
   TopoDS_Edge interpolation_curve(std::vector<Point<3> > &curve_points,
 				   const Point<3> direction,