From d6f3c7cbd7f089ffecaf342496e58aead162bf9b Mon Sep 17 00:00:00 2001 From: amola Date: Thu, 2 Oct 2014 12:41:10 +0200 Subject: [PATCH] an attempt to improve opencascade directional projection with informations on local normal 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 --- examples/step-54/parameters.prm | 2 +- examples/step-54/step-54.cc | 2 +- include/deal.II/opencascade/utilities.h | 16 +++ source/opencascade/boundary_lib.cc | 53 ++++++++-- source/opencascade/utilities.cc | 125 +++++++++++++++++++----- 5 files changed, 165 insertions(+), 33 deletions(-) diff --git a/examples/step-54/parameters.prm b/examples/step-54/parameters.prm index 6f276eabe2..534c102828 100644 --- 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 ba1c919973..7fecec4a9e 100644 --- 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 629a193311..b3d033fe4e 100644 --- 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 585ec527ed..f8110886be 100644 --- 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; i0 ? - 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 = ["< 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 ebeb452bc8..348e4ee0c9 100644 --- a/source/opencascade/utilities.cc +++ b/source/opencascade/utilities.cc @@ -50,6 +50,7 @@ #include #include #include +#include #include #include @@ -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 > &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) { -- 2.39.5