From: Andrea Mola Date: Tue, 14 Oct 2014 15:45:03 +0000 (+0200) Subject: step-54 is now completed (it only needs comments). three different projectors have... X-Git-Tag: v8.2.0-rc1~81^2~4 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=77aa8ff71d8ec7abb9a42d14d977a012cf517677;p=dealii.git step-54 is now completed (it only needs comments). three different projectors have now been created in boundary_lib: normal, directional and normal_to_mesh. to adjust tolerances, a new method retrieving tolerances from a shape has been added 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. --- diff --git a/examples/step-54/coarse_circle.inp b/examples/step-54/coarse_circle.inp deleted file mode 100644 index dda97cde9e..0000000000 --- a/examples/step-54/coarse_circle.inp +++ /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 index 5ea249c4dc..0000000000 --- a/examples/step-54/coarse_sphere.inp +++ /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 index 8080f6733f..0000000000 --- a/examples/step-54/initial_mesh_3d.inp +++ /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 index 0000000000..ddf377e281 --- /dev/null +++ b/examples/step-54/initial_mesh_3d.vtk @@ -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 index 534c102828..0000000000 --- a/examples/step-54/parameters.prm +++ /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 - - - diff --git a/examples/step-54/step-54.cc b/examples/step-54/step-54.cc index 7fecec4a9e..2c030087ea 100644 --- a/examples/step-54/step-54.cc +++ b/examples/step-54/step-54.cc @@ -43,6 +43,7 @@ #include #include +// And here are the headers of the opencascade support classes and functions: #include #include @@ -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; - - 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 > - (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::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 compounds; std::vector compsolids; std::vector solids; std::vector shells; - std::vector wires; + std::vector 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::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 + class NormalToMeshProjectionBoundary : public Boundary + { + 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 + project_to_manifold (const std::vector > &surrounding_points, + const Point &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 diff --git a/include/deal.II/opencascade/utilities.h b/include/deal.II/opencascade/utilities.h index b3d033fe4e..8b6c27c3d0 100644 --- a/include/deal.II/opencascade/utilities.h +++ b/include/deal.II/opencascade/utilities.h @@ -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 diff --git a/source/opencascade/boundary_lib.cc b/source/opencascade/boundary_lib.cc index f8110886be..5a18647530 100644 --- a/source/opencascade/boundary_lib.cc +++ b/source/opencascade/boundary_lib.cc @@ -79,6 +79,7 @@ namespace OpenCASCADE return closest_point(sh, candidate, out_shape, u, v); } + /*============================== DirectionalProjectionBoundary ==============================*/ template DirectionalProjectionBoundary::DirectionalProjectionBoundary(const TopoDS_Shape &sh, @@ -94,6 +95,46 @@ namespace OpenCASCADE template Point DirectionalProjectionBoundary:: + project_to_manifold (const std::vector > &surrounding_points, + const Point &candidate) const + { + TopoDS_Shape out_shape; + double u=0, v=0; + for (unsigned int i=0; i + NormalToMeshProjectionBoundary::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 + Point NormalToMeshProjectionBoundary:: project_to_manifold (const std::vector > &surrounding_points, const Point &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 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 = ["< 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 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 ArclengthProjectionLineManifold::ArclengthProjectionLineManifold(const TopoDS_Shape &sh, const double tolerance): diff --git a/source/opencascade/boundary_lib.inst.in b/source/opencascade/boundary_lib.inst.in index 8194c9a707..0b8bb6e822 100644 --- a/source/opencascade/boundary_lib.inst.in +++ b/source/opencascade/boundary_lib.inst.in @@ -20,6 +20,7 @@ for (deal_II_dimension : DIMENSIONS) { template class NormalProjectionBoundary; template class DirectionalProjectionBoundary; + template class NormalToMeshProjectionBoundary; template class ArclengthProjectionLineManifold; } diff --git a/source/opencascade/utilities.cc b/source/opencascade/utilities.cc index 348e4ee0c9..c44b63ade4 100644 --- a/source/opencascade/utilities.cc +++ b/source/opencascade/utilities.cc @@ -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 faces; + std::vector edges; + std::vector vertices; + + extract_geometrical_shapes(shape, + faces, + edges, + vertices); + + for (unsigned int i=0; i result; for (int i=0; i 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());