an attempt to improve opencascade directional projection with informations on local...

an other attempt to improve opencascade directional projection with informations on local normal

version with directional projection depending on normals of surrounding points is now working. to this end, a method has been added to utilities, to return differential forms (normal and mean curvature) along with closest point on the CAD surface

trying to improve the directional projection

improved cell refinement on CAD with directional projector, by adding a separated treatment for case with 2 or 8 surrounding points (edge or cell refinement). the 2 points case uses some infos about the surface normal, the 8 case comes up with an approximated cell normal

diff --git a/examples/step-54/parameters.prm b/examples/step-54/parameters.prm
index 6f276eabe2390ffd4ec26c3ae7d3698b481519c3..534c102828443abd351ffa94c50947e995c6e3ae 100644 (file)
--- a/examples/step-54/parameters.prm
+++ b/examples/step-54/parameters.prm
@@ -1,6 +1,6 @@
 # Listing of Parameters
 # ---------------------
-set Number of cycles                = 6
+set Number of cycles                = 7
 
 
 subsection Quadrature rules

diff --git a/examples/step-54/step-54.cc b/examples/step-54/step-54.cc
index ba1c919973911b48c0801fff2236f3fd85d9a6c2..7fecec4a9e3a52178f055dcc960609eae74179c4 100644 (file)
--- a/examples/step-54/step-54.cc
+++ b/examples/step-54/step-54.cc
@@ -289,7 +289,7 @@ namespace Step54
   {
 
     std::string filename = ( Utilities::int_to_string(3) +
-                             "d_mesh_" +
+                             "d_meshhh_" +
                              Utilities::int_to_string(cycle) +
                              ".inp" );
   std::ofstream logfile(filename.c_str());

diff --git a/include/deal.II/opencascade/utilities.h b/include/deal.II/opencascade/utilities.h
index 629a19331100da0d0426f4f31a18186ccaac66fa..b3d033fe4e5bbdb9eaf17227051120c92030fe5d 100644 (file)
--- a/include/deal.II/opencascade/utilities.h
+++ b/include/deal.II/opencascade/utilities.h
@@ -223,6 +223,22 @@ namespace OpenCASCADE
                          double &v,
                          const double tolerance=1e-7);
 
+
+  /**
+   * Get the closest point to the given topological shape. If the
+   * shape is not elementary, all its subshapes are iterated, faces
+   * first, then edges, and the closest point is returned together
+   * with the shape which contains it, and the u v coordinates of the
+   * point. If the returned shape is an edge, then only the u
+   * coordinate is filled with sensible information, and the v
+   * coordinate is set to zero.
+   */
+  Point<3>  closest_point_and_differential_forms(const TopoDS_Shape in_shape,
+                                                 const Point<3> origin,
+                                                 Point<3> &surface_normal,
+                                                 double &mean_curvature,
+                                                 const double tolerance=1e-7);
+
   /**
    * Intersect a line passing through the given @p origin point along
    * @p direction and the given topological shape. If there is more than

diff --git a/source/opencascade/boundary_lib.cc b/source/opencascade/boundary_lib.cc
index 585ec527ed055fc5f3e25c5e7badd947afc64eef..f8110886be32e94757753751ed6b804be7154540 100644 (file)
--- a/source/opencascade/boundary_lib.cc
+++ b/source/opencascade/boundary_lib.cc
@@ -99,14 +99,55 @@ namespace OpenCASCADE
   {
     TopoDS_Shape out_shape;
     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)
-      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),
+        {
+        Point<3> surface_normal;
+        double mean_curvature;
+        Assert(closest_point_and_differential_forms(sh, surrounding_points[i], surface_normal, mean_curvature)
+               .distance(surrounding_points[i]) <
+             1e3*std::max(tolerance*surrounding_points[i].norm(), tolerance),
              ExcPointNotOnManifold(surrounding_points[i]));
-
-    return line_intersection(sh, candidate, direction, tolerance);
+        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);
   }
 
   /*============================== ArclengthProjectionLineManifold ==============================*/

diff --git a/source/opencascade/utilities.cc b/source/opencascade/utilities.cc
index ebeb452bc82ea02f8843fec0eab78422106a9a23..348e4ee0c99b25d599fe581f6a9501d23c61f8bc 100644 (file)
--- a/source/opencascade/utilities.cc
+++ b/source/opencascade/utilities.cc
@@ -50,6 +50,7 @@
 #include <Geom_BoundedCurve.hxx>
 #include <GeomAPI_Interpolate.hxx>
 #include <GeomConvert_CompCurveToBSplineCurve.hxx>
