Manifolds have names ending in 'Manifold': now that Boundary is gone these names
do not make sense.
For backwards compatibility, the old classes now simply inherit from the
correctly-named Manifold versions. This requires marking some things that were
private now as protected.
--- /dev/null
+Changed: The OpenCASCADE Manifold classes with names ending in Boundary (i.e.,
+NormalProjectionBoundary, DirectionalProjectionBoundary, and
+NormalToMeshProjectionBoundary) have been deprecated in favor of renamed classes
+ending in Manifold (i.e., NormalProjectionManifold,
+DirectionalProjectionManifold, and NormalToMeshProjectionManifold).
+<br>
+(David Wells, 2018/05/16)
#ifdef DEAL_II_WITH_OPENCASCADE
-#include <deal.II/opencascade/utilities.h>
-#include <deal.II/grid/manifold.h>
-
-// opencascade needs "HAVE_CONFIG_H" to be exported...
-#define HAVE_CONFIG_H
-#include <Adaptor3d_Curve.hxx>
-#include <Adaptor3d_HCurve.hxx>
-#include <BRepAdaptor_Curve.hxx>
-#undef HAVE_CONFIG_H
+#include <deal.II/opencascade/manifold_lib.h>
DEAL_II_NAMESPACE_OPEN
* not be lying on the given shape, raising an exception.
*
* @author Luca Heltai, Andrea Mola, 2011--2014.
+ *
+ * @deprecated Use NormalProjectionManifold instead, which is identical to
+ * this class but satisfies the modern Manifold-based naming convention.
*/
template <int dim, int spacedim>
- class NormalProjectionBoundary : public FlatManifold<dim,spacedim>
+ class DEAL_II_DEPRECATED NormalProjectionBoundary
+ : public NormalProjectionManifold<dim, spacedim>
{
public:
-
/**
- * The standard constructor takes a generic TopoDS_Shape @p sh, and a
- * tolerance used to compute distances internally.
- *
- * The TopoDS_Shape can be arbitrary, i.e., a collection of shapes, faces,
- * edges or a single face or edge.
+ * Inherit all constructors.
*/
- NormalProjectionBoundary(const TopoDS_Shape &sh,
- const double tolerance=1e-7);
+ using NormalProjectionManifold<dim, spacedim>::NormalProjectionManifold;
/**
* Clone the current Manifold.
*/
virtual std::unique_ptr<Manifold<dim,spacedim> > clone() const override;
-
- /**
- * 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 normal projection
- * algorithms.
- */
- virtual Point<spacedim>
- project_to_manifold (const ArrayView<const Point<spacedim>> &surrounding_points,
- const Point<spacedim> &candidate) const override;
-
-
- 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;
-
- /**
- * Relative tolerance used by this class to compute distances.
- */
- const double tolerance;
};
/**
* intersect the shape. An exception is thrown when this happens.
*
* @author Luca Heltai, Andrea Mola, 2011--2014.
+ *
+ * @deprecated Use DirectionalProjectionManifold instead, which is identical to
+ * this class but satisfies the modern Manifold-based naming convention.
*/
template <int dim, int spacedim>
- class DirectionalProjectionBoundary : public FlatManifold<dim,spacedim>
+ class DEAL_II_DEPRECATED DirectionalProjectionBoundary
+ : public DirectionalProjectionManifold<dim, spacedim>
{
public:
/**
- * Construct a Manifold object which will project points on the
- * TopoDS_Shape @p sh, along the given @p direction.
+ * Inherit all constructors.
*/
- DirectionalProjectionBoundary(const TopoDS_Shape &sh,
- const Tensor<1,spacedim> &direction,
- const double tolerance=1e-7);
+ using DirectionalProjectionManifold<dim, spacedim>::DirectionalProjectionManifold;
/**
* Clone the current Manifold.
*/
virtual std::unique_ptr<Manifold<dim,spacedim> > clone() const override;
-
- /**
- * 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 ArrayView<const Point<spacedim>> &surrounding_points,
- const Point<spacedim> &candidate) const override;
-
- 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 Tensor<1, spacedim> direction;
-
- /**
- * Relative tolerance used by this class to compute distances.
- */
- const double tolerance;
};
* happens.
*
* @author Luca Heltai, Andrea Mola, 2011--2014.
- */
- template <int dim, int spacedim>
- class NormalToMeshProjectionBoundary : public FlatManifold<dim,spacedim>
- {
- public:
- /**
- * Construct a Manifold object which will project points on the
- * TopoDS_Shape @p sh, along a direction which is approximately normal to
- * the mesh cell.
- */
- NormalToMeshProjectionBoundary(const TopoDS_Shape &sh,
- const double tolerance=1e-7);
-
- /**
- * Clone the current Manifold.
- */
- virtual std::unique_ptr<Manifold<dim,spacedim> > clone() const override;
-
- /**
- * 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.
- */
- virtual Point<spacedim>
- project_to_manifold (const ArrayView<const Point<spacedim>> &surrounding_points,
- const Point<spacedim> &candidate) const override;
-
- 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;
-
- /**
- * Relative tolerance used by this class to compute distances.
- */
- const double tolerance;
- };
-
- /**
- * A Manifold object based on OpenCASCADE TopoDS_Shape objects which have
- * topological dimension equal to one (TopoDS_Edge or TopoDS_Wire) where new
- * points are located at the arclength average of the surrounding points. If
- * the given TopoDS_Shape can be casted to a periodic (closed) curve, then
- * this information is used internally to set the periodicity of the base
- * ChartManifold class.
- *
- * This class can only work on TopoDS_Edge or TopoDS_Wire objects, and it
- * only makes sense when spacedim is three. If you use an object of
- * topological dimension different from one, an exception is throw.
*
- * In debug mode there is an additional sanity check to make sure that the
- * surrounding points 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.
+ * @deprecated Use NormalToMeshProjectionManifold instead, which is identical to
+ * this class but satisfies the modern Manifold-based naming convention.
*/
template <int dim, int spacedim>
- class ArclengthProjectionLineManifold : public ChartManifold<dim,spacedim,1>
+ class DEAL_II_DEPRECATED NormalToMeshProjectionBoundary
+ : public NormalToMeshProjectionManifold<dim, spacedim>
{
public:
/**
- * Default constructor with a TopoDS_Edge.
+ * Inherit all constructors.
*/
- ArclengthProjectionLineManifold(const TopoDS_Shape &sh,
- const double tolerance=1e-7);
+ using NormalToMeshProjectionManifold<dim, spacedim>::NormalToMeshProjectionManifold;
/**
* Clone the current Manifold.
*/
virtual std::unique_ptr<Manifold<dim,spacedim> > clone() const override;
-
- /**
- * Given a point on real space, find its arclength parameter. Throws an
- * error in debug mode, if the point is not on the TopoDS_Edge given at
- * construction time.
- */
- virtual Point<1>
- pull_back(const Point<spacedim> &space_point) const override;
-
- /**
- * Given an arclength parameter, find its image in real space.
- */
- virtual Point<spacedim>
- push_forward(const Point<1> &chart_point) const override;
-
- private:
- /**
- * The actual shape used to build this object.
- */
- const TopoDS_Shape sh;
-
- /**
- * A Curve adaptor. This is the one which is used in the computations, and
- * it points to the right one above.
- */
- Handle_Adaptor3d_HCurve curve;
-
- /**
- * Relative tolerance used in all internal computations.
- */
- const double tolerance;
-
- /**
- * The total length of the curve. This is also used as a period if the
- * edge is periodic.
- */
- const double length;
- };
-
- /**
- * Manifold description for the face of a CAD imported usign OpenCASCADE.
- *
- * @ingroup manifold
- *
- * @author Andrea Mola, Mauro Bardelloni, 2016
- */
- template <int dim, int spacedim>
- class NURBSPatchManifold : public ChartManifold<dim, spacedim, 2>
- {
- public:
- /**
- * The constructor takes an OpenCASCADE TopoDS_Face @p face and an optional
- * @p tolerance. This class uses the interval OpenCASCADE variables @var u,
- * @var v to descrive the manifold.
- */
- NURBSPatchManifold(const TopoDS_Face &face, const double tolerance = 1e-7);
-
- /**
- * Clone the current Manifold.
- */
- virtual std::unique_ptr<Manifold<dim,spacedim> > clone() const override;
-
- /**
- * Pull back the given point from the Euclidean space. Will return the uv
- * coordinates associated with the point @p space_point.
- */
- virtual Point<2>
- pull_back(const Point<spacedim> &space_point) const override;
-
- /**
- * Given a @p chart_point in the uv coordinate system, this method returns the
- * Euclidean coordinates associated.
- */
- virtual Point<spacedim>
- push_forward(const Point<2> &chart_point) const override;
-
- /**
- * Given a point in the spacedim dimensional Euclidean space, this
- * method returns the derivatives of the function $F$ that maps from
- * the uv coordinate system to the Euclidean coordinate
- * system. In other words, it is a matrix of size
- * $\text{spacedim}\times\text{chartdim}$.
- *
- * This function is used in the computations required by the
- * get_tangent_vector() function.
- *
- * Refer to the general documentation of this class for more information.
- */
- virtual
- DerivativeForm<1,2,spacedim>
- push_forward_gradient(const Point<2> &chart_point) const override;
-
- private:
- /**
- * Return a tuple representing the minimum and maximum values of u
- * and v. Precisely, it returns (u_min, u_max, v_min, v_max)
- */
- std::tuple<double, double, double, double>
- get_uv_bounds() const;
-
- /**
- * An OpenCASCADE TopoDS_Face @p face given by the CAD.
- */
- TopoDS_Face face;
-
- /**
- * Tolerance used by OpenCASCADE to identify points in each
- * operation.
- */
- double tolerance;
};
-
}
/*@}*/
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2014 - 2018 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 dealii_occ_manifold_lib_h
+#define dealii_occ_manifold_lib_h
+
+#include <deal.II/base/config.h>
+
+#ifdef DEAL_II_WITH_OPENCASCADE
+
+#include <deal.II/opencascade/utilities.h>
+#include <deal.II/grid/manifold.h>
+
+// opencascade needs "HAVE_CONFIG_H" to be exported...
+#define HAVE_CONFIG_H
+#include <Adaptor3d_Curve.hxx>
+#include <Adaptor3d_HCurve.hxx>
+#include <BRepAdaptor_Curve.hxx>
+#undef HAVE_CONFIG_H
+
+DEAL_II_NAMESPACE_OPEN
+
+/**
+ * @addtogroup OpenCASCADE
+ * @{
+ */
+
+namespace OpenCASCADE
+{
+ /**
+ * A Manifold 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 are 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 lying on the given shape, raising an exception.
+ *
+ * @author Luca Heltai, Andrea Mola, 2011--2014.
+ */
+ template <int dim, int spacedim>
+ class NormalProjectionManifold : public FlatManifold<dim,spacedim>
+ {
+ public:
+
+ /**
+ * The standard constructor takes a generic TopoDS_Shape @p sh, and a
+ * tolerance used to compute distances internally.
+ *
+ * The TopoDS_Shape can be arbitrary, i.e., a collection of shapes, faces,
+ * edges or a single face or edge.
+ */
+ NormalProjectionManifold(const TopoDS_Shape &sh,
+ const double tolerance=1e-7);
+
+ /**
+ * Clone the current Manifold.
+ */
+ virtual std::unique_ptr<Manifold<dim,spacedim> > clone() const override;
+
+ /**
+ * 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 normal projection
+ * algorithms.
+ */
+ virtual Point<spacedim>
+ project_to_manifold (const ArrayView<const Point<spacedim>> &surrounding_points,
+ const Point<spacedim> &candidate) const override;
+
+
+ protected:
+ /**
+ * 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;
+
+ /**
+ * Relative tolerance used by this class to compute distances.
+ */
+ const double tolerance;
+ };
+
+ /**
+ * A Manifold 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 them onto the manifold along the
+ * direction specified 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.
+ *
+ * Notice that this type of Manifold 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 direction you use at construction time does not
+ * intersect the shape. An exception is thrown when this happens.
+ *
+ * @author Luca Heltai, Andrea Mola, 2011--2014.
+ */
+ template <int dim, int spacedim>
+ class DirectionalProjectionManifold : public FlatManifold<dim,spacedim>
+ {
+ public:
+ /**
+ * Construct a Manifold object which will project points on the
+ * TopoDS_Shape @p sh, along the given @p direction.
+ */
+ DirectionalProjectionManifold(const TopoDS_Shape &sh,
+ const Tensor<1,spacedim> &direction,
+ const double tolerance=1e-7);
+
+ /**
+ * Clone the current Manifold.
+ */
+ virtual std::unique_ptr<Manifold<dim,spacedim> > clone() const override;
+
+ /**
+ * 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 ArrayView<const Point<spacedim>> &surrounding_points,
+ const Point<spacedim> &candidate) const override;
+
+ protected:
+ /**
+ * 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 Tensor<1, spacedim> direction;
+
+ /**
+ * Relative tolerance used by this class to compute distances.
+ */
+ const double tolerance;
+ };
+
+
+ /**
+ * A Manifold 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 them 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 one face. If that is not the case,
+ * an Exception 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 Manifold 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
+ * happens.
+ *
+ * @author Luca Heltai, Andrea Mola, 2011--2014.
+ */
+ template <int dim, int spacedim>
+ class NormalToMeshProjectionManifold : public FlatManifold<dim,spacedim>
+ {
+ public:
+ /**
+ * Construct a Manifold object which will project points on the
+ * TopoDS_Shape @p sh, along a direction which is approximately normal to
+ * the mesh cell.
+ */
+ NormalToMeshProjectionManifold(const TopoDS_Shape &sh,
+ const double tolerance=1e-7);
+
+ /**
+ * Clone the current Manifold.
+ */
+ virtual std::unique_ptr<Manifold<dim,spacedim> > clone() const override;
+
+ /**
+ * 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.
+ */
+ virtual Point<spacedim>
+ project_to_manifold (const ArrayView<const Point<spacedim>> &surrounding_points,
+ const Point<spacedim> &candidate) const override;
+
+ protected:
+ /**
+ * 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;
+
+ /**
+ * Relative tolerance used by this class to compute distances.
+ */
+ const double tolerance;
+ };
+
+ /**
+ * A Manifold object based on OpenCASCADE TopoDS_Shape objects which have
+ * topological dimension equal to one (TopoDS_Edge or TopoDS_Wire) where new
+ * points are located at the arclength average of the surrounding points. If
+ * the given TopoDS_Shape can be casted to a periodic (closed) curve, then
+ * this information is used internally to set the periodicity of the base
+ * ChartManifold class.
+ *
+ * This class can only work on TopoDS_Edge or TopoDS_Wire objects, and it
+ * only makes sense when spacedim is three. If you use an object of
+ * topological dimension different from one, an exception is throw.
+ *
+ * In debug mode there is an additional sanity check to make sure that the
+ * surrounding points 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 ArclengthProjectionLineManifold : public ChartManifold<dim,spacedim,1>
+ {
+ public:
+ /**
+ * Default constructor with a TopoDS_Edge.
+ */
+ ArclengthProjectionLineManifold(const TopoDS_Shape &sh,
+ const double tolerance=1e-7);
+
+ /**
+ * Clone the current Manifold.
+ */
+ virtual std::unique_ptr<Manifold<dim,spacedim> > clone() const override;
+
+ /**
+ * Given a point on real space, find its arclength parameter. Throws an
+ * error in debug mode, if the point is not on the TopoDS_Edge given at
+ * construction time.
+ */
+ virtual Point<1>
+ pull_back(const Point<spacedim> &space_point) const override;
+
+ /**
+ * Given an arclength parameter, find its image in real space.
+ */
+ virtual Point<spacedim>
+ push_forward(const Point<1> &chart_point) const override;
+
+ protected:
+ /**
+ * The actual shape used to build this object.
+ */
+ const TopoDS_Shape sh;
+
+ /**
+ * A Curve adaptor. This is the one which is used in the computations, and
+ * it points to the right one above.
+ */
+ Handle_Adaptor3d_HCurve curve;
+
+ /**
+ * Relative tolerance used in all internal computations.
+ */
+ const double tolerance;
+
+ /**
+ * The total length of the curve. This is also used as a period if the
+ * edge is periodic.
+ */
+ const double length;
+ };
+
+ /**
+ * Manifold description for the face of a CAD imported usign OpenCASCADE.
+ *
+ * @ingroup manifold
+ *
+ * @author Andrea Mola, Mauro Bardelloni, 2016
+ */
+ template <int dim, int spacedim>
+ class NURBSPatchManifold : public ChartManifold<dim, spacedim, 2>
+ {
+ public:
+ /**
+ * The constructor takes an OpenCASCADE TopoDS_Face @p face and an optional
+ * @p tolerance. This class uses the interval OpenCASCADE variables @var u,
+ * @var v to descrive the manifold.
+ */
+ NURBSPatchManifold(const TopoDS_Face &face, const double tolerance = 1e-7);
+
+ /**
+ * Clone the current Manifold.
+ */
+ virtual std::unique_ptr<Manifold<dim,spacedim> > clone() const override;
+
+ /**
+ * Pull back the given point from the Euclidean space. Will return the uv
+ * coordinates associated with the point @p space_point.
+ */
+ virtual Point<2>
+ pull_back(const Point<spacedim> &space_point) const override;
+
+ /**
+ * Given a @p chart_point in the uv coordinate system, this method returns the
+ * Euclidean coordinates associated.
+ */
+ virtual Point<spacedim>
+ push_forward(const Point<2> &chart_point) const override;
+
+ /**
+ * Given a point in the spacedim dimensional Euclidean space, this
+ * method returns the derivatives of the function $F$ that maps from
+ * the uv coordinate system to the Euclidean coordinate
+ * system. In other words, it is a matrix of size
+ * $\text{spacedim}\times\text{chartdim}$.
+ *
+ * This function is used in the computations required by the
+ * get_tangent_vector() function.
+ *
+ * Refer to the general documentation of this class for more information.
+ */
+ virtual
+ DerivativeForm<1,2,spacedim>
+ push_forward_gradient(const Point<2> &chart_point) const override;
+
+ protected:
+ /**
+ * Return a tuple representing the minimum and maximum values of u
+ * and v. Precisely, it returns (u_min, u_max, v_min, v_max)
+ */
+ std::tuple<double, double, double, double>
+ get_uv_bounds() const;
+
+ /**
+ * An OpenCASCADE TopoDS_Face @p face given by the CAD.
+ */
+ TopoDS_Face face;
+
+ /**
+ * Tolerance used by OpenCASCADE to identify points in each
+ * operation.
+ */
+ double tolerance;
+ };
+
+}
+
+/*@}*/
+
+DEAL_II_NAMESPACE_CLOSE
+
+
+#endif // DEAL_II_WITH_OPENCASCADE
+#endif // dealii_occ_manifold_lib_h
SET(_src
utilities.cc
boundary_lib.cc
+ manifold_lib.cc
)
SET(_inst
boundary_lib.inst.in
+ manifold_lib.inst.in
utilities.inst.in
)
#ifdef DEAL_II_WITH_OPENCASCADE
-
-#include <GCPnts_AbscissaPoint.hxx>
-#include <BRepAdaptor_Curve.hxx>
-#include <BRepAdaptor_CompCurve.hxx>
-#include <BRepAdaptor_HCurve.hxx>
-#include <BRepAdaptor_HCompCurve.hxx>
-#include <GCPnts_AbscissaPoint.hxx>
-#include <ShapeAnalysis_Curve.hxx>
-#include <BRep_Tool.hxx>
-#include <BRepTools.hxx>
-#include <ShapeAnalysis_Surface.hxx>
-#include <TopoDS.hxx>
-
-#include <Standard_Version.hxx>
-#if (OCC_VERSION_MAJOR < 7)
-# include <Handle_Adaptor3d_HCurve.hxx>
-#endif
-
-
DEAL_II_NAMESPACE_OPEN
-
namespace OpenCASCADE
{
-
-
- namespace
- {
- /**
- * Return a Geometrical curve representation for the given
- * TopoDS_Shape. This function will fail when the given shape is
- * not of topological dimension one.
- */
- Handle_Adaptor3d_HCurve curve_adaptor(const TopoDS_Shape &shape)
- {
- Assert( (shape.ShapeType() == TopAbs_WIRE) ||
- (shape.ShapeType() == TopAbs_EDGE),
- ExcUnsupportedShape());
- if (shape.ShapeType() == TopAbs_WIRE)
- return Handle(BRepAdaptor_HCompCurve)(new BRepAdaptor_HCompCurve(TopoDS::Wire(shape)));
- else if (shape.ShapeType() == TopAbs_EDGE)
- return Handle(BRepAdaptor_HCurve)(new BRepAdaptor_HCurve(TopoDS::Edge(shape)));
-
- Assert(false, ExcInternalError());
- return Handle(BRepAdaptor_HCurve)(new BRepAdaptor_HCurve());
- }
-
-
-
-// Helper internal functions.
- double shape_length(const TopoDS_Shape &sh)
- {
- Handle_Adaptor3d_HCurve adapt = curve_adaptor(sh);
- return GCPnts_AbscissaPoint::Length(adapt->GetCurve());
- }
- }
-
- /*============================== NormalProjectionBoundary ==============================*/
template <int dim, int spacedim>
- NormalProjectionBoundary<dim,spacedim>::NormalProjectionBoundary(const TopoDS_Shape &sh,
- const double tolerance) :
- sh(sh),
- tolerance(tolerance)
- {
- Assert(spacedim == 3, ExcNotImplemented());
- }
-
-
-
- template<int dim, int spacedim>
std::unique_ptr<Manifold<dim, spacedim> >
NormalProjectionBoundary<dim,spacedim>::clone() const
{
- return std::unique_ptr<Manifold<dim,spacedim> >(new NormalProjectionBoundary(sh, tolerance));
+ return std::unique_ptr<Manifold<dim,spacedim> >(
+ new NormalProjectionBoundary(this->sh, this->tolerance));
}
template <int dim, int spacedim>
- Point<spacedim> NormalProjectionBoundary<dim,spacedim>::
- project_to_manifold (const ArrayView<const Point<spacedim>> &surrounding_points,
- const Point<spacedim> &candidate) const
- {
- (void)surrounding_points;
-#ifdef DEBUG
- for (unsigned int i=0; i<surrounding_points.size(); ++i)
- Assert(closest_point(sh, surrounding_points[i], tolerance)
- .distance(surrounding_points[i]) <
- std::max(tolerance*surrounding_points[i].norm(), tolerance),
- ExcPointNotOnManifold<spacedim>(surrounding_points[i]));
-#endif
- return closest_point(sh, candidate,tolerance);
- }
-
-
- /*============================== DirectionalProjectionBoundary ==============================*/
- template <int dim, int spacedim>
- DirectionalProjectionBoundary<dim,spacedim>::DirectionalProjectionBoundary(const TopoDS_Shape &sh,
- const Tensor<1,spacedim> &direction,
- const double tolerance) :
- sh(sh),
- direction(direction),
- tolerance(tolerance)
- {
- Assert(spacedim == 3, ExcNotImplemented());
- }
-
-
-
- template<int dim, int spacedim>
std::unique_ptr<Manifold<dim, spacedim> >
DirectionalProjectionBoundary<dim,spacedim>::clone() const
{
- return std::unique_ptr<Manifold<dim,spacedim> >
- (new DirectionalProjectionBoundary(sh, direction, tolerance));
- }
-
-
-
- template <int dim, int spacedim>
- Point<spacedim> DirectionalProjectionBoundary<dim,spacedim>::
- project_to_manifold (const ArrayView<const Point<spacedim>> &surrounding_points,
- const Point<spacedim> &candidate) const
- {
- (void)surrounding_points;
-#ifdef DEBUG
- for (unsigned int i=0; i<surrounding_points.size(); ++i)
- Assert(closest_point(sh, surrounding_points[i],tolerance)
- .distance(surrounding_points[i]) <
- std::max(tolerance*surrounding_points[i].norm(), tolerance),
- ExcPointNotOnManifold<spacedim>(surrounding_points[i]));
-#endif
- return line_intersection(sh, candidate, direction, tolerance);
+ return std::unique_ptr<Manifold<dim,spacedim> >(
+ new DirectionalProjectionBoundary(this->sh, this->direction, this->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());
- Assert(std::get<0>(count_elements(sh)) > 0,
- ExcMessage("NormalToMeshProjectionBoundary needs a shape containing faces to operate."));
- }
-
- template<int dim, int spacedim>
std::unique_ptr<Manifold<dim, spacedim> >
NormalToMeshProjectionBoundary<dim,spacedim>::clone() const
{
- return std::unique_ptr<Manifold<dim, spacedim> >
- (new NormalToMeshProjectionBoundary<dim,spacedim>(sh,tolerance));
- }
-
-
- template <int dim, int spacedim>
- Point<spacedim> NormalToMeshProjectionBoundary<dim,spacedim>::
- project_to_manifold (const ArrayView<const Point<spacedim>> &surrounding_points,
- const Point<spacedim> &candidate) const
- {
- TopoDS_Shape out_shape;
- Tensor<1,3> average_normal;
-#ifdef DEBUG
- for (unsigned int i=0; i<surrounding_points.size(); ++i)
- {
- Assert(closest_point(sh, surrounding_points[i], tolerance)
- .distance(surrounding_points[i]) <
- std::max(tolerance*surrounding_points[i].norm(), tolerance),
- ExcPointNotOnManifold<spacedim>(surrounding_points[i]));
- }
-#endif
-
- switch (surrounding_points.size())
- {
- case 2:
- {
- for (unsigned int i=0; i<surrounding_points.size(); ++i)
- {
- std::tuple<Point<3>, Tensor<1,3>, double, double>
- p_and_diff_forms =
- closest_point_and_differential_forms(sh,
- surrounding_points[i],
- tolerance);
- average_normal += std::get<1>(p_and_diff_forms);
- }
-
- 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."));
-
- Tensor<1,3> T = surrounding_points[0]-surrounding_points[1];
- T /= T.norm();
- average_normal = average_normal-(average_normal*T)*T;
- average_normal /= average_normal.norm();
- break;
- }
- case 4:
- {
- Tensor<1,3> u = surrounding_points[1]-surrounding_points[0];
- Tensor<1,3> v = surrounding_points[2]-surrounding_points[0];
- const double n1_coords[3] = {u[1] *v[2]-u[2] *v[1],u[2] *v[0]-u[0] *v[2],u[0] *v[1]-u[1] *v[0]};
- Tensor<1,3> n1(n1_coords);
- n1 = n1/n1.norm();
- u = surrounding_points[2]-surrounding_points[3];
- v = surrounding_points[1]-surrounding_points[3];
- const double n2_coords[3] = {u[1] *v[2]-u[2] *v[1],u[2] *v[0]-u[0] *v[2],u[0] *v[1]-u[1] *v[0]};
- Tensor<1,3> n2(n2_coords);
- n2 = n2/n2.norm();
-
- average_normal = (n1+n2)/2.0;
-
- Assert(average_normal.norm() > tolerance,
- 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;
- }
- case 8:
- {
- Tensor<1,3> u = surrounding_points[1]-surrounding_points[0];
- Tensor<1,3> v = surrounding_points[2]-surrounding_points[0];
- const double n1_coords[3] = {u[1] *v[2]-u[2] *v[1],u[2] *v[0]-u[0] *v[2],u[0] *v[1]-u[1] *v[0]};
- Tensor<1,3> n1(n1_coords);
- n1 = n1/n1.norm();
- u = surrounding_points[2]-surrounding_points[3];
- v = surrounding_points[1]-surrounding_points[3];
- const double n2_coords[3] = {u[1] *v[2]-u[2] *v[1],u[2] *v[0]-u[0] *v[2],u[0] *v[1]-u[1] *v[0]};
- Tensor<1,3> n2(n2_coords);
- n2 = n2/n2.norm();
- u = surrounding_points[4]-surrounding_points[7];
- v = surrounding_points[6]-surrounding_points[7];
- const double n3_coords[3] = {u[1] *v[2]-u[2] *v[1],u[2] *v[0]-u[0] *v[2],u[0] *v[1]-u[1] *v[0]};
- Tensor<1,3> n3(n3_coords);
- n3 = n3/n3.norm();
- u = surrounding_points[6]-surrounding_points[7];
- v = surrounding_points[5]-surrounding_points[7];
- const double n4_coords[3] = {u[1] *v[2]-u[2] *v[1],u[2] *v[0]-u[0] *v[2],u[0] *v[1]-u[1] *v[0]};
- Tensor<1,3> n4(n4_coords);
- n4 = n4/n4.norm();
-
- average_normal = (n1+n2+n3+n4)/4.0;
-
- Assert(average_normal.norm() > tolerance,
- 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;
- }
- default:
- {
- AssertThrow(false, ExcNotImplemented());
- break;
- }
- }
-
- return line_intersection(sh, candidate, average_normal, tolerance);
- }
-
-
- /*============================== ArclengthProjectionLineManifold ==============================*/
- template <int dim, int spacedim>
- ArclengthProjectionLineManifold<dim,spacedim>::ArclengthProjectionLineManifold(const TopoDS_Shape &sh,
- const double tolerance):
-
- ChartManifold<dim,spacedim,1>(sh.Closed() ?
- Point<1>(shape_length(sh)) :
- Point<1>()),
- sh(sh),
- curve(curve_adaptor(sh)),
- tolerance(tolerance),
- length(shape_length(sh))
- {
- Assert(spacedim >= 2, ExcImpossibleInDimSpacedim(dim, spacedim));
+ return std::unique_ptr<Manifold<dim,spacedim> >(
+ new NormalToMeshProjectionBoundary(this->sh, this->tolerance));
}
-
-
-
- template<int dim, int spacedim>
- std::unique_ptr<Manifold<dim, spacedim> >
- ArclengthProjectionLineManifold<dim, spacedim>::clone() const
- {
- return std::unique_ptr<Manifold<dim, spacedim> >
- (new ArclengthProjectionLineManifold(sh,tolerance));
- }
-
-
-
- template <int dim, int spacedim>
- Point<1>
- ArclengthProjectionLineManifold<dim,spacedim>::pull_back(const Point<spacedim> &space_point) const
- {
- double t (0.0);
- ShapeAnalysis_Curve curve_analysis;
- gp_Pnt proj;
- const double dist = curve_analysis.Project(curve->GetCurve(), point(space_point), tolerance, proj, t, true);
- Assert(dist < tolerance*length, ExcPointNotOnManifold<spacedim>(space_point));
- (void)dist; // Silence compiler warning in Release mode.
- return Point<1>(GCPnts_AbscissaPoint::Length(curve->GetCurve(),curve->GetCurve().FirstParameter(),t));
- }
-
-
-
- template <int dim, int spacedim>
- Point<spacedim>
- ArclengthProjectionLineManifold<dim,spacedim>::push_forward(const Point<1> &chart_point) const
- {
- GCPnts_AbscissaPoint AP(curve->GetCurve(), chart_point[0], curve->GetCurve().FirstParameter());
- gp_Pnt P = curve->GetCurve().Value(AP.Parameter());
- return point<spacedim>(P);
- }
-
- template <int dim, int spacedim>
- NURBSPatchManifold<dim, spacedim>::
- NURBSPatchManifold( const TopoDS_Face &face,
- const double tolerance)
- :
- face(face),
- tolerance(tolerance)
- {}
-
-
-
- template<int dim, int spacedim>
- std::unique_ptr<Manifold<dim, spacedim> >
- NURBSPatchManifold<dim,spacedim>::clone() const
- {
- return std::unique_ptr<Manifold<dim, spacedim> >
- (new NURBSPatchManifold<dim,spacedim>(face,tolerance));
- }
-
-
-
- template <int dim, int spacedim> Point<2>
- NURBSPatchManifold<dim, spacedim>::
- pull_back(const Point<spacedim> &space_point) const
- {
- Handle(Geom_Surface) SurfToProj = BRep_Tool::Surface(face);
-
- ShapeAnalysis_Surface projector(SurfToProj);
- gp_Pnt2d proj_params = projector.ValueOfUV(point(space_point), tolerance);
-
- double u = proj_params.X();
- double v = proj_params.Y();
-
- return Point<2>(u,v);
- }
-
- template <int dim, int spacedim>
- Point<spacedim>
- NURBSPatchManifold<dim, spacedim>::
- push_forward(const Point<2> &chart_point) const
- {
- return ::dealii::OpenCASCADE::push_forward<spacedim>(face, chart_point[0], chart_point[1]);
- }
-
- template <int dim, int spacedim>
- DerivativeForm<1,2,spacedim>
- NURBSPatchManifold<dim,spacedim>::
- push_forward_gradient(const Point<2> &chart_point) const
- {
- DerivativeForm<1,2,spacedim> DX;
- Handle(Geom_Surface) surf = BRep_Tool::Surface(face);
-
- gp_Pnt q;
- gp_Vec Du, Dv;
- surf->D1(chart_point[0],chart_point[1], q, Du, Dv);
-
- DX[0][0] = Du.X();
- DX[1][0] = Du.Y();
- if (spacedim>2)
- DX[2][0] = Du.Z();
- else
- Assert(std::abs(Du.Z()) < tolerance,
- ExcMessage("Expecting derivative along Z to be zero! Bailing out."));
- DX[0][1] = Dv.X();
- DX[1][1] = Dv.Y();
- if (spacedim>2)
- DX[2][1] = Dv.Z();
- else
- Assert(std::abs(Dv.Z()) < tolerance,
- ExcMessage("Expecting derivative along Z to be zero! Bailing out."));
- return DX;
- }
-
- template <int dim, int spacedim>
- std::tuple<double, double, double, double>
- NURBSPatchManifold<dim,spacedim>::
- get_uv_bounds() const
- {
- Standard_Real umin, umax, vmin, vmax;
- BRepTools::UVBounds(face, umin, umax, vmin, vmax);
- return std::make_tuple(umin, umax, vmin, vmax);
- }
-
// Explicit instantiations
#include "boundary_lib.inst"
-
-} // end namespace OpenCASCADE
+}
DEAL_II_NAMESPACE_CLOSE
template class NormalProjectionBoundary<deal_II_dimension, 3>;
template class DirectionalProjectionBoundary<deal_II_dimension, 3>;
template class NormalToMeshProjectionBoundary<deal_II_dimension, 3>;
- template class ArclengthProjectionLineManifold<deal_II_dimension, 3>;
- template class NURBSPatchManifold<deal_II_dimension, 3>;
#if deal_II_dimension <= 2
template class DirectionalProjectionBoundary<deal_II_dimension, 2>;
- template class ArclengthProjectionLineManifold<deal_II_dimension, 2>;
- template class NURBSPatchManifold<deal_II_dimension, 2>;
#endif
}
-
-
-
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2014 - 2018 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.
+//
+// ---------------------------------------------------------------------
+
+
+#include <deal.II/opencascade/manifold_lib.h>
+
+#ifdef DEAL_II_WITH_OPENCASCADE
+
+
+#include <GCPnts_AbscissaPoint.hxx>
+#include <BRepAdaptor_Curve.hxx>
+#include <BRepAdaptor_CompCurve.hxx>
+#include <BRepAdaptor_HCurve.hxx>
+#include <BRepAdaptor_HCompCurve.hxx>
+#include <GCPnts_AbscissaPoint.hxx>
+#include <ShapeAnalysis_Curve.hxx>
+#include <BRep_Tool.hxx>
+#include <BRepTools.hxx>
+#include <ShapeAnalysis_Surface.hxx>
+#include <TopoDS.hxx>
+
+#include <Standard_Version.hxx>
+#if (OCC_VERSION_MAJOR < 7)
+# include <Handle_Adaptor3d_HCurve.hxx>
+#endif
+
+
+DEAL_II_NAMESPACE_OPEN
+
+
+namespace OpenCASCADE
+{
+
+
+ namespace
+ {
+ /**
+ * Return a Geometrical curve representation for the given
+ * TopoDS_Shape. This function will fail when the given shape is
+ * not of topological dimension one.
+ */
+ Handle_Adaptor3d_HCurve curve_adaptor(const TopoDS_Shape &shape)
+ {
+ Assert( (shape.ShapeType() == TopAbs_WIRE) ||
+ (shape.ShapeType() == TopAbs_EDGE),
+ ExcUnsupportedShape());
+ if (shape.ShapeType() == TopAbs_WIRE)
+ return Handle(BRepAdaptor_HCompCurve)(new BRepAdaptor_HCompCurve(TopoDS::Wire(shape)));
+ else if (shape.ShapeType() == TopAbs_EDGE)
+ return Handle(BRepAdaptor_HCurve)(new BRepAdaptor_HCurve(TopoDS::Edge(shape)));
+
+ Assert(false, ExcInternalError());
+ return Handle(BRepAdaptor_HCurve)(new BRepAdaptor_HCurve());
+ }
+
+
+
+// Helper internal functions.
+ double shape_length(const TopoDS_Shape &sh)
+ {
+ Handle_Adaptor3d_HCurve adapt = curve_adaptor(sh);
+ return GCPnts_AbscissaPoint::Length(adapt->GetCurve());
+ }
+ }
+
+ /*============================== NormalProjectionManifold ==============================*/
+ template <int dim, int spacedim>
+ NormalProjectionManifold<dim,spacedim>::NormalProjectionManifold(const TopoDS_Shape &sh,
+ const double tolerance) :
+ sh(sh),
+ tolerance(tolerance)
+ {
+ Assert(spacedim == 3, ExcNotImplemented());
+ }
+
+
+
+ template<int dim, int spacedim>
+ std::unique_ptr<Manifold<dim, spacedim> >
+ NormalProjectionManifold<dim,spacedim>::clone() const
+ {
+ return std::unique_ptr<Manifold<dim,spacedim> >(new NormalProjectionManifold(sh, tolerance));
+ }
+
+
+
+ template <int dim, int spacedim>
+ Point<spacedim> NormalProjectionManifold<dim,spacedim>::
+ project_to_manifold (const ArrayView<const Point<spacedim>> &surrounding_points,
+ const Point<spacedim> &candidate) const
+ {
+ (void)surrounding_points;
+#ifdef DEBUG
+ for (unsigned int i=0; i<surrounding_points.size(); ++i)
+ Assert(closest_point(sh, surrounding_points[i], tolerance)
+ .distance(surrounding_points[i]) <
+ std::max(tolerance*surrounding_points[i].norm(), tolerance),
+ ExcPointNotOnManifold<spacedim>(surrounding_points[i]));
+#endif
+ return closest_point(sh, candidate,tolerance);
+ }
+
+
+ /*============================== DirectionalProjectionManifold ==============================*/
+ template <int dim, int spacedim>
+ DirectionalProjectionManifold<dim,spacedim>::DirectionalProjectionManifold(const TopoDS_Shape &sh,
+ const Tensor<1,spacedim> &direction,
+ const double tolerance) :
+ sh(sh),
+ direction(direction),
+ tolerance(tolerance)
+ {
+ Assert(spacedim == 3, ExcNotImplemented());
+ }
+
+
+
+ template<int dim, int spacedim>
+ std::unique_ptr<Manifold<dim, spacedim> >
+ DirectionalProjectionManifold<dim,spacedim>::clone() const
+ {
+ return std::unique_ptr<Manifold<dim,spacedim> >
+ (new DirectionalProjectionManifold(sh, direction, tolerance));
+ }
+
+
+
+ template <int dim, int spacedim>
+ Point<spacedim> DirectionalProjectionManifold<dim,spacedim>::
+ project_to_manifold (const ArrayView<const Point<spacedim>> &surrounding_points,
+ const Point<spacedim> &candidate) const
+ {
+ (void)surrounding_points;
+#ifdef DEBUG
+ for (unsigned int i=0; i<surrounding_points.size(); ++i)
+ Assert(closest_point(sh, surrounding_points[i],tolerance)
+ .distance(surrounding_points[i]) <
+ std::max(tolerance*surrounding_points[i].norm(), tolerance),
+ ExcPointNotOnManifold<spacedim>(surrounding_points[i]));
+#endif
+ return line_intersection(sh, candidate, direction, tolerance);
+ }
+
+
+
+ /*============================== NormalToMeshProjectionManifold ==============================*/
+ template <int dim, int spacedim>
+ NormalToMeshProjectionManifold<dim,spacedim>::NormalToMeshProjectionManifold(const TopoDS_Shape &sh,
+ const double tolerance) :
+ sh(sh),
+ tolerance(tolerance)
+ {
+ Assert(spacedim == 3, ExcNotImplemented());
+ Assert(std::get<0>(count_elements(sh)) > 0,
+ ExcMessage("NormalToMeshProjectionManifold needs a shape containing faces to operate."));
+ }
+
+ template<int dim, int spacedim>
+ std::unique_ptr<Manifold<dim, spacedim> >
+ NormalToMeshProjectionManifold<dim,spacedim>::clone() const
+ {
+ return std::unique_ptr<Manifold<dim, spacedim> >
+ (new NormalToMeshProjectionManifold<dim,spacedim>(sh,tolerance));
+ }
+
+
+ template <int dim, int spacedim>
+ Point<spacedim> NormalToMeshProjectionManifold<dim,spacedim>::
+ project_to_manifold (const ArrayView<const Point<spacedim>> &surrounding_points,
+ const Point<spacedim> &candidate) const
+ {
+ TopoDS_Shape out_shape;
+ Tensor<1,3> average_normal;
+#ifdef DEBUG
+ for (unsigned int i=0; i<surrounding_points.size(); ++i)
+ {
+ Assert(closest_point(sh, surrounding_points[i], tolerance)
+ .distance(surrounding_points[i]) <
+ std::max(tolerance*surrounding_points[i].norm(), tolerance),
+ ExcPointNotOnManifold<spacedim>(surrounding_points[i]));
+ }
+#endif
+
+ switch (surrounding_points.size())
+ {
+ case 2:
+ {
+ for (unsigned int i=0; i<surrounding_points.size(); ++i)
+ {
+ std::tuple<Point<3>, Tensor<1,3>, double, double>
+ p_and_diff_forms =
+ closest_point_and_differential_forms(sh,
+ surrounding_points[i],
+ tolerance);
+ average_normal += std::get<1>(p_and_diff_forms);
+ }
+
+ 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."));
+
+ Tensor<1,3> T = surrounding_points[0]-surrounding_points[1];
+ T /= T.norm();
+ average_normal = average_normal-(average_normal*T)*T;
+ average_normal /= average_normal.norm();
+ break;
+ }
+ case 4:
+ {
+ Tensor<1,3> u = surrounding_points[1]-surrounding_points[0];
+ Tensor<1,3> v = surrounding_points[2]-surrounding_points[0];
+ const double n1_coords[3] = {u[1] *v[2]-u[2] *v[1],u[2] *v[0]-u[0] *v[2],u[0] *v[1]-u[1] *v[0]};
+ Tensor<1,3> n1(n1_coords);
+ n1 = n1/n1.norm();
+ u = surrounding_points[2]-surrounding_points[3];
+ v = surrounding_points[1]-surrounding_points[3];
+ const double n2_coords[3] = {u[1] *v[2]-u[2] *v[1],u[2] *v[0]-u[0] *v[2],u[0] *v[1]-u[1] *v[0]};
+ Tensor<1,3> n2(n2_coords);
+ n2 = n2/n2.norm();
+
+ average_normal = (n1+n2)/2.0;
+
+ Assert(average_normal.norm() > tolerance,
+ 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;
+ }
+ case 8:
+ {
+ Tensor<1,3> u = surrounding_points[1]-surrounding_points[0];
+ Tensor<1,3> v = surrounding_points[2]-surrounding_points[0];
+ const double n1_coords[3] = {u[1] *v[2]-u[2] *v[1],u[2] *v[0]-u[0] *v[2],u[0] *v[1]-u[1] *v[0]};
+ Tensor<1,3> n1(n1_coords);
+ n1 = n1/n1.norm();
+ u = surrounding_points[2]-surrounding_points[3];
+ v = surrounding_points[1]-surrounding_points[3];
+ const double n2_coords[3] = {u[1] *v[2]-u[2] *v[1],u[2] *v[0]-u[0] *v[2],u[0] *v[1]-u[1] *v[0]};
+ Tensor<1,3> n2(n2_coords);
+ n2 = n2/n2.norm();
+ u = surrounding_points[4]-surrounding_points[7];
+ v = surrounding_points[6]-surrounding_points[7];
+ const double n3_coords[3] = {u[1] *v[2]-u[2] *v[1],u[2] *v[0]-u[0] *v[2],u[0] *v[1]-u[1] *v[0]};
+ Tensor<1,3> n3(n3_coords);
+ n3 = n3/n3.norm();
+ u = surrounding_points[6]-surrounding_points[7];
+ v = surrounding_points[5]-surrounding_points[7];
+ const double n4_coords[3] = {u[1] *v[2]-u[2] *v[1],u[2] *v[0]-u[0] *v[2],u[0] *v[1]-u[1] *v[0]};
+ Tensor<1,3> n4(n4_coords);
+ n4 = n4/n4.norm();
+
+ average_normal = (n1+n2+n3+n4)/4.0;
+
+ Assert(average_normal.norm() > tolerance,
+ 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;
+ }
+ default:
+ {
+ AssertThrow(false, ExcNotImplemented());
+ break;
+ }
+ }
+
+ return line_intersection(sh, candidate, average_normal, tolerance);
+ }
+
+
+ /*============================== ArclengthProjectionLineManifold ==============================*/
+ template <int dim, int spacedim>
+ ArclengthProjectionLineManifold<dim,spacedim>::ArclengthProjectionLineManifold(const TopoDS_Shape &sh,
+ const double tolerance):
+
+ ChartManifold<dim,spacedim,1>(sh.Closed() ?
+ Point<1>(shape_length(sh)) :
+ Point<1>()),
+ sh(sh),
+ curve(curve_adaptor(sh)),
+ tolerance(tolerance),
+ length(shape_length(sh))
+ {
+ Assert(spacedim >= 2, ExcImpossibleInDimSpacedim(dim, spacedim));
+ }
+
+
+
+ template<int dim, int spacedim>
+ std::unique_ptr<Manifold<dim, spacedim> >
+ ArclengthProjectionLineManifold<dim, spacedim>::clone() const
+ {
+ return std::unique_ptr<Manifold<dim, spacedim> >
+ (new ArclengthProjectionLineManifold(sh,tolerance));
+ }
+
+
+
+ template <int dim, int spacedim>
+ Point<1>
+ ArclengthProjectionLineManifold<dim,spacedim>::pull_back(const Point<spacedim> &space_point) const
+ {
+ double t (0.0);
+ ShapeAnalysis_Curve curve_analysis;
+ gp_Pnt proj;
+ const double dist = curve_analysis.Project(curve->GetCurve(), point(space_point), tolerance, proj, t, true);
+ Assert(dist < tolerance*length, ExcPointNotOnManifold<spacedim>(space_point));
+ (void)dist; // Silence compiler warning in Release mode.
+ return Point<1>(GCPnts_AbscissaPoint::Length(curve->GetCurve(),curve->GetCurve().FirstParameter(),t));
+ }
+
+
+
+ template <int dim, int spacedim>
+ Point<spacedim>
+ ArclengthProjectionLineManifold<dim,spacedim>::push_forward(const Point<1> &chart_point) const
+ {
+ GCPnts_AbscissaPoint AP(curve->GetCurve(), chart_point[0], curve->GetCurve().FirstParameter());
+ gp_Pnt P = curve->GetCurve().Value(AP.Parameter());
+ return point<spacedim>(P);
+ }
+
+ template <int dim, int spacedim>
+ NURBSPatchManifold<dim, spacedim>::
+ NURBSPatchManifold( const TopoDS_Face &face,
+ const double tolerance)
+ :
+ face(face),
+ tolerance(tolerance)
+ {}
+
+
+
+ template<int dim, int spacedim>
+ std::unique_ptr<Manifold<dim, spacedim> >
+ NURBSPatchManifold<dim,spacedim>::clone() const
+ {
+ return std::unique_ptr<Manifold<dim, spacedim> >
+ (new NURBSPatchManifold<dim,spacedim>(face,tolerance));
+ }
+
+
+
+ template <int dim, int spacedim> Point<2>
+ NURBSPatchManifold<dim, spacedim>::
+ pull_back(const Point<spacedim> &space_point) const
+ {
+ Handle(Geom_Surface) SurfToProj = BRep_Tool::Surface(face);
+
+ ShapeAnalysis_Surface projector(SurfToProj);
+ gp_Pnt2d proj_params = projector.ValueOfUV(point(space_point), tolerance);
+
+ double u = proj_params.X();
+ double v = proj_params.Y();
+
+ return Point<2>(u,v);
+ }
+
+ template <int dim, int spacedim>
+ Point<spacedim>
+ NURBSPatchManifold<dim, spacedim>::
+ push_forward(const Point<2> &chart_point) const
+ {
+ return ::dealii::OpenCASCADE::push_forward<spacedim>(face, chart_point[0], chart_point[1]);
+ }
+
+ template <int dim, int spacedim>
+ DerivativeForm<1,2,spacedim>
+ NURBSPatchManifold<dim,spacedim>::
+ push_forward_gradient(const Point<2> &chart_point) const
+ {
+ DerivativeForm<1,2,spacedim> DX;
+ Handle(Geom_Surface) surf = BRep_Tool::Surface(face);
+
+ gp_Pnt q;
+ gp_Vec Du, Dv;
+ surf->D1(chart_point[0],chart_point[1], q, Du, Dv);
+
+ DX[0][0] = Du.X();
+ DX[1][0] = Du.Y();
+ if (spacedim>2)
+ DX[2][0] = Du.Z();
+ else
+ Assert(std::abs(Du.Z()) < tolerance,
+ ExcMessage("Expecting derivative along Z to be zero! Bailing out."));
+ DX[0][1] = Dv.X();
+ DX[1][1] = Dv.Y();
+ if (spacedim>2)
+ DX[2][1] = Dv.Z();
+ else
+ Assert(std::abs(Dv.Z()) < tolerance,
+ ExcMessage("Expecting derivative along Z to be zero! Bailing out."));
+ return DX;
+ }
+
+ template <int dim, int spacedim>
+ std::tuple<double, double, double, double>
+ NURBSPatchManifold<dim,spacedim>::
+ get_uv_bounds() const
+ {
+ Standard_Real umin, umax, vmin, vmax;
+ BRepTools::UVBounds(face, umin, umax, vmin, vmax);
+ return std::make_tuple(umin, umax, vmin, vmax);
+ }
+
+// Explicit instantiations
+#include "manifold_lib.inst"
+} // end namespace OpenCASCADE
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2014 - 2017 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 NormalProjectionManifold<deal_II_dimension, 3>;
+ template class DirectionalProjectionManifold<deal_II_dimension, 3>;
+ template class NormalToMeshProjectionManifold<deal_II_dimension, 3>;
+ template class ArclengthProjectionLineManifold<deal_II_dimension, 3>;
+ template class NURBSPatchManifold<deal_II_dimension, 3>;
+#if deal_II_dimension <= 2
+ template class DirectionalProjectionManifold<deal_II_dimension, 2>;
+ template class ArclengthProjectionLineManifold<deal_II_dimension, 2>;
+ template class NURBSPatchManifold<deal_II_dimension, 2>;
+#endif
+}