Fixed indentation.
Attempt to fix create tria test.
Second attempt. This is the correct output.
Added create_tria_00.
Added supersimple test.
attempt to fix the normal projection tests that were failing by switching to wigley.igs surface (instead of side.igs which seemed to have something weird)
switching to wigley.iges in create_tria_01.cc too
added test normal_to_mesh_projection_02.cc on surface mesh refinement with projection normal to mesh direction. in the process, this highlited a bug in the boundary_lib.cc implementation of normalToMeshBoundary. Now also the results of step-54 look even better
Fixed output for iges write functions.
Added two simple tests.
Fixed a little the documentation, and improved slightly step-54.
Fixed indent.
By a mathematical perspective, the engine lying under the hood of CAD modelers is represented
by analytical geometry, and in particular by parametric curves and surfaces such as B-splines
and NURBS.
-Thus, once a virtual model is ready, all the geometrical features of the desired object
+Once a virtual model is ready, all the geometrical features of the desired object
are stored in files which materially contain the coefficients of the parametric surfaces and
curves composing the object. Depending on the specific CAD tool used to define the geometrical model,
there are of course several different file formats in which the information of a CAD model can
<h3> The boundary projector classes </h3>
-To import and interrogate CAD models, in the current deal.II release we developed a series of
-classes that wrap functions of the OpenCASCADE open source library for CAD modeling. Once imported
-from an IGES file, the model is stored in a TopoDS_Shape, which is the generic topological entity
+To import and interrogate CAD models, the deal.II library contains a series of
+wrapper functions for the OpenCASCADE open source library for CAD modeling. These functions allow to import
+IGES files into OpenCASCADE native objects, and wrap them inside a series of Manifold classes.
+
+Once imported from an IGES file, the model is stored in a TopoDS_Shape, which is the generic topological entity
defined into the OpenCASCADE framework. From a TopoDS_Shape, it is then possible to access to all the
sub-shapes (such as vertices, edges and faces) composing it, along with their geometrical description.
-In the deal.II framework, the topological entities composing a shape are used to create an object
-of the Boundary class. In Step-6 we saw how to build a HyperBallBoundary and assign
+In the deal.II framework, the topological entities composing a shape are used to create objects
+of the Manifold or Boundary classes. In Step-6 we saw how to build a HyperBallBoundary and assign
it to a set of faces (or cells, for co-dimension 1) of a Triangulation, to have cells and faces
refined on a sphere or circle.
-The new functions of the CAD modeling interface have been designed to retain the same structure. Hence, we build a
-projector object using the imported CAD shapes, and using the
-very same procedure we were previously using with HyperBallBoundary we assign such projector object to cells,
+The functions of the CAD modeling interface have been designed to retain the same structure, allowing the user to build a
+projector object using the imported CAD shapes, maintaining the
+very same procedure we use with HyperBallBoundary, i.e., assigning such projector object to cells,
faces or edges of a coarse mesh. At each refinement cycle, the new mesh nodes will be then automatically generated
-on the specified geometry. But differently from a spherical or circular boundary, a boundary with a complex geometry
-poses proplems as to where exactly place on the prescribed shape the new nodes created upon refinement.
+on the specified geometry.
+
+Differently from a spherical or circular boundary, a boundary with a complex geometry
+poses proplems as to where it is best to place the new nodes created upon refinement on the prescribed shape.
HyperBallBoundary first creates the new nodes on the cell to be refined by averaging the
surrounding points in the same way as FlatManifold does. Successively, it goes on to project such nodes on the circle or
-sphere along the radial direction. On such geometry, the radial direction ensures that the newly generated nodes remain evenly
-spaced across the refinement levels required.
-In the case of an arbitrary and complex shape though, the correct direction of projection cannot be identified that easily.
-That is why we have implemented several different projector classes that employ different projection strategies. A first
-projector, the ArclengthProjectionLineManifold, is to be used only for edge refinement. It is built assigning it a TopoDS_Edge,
+sphere along the radial direction. On such a geometry, the radial direction ensures that the newly generated nodes remain evenly
+spaced when remaining on a given refinement level.
+
+In the case of an arbitrary and complex shape though, the best direction of the projection cannot be identified that easily.
+The library provides several projector classes that employ different projection strategies. A first
+projector, the ArclengthProjectionLineManifold, is to be used only for edge refinement. It is built assigning it a topological shape
+of dimension one, either a TopoDS_Edge or a TopoDS_Wire (which is a compound shape, made of several connected TopoDS_Edges)
and refines a mesh edge finding the new vertex as the point splitting in two even parts the curvilinear length of the CAD curve
portion that lies between the vertices of the original edge.
<img src="http://www.dealii.org/images/steps/developer/step-54.CurveSplit.png" alt="" width="500">
-A different projection stategy has been implemented in the NormalProjectionBoundary. The TopoDS_Shape assigned in construction can be arbitrary
-(a collection of shapes, faces, edges or a single face or edge willall work). The new cell nodes are first computed by averaging the
+A different projection stategy has been implemented in the NormalProjectionBoundary. The TopoDS_Shape assigned at construction time can be arbitrary
+(a collection of shapes, faces, edges or a single face or edge will all work). The new cell nodes are first computed by averaging the
surrounding points in the same way as FlatManifold does. In a second moment all the new nodes will be projected onto the TopoDS_Shape
-along the direction normal to the shape. If no normal projection is available, the closest point ---typically lying on the shape boundary---
-to the original one is selected. In cases instead where more than one normal projection is available, the closest of them to the original
-point is again selected. If the shape is composed by several sub-shapes, the projection is carried out onto every single
+along the direction normal to the shape. If no normal projection is available, the point which is closest to the shape ---typically lying on the shape boundary---
+is selected. If the shape is composed by several sub-shapes, the projection is carried out onto every single
sub-shape, and the closest projection point point is selected.
<img src="http://www.dealii.org/images/steps/developer/step-54.NormalProjectionEdge.png" alt="" width="500">
<img src="http://www.dealii.org/images/steps/developer/step-54.NormalProjection.png" alt="" width="500">
-As we are about to experience, in presence of arbitrary shapes setting the projection direction as that normal to the CAD surface will not
-lead to surface mesh elements of suitable quality. This because the direction normal to the CAD surface has in principle nothing to do with the
-direction along which the mesh needs the new nodes to be located. For this reason we implemented the DirectionalProjectionBoundary, which
-is constructed assigning a TopoDS_Shape (containing at least a face) and a direction along which all the projections will be carried out.
+As we are about to experience, for some shapes, setting the projection direction as that normal to the CAD surface will not
+lead to surface mesh elements of suitable quality. This is because the direction normal to the CAD surface has in principle nothing to do with the
+direction along which the mesh needs the new nodes to be located. The DirectionalProjectionBoundary, in this case, can help. This class
+is constructed assigning a TopoDS_Shape (containing at least a face) and a direction along which all the projections will be carried out. New points will be
+computed by first averaging the surrounding points (as in the FlatManifold case), and then taking the closest intersection between the topological shape
+and the line passing throught the resulting point, along the direction used at construction time.
In this way, the user will have a higher control on the projection direction to be enforced to ensure good mesh quality.
<img src="http://www.dealii.org/images/steps/developer/step-54.DirectionalProjection.png" alt="" width="500">
Of course the latter approach results effective only when the orientation of the surface is rather uniform, so that a single projection
direction can be identified. In cases in which the surface direction is approaching the projection direction, it is even possible that the
-directional projection is not found. To overcome these problems, in the NormalToMeshProjectionBoundary class we also implemented a third
-projection algorithm. The NormalToMeshProjectionBoundary is built assigning a TopoDS_Shape (containing at least a face) to the contructor,
+directional projection is not found. To overcome these problems, the NormalToMeshProjectionBoundary class implements a third
+projection algorithm. The NormalToMeshProjectionBoundary is built assigning a TopoDS_Shape (containing at least one face) to the contructor,
and works exactly like a DirectionalProjection. But, as the name of the class suggests, NormalToMeshProjectionBoundary tries to come up
with a suitable estimate of the direction normal to the mesh elements to be refined, and uses it for the projection of the new nodes
onto the CAD surface. If we consider a mesh edge in a 2D space, the direction of its axis is a direction along which to split it in order to
the cell normal.
In the next figure, which is inpired to the geometry considered in this tutorial, we make an attempt to compare the behavior of the
-three projectors considered. As can be seen on the left, given the original cell (in blue), the new point found with the normal projection is found
+three projectors considered. As can be seen on the left, given the original cell (in blue), the new point found with the normal projection is
in a position which does not allow for the generation of evenly spaced new elements (in red). The situation will get worse in further refinement steps.
-Since the geometry here considered is somehow perpendicular to the horizontal direction, we notice that the directional projection (central image) defined
+Since the geometry we considered is somehow perpendicular to the horizontal direction, the directional projection (central image) defined
with horizontal direction as the projection direction, does a rather good job in getting the new mesh point. Yet, since the surface is almost horizontal
-at the bottom of the picture, we can expect problems in such region with a directional projection, if further refinement steps are carried out. Finally,
+at the bottom of the picture, we can expect problems in those regions when further refinement steps are carried out. Finally,
the picture on the right shows that a node located on the cell axis will result in two new cells having the same length. Of course the situation
in 3D gets a little more complicated than that described in this simple 2D case. Nevertheless, the results of this test confirm that the normal to the
mesh direction is the best approach among the three tested, when arbitrarily shaped surfaces are considered.
// @sect3{Include files}
-// We start with including a bunch of files that we will use
-// in the various parts of the program. Most of them have been discussed in
+// We start with including a bunch of files that we will use in the
+// various parts of the program. Most of them have been discussed in
// previous tutorials already:
#include <deal.II/grid/tria.h>
#include <deal.II/grid/tria_iterator.h>
#include <deal.II/numerics/data_out.h>
#include <deal.II/numerics/vector_tools.h>
-// And here are the headers of the opencascade support classes and functions:
+// These are the headers of the opencascade support classes and
+// functions. Notice that these will contain sensible data only if you
+// compiled your deal.II library with support for opencascade, i.e.,
+// specifying -DDEAL_II_WITH_OPENCASCADE=ON and
+// -DOPENCASCADE_DIR=/path/to/your/opencascade/installation
#include <deal.II/opencascade/boundary_lib.h>
#include <deal.II/opencascade/utilities.h>
-// Finally, a few C++ standard header files that we will need:
+// Finally, a few C++ standard header files
#include <cmath>
#include <iostream>
#include <fstream>
#include <string>
-// The last part of this preamble is to import everything in the dealii
-// namespace into the one into which everything in this program will go:
+// We isolate the rest of the program in its own namespace
namespace Step54
{
using namespace dealii;
// @sect3{The TriangulationOnCAD class}
- // 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 a flag indicating the surface projection strategy we want to
- // test. Along with the input argument, the only other member of the class
- // is the triangulation we want to play with.
- // The member functions of this class are similar to those that in most of the
- // other tutorial programs deal with the setup of the grid for the
- // simulations.
+ // 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 a flag indicating the surface projection strategy
+ // we want to use. Along with the input arguments, the only other
+ // member of the class is the triangulation we want to play with.
+ // The member functions of this class are similar to what you can
+ // find in most of the other tutorial programs in the setup stage of
+ // the grid for the simulations.
class TriangulationOnCAD
{
public:
- TriangulationOnCAD(const std::string &initial_mesh_filename,
- const std::string &output_filename,
- const unsigned int &surface_projection_kind = 0);
+ // An enumeration type, to specify the type of projection to use
+ enum ProjectionType { NormalProjection = 0,
+ DirectionalProjection = 1,
+ NormalToMeshProjection = 2
+ };
+
+
+ TriangulationOnCAD(const std::string initial_mesh_filename,
+ const std::string cad_file_name,
+ const std::string output_filename,
+ const ProjectionType surface_projection_kind = NormalProjection);
~TriangulationOnCAD();
Triangulation<2, 3> tria;
- const std::string &initial_mesh_filename;
- const std::string &output_filename;
- const unsigned int &surface_projection_kind;
+ const std::string initial_mesh_filename;
+ const std::string cad_file_name;
+ const std::string output_filename;
+ const ProjectionType surface_projection_kind;
};
// 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(const std::string &initial_mesh_filename,
- const std::string &output_filename,
- const unsigned int &surface_projection_kind)
+ // names, and the enumeration type that etermines which kind of
+ // surface projector is used in the mesh refinement cycles (see
+ // below for details).
+
+ TriangulationOnCAD::TriangulationOnCAD(const std::string initial_mesh_filename,
+ const std::string cad_file_name,
+ const std::string output_filename,
+ const ProjectionType surface_projection_kind)
:
initial_mesh_filename(initial_mesh_filename),
+ cad_file_name(cad_file_name),
output_filename(output_filename),
surface_projection_kind(surface_projection_kind)
{
// @sect4{TriangulationOnCAD::read_domain}
- // The following function represents the core of thhis 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.
+ // The following function represents the core of this program. In
+ // this function we import the CAD shape upon which we want to
+ // generate and refine our triangulation. We assume that the CAD
+ // surface is contained in the @p cad_file_name file (we provide an
+ // example IGES file in this directory called
+ // "DTMB-5415_bulbous_bow.iges" that represents the bulbous bow of a
+ // ship). The presence of several convex and concave high curvature
+ // regions makes the geometry we provided a particularly meaningful
+ // example.
//
- // So, after importing the hull bow surface, we extract some of the curves and surfaces
- // composing it, and use them to generate a set of projectors. Such projectors substantially
- // define the rules the Triangulation has to follow to position each new node during the cell
+ // After importing the hull bow surface, we extract some of the
+ // curves and surfaces composing it, and use them to generate a set
+ // of projectors. Such projectors define the rules the Triangulation
+ // has to follow to position each new node during the cell
// refinement.
//
- // To initialize 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 placed on the CAD shape.
+ // To initialize the Triangulation, as done in previous tutorial
+ // programs, we import a pre-existing grid saved in VTK format. It
+ // is assumed here that the user has generated a coarse mesh
+ // externally, which matches the IGES geometry. At the moment, the
+ // deal.II library does not support automatic generation of such
+ // meshes, but there are several tools which can provide you with
+ // reasonable initial meshes starting from CAD files.
//
- // So, after importing both the IGES geometry and the initial mesh, we assign the projectors
- // previously discussed to each of the edges and cells which will have to be
- // refined on the CAD surface.
+ // In our example, the imported mesh is composed of a single
+ // quadrilateral cell whose vertices have been placed on the CAD
+ // shape.
//
- // In this tutorial, we will test the three different CAD surface projectors described in
- // the introduction, and will analyze the results obtained with each of them.
- // As mentioned, each of such projection strategiy has been
- // implemented in a different class, which can be assigned to the set_manifold method
- // of a Triangulation class.
+ // After importing both the IGES geometry and the initial mesh, we
+ // assign the projectors previously discussed to each of the edges
+ // and cells which will have to be refined on the CAD surface.
//
-
-
+ // In this tutorial, we will test the three different CAD surface
+ // projectors described in the introduction, and will analyze the
+ // results obtained with each of them. As mentioned, each of such
+ // projection strategiy has been implemented in a different class,
+ // which can be assigned to the set_manifold method of the
+ // Triangulation object.
void TriangulationOnCAD::read_domain()
{
- // The following function allows for the CAD file of interest (in IGES format) to be imported.
- // The function arguments 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.
-
- TopoDS_Shape bow_surface = OpenCASCADE::read_IGES("DTMB-5415_bulbous_bow.iges",1e-3);
-
- // 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
-
+ // The following function is used to import the given CAD file.
+ // The function arguments 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
+ // (which is the typical unit of measure for most IGES files),
+ // 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.
+
+ TopoDS_Shape bow_surface = OpenCASCADE::read_IGES(cad_file_name,1e-3);
+
+ // 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.
+
+ // The following method extracts the tolerance of the given shape
double tolerance = OpenCASCADE::get_shape_tolerance(bow_surface);
// 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.
+ // We now want to extract from the generic shape, a set of
+ // composite sub-shapes. In particular, each face of the CAD file
+ // is composed of a trimming curve of type TopoDS_Wire, which is
+ // the collection of TopoDS_Edges that compose the boundary of a
+ // surface, and a NURBS description of the surface itself. We will
+ // use a line projector to associate the boundary of our
+ // Triangulation to the wire delimiting the surface. To extract
+ // all compound sub-shapes, like wires, shells, or solids, 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, which will be filled with all
+ // compound shapes found in the given topological shape:
std::vector<TopoDS_Compound> compounds;
std::vector<TopoDS_CompSolid> compsolids;
std::vector<TopoDS_Solid> solids;
cell->set_manifold_id(1);
// We also get an iterator to its faces, and assign each of them to manifold_id 2.
-
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.
+ // 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.
+ // 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. We first check
+ // that the wires vector contains at least one element:
+ Assert(wires.size() > 0,
+ ExcMessage("I could not find any wire in the CAD file you gave me. Bailing out."));
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
+ // 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.
+ // The surface projector is created according to what is specified
+ // with the surface_projection_kind option of the constructor.
switch (surface_projection_kind)
{
- case 0:
- // If surface_projection_kind value is 0, we select the NormalProjectionBoundary. The new mesh points will initially
- // generated at the baricenter of the cell/edge considered, 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);
+ case NormalProjection:
+ // If surface_projection_kind value is 0, we select the
+ // NormalProjectionBoundary. The new mesh points will
+ // initially generated at the baricenter of the cell/edge
+ // considered, 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);
// Once created, the normal projector is assigned to the manifold having id 1.
tria.set_manifold(1,normal_projector);
break;
- case 1:
- // If surface_projection_kind 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);
+ case DirectionalProjection:
+ // If surface_projection_kind 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 surface_projection_kind 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);
+ case NormalToMeshProjection:
+ // If surface_projection_kind 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."));
+ // We should never get here
+ AssertThrow(false, ExcInternalError());
break;
}
// @sect4{TriangulationOnCAD::refine_and_resize}
- // 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 Triangulation produced.
-
+ // 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 Triangulation produced.
void TriangulationOnCAD::refine_and_resize()
{
// @sect4{TriangulationOnCAD::output_results}
// Outputting the results of our computations is a rather mechanical
- // tasks. All the components of this function have been discussed before.
+ // tasks. All the components of this function have been discussed
+ // before.
void TriangulationOnCAD::output_results(const unsigned int cycle)
{
deallog.depth_console (3);
std::string in_mesh_filename = "initial_mesh_3d.vtk";
+ std::string cad_file_name = "DTMB-5415_bulbous_bow.iges";
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);
+ TriangulationOnCAD tria_on_cad_norm(in_mesh_filename,
+ cad_file_name,
+ out_mesh_filename,
+ TriangulationOnCAD::NormalProjection);
tria_on_cad_norm.run();
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);
+ TriangulationOnCAD tria_on_cad_dir(in_mesh_filename,
+ cad_file_name,
+ out_mesh_filename,
+ TriangulationOnCAD::DirectionalProjection);
tria_on_cad_dir.run();
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);
+ TriangulationOnCAD tria_on_cad_norm_to_mesh(in_mesh_filename,
+ cad_file_name,
+ out_mesh_filename,
+ TriangulationOnCAD::NormalToMeshProjection);
tria_on_cad_norm_to_mesh.run();
cout<<"----------------------------------------------------------"<<endl;
cout<<endl;
cout<<endl;
-
-
-
}
catch (std::exception &exc)
{
const Point<spacedim> &candidate) const
{
TopoDS_Shape out_shape;
- double u=0, v=0;
Point<3> average_normal(0.0,0.0,0.0);
#ifdef DEBUG
for (unsigned int i=0; i<surrounding_points.size(); ++i)
{
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));
+ Point<3> n1(u(1)*v(2)-u(2)*v(1),u(2)*v(0)-u(0)*v(2),u(0)*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));
+ Point<3> n2(u(1)*v(2)-u(2)*v(1),u(2)*v(0)-u(0)*v(2),u(0)*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));
+ Point<3> n3(u(1)*v(2)-u(2)*v(1),u(2)*v(0)-u(0)*v(2),u(0)*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));
+ Point<3> n4(u(1)*v(2)-u(2)*v(1),u(2)*v(0)-u(0)*v(2),u(0)*v(1)-u(1)*v(0));
n4 = n4/n4.norm();
+ //for (unsigned int i=0; i<surrounding_points.size(); ++i)
+ // cout<<surrounding_points[i]<<endl;
+ //cout<<"-"<<endl;
+ //cout<<n1<<endl;cout<<n2<<endl;cout<<n3<<endl;cout<<n4<<endl;
average_normal = (n1+n2+n3+n4)/4.0;
Assert(average_normal.norm() > tolerance,
- 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."));
+ ExcMessage("Failed to refine cell: the normal estimated via the surrounding points turns out to be a null vector, making the projection direction undetermined."));
average_normal /= average_normal.norm();
break;
--- /dev/null
+//-----------------------------------------------------------
+//
+// Copyright (C) 2014 by the deal.II authors
+//
+// This file is subject to LGPL and may not be distributed
+// without copyright and license information. Please refer
+// to the file deal.II/doc/license.html for the text and
+// further information on this license.
+//
+//-----------------------------------------------------------
+
+// Read the file side.iges, and
+// create a single-cell Triangulation
+
+#include "../tests.h"
+
+#include <deal.II/opencascade/boundary_lib.h>
+
+#include <fstream>
+#include <base/logstream.h>
+#include <grid/tria.h>
+#include <grid/tria_accessor.h>
+#include <grid/grid_out.h>
+#include <grid/tria_iterator.h>
+#include <grid/grid_generator.h>
+#include <grid/tria_boundary_lib.h>
+
+#include <gp_Pnt.hxx>
+#include <gp_Dir.hxx>
+#include <gp_Ax2.hxx>
+#include <GC_MakeCircle.hxx>
+#include <BRepBuilderAPI_MakeEdge.hxx>
+#include <BRepBuilderAPI_MakeFace.hxx>
+#include <BRepBuilderAPI_MakeWire.hxx>
+#include <TopoDS_Wire.hxx>
+#include <TopoDS_Edge.hxx>
+#include <TopoDS_Face.hxx>
+#include <BRepPrimAPI_MakeSphere.hxx>
+
+using namespace OpenCASCADE;
+
+int main ()
+{
+ std::ofstream logfile("output");
+
+ TopoDS_Shape sh = read_IGES(SOURCE_DIR "/iges_files/wigley.iges", 0.001);
+ std::vector<TopoDS_Face> faces;
+ std::vector<TopoDS_Edge> edges;
+ std::vector<TopoDS_Vertex> vertices;
+
+ extract_geometrical_shapes(sh, faces, edges, vertices);
+
+ Triangulation<2,3> tria;
+ create_triangulation(faces[0], tria);
+
+ // You can open the generated file with gmsh
+ GridOut gridout;
+ gridout.write_msh (tria, logfile);
+
+ return 0;
+}
+
--- /dev/null
+$NOD
+4
+1 -1.25 0 0.15625
+2 1.25 -3.88578e-16 0.15625
+3 -1.25 6.93889e-18 -0.15625
+4 1.25 6.93889e-18 -0.15625
+$ENDNOD
+$ELM
+1
+1 3 0 0 4 1 2 4 3
+$ENDELM
#include <TopoDS_Face.hxx>
#include <BRepPrimAPI_MakeSphere.hxx>
+
using namespace OpenCASCADE;
int main ()
{
std::ofstream logfile("output");
- TopoDS_Shape sh = read_IGES(SOURCE_DIR "/iges_files/side.iges");
+ TopoDS_Shape sh = read_IGES(SOURCE_DIR "/iges_files/wigley.iges", 0.001);
std::vector<TopoDS_Face> faces;
std::vector<TopoDS_Edge> edges;
std::vector<TopoDS_Vertex> vertices;
extract_geometrical_shapes(sh, faces, edges, vertices);
// Create a boundary projector on the first face.
- NormalProjectionBoundary<2,3> boundary(faces[0], 1e-9);
+ NormalProjectionBoundary<2,3> boundary(faces[0]);
// Create a Triangulation with a single cell
Triangulation<2,3> tria;
$NOD
25
-1 0 0 0.001
-2 -2.11667e-14 0.001 0.001
-3 0 0 0
-4 -2.11667e-14 0.001 0
-5 -0.000234467 0.000444022 0.001
-6 3.40745e-16 -4.38101e-16 0.0005
-7 -2.11537e-14 0.001 0.0005
-8 -0.000234467 0.000444022 0
-9 -0.000234467 0.000444022 0.0005
-10 -0.000139999 0.00020943 0.001
-11 -0.00021601 0.000764021 0.001
-12 5.11117e-16 -6.57151e-16 0.00075
-13 5.11117e-16 -6.57151e-16 0.00025
-14 -2.11472e-14 0.001 0.00075
-15 -2.11472e-14 0.001 0.00025
-16 -0.000139999 0.00020943 0
-17 -0.00021601 0.000764021 0
-18 -0.000234467 0.000444022 0.00075
-19 -0.000234467 0.000444022 0.00025
-20 -0.000139999 0.00020943 0.0005
-21 -0.00021601 0.000764021 0.0005
-22 -0.000139999 0.00020943 0.00075
-23 -0.00021601 0.000764021 0.00075
-24 -0.000139999 0.00020943 0.00025
-25 -0.00021601 0.000764021 0.00025
+1 -1.25 0 0.15625
+2 1.25 -3.88578e-16 0.15625
+3 -1.25 6.93889e-18 -0.15625
+4 1.25 6.93889e-18 -0.15625
+5 1.54074e-33 0.124767 0.15625
+6 -1.25 -4.16334e-17 2.77556e-17
+7 1.25 -3.747e-16 2.77556e-17
+8 -1.11022e-16 2.08167e-17 -0.15625
+9 -2.24033e-17 0.113663 -0.0466831
+10 -0.628108 0.0935705 0.15625
+11 0.628108 0.0935705 0.15625
+12 -1.25 0 0.078125
+13 -1.25 1.38778e-17 -0.078125
+14 1.25 -3.88578e-16 0.078125
+15 1.25 -2.91434e-16 -0.078125
+16 -0.625 1.73472e-17 -0.15625
+17 0.625 1.73472e-17 -0.15625
+18 -1.26424e-17 0.124766 0.0547831
+19 -1.8952e-18 0.0638253 -0.109282
+20 -0.628157 0.0898099 -0.0312632
+21 0.628157 0.0898099 -0.0312632
+22 -0.628298 0.0935509 0.0654633
+23 0.628298 0.0935509 0.0654633
+24 -0.626331 0.0507074 -0.105812
+25 0.626331 0.0507074 -0.105812
$ENDNOD
$ELM
16
pts.push_back(Point<3>(1,0,0));
TopoDS_Edge edge = interpolation_curve(pts);
- write_IGES(edge, "output");
-
+ write_IGES(edge, "tmp.iges");
+ std::ifstream in("tmp.iges");
+ std::ofstream out("output");
+ std::string line;
+ unsigned int counter = 5;
+ while(counter--) std::getline(in, line);
+ while(std::getline(in, line))
+ out << line << std::endl;
+ in.close();
+ out.close();
return 0;
}
- S0000001
-,,31HOpen CASCADE IGES processor 6.7,13HFilename.iges, G0000001
-16HOpen CASCADE 6.7,31HOpen CASCADE IGES processor 6.7,32,308,15,308,15,G0000002
-,1.,2,2HMM,1,0.01,15H20140813.171933,1E-07,1.309333,6Hheltai,,11,0, G0000003
-15H20140813.171933,; G0000004
126 1 0 0 0 0 0 000000000D0000001
126 0 0 4 0 0D0000002
126,5,3,1,0,1,0,0.,0.,0.,0.,1.,2.,3.,3.,3.,3.,1.,1.,1.,1.,1.,1., 0000001P0000001
int main ()
{
TopoDS_Shape sh = read_IGES(SOURCE_DIR "/iges_files/wigley.iges");
- unsigned int nf=0, ne=0, nv=0;
- count_elements(sh, nf, ne, nv);
-
- write_IGES(sh, "output");
+ write_IGES(sh, "tmp.iges");
+ std::ifstream in("tmp.iges");
+ std::ofstream out("output");
+ std::string line;
+ unsigned int counter = 5;
+ while(counter--) std::getline(in,line);
+ while(std::getline(in,line))
+ out << line << std::endl;
+ in.close();
+ out.close();
return 0;
}
- S0000001
-,,31HOpen CASCADE IGES processor 6.7,13HFilename.iges, G0000001
-16HOpen CASCADE 6.7,31HOpen CASCADE IGES processor 6.7,32,308,15,308,15,G0000002
-,1.,2,2HMM,1,0.01,15H20140813.132822,1E-07,1.25,6Hheltai,,11,0, G0000003
-15H20140813.132822,; G0000004
144 1 0 0 0 0 0 000000000D0000001
144 0 0 1 0 0D0000002
128 2 0 0 0 0 0 000010000D0000003
--- /dev/null
+//-----------------------------------------------------------
+//
+// Copyright (C) 2014 by the deal.II authors
+//
+// This file is subject to LGPL and may not be distributed
+// without copyright and license information. Please refer
+// to the file deal.II/doc/license.html for the text and
+// further information on this license.
+//
+//-----------------------------------------------------------
+
+// Create a BSpline surface, and test axis projection.
+
+#include "../tests.h"
+#include <fstream>
+#include <base/logstream.h>
+
+#include <deal.II/opencascade/utilities.h>
+#include <deal.II/opencascade/boundary_lib.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/grid_out.h>
+#include <deal.II/grid/grid_generator.h>
+
+#include <TopTools.hxx>
+#include <TopoDS_Shape.hxx>
+#include <TopoDS_Edge.hxx>
+#include <TopoDS_Face.hxx>
+#include <BRepFill.hxx>
+#include <Standard_Stream.hxx>
+
+using namespace OpenCASCADE;
+
+int main ()
+{
+ std::ofstream logfile("output");
+
+ // Create a bspline passing through the points
+ std::vector<Point<3> > pts;
+ pts.push_back(Point<3>(0,0,0));
+ pts.push_back(Point<3>(0,1,0));
+ pts.push_back(Point<3>(1,1,0));
+ pts.push_back(Point<3>(1,0,0));
+
+ TopoDS_Edge edge1 = interpolation_curve(pts);
+ for(unsigned int i=0; i<pts.size(); ++i)
+ pts[i] += Point<3>(0,0,1);
+ TopoDS_Edge edge2 = interpolation_curve(pts);
+
+ TopoDS_Face face = BRepFill::Face (edge1, edge2);
+
+
+ NormalToMeshProjectionBoundary<1,3> manifold(face);
+
+ Triangulation<1,3> tria;
+ GridGenerator::hyper_cube(tria);
+
+ tria.begin_active()->set_all_manifold_ids(1);
+ tria.set_manifold(1, manifold);
+
+ tria.refine_global(4);
+ GridOut gridout;
+ gridout.write_msh(tria, logfile);
+
+ return 0;
+}
--- /dev/null
+$NOD
+17
+1 0 0 0
+2 1 0 0
+3 0.5 1.125 0
+4 -0.213981 0.768714 0
+5 1.21398 0.768714 0
+6 -0.207824 0.356289 0
+7 0.0929699 1.04713 0
+8 0.90703 1.04713 0
+9 1.20782 0.356289 0
+10 -0.117453 0.170246 0
+11 -0.251277 0.561899 0
+12 -0.0871235 0.937269 0
+13 0.292624 1.10625 0
+14 0.707376 1.10625 0
+15 1.08712 0.937269 0
+16 1.25128 0.561899 0
+17 1.11745 0.170246 0
+$ENDNOD
+$ELM
+16
+1 1 0 0 2 1 10
+2 1 0 0 2 10 6
+3 1 0 0 2 6 11
+4 1 0 0 2 11 4
+5 1 0 0 2 4 12
+6 1 0 0 2 12 7
+7 1 0 0 2 7 13
+8 1 0 0 2 13 3
+9 1 0 0 2 3 14
+10 1 0 0 2 14 8
+11 1 0 0 2 8 15
+12 1 0 0 2 15 5
+13 1 0 0 2 5 16
+14 1 0 0 2 16 9
+15 1 0 0 2 9 17
+16 1 0 0 2 17 2
+$ENDELM
--- /dev/null
+//-----------------------------------------------------------
+//
+// Copyright (C) 2014 by the deal.II authors
+//
+// This file is subject to LGPL and may not be distributed
+// without copyright and license information. Please refer
+// to the file deal.II/doc/license.html for the text and
+// further information on this license.
+//
+//-----------------------------------------------------------
+
+// Create a BSpline surface, and test axis projection.
+
+#include "../tests.h"
+#include <fstream>
+#include <base/logstream.h>
+
+#include <deal.II/opencascade/utilities.h>
+#include <deal.II/opencascade/boundary_lib.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/grid_out.h>
+#include <deal.II/grid/grid_generator.h>
+
+#include <TopTools.hxx>
+#include <TopoDS_Shape.hxx>
+#include <TopoDS_Edge.hxx>
+#include <TopoDS_Face.hxx>
+#include <BRepFill.hxx>
+#include <Standard_Stream.hxx>
+
+using namespace OpenCASCADE;
+
+int main ()
+{
+ std::ofstream logfile("output");
+
+ // Create a bspline passing through the points
+ std::vector<Point<3> > pts;
+ pts.push_back(Point<3>(0,0,0));
+ pts.push_back(Point<3>(1,0,0));
+ pts.push_back(Point<3>(1,0,-1));
+ pts.push_back(Point<3>(0,0,-1));
+
+ TopoDS_Edge edge1 = interpolation_curve(pts);
+ for(unsigned int i=0; i<pts.size(); ++i)
+ pts[i] += Point<3>(0,1,0);
+ TopoDS_Edge edge2 = interpolation_curve(pts);
+
+ TopoDS_Face face = BRepFill::Face (edge1, edge2);
+
+
+ NormalToMeshProjectionBoundary<2,3> manifold(face);
+
+ Triangulation<2,3> tria;
+ GridGenerator::hyper_cube(tria);
+
+ tria.begin_active()->set_all_manifold_ids(1);
+ tria.set_manifold(1, manifold);
+
+ tria.refine_global(2);
+ GridOut gridout;
+ gridout.write_msh(tria, logfile);
+
+ return 0;
+}
--- /dev/null
+$NOD
+25
+1 0 0 0
+2 1 0 0
+3 0 1 0
+4 1 1 0
+5 0.5 0 0.245356
+6 0 0.5 0
+7 1 0.5 -4.50756e-17
+8 0.5 1 0.245356
+9 0.5 0.5 0.245356
+10 0.234806 0 0.15364
+11 0.789901 0 0.20399
+12 0 0.25 0
+13 0 0.75 0
+14 1 0.25 -2.75844e-17
+15 1 0.75 -3.12037e-17
+16 0.234806 1 0.15364
+17 0.789901 1 0.20399
+18 0.5 0.25 0.245356
+19 0.5 0.75 0.245356
+20 0.234806 0.5 0.15364
+21 0.789901 0.5 0.20399
+22 0.234806 0.25 0.15364
+23 0.789901 0.25 0.20399
+24 0.234806 0.75 0.15364
+25 0.789901 0.75 0.20399
+$ENDNOD
+$ELM
+16
+1 3 0 0 4 1 10 22 12
+2 3 0 0 4 10 5 18 22
+3 3 0 0 4 12 22 20 6
+4 3 0 0 4 22 18 9 20
+5 3 0 0 4 5 11 23 18
+6 3 0 0 4 11 2 14 23
+7 3 0 0 4 18 23 21 9
+8 3 0 0 4 23 14 7 21
+9 3 0 0 4 6 20 24 13
+10 3 0 0 4 20 9 19 24
+11 3 0 0 4 13 24 16 3
+12 3 0 0 4 24 19 8 16
+13 3 0 0 4 9 21 25 19
+14 3 0 0 4 21 7 15 25
+15 3 0 0 4 19 25 17 8
+16 3 0 0 4 25 15 4 17
+$ENDELM
--- /dev/null
+//-----------------------------------------------------------
+//
+// Copyright (C) 2014 by the deal.II authors
+//
+// This file is subject to LGPL and may not be distributed
+// without copyright and license information. Please refer
+// to the file deal.II/doc/license.html for the text and
+// further information on this license.
+//
+//-----------------------------------------------------------
+
+// Read the file side.iges, attach it to a projector, create a
+// single-cell Triangulation, ask the center of the single cell, and
+// then repeat by using the version which queries the manifold.
+
+#include "../tests.h"
+
+#include <deal.II/opencascade/boundary_lib.h>
+
+#include <fstream>
+#include <base/logstream.h>
+#include <grid/tria.h>
+#include <grid/tria_accessor.h>
+#include <grid/grid_out.h>
+#include <grid/tria_iterator.h>
+#include <grid/grid_generator.h>
+#include <grid/tria_boundary_lib.h>
+
+#include <gp_Pnt.hxx>
+#include <gp_Dir.hxx>
+#include <gp_Ax2.hxx>
+#include <GC_MakeCircle.hxx>
+#include <BRepBuilderAPI_MakeEdge.hxx>
+#include <BRepBuilderAPI_MakeFace.hxx>
+#include <BRepBuilderAPI_MakeWire.hxx>
+#include <TopoDS_Wire.hxx>
+#include <TopoDS_Edge.hxx>
+#include <TopoDS_Face.hxx>
+#include <BRepPrimAPI_MakeSphere.hxx>
+
+
+using namespace OpenCASCADE;
+
+int main ()
+{
+ initlog();
+
+ TopoDS_Shape sh = read_IGES(SOURCE_DIR "/iges_files/wigley.iges", 0.001);
+ std::vector<TopoDS_Face> faces;
+ std::vector<TopoDS_Edge> edges;
+ std::vector<TopoDS_Vertex> vertices;
+
+ extract_geometrical_shapes(sh, faces, edges, vertices);
+
+
+ // Create a boundary projector on the first face.
+ NormalProjectionBoundary<2,3> boundary(faces[0]);
+
+ // Create a Triangulation with a single cell
+ Triangulation<2,3> tria;
+ create_triangulation(faces[0], tria);
+
+ // Set the boundary
+ tria.set_manifold(1, boundary);
+ tria.begin()->set_all_manifold_ids(1);
+
+ deallog << "Ncells: " << tria.n_active_cells() << std::endl
+ << "Cell[0] center: " << tria.begin()->center() << std::endl
+ << "Projected center: " << tria.begin()->center(true) << std::endl;
+
+
+ return 0;
+}
+
--- /dev/null
+
+DEAL::Ncells: 1
+DEAL::Cell[0] center: 0.00000 -9.36751e-17 0.00000
+DEAL::Projected center: -2.24033e-17 0.113663 -0.0466831
--- /dev/null
+//-----------------------------------------------------------
+//
+// Copyright (C) 2014 by the deal.II authors
+//
+// This file is subject to LGPL and may not be distributed
+// without copyright and license information. Please refer
+// to the file deal.II/doc/license.html for the text and
+// further information on this license.
+//
+//-----------------------------------------------------------
+
+// Read the file side.iges, attach it to a projector, create a
+// single-cell Triangulation, ask the center of the single cell, and
+// then repeat by using the version which queries the manifold.
+
+#include "../tests.h"
+
+#include <deal.II/opencascade/boundary_lib.h>
+
+#include <fstream>
+#include <base/logstream.h>
+#include <grid/tria.h>
+#include <grid/tria_accessor.h>
+#include <grid/grid_out.h>
+#include <grid/tria_iterator.h>
+#include <grid/grid_generator.h>
+#include <grid/tria_boundary_lib.h>
+
+#include <gp_Pnt.hxx>
+#include <gp_Dir.hxx>
+#include <gp_Ax2.hxx>
+#include <GC_MakeCircle.hxx>
+#include <BRepBuilderAPI_MakeEdge.hxx>
+#include <BRepBuilderAPI_MakeFace.hxx>
+#include <BRepBuilderAPI_MakeWire.hxx>
+#include <TopoDS_Wire.hxx>
+#include <TopoDS_Edge.hxx>
+#include <TopoDS_Face.hxx>
+#include <BRepPrimAPI_MakeSphere.hxx>
+
+
+using namespace OpenCASCADE;
+
+int main ()
+{
+ initlog();
+
+ TopoDS_Shape sh = read_IGES(SOURCE_DIR "/iges_files/wigley.iges", 0.001);
+ std::vector<TopoDS_Face> faces;
+ std::vector<TopoDS_Edge> edges;
+ std::vector<TopoDS_Vertex> vertices;
+
+ extract_geometrical_shapes(sh, faces, edges, vertices);
+
+
+ // Create a boundary projector on the first face.
+ DirectionalProjectionBoundary<2,3> boundary(faces[0], Point<3>(0.,1.,0.));
+
+ // Create a Triangulation with a single cell
+ Triangulation<2,3> tria;
+ create_triangulation(faces[0], tria);
+
+ // Set the boundary
+ tria.set_manifold(1, boundary);
+ tria.begin()->set_all_manifold_ids(1);
+
+ deallog << "Ncells: " << tria.n_active_cells() << std::endl
+ << "Cell[0] center: " << tria.begin()->center() << std::endl
+ << "Projected center: " << tria.begin()->center(true) << std::endl;
+
+
+ return 0;
+}
+
--- /dev/null
+
+DEAL::Ncells: 1
+DEAL::Cell[0] center: 0.00000 -9.36751e-17 0.00000
+DEAL::Projected center: 0.00000 0.124760 0.00000
--- /dev/null
+//-----------------------------------------------------------
+//
+// Copyright (C) 2014 by the deal.II authors
+//
+// This file is subject to LGPL and may not be distributed
+// without copyright and license information. Please refer
+// to the file deal.II/doc/license.html for the text and
+// further information on this license.
+//
+//-----------------------------------------------------------
+
+// Read the file side.iges, attach it to a projector, create a
+// single-cell Triangulation, ask the center of the single cell, and
+// then repeat by using the version which queries the manifold.
+
+#include "../tests.h"
+
+#include <deal.II/opencascade/boundary_lib.h>
+
+#include <fstream>
+#include <base/logstream.h>
+#include <grid/tria.h>
+#include <grid/tria_accessor.h>
+#include <grid/grid_out.h>
+#include <grid/tria_iterator.h>
+#include <grid/grid_generator.h>
+#include <grid/tria_boundary_lib.h>
+
+#include <gp_Pnt.hxx>
+#include <gp_Dir.hxx>
+#include <gp_Ax2.hxx>
+#include <GC_MakeCircle.hxx>
+#include <BRepBuilderAPI_MakeEdge.hxx>
+#include <BRepBuilderAPI_MakeFace.hxx>
+#include <BRepBuilderAPI_MakeWire.hxx>
+#include <TopoDS_Wire.hxx>
+#include <TopoDS_Edge.hxx>
+#include <TopoDS_Face.hxx>
+#include <BRepPrimAPI_MakeSphere.hxx>
+
+
+using namespace OpenCASCADE;
+
+int main ()
+{
+ initlog();
+
+ TopoDS_Shape sh = read_IGES(SOURCE_DIR "/iges_files/wigley.iges", 0.001);
+ std::vector<TopoDS_Face> faces;
+ std::vector<TopoDS_Edge> edges;
+ std::vector<TopoDS_Vertex> vertices;
+
+ extract_geometrical_shapes(sh, faces, edges, vertices);
+
+
+ // Create a boundary projector on the first face.
+ NormalToMeshProjectionBoundary<2,3> boundary(faces[0]);
+
+ // Create a Triangulation with a single cell
+ Triangulation<2,3> tria;
+ create_triangulation(faces[0], tria);
+
+ // Set the boundary
+ tria.set_manifold(1, boundary);
+ tria.begin()->set_all_manifold_ids(1);
+
+ deallog << "Ncells: " << tria.n_active_cells() << std::endl
+ << "Cell[0] center: " << tria.begin()->center() << std::endl
+ << "Projected center: " << tria.begin()->center(true) << std::endl;
+
+
+ return 0;
+}
+
--- /dev/null
+
+DEAL::Ncells: 1
+DEAL::Cell[0] center: 0.00000 -9.36751e-17 0.00000
+DEAL::Projected center: 1.34026e-17 0.122612 -0.0205624
projects / dealii.git / commitdiff