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 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;
+
/**
- * Exception thrown when the point specified as argument does not
- * lie between #tolerance from the given TopoDS_Shape.
+ * Relative tolerance used by this class to compute distances.
*/
- DeclException1 (ExcPointNotOnManifold,
- Point<spacedim>,
- <<"The point [ "<<arg1<<" ] is not on the manifold.");
+ const double tolerance;
+ };
+
+ /**
+ * A Boundary object based on OpenCASCADE TopoDS_Shape where new
+ * points are first computed using the FlatManifold class, and then
+ * projected along the direction given at construction time, 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.
+ *
+ * @author Luca Heltai, Andrea Mola, 2011--2014.
+ */
+ template <int dim, int spacedim>
+ class AxisProjectionBoundary : public Boundary<dim,spacedim> {
+ public:
+ AxisProjectionBoundary(const TopoDS_Shape sh,
+ const Point<3> direction,
+ 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;
private:
/**
*/
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.
*/
double &v,
const double tolerance=1e-7);
+ /**
+ * Intersect a line passing through the given #origin point along
+ * #direction and the given topological shape. If there is more than
+ * one intersection, it will return the closest one.
+ *
+ * The optional #tolerance parameter is used to compute distances.
+ */
+ Point<3> axis_intersection(const TopoDS_Shape in_shape,
+ const Point<3> origin,
+ const Point<3> direction,
+ const double tolerance=1e-7);
+
+
/**
* Convert OpenCASCADE point into a Point<3>.
*/
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);
-
-}
+
+
+ /**
+ * Exception thrown when the point specified as argument does not
+ * lie between #tolerance from the given TopoDS_Shape.
+ */
+ DeclException1 (ExcPointNotOnManifold,
+ Point<3>,
+ <<"The point [ "<<arg1<<" ] is not on the manifold.");
+
+ /**
+ * Exception thrown when the point specified as argument cannot be
+ * projected to the manifold.
+ */
+ DeclException1 (ExcProjectionFailed,
+ Point<3>,
+ <<"Projection of point [ "<< arg1
+ << " ] failed.");
+
+ /**
+ * Thrown when internal OpenCASCADE utilities fail to return the OK
+ * status.
+ */
+ DeclException1 (ExcOCCError,
+ IFSelect_ReturnStatus,
+ <<"An OpenCASCADE routine failed with return status "
+ <<arg1);
+}
/*@}*/
DEAL_II_NAMESPACE_CLOSE
namespace OpenCASCADE
{
+
+ /*============================== NormalProjectionBoundary ==============================*/
template <int dim, int spacedim>
NormalProjectionBoundary<dim,spacedim>::NormalProjectionBoundary(const TopoDS_Shape sh,
const double tolerance) :
.distance(surrounding_points[i]) < (surrounding_points[i].norm()>0 ?
tolerance*surrounding_points[i].norm() :
tolerance),
- ExcPointNotOnManifold(surrounding_points[i]));
+ // ExcPointNotOnManifold(surrounding_points[i]));
+ ExcMessage("Points are not on manifold"));
return closest_point(sh, candidate, out_shape, u, v);
}
+
+ /*============================== AxisProjectionBoundary ==============================*/
+ template <int dim, int spacedim>
+ AxisProjectionBoundary<dim,spacedim>::AxisProjectionBoundary(const TopoDS_Shape sh,
+ const Point<3> direction,
+ const double tolerance) :
+ sh(sh),
+ direction(direction),
+ tolerance(tolerance)
+ {
+ Assert(spacedim == 3, ExcImpossibleInDim(spacedim));
+ }
+
+
+ template <int dim, int spacedim>
+ Point<spacedim> AxisProjectionBoundary<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 axis_intersection(sh, candidate, direction, tolerance);
+ }
+
// Explicit instantiations
#include "occ_boundary_lib.inst"
for (deal_II_dimension : DIMENSIONS)
{
template class NormalProjectionBoundary<deal_II_dimension, 3>;
+ template class AxisProjectionBoundary<deal_II_dimension, 3>;
}
#include <TColStd_SequenceOfTransient.hxx>
#include <TColStd_HSequenceOfTransient.hxx>
#include <TopExp_Explorer.hxx>
-// #include <gp_Pnt.hxx>
-// #include <gp_Vec.hxx>
+#include <gp_Pnt.hxx>
+#include <gp_Lin.hxx>
+#include <gp_Vec.hxx>
#include <GeomAPI_ProjectPointOnSurf.hxx>
#include <GeomAPI_ProjectPointOnCurve.hxx>
+#include <IntCurvesFace_ShapeIntersector.hxx>
// #include <Standard_Real.hxx>
// #include <Standard_Integer.hxx>
// #include <BRep_Tool.hxx>
}
+ Point<3> axis_intersection(const TopoDS_Shape in_shape,
+ const Point<3> origin,
+ const Point<3> direction,
+ const double tolerance)
+ {
+ // translating original Point<dim> to gp point
+
+ gp_Pnt P0 = Pnt(origin);
+ gp_Ax1 gpaxis(P0, gp_Dir(direction(0), direction(1), direction(2)));
+ gp_Lin line(gpaxis);
+
+ // destination point
+ gp_Pnt Pproj(0.0,0.0,0.0);
+
+ // we prepare now the surface for the projection we get the whole
+ // shape from the iges model
+ IntCurvesFace_ShapeIntersector Inters;
+ Inters.Load(in_shape,tolerance);
+
+ Inters.PerformNearest(line,-RealLast(),+RealLast());
+ Assert(Inters.IsDone(), ExcMessage("Could not project point."));
+ return Pnt(Inters.Pnt(1));
+ }
TopoDS_Edge interpolation_curve(std::vector<Point<3> > &curve_points,
const Point<3> direction,
projects / dealii.git / commitdiff