#include <deal.II/grid/grid_out.h>
#include <deal.II/grid/tria_boundary_lib.h>
+// And here are the headers of the opencascade support classes and functions:
#include <deal.II/opencascade/boundary_lib.h>
#include <deal.II/opencascade/utilities.h>
// @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();
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<2> > 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<Quadrature<2> >
- (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<GeometryInfo<2>::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<TopoDS_Compound> compounds;
std::vector<TopoDS_CompSolid> compsolids;
std::vector<TopoDS_Solid> solids;
std::vector<TopoDS_Shell> shells;
- std::vector<TopoDS_Wire> wires;
+ std::vector<TopoDS_Wire> wires;
OpenCASCADE::extract_compound_shapes(bow_surface,
compounds,
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<GeometryInfo<2>::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();
-
}
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);
}
void TriangulationOnCAD::run()
{
- read_parameters("parameters.prm");
read_domain();
-
+ unsigned int n_cycles = 5;
for (unsigned int cycle=0; cycle<n_cycles; ++cycle)
{
refine_and_resize();
- output_results(cycle);
+ output_results(cycle+1);
}
}
using namespace dealii;
using namespace Step54;
- const unsigned int degree = 1;
- const unsigned int mapping_degree = 1;
deallog.depth_console (3);
- TriangulationOnCAD triangulation_on_cad(degree, mapping_degree);
- triangulation_on_cad.run();
+
+ std::string in_mesh_filename = "initial_mesh_3d.vtk";
+
+ cout<<"----------------------------------------------------------"<<endl;
+ cout<<"Testing projection in direction normal to CAD surface"<<endl;
+ cout<<"----------------------------------------------------------"<<endl;
+ std::string out_mesh_filename = ( "3d_mesh_normal_projection" );
+ TriangulationOnCAD tria_on_cad_norm(in_mesh_filename,out_mesh_filename,0);
+ tria_on_cad_norm.run();
+ cout<<"----------------------------------------------------------"<<endl;
+ cout<<endl;
+ cout<<endl;
+
+ cout<<"----------------------------------------------------------"<<endl;
+ cout<<"Testing projection in y-axis direction"<<endl;
+ cout<<"----------------------------------------------------------"<<endl;
+ out_mesh_filename = ( "3d_mesh_directional_projection" );
+ TriangulationOnCAD tria_on_cad_dir(in_mesh_filename,out_mesh_filename,1);
+ tria_on_cad_dir.run();
+ cout<<"----------------------------------------------------------"<<endl;
+ cout<<endl;
+ cout<<endl;
+
+ cout<<"----------------------------------------------------------"<<endl;
+ cout<<"Testing projection in direction normal to mesh elements"<<endl;
+ cout<<"----------------------------------------------------------"<<endl;
+ out_mesh_filename = ( "3d_mesh_normal_to_mesh_projection" );
+ TriangulationOnCAD tria_on_cad_norm_to_mesh(in_mesh_filename,out_mesh_filename,2);
+ tria_on_cad_norm_to_mesh.run();
+ cout<<"----------------------------------------------------------"<<endl;
+ cout<<endl;
+ cout<<endl;
+
+
}
catch (std::exception &exc)
/**
* A Boundary object based on OpenCASCADE TopoDS_Shape where where
* new points are first computed by averaging the surrounding points
- * in the same way as FlatManifold does, and then projecting in the
+ * in the same way as FlatManifold does, and are then projected in the
* normal direction using OpenCASCADE utilities.
*
* This class makes no assumptions on the shape you pass to it, and
/**
* A Boundary object based on OpenCASCADE TopoDS_Shape where new
* points are first computed by averaging the surrounding points in
- * the same way as FlatManifold does, and then projecting in onto
+ * the same way as FlatManifold does, and then projecting them onto
* the manifold along the direction specified at construction time
* using OpenCASCADE utilities.
*
const double tolerance;
};
+
+ /**
+ * A Boundary object based on OpenCASCADE TopoDS_Shape where new
+ * points are first computed by averaging the surrounding points in
+ * the same way as FlatManifold does, and then projecting it using
+ * OpenCASCADE utilities onto the manifold along a direction which
+ * is an estimation of the surrounding points (hence mesh cell) normal.
+ * The direction normal to the mesh is particularly useful because
+ * it is the direction in which the mesh is missing nodes. For
+ * instance, during the refinement of a cell a new node is initially
+ * created around the baricenter of the cell. This location somehow
+ * ensures a uniform distance from the nodes of the old cell.
+ * Projecting such cell baricenter onto the CAD surface in the direction
+ * normal to the original cell will then retain uniform distance from
+ * the points of the original cell. Of course, at the stage of mesh
+ * generation, no dof handler nor finite element are defined, and
+ * such direction has to be estimated.
+ * For the case in which 8 surrounding points are present, 4 different
+ * triangles are identified with the points assigned, and the normals
+ * of such triangles are averaged to obtain the approximation of
+ * the normal to the cell.
+ * The case in which 2 surrounding points are present (i.e.:a cell
+ * edge is being refined) is of course more tricky. The average of
+ * the CAD surface normals at the 2 surrounding points is first
+ * computed, and then projected onto the plane normal to the
+ * segment linking the surrounding points. This again is an attempt
+ * to have the new point with equal distance with respect to the
+ * surrounding points
+ *
+ * This class only operates with CAD faces and makes the
+ * assumption that the shape you pass to it
+ * contains at least a face. If that is not the case, an Exeption
+ * is thrown. 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.
+ *
+ *
+ * Notice that this type of Boundary descriptor may fail to give
+ * results if the triangulation to be refined is close to the
+ * boundary of the given TopoDS_Shape, or when the normal direction
+ * estimated from the surrounding points does not intersect the shape.
+ * An exception
+ * is thrown when this appens.
+ *
+ * @author Luca Heltai, Andrea Mola, 2011--2014.
+ */
+ template <int dim, int spacedim>
+ class NormalToMeshProjectionBoundary : public Boundary<dim,spacedim>
+ {
+ 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<spacedim>
+ 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 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
return closest_point(sh, candidate, out_shape, u, v);
}
+
/*============================== DirectionalProjectionBoundary ==============================*/
template <int dim, int spacedim>
DirectionalProjectionBoundary<dim,spacedim>::DirectionalProjectionBoundary(const TopoDS_Shape &sh,
template <int dim, int spacedim>
Point<spacedim> DirectionalProjectionBoundary<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]) <
+ std::max(tolerance*surrounding_points[i].norm(), tolerance),
+ ExcPointNotOnManifold(surrounding_points[i]));
+
+ return line_intersection(sh, candidate, direction, tolerance);
+
+ }
+
+
+
+ /*============================== NormalToMeshProjectionBoundary ==============================*/
+ template <int dim, int spacedim>
+ NormalToMeshProjectionBoundary<dim,spacedim>::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 <int dim, int spacedim>
+ Point<spacedim> NormalToMeshProjectionBoundary<dim,spacedim>::
project_to_manifold (const std::vector<Point<spacedim> > &surrounding_points,
const Point<spacedim> &candidate) const
{
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)
- {
- Point<3> 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 = ["<<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);
+ std::max(tolerance*surrounding_points[i].norm(), tolerance),
+ ExcPointNotOnManifold(surrounding_points[i]));
+ }
+
+ switch (surrounding_points.size())
+ {
+ case 2:
+ {
+ for (unsigned int i=0; i<surrounding_points.size(); ++i)
+ {
+ Point<3> 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 <int dim, int spacedim>
ArclengthProjectionLineManifold<dim,spacedim>::ArclengthProjectionLineManifold(const TopoDS_Shape &sh,
const double tolerance):
AssertThrow(OK, ExcMessage("Failed to write IGES file."));
}
+ double get_shape_tolerance(const TopoDS_Shape &shape)
+ {
+ double tolerance = 0.0;
+
+ std::vector<TopoDS_Face> faces;
+ std::vector<TopoDS_Edge> edges;
+ std::vector<TopoDS_Vertex> vertices;
+
+ extract_geometrical_shapes(shape,
+ faces,
+ edges,
+ vertices);
+
+ for (unsigned int i=0; i<vertices.size(); ++i)
+ tolerance = fmax(tolerance,BRep_Tool::Tolerance(vertices[i]));
+
+ for (unsigned int i=0; i<edges.size(); ++i)
+ tolerance = fmax(tolerance,BRep_Tool::Tolerance(edges[i]));
+
+ for (unsigned int i=0; i<faces.size(); ++i)
+ tolerance = fmax(tolerance,BRep_Tool::Tolerance(faces[i]));
+
+
+ return tolerance;
+ }
TopoDS_Shape intersect_plane(const TopoDS_Shape &in_shape,
const double c_x,
Assert(Inters.IsDone(), ExcMessage("Could not project point."));
double minDistance = 1e7;
- double distance;
+ 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;
+ {
+ distance = point(origin).Distance(Inters.Pnt(i+1));
+ //cout<<"Point "<<i<<": "<<point(Inters.Pnt(i+1))<<" distance: "<<distance<<endl;
if (distance < minDistance)
- {
- minDistance = distance;
- result = Pnt(Inters.Pnt(i+1));
- lowest_dist_int = i+1;
- }
- }
+ {
+ minDistance = distance;
+ result = point(Inters.Pnt(i+1));
+ lowest_dist_int = i+1;
+ }
+ }
return result;
}
// 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))
- {
+ 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))
+ {
TopLoc_Location L;
Standard_Real First;
Standard_Real Last;
// curve upon which the edge is defined
Handle(Geom_Curve) CurveToProj = BRep_Tool::Curve(edge,L,First,Last);
- GeomAPI_ProjectPointOnCurve Proj(Pnt(origin),CurveToProj);
+ GeomAPI_ProjectPointOnCurve Proj(point(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;
- }
- }
- }
+ {
+ 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);
// 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);
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());
projects / dealii.git / commitdiff