]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Rename the opencascade Manifold classes.
authorDavid Wells <wellsd2@rpi.edu>
Sat, 12 May 2018 18:15:28 +0000 (14:15 -0400)
committerDavid Wells <drwells@vt.edu>
Wed, 16 May 2018 22:18:45 +0000 (18:18 -0400)
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.

doc/news/changes/incompatibilities/20180516Wells [new file with mode: 0644]
include/deal.II/opencascade/boundary_lib.h
include/deal.II/opencascade/manifold_lib.h [new file with mode: 0644]
source/opencascade/CMakeLists.txt
source/opencascade/boundary_lib.cc
source/opencascade/boundary_lib.inst.in
source/opencascade/manifold_lib.cc [new file with mode: 0644]
source/opencascade/manifold_lib.inst.in [new file with mode: 0644]

diff --git a/doc/news/changes/incompatibilities/20180516Wells b/doc/news/changes/incompatibilities/20180516Wells
new file mode 100644 (file)
index 0000000..c5f9b69
--- /dev/null
@@ -0,0 +1,7 @@
+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)
index 76b03fea1c58fa2c91c9d2122da3d53cee4ad077..6e591b4cecff0b0c55b3d2f5df6bb83bb5b3791e 100644 (file)
 
 #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
 
@@ -60,54 +52,24 @@ namespace OpenCASCADE
    * 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;
   };
 
   /**
@@ -130,55 +92,24 @@ namespace OpenCASCADE
    * 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;
   };
 
 
@@ -225,194 +156,25 @@ namespace OpenCASCADE
    * 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;
   };
-
 }
 
 /*@}*/
diff --git a/include/deal.II/opencascade/manifold_lib.h b/include/deal.II/opencascade/manifold_lib.h
new file mode 100644 (file)
index 0000000..f60ce39
--- /dev/null
@@ -0,0 +1,424 @@
+// ---------------------------------------------------------------------
+//
+// 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
index 0d9905db344c587567f2e7a8b5de2f50238b407c..7a78dbbfe1b9cb05e40c1ce9f9f9fb6cb34c0eb2 100644 (file)
@@ -18,10 +18,12 @@ INCLUDE_DIRECTORIES(BEFORE ${CMAKE_CURRENT_BINARY_DIR})
 SET(_src
   utilities.cc
   boundary_lib.cc
+  manifold_lib.cc
   )
 
 SET(_inst
   boundary_lib.inst.in
+  manifold_lib.inst.in
   utilities.inst.in
   )
 
index ab936f37d43dc01675d934889173ad21341677f2..3071532bea60adca59a16c5c3570b91f371521ed 100644 (file)
 
 #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
 
index 9e8f49ad097ebaaad7e6761d7bf678d06abd2b8a..46a478c34a8fa6a00548d855196a028975d4f7a9 100644 (file)
@@ -20,14 +20,7 @@ for (deal_II_dimension : DIMENSIONS)
     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
 }
-
-
-
diff --git a/source/opencascade/manifold_lib.cc b/source/opencascade/manifold_lib.cc
new file mode 100644 (file)
index 0000000..02462de
--- /dev/null
@@ -0,0 +1,425 @@
+// ---------------------------------------------------------------------
+//
+// 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
diff --git a/source/opencascade/manifold_lib.inst.in b/source/opencascade/manifold_lib.inst.in
new file mode 100644 (file)
index 0000000..c227a7f
--- /dev/null
@@ -0,0 +1,30 @@
+// ---------------------------------------------------------------------
+//
+// 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
+}

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.