+#include <GeomLProp_SLProps.hxx>
 
 #include <GCPnts_AbscissaPoint.hxx>
 #include <ShapeAnalysis_Surface.hxx>
@@ -314,9 +315,28 @@ namespace OpenCASCADE
     IntCurvesFace_ShapeIntersector Inters;
     Inters.Load(in_shape,tolerance);
 
-    Inters.PerformNearest(line,-RealLast(),+RealLast());
+    // Keep in mind: PerformNearest sounds pretty but DOESN'T WORK!!!
+    // The closest point must be found by hand
+    Inters.Perform(line,-RealLast(),+RealLast());
     Assert(Inters.IsDone(), ExcMessage("Could not project point."));
-    return point(Inters.Pnt(1));
+
+    double minDistance = 1e7;
+    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;
+        if (distance < minDistance)
+           {
+           minDistance = distance;
+           result = Pnt(Inters.Pnt(i+1));
+           lowest_dist_int = i+1;
+           }
+        }
+
+    return result;
   }
 
   TopoDS_Edge interpolation_curve(std::vector<Point<3> > &curve_points,
@@ -366,6 +386,7 @@ namespace OpenCASCADE
     gp_Pnt tmp_proj(0.0,0.0,0.0);
 
     unsigned int counter = 0;
+    unsigned int face_counter = 0;
     u=0;
     v=0;
 
@@ -392,38 +413,92 @@ namespace OpenCASCADE
             v=proj_params.Y();
             ++counter;
           }
+        ++face_counter;
       }
 
-    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;
-
-            // the projection function needs a Curve, so we obtain the
-            // curve upon which the edge is defined
-            Handle(Geom_Curve) CurveToProj = BRep_Tool::Curve(edge,L,First,Last);
-
-            GeomAPI_ProjectPointOnCurve Proj(point(origin),CurveToProj);
-            unsigned int num_proj_points = Proj.NbPoints();
-            if ((num_proj_points > 0) && (Proj.LowerDistance() < minDistance))
+    // face counter tells us if the shape contained faces: if it does, there is no need
+    // 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))
               {
-                minDistance = Proj.LowerDistance();
-                Pproj = Proj.NearestPoint();
-                out_shape = edge;
-                u=Proj.LowerDistanceParameter();
-                ++counter;
+              TopLoc_Location L;
+              Standard_Real First;
+              Standard_Real Last;
+
+              // the projection function needs a Curve, so we obtain the
+              // curve upon which the edge is defined
+              Handle(Geom_Curve) CurveToProj = BRep_Tool::Curve(edge,L,First,Last);
+
+              GeomAPI_ProjectPointOnCurve Proj(Pnt(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;
+                 }
               }
-          }
-      }
+         }
 
     Assert(counter > 0, ExcMessage("Could not find projection points."));
     return point(Pproj);
   }
 
+
+  Point<3> closest_point_and_differential_forms(const TopoDS_Shape in_shape,
+                                                const Point<3> origin,
+                                                Point<3> &surface_normal,
+                                                double &mean_curvature,
+                                                const double tolerance)
+
+  {
+
+    TopoDS_Shape out_shape;
+    double u, v;
+
+    Point<3> result = closest_point(in_shape,
+                                    origin,
+                                    out_shape,
+                                    u,
+                                    v,
+                                    tolerance);
+
+    // 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;   
+    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);
+    TopoDS_Face face = TopoDS::Face(exp.Current());
+
+    Handle(Geom_Surface) SurfToProj = BRep_Tool::Surface(face);
+    GeomLProp_SLProps props(SurfToProj, u, v, 1, tolerance);
+    gp_Dir Normal = props.Normal(); 
+    Standard_Real Mean_Curvature = props.MeanCurvature();
+
+    surface_normal = Point<3>(Normal.X(),Normal.Y(),Normal.Z());
+    mean_curvature = double(Mean_Curvature);
+
+    return result;
+  }
+
+
+
   void create_triangulation(const TopoDS_Face &face,
                             Triangulation<2,3> &tria)
   {