Added interpolation curve boundary test.
Updated documentation.
--- /dev/null
+// ---------------------------------------------------------------------
+// $Id$
+//
+// Copyright (C) 2014 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+#ifndef __deal2__occ_boundary_lib_h
+#define __deal2__occ_boundary_lib_h
+
+#include <deal.II/base/config.h>
+
+#ifdef DEAL_II_WITH_OPENCASCADE
+
+#include <deal.II/grid/occ_utilities.h>
+#include <deal.II/grid/tria_boundary.h>
+
+DEAL_II_NAMESPACE_OPEN
+
+namespace OpenCASCADE
+{
+ /**
+ * @addtogroup OpenCASCADE
+ * @{
+ *
+ * A Boundary object based on OpenCASCADE TopoDS_Shape where new
+ * points are first computed using the FlatManifold class, and then
+ * projected in the normal direction using OpenCASCADE utilities.
+ *
+ * This class makes no assumptions on the shape you pass to it, and
+ * the topological dimension of the Manifold is inferred from the
+ * TopoDS_Shape itself. 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.
+ *
+ * This could happen, for example, if you are trying to use a shape
+ * of type TopoDS_Edge when projecting on a face. In this case, the
+ * vertices of the face would be collapsed to the edge, and your
+ * surrounding points would not be on lying on the given shape,
+ * raising an exception.
+ *
+ * @author Luca Heltai, Andrea Mola, 2011--2014.
+ */
+ template <int dim, int spacedim>
+ class NormalProjectionBoundary : public Boundary<dim,spacedim> {
+ public:
+ NormalProjectionBoundary(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 #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 normal
+ * projection algorithms.
+ */
+ virtual Point<spacedim>
+ project_to_manifold (const std::vector<Point<spacedim> > &surrounding_points,
+ const Point<spacedim> &candidate) const;
+
+ /**
+ * Exception thrown when the point specified as argument does not
+ * lie between #tolerance from the given TopoDS_Shape.
+ */
+ DeclException1 (ExcPointNotOnManifold,
+ Point<spacedim>,
+ <<"The point [ "<<arg1<<" ] is not on the manifold.");
+
+ private:
+ /**
+ * The topological shape which is used internally to project
+ * points. You can construct one 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;
+
+ /**
+ * Relative tolerance used by this class to compute distances.
+ */
+ const double tolerance;
+ };
+
+}
+
+/*@}*/
+
+DEAL_II_NAMESPACE_CLOSE
+
+
+#endif // DEAL_II_WITH_OPENCASCADE
+
+/*------------------------------ occ_boundary_lib.h ------------------------------*/
+#endif
+/*------------------------------ occ_boundary_lib.h ------------------------------*/
#ifndef __deal2__occ_utilities_h
#define __deal2__occ_utilities_h
+#include <deal.II/base/config.h>
+
#ifdef DEAL_II_WITH_OPENCASCADE
#include <string>
#include <Geom_Curve.hxx>
#include <gp_Pnt.hxx>
-#include <base/point.h>
+#include <deal.II/base/point.h>
DEAL_II_NAMESPACE_OPEN
/**
* Creates a 3D smooth BSpline curve passing through the points in
- * the assigned vector. The points are reordered internally
- * according to their scalar product with the direction, if
- * direction is different from zero, otherwise they are used as
+ * the assigned vector, and store it in the returned TopoDS_Shape
+ * (which is of type TopoDS_Edge). The points are reordered
+ * internally according to their scalar product with the direction,
+ * if direction is different from zero, otherwise they are used as
* passed. Notice that this function changes the input points if
* required by the algorithm.
*
* This class is used to interpolate a BsplineCurve passing through
- * an array of points, with a C2 Continuity.
+ * an array of points, with a C2 Continuity. If the optional
+ * parameter #closed is set to true, then the curve will be C2 at
+ * all points execpt the first (where only C1 continuity will be
+ * given), and it will be a closed curve.
+ *
+ * The curve is garanteed to be at distance #tolerance from the
+ * input points. If the algorithm fails in generating such a curve,
+ * an exception is thrown.
*/
- TopoDS_Shape interpolation_curve(std::vector<dealii::Point<3> > &curve_points,
- const dealii::Point<3> direction=dealii::Point<3>(),
- const double tolerance=1e-7);
+ TopoDS_Edge interpolation_curve(std::vector<dealii::Point<3> > &curve_points,
+ const dealii::Point<3> direction=dealii::Point<3>(),
+ const bool closed=false,
+ const double tolerance=1e-7);
+
/**
* Get the closest point to the given topological shape. If the
/**
* Convert OpenCASCADE point into a Point<3>.
*/
- Point<3> Pnt(const gp_Pnt &p);
+ inline Point<3> Pnt(const gp_Pnt &p);
/**
* Convert Point<3> into OpenCASCADE point.
*/
- gp_Pnt Pnt(const Point<3> &p);
+ inline gp_Pnt Pnt(const Point<3> &p);
/**
* Sort two points according to their scalar product with
* direction. If the norm of the direction is zero, then use
- * lexycographical ordering.
+ * lexycographical ordering. The optional parameter is used as a
+ * relative tolerance when comparing objects.
*/
inline bool point_compare(const dealii::Point<3> &p1, const dealii::Point<3> &p2,
const dealii::Point<3> direction=Point<3>(),
manifold.cc
manifold_lib.cc
occ_utilities.cc
+ occ_boundary_lib.cc
persistent_tria.cc
tria_accessor.cc
tria_boundary.cc
intergrid_map.inst.in
manifold.inst.in
manifold_lib.inst.in
+ occ_boundary_lib.inst.in
tria_accessor.inst.in
tria_boundary.inst.in
tria_boundary_lib.inst.in
--- /dev/null
+#include <deal.II/grid/occ_boundary_lib.h>
+
+#ifdef DEAL_II_WITH_OPENCASCADE
+
+DEAL_II_NAMESPACE_OPEN
+
+namespace OpenCASCADE
+{
+ template <int dim, int spacedim>
+ NormalProjectionBoundary<dim,spacedim>::NormalProjectionBoundary(const TopoDS_Shape sh,
+ const double tolerance) :
+ sh(sh),
+ tolerance(tolerance)
+ {
+ Assert(spacedim == 3, ExcImpossibleInDim(spacedim));
+ }
+
+
+ template <int dim, int spacedim>
+ Point<spacedim> NormalProjectionBoundary<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]) < (surrounding_points[i].norm()>0 ?
+ tolerance*surrounding_points[i].norm() :
+ tolerance),
+ ExcPointNotOnManifold(surrounding_points[i]));
+
+ return closest_point(sh, candidate, out_shape, u, v);
+ }
+
+ // Explicit instantiations
+#include "occ_boundary_lib.inst"
+
+} // end namespace OpenCASCADE
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif
--- /dev/null
+// ---------------------------------------------------------------------
+// $Id$
+//
+// Copyright (C) 1998 - 2013 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+
+for (deal_II_dimension : DIMENSIONS)
+ {
+ template class NormalProjectionBoundary<deal_II_dimension, 3>;
+ }
+
+
+
#include <deal.II/grid/occ_utilities.h>
-#include <deal.II/base/config.h>
#ifdef DEAL_II_WITH_OPENCASCADE
}
inline bool point_compare(const dealii::Point<3> &p1, const dealii::Point<3> &p2,
- const dealii::Point<3> direction=Point<3>(),
- const double tolerance=1e-10)
+ const dealii::Point<3> direction,
+ const double tolerance)
{
- const double rel_tol=std::max(p1.norm(), p2.norm())*tolerance;
+ const double rel_tol=std::max(tolerance, std::max(p1.norm(), p2.norm())*tolerance);
if(direction.norm() > 0.0)
return (p1*direction < p2*direction-rel_tol);
else
- TopoDS_
Shape interpolation_curve(std::vector<Point<3> > &curve_points,
+ TopoDS_
Edge interpolation_curve(std::vector<Point<3> > &curve_points,
const Point<3> direction,
+ const bool closed,
const double tolerance)
{
}
- GeomAPI_Interpolate bspline_generator(vertices, false, tolerance);
+ GeomAPI_Interpolate bspline_generator(vertices, closed, tolerance);
bspline_generator.Perform();
Assert( (bspline_generator.IsDone()), ExcMessage("Interpolated bspline generation failed"));
Handle(Geom_BSplineCurve) bspline = bspline_generator.Curve();
- TopoDS_Shape out_shape = BRepBuilderAPI_MakeEdge(bspline);
+ TopoDS_Edge out_shape = BRepBuilderAPI_MakeEdge(bspline);
return out_shape;
}
--- /dev/null
+
+//---------------------------- occ_circle_projection.cc ---------------------------
+// $Id: testsuite.html 13373 2006-07-13 13:12:08Z kanschat $
+// Version: $Name$
+//
+// Copyright (C) 2005 by the deal.II authors
+//
+// This file is subject to QPL 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.
+//
+//---------------------------- occ_circle_projection.cc ---------------------------
+
+
+// Create a circle, a Triangulation, and try to project normally on
+// it.
+
+#include "../tests.h"
+
+#include <deal.II/grid/occ_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>
+
+using namespace OpenCASCADE;
+
+int main ()
+{
+ std::ofstream logfile("output");
+ deallog.attach(logfile);
+ deallog.depth_console(0);
+
+ // The circle passing through the
+ // vertices of the unit square
+ gp_Dir z_axis(0.,0.,1.);
+ gp_Pnt center(.5,.5,0.);
+ gp_Ax2 axis(center, z_axis);
+ Standard_Real radius(std::sqrt(2.)/2.);
+
+ Handle(Geom_Curve) circle = GC_MakeCircle(axis, radius);
+ TopoDS_Edge edge = BRepBuilderAPI_MakeEdge(circle);
+
+
+ // Create a boundary projector.
+ NormalProjectionBoundary<2,3> boundary_line(edge);
+
+ // This one is for checking: This
+ // is what deal.II would do for a
+ // circle.
+ HyperBallBoundary<2,3> boundary_line_deal (Point<3>(.5,.5,0),
+ std::sqrt(2.)/2.);
+
+
+
+ // The unit square.
+ Triangulation<2,3> tria;
+ GridGenerator::hyper_cube(tria);
+
+ // Set the exterior boundary
+ tria.set_boundary(0, boundary_line);
+ // Uncomment the following line to
+ // see the difference.
+// tria.set_boundary(0, boundary_line_deal);
+
+
+ // This is here to ignore the
+ // points created in the interior
+ // of the face.
+ tria.begin()->set_material_id(1);
+
+ // We refine twice, and expect the
+ // outer points to end up on the
+ // circle.
+ tria.refine_global(2);
+
+
+ // You can open the generated file
+ // with paraview.
+ GridOut gridout;
+ gridout.write_ucd (tria, logfile);
+
+ return 0;
+}
+
--- /dev/null
+
+25 16 0 0 0
+1 0.00000 0.00000 0.00000
+2 1.00000 0.00000 0.00000
+3 0.00000 1.00000 0.00000
+4 1.00000 1.00000 0.00000
+5 0.500000 -0.207107 0.00000
+6 -0.207107 0.500000 0.00000
+7 1.20711 0.500000 0.00000
+8 0.500000 1.20711 0.00000
+9 0.500000 0.500000 0.00000
+10 0.229402 -0.153281 0.00000
+11 0.770598 -0.153281 0.00000
+12 -0.153281 0.229402 0.00000
+13 -0.153281 0.770598 0.00000
+14 1.15328 0.229402 0.00000
+15 1.15328 0.770598 0.00000
+16 0.229402 1.15328 0.00000
+17 0.770598 1.15328 0.00000
+18 0.500000 0.146447 0.00000
+19 0.500000 0.853553 0.00000
+20 0.146447 0.500000 0.00000
+21 0.853553 0.500000 0.00000
+22 0.189433 0.189433 0.00000
+23 0.810567 0.189433 0.00000
+24 0.189433 0.810567 0.00000
+25 0.810567 0.810567 0.00000
+1 1 quad 1 10 22 12
+2 1 quad 10 5 18 22
+3 1 quad 12 22 20 6
+4 1 quad 22 18 9 20
+5 1 quad 5 11 23 18
+6 1 quad 11 2 14 23
+7 1 quad 18 23 21 9
+8 1 quad 23 14 7 21
+9 1 quad 6 20 24 13
+10 1 quad 20 9 19 24
+11 1 quad 13 24 16 3
+12 1 quad 24 19 8 16
+13 1 quad 9 21 25 19
+14 1 quad 21 7 15 25
+15 1 quad 19 25 17 8
+16 1 quad 25 15 4 17
--- /dev/null
+
+//---------------------------- occ_circle_projection.cc ---------------------------
+// $Id: testsuite.html 13373 2006-07-13 13:12:08Z kanschat $
+// Version: $Name$
+//
+// Copyright (C) 2005 by the deal.II authors
+//
+// This file is subject to QPL 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.
+//
+//---------------------------- occ_circle_projection.cc ---------------------------
+
+
+// Create a Triangulation, interpolate its boundary points to a close
+// smooth BSpline, and use that as a Boundary Descriptor.
+
+#include "../tests.h"
+
+#include <deal.II/grid/occ_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>
+
+using namespace OpenCASCADE;
+
+int main ()
+{
+ std::ofstream logfile("output");
+ deallog.attach(logfile);
+ deallog.depth_console(0);
+
+ // The curve passing through the vertices of the unit square.
+ std::vector<Point<3> > vertices;
+ vertices.push_back(Point<3>(0,0,0));
+ vertices.push_back(Point<3>(1,0,0));
+ vertices.push_back(Point<3>(1,1,0));
+ vertices.push_back(Point<3>(0,1,0));
+ TopoDS_Shape shape = interpolation_curve(vertices, Point<3>(), true);
+
+ // Create a boundary projector.
+ NormalProjectionBoundary<2,3> boundary_line(shape);
+
+ // The unit square.
+ Triangulation<2,3> tria;
+ GridGenerator::hyper_cube(tria);
+
+ // Set the exterior boundary
+ tria.set_boundary(0, boundary_line);
+
+ // This is here to ignore the points created in the interior of the
+ // face.
+ tria.begin()->set_material_id(1);
+
+ // We refine twice, and expect the outer points to end up on a
+ // smooth curve interpolating the square vertices.
+ tria.refine_global(2);
+
+
+ // You can open the generated file with paraview.
+ GridOut gridout;
+ gridout.write_msh (tria, logfile);
+
+ std::ofstream of("/tmp/grid.msh");
+ gridout.write_msh(tria, of);
+
+ return 0;
+}
+
--- /dev/null
+
+$NOD
+25
+1 0.00000 0.00000 0.00000
+2 1.00000 0.00000 0.00000
+3 0.00000 1.00000 0.00000
+4 1.00000 1.00000 0.00000
+5 0.439917 -0.236739 0.00000
+6 -0.310022 0.473832 0.00000
+7 1.14617 0.492348 0.00000
+8 0.488859 1.16307 0.00000
+9 0.470616 0.486564 0.00000
+10 0.207735 -0.140071 0.00000
+11 0.764326 -0.218873 0.00000
+12 -0.202060 0.202837 0.00000
+13 -0.224611 0.783020 0.00000
+14 1.10302 0.237503 0.00000
+15 1.12328 0.758622 0.00000
+16 0.234185 1.11236 0.00000
+17 0.760208 1.13858 0.00000
+18 0.455267 0.124913 0.00000
+19 0.479738 0.824816 0.00000
+20 0.0802967 0.480198 0.00000
+21 0.808395 0.489456 0.00000
+22 0.142719 0.173942 0.00000
+23 0.773464 0.171896 0.00000
+24 0.152383 0.790482 0.00000
+25 0.784659 0.794181 0.00000
+$ENDNOD
+$ELM
+16
+1 3 1 0 4 1 10 22 12
+2 3 1 0 4 10 5 18 22
+3 3 1 0 4 12 22 20 6
+4 3 1 0 4 22 18 9 20
+5 3 1 0 4 5 11 23 18
+6 3 1 0 4 11 2 14 23
+7 3 1 0 4 18 23 21 9
+8 3 1 0 4 23 14 7 21
+9 3 1 0 4 6 20 24 13
+10 3 1 0 4 20 9 19 24
+11 3 1 0 4 13 24 16 3
+12 3 1 0 4 24 19 8 16
+13 3 1 0 4 9 21 25 19
+14 3 1 0 4 21 7 15 25
+15 3 1 0 4 19 25 17 8
+16 3 1 0 4 25 15 4 17
+$ENDELM
projects / dealii.git / commitdiff