]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Added ArclengthProjectionLineManifold.
authorLuca Heltai <luca.heltai@sissa.it>
Wed, 13 Aug 2014 15:22:44 +0000 (17:22 +0200)
committerLuca Heltai <luca.heltai@sissa.it>
Fri, 31 Oct 2014 10:06:49 +0000 (11:06 +0100)
Added first arclength test

Added axis boundary 1.

Added arclength boundary 2.

Fixed loop in projection for non degenerate edges.

Added axis boundary 2.

Added Degenerate Edge Checks and Assertions.

Added two functions to split opencascade objects and a test.

Added test for extraction of geometrical entities.

Added test for extraction of compound entities.

Added creation of a single cell triangulation, on top of a surface.

Fixed possibility to use also Wires.

Fixed indentation.

Removed creation of local files.

Fixed typos in the documentation.

Fixed typos in the documentation.

Addressed @tamiko and @guidokanschat comments on configs.

Renamed OPENCASCADE_LIBRARIES to _opencascade_libraries, and removed old TRILINOS cut and paste magic.

Removed unnecessary include files, removed commented lines, fixed header for test files.

Attempt to fix tolerance.

Fixed indentation issue.

Fixed WB comments.

Fixed indentation.

Missed one 'one'.

Replaced ImpossibleInDim with ExcNotImplemented

32 files changed:
cmake/config/Make.global_options.in
cmake/modules/FindOPENCASCADE.cmake
include/deal.II/grid/occ_boundary_lib.h
include/deal.II/grid/occ_utilities.h
source/grid/occ_boundary_lib.cc
source/grid/occ_boundary_lib.inst.in
source/grid/occ_utilities.cc
tests/opencascade/arclength_boundary_01.cc [new file with mode: 0644]
tests/opencascade/arclength_boundary_01.output [new file with mode: 0644]
tests/opencascade/arclength_boundary_02.cc [new file with mode: 0644]
tests/opencascade/arclength_boundary_02.output [new file with mode: 0644]
tests/opencascade/axis_boundary_01.cc [new file with mode: 0644]
tests/opencascade/axis_boundary_01.output [new file with mode: 0644]
tests/opencascade/axis_boundary_02.cc [new file with mode: 0644]
tests/opencascade/axis_boundary_02.output [new file with mode: 0644]
tests/opencascade/circle_normal_projection.cc
tests/opencascade/closest_point_to_shape.cc
tests/opencascade/create_tria_01.cc [new file with mode: 0644]
tests/opencascade/create_tria_01.output [new file with mode: 0644]
tests/opencascade/iges_create.cc
tests/opencascade/iges_describe.cc [new file with mode: 0644]
tests/opencascade/iges_describe.output [new file with mode: 0644]
tests/opencascade/iges_describe_02.cc [new file with mode: 0644]
tests/opencascade/iges_describe_02.output [new file with mode: 0644]
tests/opencascade/iges_files/side.iges [new file with mode: 0644]
tests/opencascade/iges_read.cc
tests/opencascade/iges_write.cc
tests/opencascade/interpolation_curve_boundary.cc
tests/opencascade/normal_projection_01.cc [new file with mode: 0644]
tests/opencascade/normal_projection_01.output [new file with mode: 0644]
tests/opencascade/plane_intersection_01.cc [new file with mode: 0644]
tests/opencascade/plane_intersection_01.output [new file with mode: 0644]

index d49b920f6a0a7e1c67149364d2596d46d9009a6e..6c6052bbae65f0bb52ff3d7c6f81a408f7ba8821 100644 (file)
@@ -50,19 +50,19 @@ DEAL_II_MINOR        = @DEAL_II_VERSION_MINOR@
 PERL                 = perl
 
 
-USE_CONTRIB_PETSC       = @MAKEFILE_PETSC@
-USE_CONTRIB_SLEPC       = @MAKEFILE_SLEPC@
-USE_CONTRIB_TRILINOS    = @MAKEFILE_TRILINOS@
-USE_CONTRIB_HDF5        = @MAKEFILE_HDF5@
-USE_CONTRIB_BLAS        = @MAKEFILE_BLAS@
-USE_CONTRIB_LAPACK      = @MAKEFILE_LAPACK@
-USE_CONTRIB_MUMPS       = @MAKEFILE_MUMPS@
-USE_CONTRIB_ARPACK      = @MAKEFILE_ARPACK@
-USE_CONTRIB_METIS       = @MAKEFILE_METIS@
-USE_CONTRIB_UMFPACK     = @MAKEFILE_UMFPACK@
-USE_CONTRIB_P4EST       = @MAKEFILE_P4EST@
-USE_CONTRIB_OPENCASCADE = @MAKEFILE_OPENCASCADE@
-DEAL_II_USE_MPI         = @MAKEFILE_MPI@
+USE_CONTRIB_PETSC                 = @MAKEFILE_PETSC@
+USE_CONTRIB_SLEPC                 = @MAKEFILE_SLEPC@
+USE_CONTRIB_TRILINOS              = @MAKEFILE_TRILINOS@
+USE_CONTRIB_HDF5                  = @MAKEFILE_HDF5@
+USE_CONTRIB_BLAS                  = @MAKEFILE_BLAS@
+USE_CONTRIB_LAPACK                = @MAKEFILE_LAPACK@
+USE_CONTRIB_MUMPS                 = @MAKEFILE_MUMPS@
+USE_CONTRIB_ARPACK                = @MAKEFILE_ARPACK@
+USE_CONTRIB_METIS                 = @MAKEFILE_METIS@
+USE_CONTRIB_UMFPACK               = @MAKEFILE_UMFPACK@
+USE_CONTRIB_P4EST                 = @MAKEFILE_P4EST@
+USE_CONTRIB_OPENCASCADE           = @MAKEFILE_OPENCASCADE@
+DEAL_II_USE_MPI                   = @MAKEFILE_MPI@
 
 OBJEXT     = o
 EXEEXT     =
index f6b03bce793038ddaceda8ce3e065b9f57d098d4..8e402ade6b23885cb25420cbf876fc9542147746 100644 (file)
@@ -48,7 +48,7 @@ IF(EXISTS ${OPENCASCADE_INCLUDE_DIR}/Standard_Version.hxx)
 ENDIF()
 
 # These seem to be pretty much the only required ones.
-SET(OPENCASCADE_LIBRARIES  
+SET(_opencascade_libraries
     TKFillet
     TKMesh
     TKernel
@@ -75,9 +75,9 @@ SET(OPENCASCADE_LIBRARIES
   )
 
 SET(_libraries "")
-FOREACH(_library ${OPENCASCADE_LIBRARIES})
-  LIST(APPEND _libraries ${_library})
-  DEAL_II_FIND_LIBRARY(${_library}
+FOREACH(_library ${_opencascade_libraries})
+  LIST(APPEND _libraries OPENCASCADE_${_library})
+  DEAL_II_FIND_LIBRARY(OPENCASCADE_${_library}
     NAMES ${_library}
     HINTS ${OPENCASCADE_DIR}
     PATH_SUFFIXES lib${LIB_SUFFIX} lib64 lib
@@ -99,6 +99,5 @@ DEAL_II_PACKAGE_HANDLE(OPENCASCADE
   USER_INCLUDE_DIRS
     REQUIRED OPENCASCADE_INCLUDE_DIR
   CLEAR
-    OPENCASCADE_LIBRARIES ${_libraries}
-    TRILINOS_SUPPORTS_CPP11 TRILINOS_HAS_C99_TR1_WORKAROUND
+    _opencascade_libraries ${_libraries}
   )
index eb53542afafb4660dc72334a9077156ff5906cc7..c52ee5a60a6ceadf2645f21c24398d3774430ebe 100644 (file)
 #include <deal.II/grid/occ_utilities.h>
 #include <deal.II/grid/tria_boundary.h>
 
+#include <BRepAdaptor_Curve.hxx>
+#include <Adaptor3d_Curve.hxx>
+
 DEAL_II_NAMESPACE_OPEN
 
-namespace OpenCASCADE 
+/**
+ * @addtogroup OpenCASCADE
+ * @{
+ */
+
+namespace OpenCASCADE
 {
   /**
-   * @addtogroup OpenCASCADE
-   * @{
-   * 
-   * A Boundary object based on OpenCASCADE TopoDS_Shape where new
-   * points are first computed using the FlatManifold class, and then
-   * projected in the normal direction using OpenCASCADE utilities.
+   * A Boundary object based on OpenCASCADE TopoDS_Shape where where
+   * new points are first computed by averaging the surrounding points
+   * in the same way as FlatManifold does, and then projecting 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
@@ -49,20 +55,29 @@ namespace OpenCASCADE
    * This could happen, for example, if you are trying to use a shape
    * of type TopoDS_Edge when projecting on a face. In this case, the
    * vertices of the face would be collapsed to the edge, and your
-   * surrounding points would not be on lying on the given shape,
-   * raising an exception.
-   * 
+   * 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 NormalProjectionBoundary : public Boundary<dim,spacedim> {
-    public:
-    NormalProjectionBoundary(const TopoDS_Shape sh, 
-                            const double tolerance=1e-7);
-    
+  class NormalProjectionBoundary : public Boundary<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 of arbitrary, i.e., a collection of
+     * shapes, faces, edges or a single face or edge.
+     */
+    NormalProjectionBoundary(const TopoDS_Shape &sh,
+                             const double tolerance=1e-7);
+
     /**
      * Perform the actual projection onto the manifold. This function,
-     * in debug mode, checks that each of the #surrounding_points is
+     * 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.
      *
@@ -71,16 +86,16 @@ namespace OpenCASCADE
      */
     virtual Point<spacedim>
     project_to_manifold (const std::vector<Point<spacedim> > &surrounding_points,
-                        const Point<spacedim> &candidate) const;
+                         const Point<spacedim> &candidate) const;
 
 
   private:
     /**
      * The topological shape which is used internally to project
-     * points. You can construct one such a shape by calling the
+     * 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;
 
     /**
@@ -91,9 +106,10 @@ namespace OpenCASCADE
 
   /**
    * A Boundary object based on OpenCASCADE TopoDS_Shape where new
-   * points are first computed using the FlatManifold class, and then
-   * projected along the direction given at construction time, using
-   * OpenCASCADE utilities.
+   * points are first computed by averaging the surrounding points in
+   * the same way as FlatManifold does, and then projecting in 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
@@ -103,40 +119,51 @@ namespace OpenCASCADE
    * 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 Boundary 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 appens.
+   *
    * @author Luca Heltai, Andrea Mola, 2011--2014.
    */
   template <int dim, int spacedim>
-  class AxisProjectionBoundary : public Boundary<dim,spacedim> {
-    public:
-    AxisProjectionBoundary(const TopoDS_Shape sh, 
-                          const Point<3> direction, 
-                          const double tolerance=1e-7);
-    
+  class DirectionalProjectionBoundary : public Boundary<dim,spacedim>
+  {
+  public:
+    /**
+     * Construct a Boundary object which will project points on the
+     * TopoDS_Shape @p sh, along the given @p direction.
+     */
+    DirectionalProjectionBoundary(const TopoDS_Shape &sh,
+                                  const Tensor<1,spacedim> &direction,
+                                  const double tolerance=1e-7);
+
     /**
      * Perform the actual projection onto the manifold. This function,
-     * in debug mode, checks that each of the #surrounding_points is
+     * 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
+     * The projected point is computed using OpenCASCADE directional
      * projection algorithms.
      */
     virtual Point<spacedim>
     project_to_manifold (const std::vector<Point<spacedim> > &surrounding_points,
-                        const Point<spacedim> &candidate) const;
+                         const Point<spacedim> &candidate) const;
 
   private:
     /**
      * The topological shape which is used internally to project
-     * points. You can construct one such a shape by calling the
+     * 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. 
+     * Direction used to project new points on the shape.
      */
     const Point<3> direction;
 
@@ -146,6 +173,70 @@ namespace OpenCASCADE
     const double tolerance;
   };
 
+  /**
+   * A Boundary 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);
+
+    /**
+     * 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;
+
+    /**
+     * Given an arclength parameter, find its image in real space.
+     */
+    virtual Point<spacedim>
+    push_forward(const Point<1> &chart_point) const;
+
+  private:
+    /**
+     * 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;
+  };
 }
 
 /*@}*/
index 074ca7ea1670afbf0f68c454db2966bf51307209..e45bf4efa7868c9586152dbef9a23cd1fd98f0e8 100644 (file)
 
 #ifdef DEAL_II_WITH_OPENCASCADE
 
+#include <deal.II/grid/tria.h>
+
 #include <string>
 #include <TopoDS_Shape.hxx>
+#include <TopoDS_Face.hxx>
 #include <TopoDS_Edge.hxx>
+#include <TopoDS_Vertex.hxx>
+#include <TopoDS_Compound.hxx>
+#include <TopoDS_CompSolid.hxx>
+#include <TopoDS_Solid.hxx>
+#include <TopoDS_Shell.hxx>
+#include <TopoDS_Wire.hxx>
 #include <IFSelect_ReturnStatus.hxx>
-#include <Geom_Plane.hxx>
-#include <Geom_Curve.hxx>
 #include <gp_Pnt.hxx>
 
 #include <deal.II/base/point.h>
 
+
 DEAL_II_NAMESPACE_OPEN
 
 /**
@@ -52,10 +60,9 @@ DEAL_II_NAMESPACE_OPEN
  * contains 3 fields: location, orientation and a myTShape handle (of
  * the TopoDS_TShape type). According to OpenCASCADE documentation,
  * myTShape and Location are used to share data between various shapes
- * and thus save huge amounts of memory. For example, an edge
- * belonging to two faces has equal Locations and myTShape fields but
- * different Orientations (Forward in context of one face and Reversed
- * in one of the other).
+ * to save memory. For example, an edge belonging to two faces has
+ * equal Locations and myTShape fields but different Orientations
+ * (Forward in context of one face and Reversed in one of the other).
  *
  * Valid shapes include collection of other shapes, solids, faces,
  * edges, vertices, etc.
@@ -85,38 +92,38 @@ DEAL_II_NAMESPACE_OPEN
  *
  * @author Luca Heltai, Andrea Mola, 2011--2014.
  */
-namespace OpenCASCADE 
+namespace OpenCASCADE
 {
   /**
-   * Count the subobjects of a shape. This function just outputs some
-   * information about the TopoDS_Shape passed as argument. It counts
-   * the number of faces, edges and vertices (the only topological
-   * entities associated with actual geometries) which are contained
-   * in the given shape.
+   * Count the subobjects of a shape. This function is useful to
+   * gather information about the TopoDS_Shape passed as argument. It
+   * counts the number of faces, edges and vertices (the only
+   * topological entities associated with actual geometries) which are
+   * contained in the given shape.
    */
   void count_elements(const TopoDS_Shape &shape,
-                     unsigned int &n_faces,
-                     unsigned int &n_edges,
-                     unsigned int &n_vertices);
-  
+                      unsigned int &n_faces,
+                      unsigned int &n_edges,
+                      unsigned int &n_vertices);
+
   /**
    * Read IGES files and translate their content into openCascade
    * topological entities. The option scale_factor is used to
    * compensate for different units being used in the IGES files and
    * in the target application. The standard unit for IGES files is
    * millimiters. The return object is a TopoDS_Shape which contains
-   * all objects from the file. 
-   */
-  TopoDS_Shape read_IGES(const std::string &filename, 
-                        const double scale_factor=1e-3);
-  
- /**
-   * Write the given topological shape into an IGES file.
+   * all objects from the file.
    */
-  void write_IGES(const TopoDS_Shape &shape, 
-                 const std::string &filename);
-                 
-  
+  TopoDS_Shape read_IGES(const std::string &filename,
+                         const double scale_factor=1e-3);
+
+  /**
+    * Write the given topological shape into an IGES file.
+    */
+  void write_IGES(const TopoDS_Shape &shape,
+                  const std::string &filename);
+
+
 
   /**
    * Perform the intersection of the given topological shape with the
@@ -125,13 +132,23 @@ namespace OpenCASCADE
    * thrown if the intersection produces an empty shape.
    */
   TopoDS_Shape  intersect_plane(const TopoDS_Shape &in_shape,
-                               const double c_x,
-                               const double c_y,
-                               const double c_z,
-                               const double c,
-                               const double tolerance=1e-7);
+                                const double c_x,
+                                const double c_y,
+                                const double c_z,
+                                const double c,
+                                const double tolerance=1e-7);
+
+  /**
+   * Try to join all edges contained in the given TopoDS_Shape into a
+   * single TopoDS_Edge, containing as few BSPlines as possible. If
+   * the input shape contains faces, they will be ignored by this
+   * function. If the contained edges cannot be joined into a single
+   * one, i.e., they form disconnected curves, an exception will be
+   * thrown.
+   */
+  TopoDS_Edge join_edges(const TopoDS_Shape &in_shape,
+                         const double tolerance=1e-7);
 
-  
   /**
    * Creates a 3D smooth BSpline curve passing through the points in
    * the assigned vector, and store it in the returned TopoDS_Shape
@@ -143,18 +160,53 @@ namespace OpenCASCADE
    *
    * This class is used to interpolate a BsplineCurve passing through
    * an array of points, with a C2 Continuity. If the optional
-   * parameter #closed is set to true, then the curve will be C2 at
+   * parameter @p closed is set to true, then the curve will be C2 at
    * all points execpt the first (where only C1 continuity will be
    * given), and it will be a closed curve.
    *
-   * The curve is garanteed to be at distance #tolerance from the
+   * The curve is garanteed to be at distance @p tolerance from the
    * input points. If the algorithm fails in generating such a curve,
    * an exception is thrown.
    */
   TopoDS_Edge interpolation_curve(std::vector<dealii::Point<3> >  &curve_points,
-                                 const dealii::Point<3> direction=dealii::Point<3>(), 
-                                 const bool closed=false,
-                                 const double tolerance=1e-7);
+                                  const dealii::Point<3> direction=dealii::Point<3>(),
+                                  const bool closed=false,
+                                  const double tolerance=1e-7);
+
+  /**
+   * Extract all subshapes from a TopoDS_Shape, and store the results
+   * into standard containers. If the shape does not contain a certain
+   * type of shape, the respective container will be empty.
+   */
+  void extract_geometrical_shapes(const TopoDS_Shape &shape,
+                                  std::vector<TopoDS_Face> &faces,
+                                  std::vector<TopoDS_Edge> &edges,
+                                  std::vector<TopoDS_Vertex> &vertices);
+
+  /**
+   * Create a triangulation from a single face. This class extract the
+   * first u and v parameter of the parametric surface making up this
+   * face, and creates a Triangulation<2,3> containing a single coarse
+   * cell reflecting this face. If the surface is not a trimmed
+   * surface, the vertices of this cell will coincide with the
+   * TopoDS_Vertex vertices of the original TopoDS_Face. This,
+   * however, is often not the case, and the user should be careful on
+   * how this mesh is used.
+   */
+  void create_triangulation(const TopoDS_Face &face,
+                            Triangulation<2,3> &tria);
+
+  /**
+   * Extract all compound shapes from a TopoDS_Shape, and store the
+   * results into standard containers. If the shape does not contain a
+   * certain type of compound, the respective container will be empty.
+   */
+  void extract_compound_shapes(const TopoDS_Shape &shape,
+                               std::vector<TopoDS_Compound> &compounds,
+                               std::vector<TopoDS_CompSolid> &compsolids,
+                               std::vector<TopoDS_Solid> &solids,
+                               std::vector<TopoDS_Shell> &shells,
+                               std::vector<TopoDS_Wire> &wires);
 
 
   /**
@@ -166,75 +218,86 @@ namespace OpenCASCADE
    * coordinate is filled with sensible information, and the v
    * coordinate is set to zero.
    */
-  Point<3> closest_point(const TopoDS_Shape in_shape, 
-                        const Point<3> origin,
-                        TopoDS_Shape &out_shape,
-                        double &u, 
-                        double &v, 
-                        const double tolerance=1e-7);
+  Point<3> closest_point(const TopoDS_Shape in_shape,
+                         const Point<3> origin,
+                         TopoDS_Shape &out_shape,
+                         double &u,
+                         double &v,
+                         const double tolerance=1e-7);
 
   /**
-   * Intersect a line passing through the given #origin point along
-   * #direction and the given topological shape. If there is more than
+   * Intersect a line passing through the given @p origin point along
+   * @p direction and the given topological shape. If there is more than
    * one intersection, it will return the closest one.
    *
-   * The optional #tolerance parameter is used to compute distances.
+   * The optional @p tolerance parameter is used to compute distances.
    */
-  Point<3> axis_intersection(const TopoDS_Shape in_shape, 
-                            const Point<3> origin, 
-                            const Point<3> direction, 
-                            const double tolerance=1e-7);
-  
+  Point<3> axis_intersection(const TopoDS_Shape in_shape,
+                             const Point<3> origin,
+                             const Point<3> direction,
+                             const double tolerance=1e-7);
+
 
   /**
    * Convert OpenCASCADE point into a Point<3>.
    */
-  inline Point<3> Pnt(const gp_Pnt &p);
+  Point<3> Pnt(const gp_Pnt &p);
 
 
   /**
    * Convert Point<3> into OpenCASCADE point.
    */
-  inline gp_Pnt Pnt(const Point<3> &p);
+  gp_Pnt Pnt(const Point<3> &p);
+
 
-  
   /**
    * Sort two points according to their scalar product with
    * direction. If the norm of the direction is zero, then use
    * lexycographical ordering. The optional parameter is used as a
    * relative tolerance when comparing objects.
    */
-  inline bool point_compare(const dealii::Point<3> &p1, const dealii::Point<3> &p2,
-                           const dealii::Point<3> direction=Point<3>(),
-                           const double tolerance=1e-10);
+  bool point_compare(const dealii::Point<3> &p1, const dealii::Point<3> &p2,
+                     const dealii::Point<3> direction=Point<3>(),
+                     const double tolerance=1e-10);
 
 
   /**
    * Exception thrown when the point specified as argument does not
-   * lie between #tolerance from the given TopoDS_Shape.
+   * lie between @p tolerance from the given TopoDS_Shape.
    */
   DeclException1 (ExcPointNotOnManifold,
-                 Point<3>,
-                 <<"The point [ "<<arg1<<" ] is not on the manifold.");
+                  Point<3>,
+                  <<"The point [ "<<arg1<<" ] is not on the manifold.");
 
   /**
    * Exception thrown when the point specified as argument cannot be
    * projected to the manifold.
-   */                        
-  DeclException1 (ExcProjectionFailed, 
-                 Point<3>,
-                 <<"Projection of point [ "<< arg1
-                 << " ] failed.");
+   */
+  DeclException1 (ExcProjectionFailed,
+                  Point<3>,
+                  <<"Projection of point [ "<< arg1
+                  << " ] failed.");
 
   /**
    * Thrown when internal OpenCASCADE utilities fail to return the OK
    * status.
    */
-  DeclException1 (ExcOCCError, 
-                 IFSelect_ReturnStatus, 
-                 <<"An OpenCASCADE routine failed with return status "
-                 <<arg1);
-} 
+  DeclException1 (ExcOCCError,
+                  IFSelect_ReturnStatus,
+                  <<"An OpenCASCADE routine failed with return status "
+                  <<arg1);
+
+  /**
+   * Trying to make curve operations on a degenerate edge.
+   */
+  DeclException0(ExcEdgeIsDegenerate);
+
+  /**
+   * Trying to make operations on the wrong type of shapes.
+   */
+  DeclException0(ExcUnsupportedShape);
+
+}
 /*@}*/
 
 DEAL_II_NAMESPACE_CLOSE
index 64724e91480c9f5b28c3800f86bd756dff2f8e05..fb95599835135c3f3b95907de699399e3b86631b 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 <TopoDS.hxx>
+#include <Adaptor3d_HCurve.hxx>
+#include <Handle_Adaptor3d_HCurve.hxx>
+
 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) :
+  NormalProjectionBoundary<dim,spacedim>::NormalProjectionBoundary(const TopoDS_Shape &sh,
+      const double tolerance) :
     sh(sh),
-    tolerance(tolerance) 
+    tolerance(tolerance)
   {
-    Assert(spacedim == 3, ExcImpossibleInDim(spacedim));
+    Assert(spacedim == 3, ExcNotImplemented());
   }
-  
-  
+
+
   template <int dim, int spacedim>
   Point<spacedim>  NormalProjectionBoundary<dim,spacedim>::
   project_to_manifold (const std::vector<Point<spacedim> > &surrounding_points,
-                      const Point<spacedim> &candidate) const {
+                       const Point<spacedim> &candidate) const
+  {
     TopoDS_Shape out_shape;
     double u=0, v=0;
-    for(unsigned int i=0; i<surrounding_points.size(); ++i)
+    for (unsigned int i=0; i<surrounding_points.size(); ++i)
       Assert(closest_point(sh, surrounding_points[i], out_shape, u, v)
-            .distance(surrounding_points[i]) < (surrounding_points[i].norm()>0 ? 
-                                                tolerance*surrounding_points[i].norm() :
-                                                tolerance),
-            // ExcPointNotOnManifold(surrounding_points[i]));
-            ExcMessage("Points are not on manifold"));
-    
+             .distance(surrounding_points[i]) <
+             std::max(tolerance*surrounding_points[i].norm(), tolerance),
+             ExcPointNotOnManifold(surrounding_points[i]));
+
     return closest_point(sh, candidate, out_shape, u, v);
   }
 
-  /*============================== AxisProjectionBoundary ==============================*/
+  /*============================== DirectionalProjectionBoundary ==============================*/
   template <int dim, int spacedim>
-  AxisProjectionBoundary<dim,spacedim>::AxisProjectionBoundary(const TopoDS_Shape sh, 
-                                                              const Point<3> direction,
-                                                              const double tolerance) :
+  DirectionalProjectionBoundary<dim,spacedim>::DirectionalProjectionBoundary(const TopoDS_Shape &sh,
+      const Tensor<1,spacedim> &direction,
+      const double tolerance) :
     sh(sh),
     direction(direction),
-    tolerance(tolerance) 
+    tolerance(tolerance)
   {
-    Assert(spacedim == 3, ExcImpossibleInDim(spacedim));
+    Assert(spacedim == 3, ExcNotImplemented());
   }
-  
-  
+
+
   template <int dim, int spacedim>
-  Point<spacedim>  AxisProjectionBoundary<dim,spacedim>::
+  Point<spacedim>  DirectionalProjectionBoundary<dim,spacedim>::
   project_to_manifold (const std::vector<Point<spacedim> > &surrounding_points,
-                      const Point<spacedim> &candidate) const {
+                       const Point<spacedim> &candidate) const
+  {
     TopoDS_Shape out_shape;
     double u=0, v=0;
-    for(unsigned int i=0; i<surrounding_points.size(); ++i)
+    for (unsigned int i=0; i<surrounding_points.size(); ++i)
       Assert(closest_point(sh, surrounding_points[i], out_shape, u, v)
-            .distance(surrounding_points[i]) < (surrounding_points[i].norm()>0 ? 
-                                                tolerance*surrounding_points[i].norm() :
-                                                tolerance),
-            ExcPointNotOnManifold(surrounding_points[i]));
-    
+             .distance(surrounding_points[i]) < (surrounding_points[i].norm()>0 ?
+                                                 tolerance*surrounding_points[i].norm() :
+                                                 tolerance),
+             ExcPointNotOnManifold(surrounding_points[i]));
+
     return axis_intersection(sh, candidate, direction, tolerance);
   }
 
-  
+  /*============================== 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>()),
+    curve(curve_adaptor(sh)),
+    tolerance(tolerance),
+    length(shape_length(sh))
+  {
+    Assert(spacedim == 3, ExcNotImplemented());
+  }
+
+
+  template <int dim, int spacedim>
+  Point<1>
+  ArclengthProjectionLineManifold<dim,spacedim>::pull_back(const Point<spacedim> &space_point) const
+  {
+    ShapeAnalysis_Curve curve_analysis;
+    gp_Pnt proj;
+    double t;
+    double dist = curve_analysis.Project(curve->GetCurve(), Pnt(space_point), tolerance, proj, t, true);
+    Assert(dist < tolerance*length, ExcPointNotOnManifold(space_point));
+    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 Pnt(P);
+  }
+
+
   // Explicit instantiations
-#include "occ_boundary_lib.inst"  
-  
+#include "occ_boundary_lib.inst"
+
 } // end namespace OpenCASCADE
 
 DEAL_II_NAMESPACE_CLOSE
index 3471c0aea80e0a5c1ed1b484f55491622cf14568..8194c9a707f35a089a0227646e8d8c4fbf510dda 100644 (file)
@@ -19,7 +19,8 @@
 for (deal_II_dimension : DIMENSIONS)
   {
        template class NormalProjectionBoundary<deal_II_dimension, 3>;
-       template class AxisProjectionBoundary<deal_II_dimension, 3>;
+       template class DirectionalProjectionBoundary<deal_II_dimension, 3>;
+       template class ArclengthProjectionLineManifold<deal_II_dimension, 3>; 
   }
 
 
index 84db8fbf3b5fc9896835deadbd8b686a6d954d66..238bba35ca499cfe38c2db80250e24b8defc14e8 100644 (file)
@@ -3,6 +3,10 @@
 #ifdef DEAL_II_WITH_OPENCASCADE
 
 #include <deal.II/base/point.h>
+#include <deal.II/base/utilities.h>
+#include <deal.II/base/exceptions.h>
+
+#include <boost/bind.hpp>
 
 #include <stdio.h>
 #include <stdlib.h>
 
 #include <TopoDS.hxx>
 #include <TopoDS_Shape.hxx>
-// #include <TopoDS_Shell.hxx>
 #include <TopoDS_Face.hxx>
 #include <TopoDS_Edge.hxx>
-#include <BRepTools.hxx>
+#include <TopExp_Explorer.hxx>
 
-// #include <XSControl_Reader.hxx>
-// #include <TopTools_SequenceOfShape.hxx>
 #include <Handle_Standard_Transient.hxx>
+
 #include <TColStd_SequenceOfTransient.hxx>
 #include <TColStd_HSequenceOfTransient.hxx>
-#include <TopExp_Explorer.hxx>
+#include <TColgp_HArray1OfPnt.hxx>
+
 #include <gp_Pnt.hxx>
 #include <gp_Lin.hxx>
 #include <gp_Vec.hxx>
 #include <GeomAPI_ProjectPointOnSurf.hxx>
 #include <GeomAPI_ProjectPointOnCurve.hxx>
 #include <IntCurvesFace_ShapeIntersector.hxx>
-// #include <Standard_Real.hxx>
-// #include <Standard_Integer.hxx>
-// #include <BRep_Tool.hxx>
-// #include <Geom_Surface.hxx>
-#include <Geom_Plane.hxx>
-// #include <Prs3d_ShapeTool.hxx>
-// #include <GeomAPI_IntSS.hxx>
-// #include <Bnd_Box.hxx>
-// #include <gp_Trsf.hxx>
-// #include <gp_Ax3.hxx>
-// #include <gp_Pln.hxx>
+
+#include <BRepTools.hxx>
+#include <BRepAdaptor_Curve.hxx>
+#include <BRepAdaptor_HCurve.hxx>
+#include <BRepAdaptor_HCompCurve.hxx>
+#include <BRepAdaptor_Surface.hxx>
+#include <BRep_Builder.hxx>
 #include <BRepBuilderAPI_Transform.hxx>
-#include <GeomConvert_CompCurveToBSplineCurve.hxx>
 #include <BRepBuilderAPI_MakeEdge.hxx>
-// #include <TColGeom_Array1OfCurve.hxx>
-#include <TColgp_HArray1OfPnt.hxx>
-// #include <Geom_Curve.hxx>
-#include <Geom_BoundedCurve.hxx>
-// #include <Geom_TrimmedCurve.hxx>
-// #include <Geom_BSplineCurve.hxx>
-#include <GeomAPI_Interpolate.hxx>
 #include <BRepAlgo_Section.hxx>
-// #include <GeomLib_Tool.hxx>
-// #include <TColGeom_Array2OfBezierSurface.hxx>
-// #include <ProjLib_ProjectOnPlane.hxx>
-// #include <Adaptor3d_HCurve.hxx>
-// #include <GeomAdaptor_HCurve.hxx>
-#include <ShapeAnalysis_Surface.hxx>
-// #include <GeomLProp_SLProps.hxx>
-// #include <BRepExtrema_DistShapeShape.hxx>
-// #include <BRepBuilderAPI_MakeVertex.hxx>
-// #include <GCPnts_AbscissaPoint.hxx>
-// #include <GeomLProp_CLProps.hxx>
 
-#include <deal.II/base/utilities.h>
-#include <deal.II/base/exceptions.h>
-
-#include <boost/bind.hpp>
-
-#include <BRep_Builder.hxx>
+#include <Geom_Plane.hxx>
+#include <Geom_BoundedCurve.hxx>
+#include <GeomAPI_Interpolate.hxx>
+#include <GeomConvert_CompCurveToBSplineCurve.hxx>
 
+#include <GCPnts_AbscissaPoint.hxx>
+#include <ShapeAnalysis_Surface.hxx>
 
 #include <vector>
 #include <algorithm>
@@ -82,17 +62,82 @@ DEAL_II_NAMESPACE_OPEN
 namespace OpenCASCADE
 {
   void count_elements(const TopoDS_Shape &shape,
-                     unsigned int &n_faces,
-                     unsigned int &n_edges,
-                     unsigned int &n_vertices) {
+                      unsigned int &n_faces,
+                      unsigned int &n_edges,
+                      unsigned int &n_vertices)
+  {
     TopExp_Explorer exp;
     n_faces=0, n_edges=0, n_vertices=0;
-    for(exp.Init(shape, TopAbs_FACE); exp.More(); exp.Next(), ++n_faces) {
-    }
-    for(exp.Init(shape, TopAbs_EDGE); exp.More(); exp.Next(), ++n_edges) {
-    }
-    for(exp.Init(shape, TopAbs_VERTEX); exp.More(); exp.Next(), ++n_vertices) {
-    }    
+    for (exp.Init(shape, TopAbs_FACE);
+         exp.More(); exp.Next(), ++n_faces)
+      {}
+    for (exp.Init(shape, TopAbs_EDGE);
+         exp.More(); exp.Next(), ++n_edges)
+      {}
+    for (exp.Init(shape, TopAbs_VERTEX);
+         exp.More(); exp.Next(), ++n_vertices)
+      {}
+  }
+
+  void extract_geometrical_shapes(const TopoDS_Shape &shape,
+                                  std::vector<TopoDS_Face> &faces,
+                                  std::vector<TopoDS_Edge> &edges,
+                                  std::vector<TopoDS_Vertex> &vertices)
+  {
+    faces.resize(0);
+    edges.resize(0);
+    vertices.resize(0);
+
+    TopExp_Explorer exp;
+    for (exp.Init(shape, TopAbs_FACE); exp.More(); exp.Next())
+      {
+        faces.push_back(TopoDS::Face(exp.Current()));
+      }
+    for (exp.Init(shape, TopAbs_EDGE); exp.More(); exp.Next())
+      {
+        edges.push_back(TopoDS::Edge(exp.Current()));
+      }
+    for (exp.Init(shape, TopAbs_VERTEX); exp.More(); exp.Next())
+      {
+        vertices.push_back(TopoDS::Vertex(exp.Current()));
+      }
+  }
+
+
+  void extract_compound_shapes(const TopoDS_Shape &shape,
+                               std::vector<TopoDS_Compound> &compounds,
+                               std::vector<TopoDS_CompSolid> &compsolids,
+                               std::vector<TopoDS_Solid> &solids,
+                               std::vector<TopoDS_Shell> &shells,
+                               std::vector<TopoDS_Wire> &wires)
+  {
+    compounds.resize(0);
+    compsolids.resize(0);
+    solids.resize(0);
+    shells.resize(0);
+    wires.resize(0);
+
+    TopExp_Explorer exp;
+    for (exp.Init(shape, TopAbs_COMPOUND); exp.More(); exp.Next())
+      {
+        compounds.push_back(TopoDS::Compound(exp.Current()));
+      }
+    for (exp.Init(shape, TopAbs_COMPSOLID); exp.More(); exp.Next())
+      {
+        compsolids.push_back(TopoDS::CompSolid(exp.Current()));
+      }
+    for (exp.Init(shape, TopAbs_SOLID); exp.More(); exp.Next())
+      {
+        solids.push_back(TopoDS::Solid(exp.Current()));
+      }
+    for (exp.Init(shape, TopAbs_SHELL); exp.More(); exp.Next())
+      {
+        shells.push_back(TopoDS::Shell(exp.Current()));
+      }
+    for (exp.Init(shape, TopAbs_WIRE); exp.More(); exp.Next())
+      {
+        wires.push_back(TopoDS::Wire(exp.Current()));
+      }
   }
 
   gp_Pnt Pnt(const Point<3> &p)
@@ -100,69 +145,69 @@ namespace OpenCASCADE
     return gp_Pnt(p(0), p(1), p(2));
   }
 
-  
+
   Point<3> Pnt(const gp_Pnt &p)
   {
     return Point<3>(p.X(), p.Y(), p.Z());
   }
-  
-  inline bool point_compare(const dealii::Point<3> &p1, const dealii::Point<3> &p2,
-                           const dealii::Point<3> direction,
-                           const double tolerance)
+
+  bool point_compare(const dealii::Point<3> &p1, const dealii::Point<3> &p2,
+                     const dealii::Point<3> direction,
+                     const double tolerance)
   {
     const double rel_tol=std::max(tolerance, std::max(p1.norm(), p2.norm())*tolerance);
-    if(direction.norm() > 0.0)
+    if (direction.norm() > 0.0)
       return (p1*direction < p2*direction-rel_tol);
-    else 
-      for(int d=2; d>=0; --d) 
-       if(p1[d] < p2[d]-rel_tol)
-         return true;
-       else if(p2[d] < p1[d]-rel_tol)
-         return false;
-         
+    else
+      for (int d=2; d>=0; --d)
+        if (p1[d] < p2[d]-rel_tol)
+          return true;
+        else if (p2[d] < p1[d]-rel_tol)
+          return false;
+
     // If we got here, for all d, none of the conditions above was
     // satisfied. The two points are equal up to tolerance
     return false;
   }
-  
 
-  TopoDS_Shape read_IGES(const std::string &filename, 
-                        const double scale_factor)
+
+  TopoDS_Shape read_IGES(const std::string &filename,
+                         const double scale_factor)
   {
     IGESControl_Reader reader;
     IFSelect_ReturnStatus stat;
-    stat = reader.ReadFile(filename.c_str()); 
-    Assert(stat == IFSelect_RetDone, 
-          ExcMessage("Error in reading file!"));
-          
+    stat = reader.ReadFile(filename.c_str());
+    Assert(stat == IFSelect_RetDone,
+           ExcMessage("Error in reading file!"));
+
     Standard_Boolean failsonly = Standard_False;
     IFSelect_PrintCount mode = IFSelect_ItemsByEntity;
     reader.PrintCheckLoad (failsonly, mode);
 
     Handle(TColStd_HSequenceOfTransient) myList = reader.GiveList("iges-faces");
-                                    //selects all IGES faces in the
-                                    //file and puts them into a list
-                                    //called MyList,
-    Standard_Integer // nIgesFaces = myList->Length(), 
-      nTransFaces = reader.TransferList(myList);
+    //selects all IGES faces in the
+    //file and puts them into a list
+    //called MyList,
+    Standard_Integer nTransFaces = reader.TransferList(myList);
 
-    AssertThrow(nTransFaces > 0, 
-               ExcMessage("Read nothing from file."));
+    AssertThrow(nTransFaces > 0,
+                ExcMessage("Read nothing from file."));
 
-                                    // Handle IGES Scale here.
+    // Handle IGES Scale here.
     gp_Pnt Origin;
     gp_Trsf scale;
     scale.SetScale (Origin, scale_factor);
 
     TopoDS_Shape sh = reader.OneShape();
     BRepBuilderAPI_Transform trans(sh, scale);
-    
+
     return trans.Shape();   // this is the actual translation
   }
-  
-  void write_IGES(const TopoDS_Shape &shape, 
-                 const std::string &filename) {
-    
+
+  void write_IGES(const TopoDS_Shape &shape,
+                  const std::string &filename)
+  {
+
     IGESControl_Controller::Init();
     IGESControl_Writer ICW ("MM", 0);
     Standard_Boolean ok = ICW.AddShape (shape);
@@ -171,96 +216,89 @@ namespace OpenCASCADE
     Standard_Boolean OK = ICW.Write (filename.c_str());
     Assert(OK, ExcMessage("Failed to write IGES file."));
   }
-  
-
-   TopoDS_Shape intersect_plane(const TopoDS_Shape &in_shape,
-                               const double c_x,
-                               const double c_y,
-                               const double c_z,
-                               const double c,
-                               const double tolerance) 
-   {
-     TopoDS_Shape out_shape;
-     Handle(Geom_Plane) plane = new Geom_Plane(c_x,c_y,c_z,c);
-
-     // This is to loop on the faces,
-     // that extracts all
-     // intersections.
-   
-     //TopExp_Explorer faceExplorer(in_shape , TopAbs_FACE);
-
-     std::vector<Handle_Geom_BoundedCurve> intersections;
-     BRepAlgo_Section section(in_shape, plane);
-     TopoDS_Shape edges = section.Shape();
-
-     TopoDS_Edge edge;
-     TopLoc_Location L;
-     Standard_Real First;
-     Standard_Real Last;   
-     gp_Pnt PIn(0.0,0.0,0.0);
-     gp_Pnt PFin(0.0,0.0,0.0);
-     gp_Pnt PMid(0.0,0.0,0.0); 
-     TopExp_Explorer edgeExplorer(edges , TopAbs_EDGE);
-     while (edgeExplorer.More())
-       {
-        edge = TopoDS::Edge(edgeExplorer.Current());
-        Handle(Geom_Curve) curve = BRep_Tool::Curve(edge,L,First,Last);
-        intersections.push_back(Handle(Geom_BoundedCurve)::DownCast(curve));
-        edgeExplorer.Next();
-       }
-  
-     // Now we build a single bspline out of all the geometrical
-     // curves
-     unsigned int numIntersEdges = intersections.size();  
-     for (unsigned int i = 0;i<intersections.size();++i)
-       {
-        if (intersections[i]->Value(intersections[i]->FirstParameter()).X() > 
-            intersections[i]->Value(intersections[i]->LastParameter()).X()  )
-          intersections[i]->Reverse();
-       }
-
-     GeomConvert_CompCurveToBSplineCurve
-       convert_bspline(intersections[0], Convert_TgtThetaOver2);
-     bool check = false, one_added = true, one_failed=true;
-     std::vector<bool> added(numIntersEdges, false);
-     added[0] = true;
-     while(one_added == true) 
-       {
-        one_added = false;
-        one_failed = false;
-        for (unsigned int i=1; i<numIntersEdges; ++i)
-          if(added[i] == false) 
-            { 
-              Handle(Geom_Curve) curve = intersections[i];
-              Handle(Geom_BoundedCurve) bcurve = Handle(Geom_BoundedCurve)::DownCast(curve);
-              check = convert_bspline.Add(bcurve,tolerance,0,1,0);
-              one_failed = one_failed || (check == false);
-              one_added = one_added || (check == true);
-              added[i] = check;
-              //cout<<i<<" -->  "<<added[i]<<"  "<<false<<endl;
-            }
-       }
-     Assert(one_failed == false,
-           ExcMessage("Bspline convertion of intersection with plane has failed."));
-
-     Handle(Geom_Curve) bspline = convert_bspline.BSplineCurve();
-
-     out_shape = BRepBuilderAPI_MakeEdge(bspline);
-
-     if (bspline->IsCN(1))
-       cout<<"Intersection with plane is at least a C1 curve"<<endl;
-     else
-       cout<<"Intersection with plane is not a C1 curve "<<endl;
-     return out_shape;
-   }
-
-
-  Point<3> axis_intersection(const TopoDS_Shape in_shape, 
-                            const Point<3> origin, 
-                            const Point<3> direction, 
-                            const double tolerance)
+
+
+  TopoDS_Shape intersect_plane(const TopoDS_Shape &in_shape,
+                               const double c_x,
+                               const double c_y,
+                               const double c_z,
+                               const double c,
+                               const double tolerance)
+  {
+    Handle(Geom_Plane) plane = new Geom_Plane(c_x,c_y,c_z,c);
+    BRepAlgo_Section section(in_shape, plane);
+    TopoDS_Shape edges = section.Shape();
+    return edges;
+  }
+
+  TopoDS_Edge join_edges(const TopoDS_Shape &in_shape,
+                         const double tolerance)
+  {
+    TopoDS_Edge out_shape;
+    TopoDS_Shape edges = in_shape;
+    std::vector<Handle_Geom_BoundedCurve> intersections;
+    TopLoc_Location L;
+    Standard_Real First;
+    Standard_Real Last;
+    gp_Pnt PIn(0.0,0.0,0.0);
+    gp_Pnt PFin(0.0,0.0,0.0);
+    gp_Pnt PMid(0.0,0.0,0.0);
+    TopExp_Explorer edgeExplorer(edges , TopAbs_EDGE);
+    TopoDS_Edge edge;
+    while (edgeExplorer.More())
+      {
+        edge = TopoDS::Edge(edgeExplorer.Current());
+        Handle(Geom_Curve) curve = BRep_Tool::Curve(edge,L,First,Last);
+        intersections.push_back(Handle(Geom_BoundedCurve)::DownCast(curve));
+        edgeExplorer.Next();
+      }
+
+    // Now we build a single bspline out of all the geometrical
+    // curves, in Lexycographical order
+    unsigned int numIntersEdges = intersections.size();
+    Assert(numIntersEdges>0, ExcMessage("No curves to process!"));
+
+    GeomConvert_CompCurveToBSplineCurve convert_bspline(intersections[0]);
+
+    bool check = false, one_added = true, one_failed=true;
+    std::vector<bool> added(numIntersEdges, false);
+    added[0] = true;
+    while (one_added == true)
+      {
+        one_added = false;
+        one_failed = false;
+        for (unsigned int i=1; i<numIntersEdges; ++i)
+          if (added[i] == false)
+            {
+              Handle(Geom_Curve) curve = intersections[i];
+              Handle(Geom_BoundedCurve) bcurve = Handle(Geom_BoundedCurve)::DownCast(curve);
+              check = convert_bspline.Add(bcurve,tolerance,0,1,0);
+              if (check == false)  // If we failed, try again with the reversed curve
+                {
+                  curve->Reverse();
+                  Handle(Geom_BoundedCurve) bcurve = Handle(Geom_BoundedCurve)::DownCast(curve);
+                  check = convert_bspline.Add(bcurve,tolerance,0,1,0);
+                }
+              one_failed = one_failed || (check == false);
+              one_added = one_added || (check == true);
+              added[i] = check;
+            }
+      }
+
+    Assert(one_failed == false,
+           ExcMessage("Joining some of the Edges failed."));
+
+    Handle(Geom_Curve) bspline = convert_bspline.BSplineCurve();
+
+    out_shape = BRepBuilderAPI_MakeEdge(bspline);
+    return out_shape;
+  }
+
+
+  Point<3> axis_intersection(const TopoDS_Shape in_shape,
+                             const Point<3> origin,
+                             const Point<3> direction,
+                             const double tolerance)
   {
     // translating original Point<dim> to gp point
 
@@ -275,111 +313,143 @@ namespace OpenCASCADE
     // shape from the iges model
     IntCurvesFace_ShapeIntersector Inters;
     Inters.Load(in_shape,tolerance);
-       
+
     Inters.PerformNearest(line,-RealLast(),+RealLast());
     Assert(Inters.IsDone(), ExcMessage("Could not project point."));
     return Pnt(Inters.Pnt(1));
   }
 
   TopoDS_Edge interpolation_curve(std::vector<Point<3> > &curve_points,
-                                  const Point<3> direction,
-                                  const bool closed,
-                                  const double tolerance)
+                                  const Point<3> direction,
+                                  const bool closed,
+                                  const double tolerance)
   {
 
     unsigned int n_vertices = curve_points.size();
 
-    if(direction*direction > 0) 
+    if (direction*direction > 0)
       {
-       std::sort(curve_points.begin(), curve_points.end(),
-                 boost::bind(&OpenCASCADE::point_compare, _1, _2, direction, tolerance));
+        std::sort(curve_points.begin(), curve_points.end(),
+                  boost::bind(&OpenCASCADE::point_compare, _1, _2, direction, tolerance));
       }
 
-                                    // set up array of vertices
+    // set up array of vertices
     Handle(TColgp_HArray1OfPnt) vertices = new TColgp_HArray1OfPnt(1,n_vertices);
     for (unsigned int vertex=0; vertex<n_vertices; ++vertex)
       {
-       vertices->SetValue(vertex+1,Pnt(curve_points[vertex]));
+        vertices->SetValue(vertex+1,Pnt(curve_points[vertex]));
       }
 
 
     GeomAPI_Interpolate bspline_generator(vertices, closed, tolerance);
     bspline_generator.Perform();
     Assert( (bspline_generator.IsDone()), ExcMessage("Interpolated bspline generation failed"));
-    
+
     Handle(Geom_BSplineCurve) bspline = bspline_generator.Curve();
     TopoDS_Edge out_shape = BRepBuilderAPI_MakeEdge(bspline);
     return out_shape;
   }
 
 
-  Point<3> closest_point(const TopoDS_Shape in_shape, 
-                        const Point<3> origin,
-                        TopoDS_Shape &out_shape,
-                        double &u, 
-                        double &v, 
-                        const double tolerance) {
-    
+  Point<3> closest_point(const TopoDS_Shape in_shape,
+                         const Point<3> origin,
+                         TopoDS_Shape &out_shape,
+                         double &u,
+                         double &v,
+                         const double tolerance)
+  {
+
     TopExp_Explorer exp;
     gp_Pnt Pproj = Pnt(origin);
 
     double minDistance = 1e7;
     gp_Pnt tmp_proj(0.0,0.0,0.0);
-    
+
     unsigned int counter = 0;
-    u=0; v=0;
-    
-    for(exp.Init(in_shape, TopAbs_FACE); exp.More(); exp.Next()) {
-      TopoDS_Face face = TopoDS::Face(exp.Current());
-      
-      // the projection function needs a surface, so we obtain the
-      // surface upon which the face is defined
-      Handle(Geom_Surface) SurfToProj = BRep_Tool::Surface(face);
-      
-      ShapeAnalysis_Surface projector(SurfToProj);
-      gp_Pnt2d proj_params = projector.ValueOfUV(Pnt(origin), tolerance);
-      
-      SurfToProj->D0(proj_params.X(),proj_params.Y(),tmp_proj);
-      
-      double distance = Pnt(tmp_proj).distance(origin);
-      if (distance < minDistance)
-       {
-         minDistance = distance;
-         Pproj = tmp_proj;
-         out_shape = face;
-         u=proj_params.X();
-         v=proj_params.Y();
-         ++counter;
-       }
-    }
-
-    for(exp.Init(in_shape, TopAbs_EDGE); exp.More(); exp.Next()) {
-      TopoDS_Edge edge = TopoDS::Edge(exp.Current());
-
-      TopLoc_Location L;
-      Standard_Real First;
-      Standard_Real Last;
-      
-      // the projection function needs a Curve, so we obtain the
-      // curve upon which the edge is defined
-      Handle(Geom_Curve) CurveToProj = BRep_Tool::Curve(edge,L,First,Last);
-
-      GeomAPI_ProjectPointOnCurve Proj(Pnt(origin),CurveToProj);
-      unsigned int num_proj_points = Proj.NbPoints();
-      if ((num_proj_points > 0) && (Proj.LowerDistance() < minDistance))
-       {
-         minDistance = Proj.LowerDistance();
-         Pproj = Proj.NearestPoint();
-         out_shape = edge;
-         u=Proj.LowerDistanceParameter();
-         ++counter;
-       }
-    }
-    
+    u=0;
+    v=0;
+
+    for (exp.Init(in_shape, TopAbs_FACE); exp.More(); exp.Next())
+      {
+        TopoDS_Face face = TopoDS::Face(exp.Current());
+
+        // the projection function needs a surface, so we obtain the
+        // surface upon which the face is defined
+        Handle(Geom_Surface) SurfToProj = BRep_Tool::Surface(face);
+
+        ShapeAnalysis_Surface projector(SurfToProj);
+        gp_Pnt2d proj_params = projector.ValueOfUV(Pnt(origin), tolerance);
+
+        SurfToProj->D0(proj_params.X(),proj_params.Y(),tmp_proj);
+
+        double distance = Pnt(tmp_proj).distance(origin);
+        if (distance < minDistance)
+          {
+            minDistance = distance;
+            Pproj = tmp_proj;
+            out_shape = face;
+            u=proj_params.X();
+            v=proj_params.Y();
+            ++counter;
+          }
+      }
+
+    for (exp.Init(in_shape, TopAbs_EDGE); exp.More(); exp.Next())
+      {
+        TopoDS_Edge edge = TopoDS::Edge(exp.Current());
+        if (!BRep_Tool::Degenerated(edge))
+          {
+            TopLoc_Location L;
+            Standard_Real First;
+            Standard_Real Last;
+
+            // the projection function needs a Curve, so we obtain the
+            // curve upon which the edge is defined
+            Handle(Geom_Curve) CurveToProj = BRep_Tool::Curve(edge,L,First,Last);
+
+            GeomAPI_ProjectPointOnCurve Proj(Pnt(origin),CurveToProj);
+            unsigned int num_proj_points = Proj.NbPoints();
+            if ((num_proj_points > 0) && (Proj.LowerDistance() < minDistance))
+              {
+                minDistance = Proj.LowerDistance();
+                Pproj = Proj.NearestPoint();
+                out_shape = edge;
+                u=Proj.LowerDistanceParameter();
+                ++counter;
+              }
+          }
+      }
+
     Assert(counter > 0, ExcMessage("Could not find projection points."));
     return Pnt(Pproj);
   }
 
+  void create_triangulation(const TopoDS_Face &face,
+                            Triangulation<2,3> &tria)
+  {
+    BRepAdaptor_Surface surf(face);
+    const double u0 = surf.FirstUParameter();
+    const double u1 = surf.LastUParameter();
+    const double v0 = surf.FirstVParameter();
+    const double v1 = surf.LastVParameter();
+
+    std::vector<CellData<2> > cells;
+    std::vector<Point<3> > vertices;
+    SubCellData t;
+
+    vertices.push_back(Pnt(surf.Value(u0,v0)));
+    vertices.push_back(Pnt(surf.Value(u1,v0)));
+    vertices.push_back(Pnt(surf.Value(u0,v1)));
+    vertices.push_back(Pnt(surf.Value(u1,v1)));
+
+    CellData<2> cell;
+    for (unsigned int i=0; i<4; ++i)
+      cell.vertices[i] = i;
+
+    cells.push_back(cell);
+    tria.create_triangulation(vertices, cells, t);
+  }
+
 } // end namespace
 
 DEAL_II_NAMESPACE_CLOSE
diff --git a/tests/opencascade/arclength_boundary_01.cc b/tests/opencascade/arclength_boundary_01.cc
new file mode 100644 (file)
index 0000000..7d7ad2c
--- /dev/null
@@ -0,0 +1,57 @@
+//-----------------------------------------------------------
+//
+//    Copyright (C) 2014 by the deal.II authors 
+//
+//    This file is subject to LGPL and may not be distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//-----------------------------------------------------------
+
+#include "../tests.h"
+#include <fstream>
+#include <base/logstream.h>
+
+#include <deal.II/grid/occ_utilities.h>
+#include <deal.II/grid/occ_boundary_lib.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/grid_out.h>
+#include <deal.II/grid/grid_generator.h>
+
+#include <TopTools.hxx>
+#include <TopoDS_Shape.hxx>
+#include <Standard_Stream.hxx>
+
+// Create a Triangulation, interpolate its boundary points to a smooth
+// BSpline, and use that as an arlength Boundary Descriptor.
+
+using namespace OpenCASCADE;
+
+int main () 
+{
+  std::ofstream logfile("output");
+  
+  // Create a bspline passign through the points
+  std::vector<Point<3> > pts;
+  pts.push_back(Point<3>(0,0,0));
+  pts.push_back(Point<3>(0,1,0));
+  pts.push_back(Point<3>(1,1,0));
+  pts.push_back(Point<3>(1,0,0));
+
+  TopoDS_Edge edge = interpolation_curve(pts);
+  ArclengthProjectionLineManifold<1,3> manifold(edge);
+  
+  Triangulation<1,3> tria;
+  GridGenerator::hyper_cube(tria);
+
+  tria.begin_active()->set_all_manifold_ids(1);
+  tria.set_manifold(1, manifold);
+  
+  tria.refine_global(4);
+  GridOut gridout;
+  gridout.write_msh(tria, logfile);
+
+  return 0;
+}
+                  
diff --git a/tests/opencascade/arclength_boundary_01.output b/tests/opencascade/arclength_boundary_01.output
new file mode 100644 (file)
index 0000000..94fa0fe
--- /dev/null
@@ -0,0 +1,39 @@
+$NOD
+17
+1  0 0 0
+2  1 0 0
+3  0.5 1.125 0
+4  -0.213053 0.770812 0
+5  1.21305 0.770812 0
+6  -0.209733 0.361552 0
+7  0.090632 1.04613 0
+8  0.909368 1.04613 0
+9  1.20973 0.361552 0
+10  -0.118931 0.172725 0
+11  -0.251438 0.566341 0
+12  -0.0875601 0.936891 0
+13  0.291313 1.106 0
+14  0.708687 1.106 0
+15  1.08756 0.936891 0
+16  1.25144 0.566341 0
+17  1.11893 0.172725 0
+$ENDNOD
+$ELM
+16
+1 1 0 0 2 1 10 
+2 1 0 0 2 10 6 
+3 1 0 0 2 6 11 
+4 1 0 0 2 11 4 
+5 1 0 0 2 4 12 
+6 1 0 0 2 12 7 
+7 1 0 0 2 7 13 
+8 1 0 0 2 13 3 
+9 1 0 0 2 3 14 
+10 1 0 0 2 14 8 
+11 1 0 0 2 8 15 
+12 1 0 0 2 15 5 
+13 1 0 0 2 5 16 
+14 1 0 0 2 16 9 
+15 1 0 0 2 9 17 
+16 1 0 0 2 17 2 
+$ENDELM
diff --git a/tests/opencascade/arclength_boundary_02.cc b/tests/opencascade/arclength_boundary_02.cc
new file mode 100644 (file)
index 0000000..c8ad93c
--- /dev/null
@@ -0,0 +1,80 @@
+//-----------------------------------------------------------
+//
+//    Copyright (C) 2014 by the deal.II authors
+//
+//    This file is subject to LGPL and may not be distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//-----------------------------------------------------------
+
+// Create a Triangulation, interpolate its boundary points to a close
+// smooth BSpline, and use that as a Boundary Descriptor.
+
+#include "../tests.h"
+
+#include <deal.II/grid/occ_utilities.h>
+#include <deal.II/grid/occ_boundary_lib.h>
+
+#include <fstream>
+#include <base/logstream.h>
+#include <grid/tria.h>
+#include <grid/tria_accessor.h>
+#include <grid/grid_out.h>
+#include <grid/tria_iterator.h>
+#include <grid/grid_generator.h>
+#include <grid/tria_boundary_lib.h>
+
+#include <gp_Pnt.hxx>
+#include <gp_Dir.hxx>
+#include <gp_Ax2.hxx>
+#include <GC_MakeCircle.hxx>
+#include <BRepBuilderAPI_MakeEdge.hxx>
+#include <BRepBuilderAPI_MakeFace.hxx>
+#include <BRepBuilderAPI_MakeWire.hxx>
+#include <TopoDS_Wire.hxx>
+#include <TopoDS_Edge.hxx>
+#include <TopoDS_Face.hxx>
+
+using namespace OpenCASCADE;
+
+int main()
+{
+  std::ofstream logfile("output");
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+
+  // A curve passing through the vertices of the unit square.
+  std::vector<Point<3> > vertices;
+  vertices.push_back(Point<3>(0.5,.5-std::sqrt(2.0)/2.0,0));
+  vertices.push_back(Point<3>(1,0,0));
+  vertices.push_back(Point<3>(1,1,0));
+  vertices.push_back(Point<3>(0,1,0));
+  vertices.push_back(Point<3>(0,0,0));
+  TopoDS_Edge edge = interpolation_curve(vertices, Point<3>(), true);
+  
+  // Create a boundary projector.
+  ArclengthProjectionLineManifold<2,3> boundary_line(edge);
+
+  // The unit square.
+  Triangulation<2,3> tria;
+  GridGenerator::hyper_cube(tria);
+
+  // Set the exterior boundary
+  tria.set_manifold(1, boundary_line);
+
+  // This is here to ignore the points created in the interior of the
+  // face.
+  tria.begin()->set_all_manifold_ids(1);
+  tria.begin()->set_manifold_id(0);
+
+  // We refine twice, and expect the outer points to end up on a
+  // smooth curve interpolating the square vertices.
+  tria.refine_global(2);
+
+
+  // You can open the generated file with gmsh
+  GridOut gridout;
+  gridout.write_msh (tria, logfile);
+  return 0;
+}
diff --git a/tests/opencascade/arclength_boundary_02.output b/tests/opencascade/arclength_boundary_02.output
new file mode 100644 (file)
index 0000000..d736f65
--- /dev/null
@@ -0,0 +1,48 @@
+
+$NOD
+25
+1  0.00000 0.00000 0.00000
+2  1.00000 0.00000 0.00000
+3  0.00000 1.00000 0.00000
+4  1.00000 1.00000 0.00000
+5  0.506410 -0.207730 0.00000
+6  -0.166677 0.501960 0.00000
+7  1.16638 0.503363 0.00000
+8  0.499605 1.18595 0.00000
+9  0.500715 0.497943 0.00000
+10  0.231188 -0.153665 0.00000
+11  0.782598 -0.172689 0.00000
+12  -0.124206 0.236965 0.00000
+13  -0.131434 0.767789 0.00000
+14  1.12068 0.239092 0.00000
+15  1.13309 0.769132 0.00000
+16  0.230639 1.14264 0.00000
+17  0.768406 1.14135 0.00000
+18  0.503562 0.145107 0.00000
+19  0.500160 0.841948 0.00000
+20  0.167019 0.499952 0.00000
+21  0.833548 0.500653 0.00000
+22  0.202251 0.190066 0.00000
+23  0.801737 0.188217 0.00000
+24  0.200003 0.804773 0.00000
+25  0.800238 0.805043 0.00000
+$ENDNOD
+$ELM
+16
+1 3 0 0 4 1 10 22 12 
+2 3 0 0 4 10 5 18 22 
+3 3 0 0 4 12 22 20 6 
+4 3 0 0 4 22 18 9 20 
+5 3 0 0 4 5 11 23 18 
+6 3 0 0 4 11 2 14 23 
+7 3 0 0 4 18 23 21 9 
+8 3 0 0 4 23 14 7 21 
+9 3 0 0 4 6 20 24 13 
+10 3 0 0 4 20 9 19 24 
+11 3 0 0 4 13 24 16 3 
+12 3 0 0 4 24 19 8 16 
+13 3 0 0 4 9 21 25 19 
+14 3 0 0 4 21 7 15 25 
+15 3 0 0 4 19 25 17 8 
+16 3 0 0 4 25 15 4 17 
+$ENDELM
diff --git a/tests/opencascade/axis_boundary_01.cc b/tests/opencascade/axis_boundary_01.cc
new file mode 100644 (file)
index 0000000..cb482c4
--- /dev/null
@@ -0,0 +1,65 @@
+//-----------------------------------------------------------
+//
+//    Copyright (C) 2014 by the deal.II authors 
+//
+//    This file is subject to LGPL and may not be distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//-----------------------------------------------------------
+
+// Create a BSpline surface, and test axis projection. 
+
+#include "../tests.h"
+#include <fstream>
+#include <base/logstream.h>
+
+#include <deal.II/grid/occ_utilities.h>
+#include <deal.II/grid/occ_boundary_lib.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/grid_out.h>
+#include <deal.II/grid/grid_generator.h>
+
+#include <TopTools.hxx>
+#include <TopoDS_Shape.hxx>
+#include <TopoDS_Edge.hxx>
+#include <TopoDS_Face.hxx>
+#include <BRepFill.hxx>
+#include <Standard_Stream.hxx>
+
+using namespace OpenCASCADE;
+
+int main () 
+{
+  std::ofstream logfile("output");
+  
+  // Create a bspline passign through the points
+  std::vector<Point<3> > pts;
+  pts.push_back(Point<3>(0,0,0));
+  pts.push_back(Point<3>(0,1,0));
+  pts.push_back(Point<3>(1,1,0));
+  pts.push_back(Point<3>(1,0,0));
+
+  TopoDS_Edge edge1 = interpolation_curve(pts);
+  for(unsigned int i=0; i<pts.size(); ++i)
+    pts[i] += Point<3>(0,0,1);
+  TopoDS_Edge edge2 = interpolation_curve(pts);
+
+  TopoDS_Face face = BRepFill::Face (edge1, edge2);
+
+
+  DirectionalProjectionBoundary<1,3> manifold(face, Point<3>(0,1,0));
+  
+  Triangulation<1,3> tria;
+  GridGenerator::hyper_cube(tria);
+
+  tria.begin_active()->set_all_manifold_ids(1);
+  tria.set_manifold(1, manifold);
+  
+  tria.refine_global(4);
+  GridOut gridout;
+  gridout.write_msh(tria, logfile);
+
+  return 0;
+}                  
diff --git a/tests/opencascade/axis_boundary_01.output b/tests/opencascade/axis_boundary_01.output
new file mode 100644 (file)
index 0000000..09b70d4
--- /dev/null
@@ -0,0 +1,39 @@
+$NOD
+17
+1  0 0 0
+2  1 0 0
+3  0.5 1.125 0
+4  0.25 1.09745 0
+5  0.75 1.09745 0
+6  0.125 1.05999 0
+7  0.375 1.11829 0
+8  0.625 1.11829 0
+9  0.875 1.05999 0
+10  0.0625 1.03342 0
+11  0.1875 1.08105 0
+12  0.3125 1.10974 0
+13  0.4375 1.12333 0
+14  0.5625 1.12333 0
+15  0.6875 1.10974 0
+16  0.8125 1.08105 0
+17  0.9375 1.03342 0
+$ENDNOD
+$ELM
+16
+1 1 0 0 2 1 10 
+2 1 0 0 2 10 6 
+3 1 0 0 2 6 11 
+4 1 0 0 2 11 4 
+5 1 0 0 2 4 12 
+6 1 0 0 2 12 7 
+7 1 0 0 2 7 13 
+8 1 0 0 2 13 3 
+9 1 0 0 2 3 14 
+10 1 0 0 2 14 8 
+11 1 0 0 2 8 15 
+12 1 0 0 2 15 5 
+13 1 0 0 2 5 16 
+14 1 0 0 2 16 9 
+15 1 0 0 2 9 17 
+16 1 0 0 2 17 2 
+$ENDELM
diff --git a/tests/opencascade/axis_boundary_02.cc b/tests/opencascade/axis_boundary_02.cc
new file mode 100644 (file)
index 0000000..67fc1a8
--- /dev/null
@@ -0,0 +1,68 @@
+
+//----------------------------  iges_read.cc  ---------------------------
+//    $Id: testsuite.html 13373 2006-07-13 13:12:08Z kanschat $
+//    Version: $Name$ 
+//
+//    Copyright (C) 2005 by the deal.II authors 
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//----------------------------  iges_read.cc  ---------------------------
+
+
+// Read goteborg.iges and dump its topological structure to the logfile.
+
+#include "../tests.h"
+#include <fstream>
+#include <base/logstream.h>
+
+#include <deal.II/grid/occ_utilities.h>
+#include <deal.II/grid/occ_boundary_lib.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/grid_out.h>
+#include <deal.II/grid/grid_generator.h>
+
+#include <TopTools.hxx>
+#include <TopoDS_Shape.hxx>
+#include <TopoDS_Edge.hxx>
+#include <TopoDS_Face.hxx>
+#include <BRepFill.hxx>
+#include <Standard_Stream.hxx>
+
+using namespace OpenCASCADE;
+
+int main () 
+{
+  std::ofstream logfile("output");
+
+  // Create a bspline passing through the points
+  std::vector<Point<3> > pts1, pts2;
+  pts1.push_back(Point<3>(0,0,0));
+  pts1.push_back(Point<3>(1,0,0));
+
+  pts2.push_back(Point<3>(0,1,0));
+  pts2.push_back(Point<3>(.5,1,1));
+  pts2.push_back(Point<3>(1,1,0));
+  
+  TopoDS_Edge edge1 = interpolation_curve(pts1);
+  TopoDS_Edge edge2 = interpolation_curve(pts2);
+  
+  TopoDS_Face face = BRepFill::Face (edge1, edge2);
+
+  DirectionalProjectionBoundary<2,3> manifold(face, Point<3>(0,0,1));
+  
+  Triangulation<2,3> tria;
+  GridGenerator::hyper_cube(tria);
+
+  tria.begin_active()->set_all_manifold_ids(1);
+  tria.set_manifold(1, manifold);
+
+  tria.refine_global(2);
+  GridOut gridout;
+  gridout.write_msh(tria, logfile);
+  
+  return 0;
+}                  
diff --git a/tests/opencascade/axis_boundary_02.output b/tests/opencascade/axis_boundary_02.output
new file mode 100644 (file)
index 0000000..a007fc2
--- /dev/null
@@ -0,0 +1,47 @@
+$NOD
+25
+1  0 0 0
+2  1 0 0
+3  0 1 0
+4  1 1 0
+5  0.5 0 0
+6  0 0.5 0
+7  1 0.5 0
+8  0.5 1 1
+9  0.5 0.5 0.5
+10  0.25 0 0
+11  0.75 0 0
+12  0 0.25 0
+13  0 0.75 0
+14  1 0.25 0
+15  1 0.75 0
+16  0.25 1 0.75
+17  0.75 1 0.75
+18  0.5 0.25 0.25
+19  0.5 0.75 0.75
+20  0.25 0.5 0.375
+21  0.75 0.5 0.375
+22  0.25 0.25 0.1875
+23  0.75 0.25 0.1875
+24  0.25 0.75 0.5625
+25  0.75 0.75 0.5625
+$ENDNOD
+$ELM
+16
+1 3 0 0 4 1 10 22 12 
+2 3 0 0 4 10 5 18 22 
+3 3 0 0 4 12 22 20 6 
+4 3 0 0 4 22 18 9 20 
+5 3 0 0 4 5 11 23 18 
+6 3 0 0 4 11 2 14 23 
+7 3 0 0 4 18 23 21 9 
+8 3 0 0 4 23 14 7 21 
+9 3 0 0 4 6 20 24 13 
+10 3 0 0 4 20 9 19 24 
+11 3 0 0 4 13 24 16 3 
+12 3 0 0 4 24 19 8 16 
+13 3 0 0 4 9 21 25 19 
+14 3 0 0 4 21 7 15 25 
+15 3 0 0 4 19 25 17 8 
+16 3 0 0 4 25 15 4 17 
+$ENDELM
index f40c3edeb7821b0db4d4f622a49392b01a38545c..6c34f047b67b9586a37524ab20365facde361def 100644 (file)
@@ -1,17 +1,13 @@
-
-//----------------------------  occ_circle_projection.cc  ---------------------------
-//    $Id: testsuite.html 13373 2006-07-13 13:12:08Z kanschat $
-//    Version: $Name$ 
+//-----------------------------------------------------------
 //
-//    Copyright (C) 2005 by the deal.II authors 
+//    Copyright (C) 2014 by the deal.II authors 
 //
-//    This file is subject to QPL and may not be  distributed
+//    This file is subject to LGPL and may not be distributed
 //    without copyright and license information. Please refer
 //    to the file deal.II/doc/license.html for the  text  and
 //    further information on this license.
 //
-//----------------------------  occ_circle_projection.cc  ---------------------------
-
+//-----------------------------------------------------------
 
 // Create a circle, a Triangulation, and try to project normally on
 // it.
@@ -76,10 +72,6 @@ int main ()
 
                                   // Set the exterior boundary
   tria.set_boundary(0, boundary_line);
-                                  // Uncomment the following line to
-                                  // see the difference.
-// tria.set_boundary(0, boundary_line_deal);
-  
 
                                   // This is here to ignore the
                                   // points created in the interior
index b1ee136e2c1471db97c3fbedc77ae02c9228467b..f87b866aa63df05a28b767b44c2a03d8e185984e 100644 (file)
@@ -1,3 +1,14 @@
+//-----------------------------------------------------------
+//
+//    Copyright (C) 2014 by the deal.II authors 
+//
+//    This file is subject to LGPL and may not be distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//-----------------------------------------------------------
+
 // Makes a OpenCASCADE cicrular arc, and project a few points onto it.
 
 
diff --git a/tests/opencascade/create_tria_01.cc b/tests/opencascade/create_tria_01.cc
new file mode 100644 (file)
index 0000000..93a0ec1
--- /dev/null
@@ -0,0 +1,73 @@
+//-----------------------------------------------------------
+//
+//    Copyright (C) 2014 by the deal.II authors 
+//
+//    This file is subject to LGPL and may not be distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//-----------------------------------------------------------
+
+// Read the first face of goteborg.iges, attach it to a projector,
+// create a single-cell Triangulation, and refine it with the new
+// projector.
+
+#include "../tests.h"
+
+#include <deal.II/grid/occ_boundary_lib.h>
+
+#include <fstream>
+#include <base/logstream.h>
+#include <grid/tria.h>
+#include <grid/tria_accessor.h>
+#include <grid/grid_out.h>
+#include <grid/tria_iterator.h>
+#include <grid/grid_generator.h>
+#include <grid/tria_boundary_lib.h>
+
+#include <gp_Pnt.hxx>
+#include <gp_Dir.hxx>
+#include <gp_Ax2.hxx>
+#include <GC_MakeCircle.hxx>
+#include <BRepBuilderAPI_MakeEdge.hxx>
+#include <BRepBuilderAPI_MakeFace.hxx>
+#include <BRepBuilderAPI_MakeWire.hxx>
+#include <TopoDS_Wire.hxx>
+#include <TopoDS_Edge.hxx>
+#include <TopoDS_Face.hxx>
+#include <BRepPrimAPI_MakeSphere.hxx>
+
+using namespace OpenCASCADE;
+
+int main () 
+{
+  std::ofstream logfile("output");
+  
+  TopoDS_Shape sh = read_IGES(SOURCE_DIR "/iges_files/goteborg.iges");
+  std::vector<TopoDS_Face> faces;
+  std::vector<TopoDS_Edge> edges;
+  std::vector<TopoDS_Vertex> vertices;
+  
+  extract_geometrical_shapes(sh, faces, edges, vertices);
+
+  // Create a boundary projector on the first face.
+  NormalProjectionBoundary<2,3> boundary(faces[0], 1e-9);
+  
+  // Create a Triangulation with a single cell
+  Triangulation<2,3> tria;
+  create_triangulation(faces[0], tria);
+
+  // Set the boundary
+  tria.set_manifold(1, boundary);
+  tria.begin()->set_all_manifold_ids(1);
+
+  tria.refine_global(2);
+
+  // You can open the generated file with gmsh
+  GridOut gridout;
+  gridout.write_msh (tria, logfile);
+  
+  return 0;
+}
+                  
diff --git a/tests/opencascade/create_tria_01.output b/tests/opencascade/create_tria_01.output
new file mode 100644 (file)
index 0000000..d358f6b
--- /dev/null
@@ -0,0 +1,47 @@
+$NOD
+25
+1  -1.98714 1.02631e-19 -0.149998
+2  -2.63758 -6.93889e-18 -0.249603
+3  -2.00886 0.314666 0.428234
+4  -3.05955 0 0.501418
+5  -2.30472 0 -0.245614
+6  -2.00073 0.167844 0.134255
+7  -2.72534 4.25119e-05 0.204053
+8  -2.54605 0.193924 0.440052
+9  -2.40412 0.0655112 0.139204
+10  -2.14241 8.03383e-19 -0.216122
+11  -2.46936 -3.70788e-20 -0.267918
+12  -1.99609 0.0965804 -0.0146004
+13  -2.00426 0.238868 0.282461
+14  -2.63221 9.34116e-07 -0.0158842
+15  -2.87756 -5.67114e-09 0.370719
+16  -2.27934 0.26212 0.42932
+17  -2.80896 0.111433 0.459585
+18  -2.35424 0.031169 -0.0529932
+19  -2.47235 0.120909 0.295137
+20  -2.20235 0.116409 0.136836
+21  -2.56575 0.0372536 0.170331
+22  -2.17413 0.0607452 -0.0340272
+23  -2.51144 0.0146272 -0.0396781
+24  -2.23928 0.183308 0.286763
+25  -2.68303 0.0677305 0.321656
+$ENDNOD
+$ELM
+16
+1 3 0 0 4 1 10 22 12 
+2 3 0 0 4 10 5 18 22 
+3 3 0 0 4 12 22 20 6 
+4 3 0 0 4 22 18 9 20 
+5 3 0 0 4 5 11 23 18 
+6 3 0 0 4 11 2 14 23 
+7 3 0 0 4 18 23 21 9 
+8 3 0 0 4 23 14 7 21 
+9 3 0 0 4 6 20 24 13 
+10 3 0 0 4 20 9 19 24 
+11 3 0 0 4 13 24 16 3 
+12 3 0 0 4 24 19 8 16 
+13 3 0 0 4 9 21 25 19 
+14 3 0 0 4 21 7 15 25 
+15 3 0 0 4 19 25 17 8 
+16 3 0 0 4 25 15 4 17 
+$ENDELM
index c8cb5bfbb17fea9adf43e4f56bb09bb1253fd0ff..3fa06155f2406b4dd9526f13b39204a8d92bd17e 100644 (file)
@@ -1,19 +1,15 @@
-
-//----------------------------  iges_read.cc  ---------------------------
-//    $Id: testsuite.html 13373 2006-07-13 13:12:08Z kanschat $
-//    Version: $Name$ 
+//-----------------------------------------------------------
 //
-//    Copyright (C) 2005 by the deal.II authors 
+//    Copyright (C) 2014 by the deal.II authors 
 //
-//    This file is subject to QPL and may not be  distributed
+//    This file is subject to LGPL and may not be distributed
 //    without copyright and license information. Please refer
 //    to the file deal.II/doc/license.html for the  text  and
 //    further information on this license.
 //
-//----------------------------  iges_read.cc  ---------------------------
-
+//-----------------------------------------------------------
 
-// Read goteborg.iges and dump its topological structure to the logfile.
+// Create a bspline, and write it as an IGES file.
 
 #include "../tests.h"
 #include <fstream>
@@ -28,10 +24,6 @@ using namespace OpenCASCADE;
 
 int main () 
 {
-  // std::ofstream logfile("output");
-  // deallog.attach(logfile);
-  // deallog.depth_console(0);
-  
   // Create a bspline passign through the points
   std::vector<Point<3> > pts;
   pts.push_back(Point<3>(0,0,0));
diff --git a/tests/opencascade/iges_describe.cc b/tests/opencascade/iges_describe.cc
new file mode 100644 (file)
index 0000000..072e3c8
--- /dev/null
@@ -0,0 +1,54 @@
+//-----------------------------------------------------------
+//
+//    Copyright (C) 2014 by the deal.II authors 
+//
+//    This file is subject to LGPL and may not be distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//-----------------------------------------------------------
+
+// Read goteborg.iges and dump its topological structure to the
+// logfile.
+
+#include "../tests.h"
+#include <fstream>
+#include <base/logstream.h>
+
+#include <deal.II/grid/occ_utilities.h>
+#include <TopTools.hxx>
+#include <TopoDS_Shape.hxx>
+#include <Standard_Stream.hxx>
+
+using namespace OpenCASCADE;
+
+int main () 
+{
+  std::ofstream logfile("output");
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+  
+  TopoDS_Shape sh = read_IGES(SOURCE_DIR "/iges_files/goteborg.iges");
+  std::vector<TopoDS_Face> faces;
+  std::vector<TopoDS_Edge> edges;
+  std::vector<TopoDS_Vertex> vertices;
+
+  extract_geometrical_shapes(sh, faces, edges, vertices);
+  
+  unsigned int nf=0, ne=0, nv=0;
+  count_elements(sh, nf, ne, nv);
+  deallog << "Shape contains " << nf << " faces, "
+         << ne << " edges, and " 
+         << nv << " vertices." << std::endl;
+  
+  if(nf != faces.size())
+    deallog << "Error!" << std::endl;
+  if(ne != edges.size())
+    deallog << "Error!" << std::endl;
+  if(nv != vertices.size())
+    deallog << "Error!" << std::endl;
+  
+  return 0;
+}
+                  
diff --git a/tests/opencascade/iges_describe.output b/tests/opencascade/iges_describe.output
new file mode 100644 (file)
index 0000000..8b1892b
--- /dev/null
@@ -0,0 +1,2 @@
+
+DEAL::Shape contains 3 faces, 12 edges, and 24 vertices.
diff --git a/tests/opencascade/iges_describe_02.cc b/tests/opencascade/iges_describe_02.cc
new file mode 100644 (file)
index 0000000..5ce78eb
--- /dev/null
@@ -0,0 +1,47 @@
+//-----------------------------------------------------------
+//
+//    Copyright (C) 2014 by the deal.II authors 
+//
+//    This file is subject to LGPL and may not be distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//-----------------------------------------------------------
+
+// Read goteborg.iges and dump its topological structure to the
+// logfile.
+
+#include "../tests.h"
+#include <fstream>
+#include <base/logstream.h>
+
+#include <deal.II/grid/occ_utilities.h>
+#include <TopTools.hxx>
+#include <TopoDS_Shape.hxx>
+#include <Standard_Stream.hxx>
+
+using namespace OpenCASCADE;
+
+int main () 
+{
+  std::ofstream logfile("output");
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+  
+  TopoDS_Shape sh = read_IGES(SOURCE_DIR "/iges_files/goteborg.iges");
+  std::vector<TopoDS_Compound> compounds;
+  std::vector<TopoDS_CompSolid> compsolids;
+  std::vector<TopoDS_Solid> solids;
+  std::vector<TopoDS_Shell> shells;
+  std::vector<TopoDS_Wire> wires;
+  extract_compound_shapes(sh, compounds, compsolids, solids, shells, wires);
+  
+  deallog << "Shape contains " << compounds.size() << " compounds, "
+         << compsolids.size() << " compound solids, "
+         << solids.size() << " solids, "
+         << shells.size() << " shells, "
+         << wires.size() << " wires." << std::endl;
+  return 0;
+}
+                  
diff --git a/tests/opencascade/iges_describe_02.output b/tests/opencascade/iges_describe_02.output
new file mode 100644 (file)
index 0000000..b760f0d
--- /dev/null
@@ -0,0 +1,2 @@
+
+DEAL::Shape contains 1 compounds, 0 compound solids, 0 solids, 0 shells, 3 wires.
diff --git a/tests/opencascade/iges_files/side.iges b/tests/opencascade/iges_files/side.iges
new file mode 100644 (file)
index 0000000..9f65cc1
--- /dev/null
@@ -0,0 +1,63 @@
+                                                                        S0000001
+,,31HOpen CASCADE IGES processor 6.7,13HFilename.iges,                  G0000001
+16HOpen CASCADE 6.7,31HOpen CASCADE IGES processor 6.7,32,308,15,308,15,G0000002
+,1.,1,4HINCH,1,0.000394,15H20140814.114632,3.93701E-09,0.049421,        G0000003
+6Hheltai,,11,0,15H20140814.114632,;                                     G0000004
+     144       1       0       0       0       0       0       000000000D0000001
+     144       0       0       1       0                               0D0000002
+     122       2       0       0       0       0       0       000010000D0000003
+     122       0       0       1       0                               0D0000004
+     126       3       0       0       0       0       0       000010000D0000005
+     126       0       0       6       0                               0D0000006
+     142       9       0       0       0       0       0       000010500D0000007
+     142       0       0       1       0                               0D0000008
+     102      10       0       0       0       0       0       000010000D0000009
+     102       0       0       1       0                               0D0000010
+     110      11       0       0       0       0       0       000010000D0000011
+     110       0       0       1       0                               0D0000012
+     110      12       0       0       0       0       0       000010000D0000013
+     110       0       0       1       0                               0D0000014
+     110      13       0       0       0       0       0       000010000D0000015
+     110       0       0       1       0                               0D0000016
+     110      14       0       0       0       0       0       000010000D0000017
+     110       0       0       1       0                               0D0000018
+     102      15       0       0       0       0       0       000010000D0000019
+     102       0       0       1       0                               0D0000020
+     110      16       0       0       0       0       0       000010000D0000021
+     110       0       0       1       0                               0D0000022
+     126      17       0       0       0       0       0       000010000D0000023
+     126       0       0       6       0                               0D0000024
+     110      23       0       0       0       0       0       000010000D0000025
+     110       0       0       2       0                               0D0000026
+     126      25       0       0       0       0       0       000010000D0000027
+     126       0       0       5       0                               0D0000028
+144,3,1,0,7;                                                     0000001P0000001
+122,5,0.,0.,0.;                                                  0000003P0000002
+126,5,3,1,0,1,0,0.,0.,0.,0.,1.,2.,3.,3.,3.,3.,1.,1.,1.,1.,1.,1., 0000005P0000003
+0.,0.,3.937007874E-02,-1.531058618E-02,1.968503937E-02,          0000005P0000004
+3.937007874E-02,-6.56167979E-03,4.593175853E-02,3.937007874E-02, 0000005P0000005
+4.593175853E-02,4.593175853E-02,3.937007874E-02,5.468066492E-02, 0000005P0000006
+1.968503937E-02,3.937007874E-02,3.937007874E-02,0.,              0000005P0000007
+3.937007874E-02,0.,3.,1.602995811E-17,3.205991622E-17,1.;        0000005P0000008
+142,0,3,9,19,3;                                                  0000007P0000009
+102,4,11,13,15,17;                                               0000009P0000010
+110,0.,1.,0.,0.,0.,0.;                                           0000011P0000011
+110,0.,0.,0.,1.,0.,0.;                                           0000013P0000012
+110,1.,0.,0.,1.,1.,0.;                                           0000015P0000013
+110,1.,1.,0.,0.,1.,0.;                                           0000017P0000014
+102,4,21,23,25,27;                                               0000019P0000015
+110,0.,0.,0.,0.,0.,3.937007874E-02;                              0000021P0000016
+126,5,3,1,0,1,0,0.,0.,0.,0.,1.,2.,3.,3.,3.,3.,1.,1.,1.,1.,1.,1., 0000023P0000017
+0.,0.,3.937007874E-02,-1.531058618E-02,1.968503937E-02,          0000023P0000018
+3.937007874E-02,-6.56167979E-03,4.593175853E-02,3.937007874E-02, 0000023P0000019
+4.593175853E-02,4.593175853E-02,3.937007874E-02,5.468066492E-02, 0000023P0000020
+1.968503937E-02,3.937007874E-02,3.937007874E-02,0.,              0000023P0000021
+3.937007874E-02,0.,3.,1.602995811E-17,3.205991622E-17,1.;        0000023P0000022
+110,3.937007874E-02,2.185478395E-18,3.937007874E-02,             0000025P0000023
+3.937007874E-02,2.185478395E-18,0.;                              0000025P0000024
+126,5,3,1,0,1,0,0.,0.,0.,0.,1.,2.,3.,3.,3.,3.,1.,1.,1.,1.,1.,1., 0000027P0000025
+3.937007874E-02,0.,0.,5.468066492E-02,1.968503937E-02,0.,        0000027P0000026
+4.593175853E-02,4.593175853E-02,0.,-6.56167979E-03,              0000027P0000027
+4.593175853E-02,0.,-1.531058618E-02,1.968503937E-02,0.,0.,0.,0., 0000027P0000028
+0.,3.,0.,0.,1.;                                                  0000027P0000029
+S      1G      4D     28P     29                                        T0000001
index 023749b95945c22f2239f34e0d1cd5a6a57e85e0..2791a880809567e08afec4d76089744329f07c72 100644 (file)
@@ -1,19 +1,16 @@
-
-//----------------------------  iges_read.cc  ---------------------------
-//    $Id: testsuite.html 13373 2006-07-13 13:12:08Z kanschat $
-//    Version: $Name$ 
+//-----------------------------------------------------------
 //
-//    Copyright (C) 2005 by the deal.II authors 
+//    Copyright (C) 2014 by the deal.II authors 
 //
-//    This file is subject to QPL and may not be  distributed
+//    This file is subject to LGPL and may not be distributed
 //    without copyright and license information. Please refer
 //    to the file deal.II/doc/license.html for the  text  and
 //    further information on this license.
 //
-//----------------------------  iges_read.cc  ---------------------------
-
+//-----------------------------------------------------------
 
-// Read goteborg.iges and dump its topological structure to the logfile.
+// Read goteborg.iges and dump its topological structure to the
+// logfile.
 
 #include "../tests.h"
 #include <fstream>
@@ -39,15 +36,6 @@ int main ()
          << ne << " edges, and " 
          << nv << " vertices." << std::endl;
 
-  // try 
-  //   {
-  //     TopoDS_Shape sh_fail = read_IGES("iges_files/nonexistent.iges");
-  //   }
-  // catch(...)
-  //   {
-  //     deallog << "Missing file correctly threw an exception." << std::endl;
-  //   }
-  
   return 0;
 }
                   
index 6430b15587aac07518ca21a6094d5ee0c4c0181e..ac66f16f8c5727ae32966d1816b07250fd6a963e 100644 (file)
@@ -1,19 +1,16 @@
-
-//----------------------------  iges_read.cc  ---------------------------
-//    $Id: testsuite.html 13373 2006-07-13 13:12:08Z kanschat $
-//    Version: $Name$ 
+//-----------------------------------------------------------
 //
-//    Copyright (C) 2005 by the deal.II authors 
+//    Copyright (C) 2014 by the deal.II authors 
 //
-//    This file is subject to QPL and may not be  distributed
+//    This file is subject to LGPL and may not be distributed
 //    without copyright and license information. Please refer
 //    to the file deal.II/doc/license.html for the  text  and
 //    further information on this license.
 //
-//----------------------------  iges_read.cc  ---------------------------
-
+//-----------------------------------------------------------
 
-// Read goteborg.iges and dump its topological structure to the logfile.
+// Read a file in iges format, and write it out again in the same
+// format.
 
 #include "../tests.h"
 #include <fstream>
@@ -28,28 +25,10 @@ using namespace OpenCASCADE;
 
 int main () 
 {
-  // std::ofstream logfile("output");
-  // deallog.attach(logfile);
-  // deallog.depth_console(0);
-  
   TopoDS_Shape sh = read_IGES(SOURCE_DIR "/iges_files/wigley.iges");
   unsigned int nf=0, ne=0, nv=0;
   count_elements(sh, nf, ne, nv);
-  // deallog << "Shape contains " << nf << " faces, "
-  //     << ne << " edges, and " 
-  //     << nv << " vertices." << std::endl;
 
   write_IGES(sh, "output");
-
-  // try 
-  //   {
-  //     TopoDS_Shape sh_fail = read_IGES("iges_files/nonexistent.iges");
-  //   }
-  // catch(...)
-  //   {
-  //     deallog << "Missing file correctly threw an exception." << std::endl;
-  //   }
-  
   return 0;
 }
-                  
index 9d43fa31cfe0dad0722384f14b829a5969bdf19d..a6de53e59a75d108e96ed4ad4bbc169331b48d0c 100644 (file)
@@ -1,17 +1,13 @@
-
-//----------------------------  occ_circle_projection.cc  ---------------------------
-//    $Id: testsuite.html 13373 2006-07-13 13:12:08Z kanschat $
-//    Version: $Name$ 
+//-----------------------------------------------------------
 //
-//    Copyright (C) 2005 by the deal.II authors 
+//    Copyright (C) 2014 by the deal.II authors 
 //
-//    This file is subject to QPL and may not be  distributed
+//    This file is subject to LGPL and may not be distributed
 //    without copyright and license information. Please refer
 //    to the file deal.II/doc/license.html for the  text  and
 //    further information on this license.
 //
-//----------------------------  occ_circle_projection.cc  ---------------------------
-
+//-----------------------------------------------------------
 
 // Create a Triangulation, interpolate its boundary points to a close
 // smooth BSpline, and use that as a Boundary Descriptor.
@@ -75,13 +71,10 @@ int main ()
   tria.refine_global(2);
 
 
-  // You can open the generated file with paraview.
+  // You can open the generated file with gmsh.
   GridOut gridout;
   gridout.write_msh (tria, logfile);
 
-  std::ofstream of("/tmp/grid.msh");
-  gridout.write_msh(tria, of);
-  
   return 0;
 }
                   
diff --git a/tests/opencascade/normal_projection_01.cc b/tests/opencascade/normal_projection_01.cc
new file mode 100644 (file)
index 0000000..31a2788
--- /dev/null
@@ -0,0 +1,78 @@
+//-----------------------------------------------------------
+//
+//    Copyright (C) 2014 by the deal.II authors 
+//
+//    This file is subject to LGPL and may not be distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//-----------------------------------------------------------
+
+// Create a circle, a Triangulation, and try to project normally on
+// it.
+
+#include "../tests.h"
+
+#include <deal.II/grid/occ_boundary_lib.h>
+
+#include <fstream>
+#include <base/logstream.h>
+#include <grid/tria.h>
+#include <grid/tria_accessor.h>
+#include <grid/grid_out.h>
+#include <grid/tria_iterator.h>
+#include <grid/grid_generator.h>
+#include <grid/tria_boundary_lib.h>
+
+#include <gp_Pnt.hxx>
+#include <gp_Dir.hxx>
+#include <gp_Ax2.hxx>
+#include <GC_MakeCircle.hxx>
+#include <BRepBuilderAPI_MakeEdge.hxx>
+#include <BRepBuilderAPI_MakeFace.hxx>
+#include <BRepBuilderAPI_MakeWire.hxx>
+#include <TopoDS_Wire.hxx>
+#include <TopoDS_Edge.hxx>
+#include <TopoDS_Face.hxx>
+#include <BRepPrimAPI_MakeSphere.hxx>
+
+using namespace OpenCASCADE;
+
+int main () 
+{
+  std::ofstream logfile("output");
+
+  gp_Pnt center(.5,.5,.5);
+  Standard_Real radius(Point<3>().distance(Pnt(center)));
+
+  
+  TopoDS_Face face = BRepPrimAPI_MakeSphere(center, radius);
+  
+                                  // Create a boundary projector.
+  NormalProjectionBoundary<3,3> sphere(face);
+  
+  
+                                  // The unit cube.
+  Triangulation<3,3> tria;
+  GridGenerator::hyper_cube(tria);
+
+                                  // Set the exterior boundary
+  tria.set_manifold(1, sphere);
+
+                                  // This is here to ignore the
+                                  // points created in the interior
+                                  // of the face.
+  tria.begin()->set_all_manifold_ids(1);
+  tria.begin()->set_manifold_id(0);
+
+  tria.refine_global(2);
+
+
+  // You can open the generated file with gmsh
+  GridOut gridout;
+  gridout.write_msh (tria, logfile);
+  
+  return 0;
+}
+                  
diff --git a/tests/opencascade/normal_projection_01.output b/tests/opencascade/normal_projection_01.output
new file mode 100644 (file)
index 0000000..5348f66
--- /dev/null
@@ -0,0 +1,195 @@
+$NOD
+125
+1  0 0 0
+2  1 0 0
+3  0 1 0
+4  1 1 0
+5  0 0 1
+6  1 0 1
+7  0 1 1
+8  1 1 1
+9  0.5 -0.112372 -0.112372
+10  -0.112372 0.5 -0.112372
+11  -0.112372 -0.112372 0.5
+12  1.11237 0.5 -0.112372
+13  1.11237 -0.112372 0.5
+14  0.5 1.11237 -0.112372
+15  -0.112372 1.11237 0.5
+16  1.11237 1.11237 0.5
+17  0.5 -0.112372 1.11237
+18  -0.112372 0.5 1.11237
+19  1.11237 0.5 1.11237
+20  0.5 1.11237 1.11237
+21  0.5 -0.366025 0.5
+22  0.5 0.5 -0.366025
+23  -0.366025 0.5 0.5
+24  1.36603 0.5 0.5
+25  0.5 1.36603 0.5
+26  0.5 0.5 1.36603
+27  0.5 0.5 0.5
+28  0.237676 -0.0836036 -0.0836036
+29  0.762324 -0.0836036 -0.0836036
+30  -0.0836036 0.237676 -0.0836036
+31  -0.0836036 0.762324 -0.0836036
+32  -0.0836036 -0.0836036 0.237676
+33  -0.0836036 -0.0836036 0.762324
+34  1.0836 0.237676 -0.0836036
+35  1.0836 0.762324 -0.0836036
+36  1.0836 -0.0836036 0.237676
+37  1.0836 -0.0836036 0.762324
+38  0.237676 1.0836 -0.0836036
+39  0.762324 1.0836 -0.0836036
+40  -0.0836036 1.0836 0.237676
+41  -0.0836036 1.0836 0.762324
+42  1.0836 1.0836 0.237676
+43  1.0836 1.0836 0.762324
+44  0.237676 -0.0836036 1.0836
+45  0.762324 -0.0836036 1.0836
+46  -0.0836036 0.237676 1.0836
+47  -0.0836036 0.762324 1.0836
+48  1.0836 0.237676 1.0836
+49  1.0836 0.762324 1.0836
+50  0.237676 1.0836 1.0836
+51  0.762324 1.0836 1.0836
+52  0.168586 -0.300103 0.5
+53  0.831414 -0.300103 0.5
+54  0.5 -0.300103 0.168586
+55  0.5 -0.300103 0.831414
+56  0.5 0.168586 -0.300103
+57  0.5 0.831414 -0.300103
+58  0.168586 0.5 -0.300103
+59  0.831414 0.5 -0.300103
+60  -0.300103 0.5 0.168586
+61  -0.300103 0.5 0.831414
+62  -0.300103 0.168586 0.5
+63  -0.300103 0.831414 0.5
+64  1.3001 0.5 0.168586
+65  1.3001 0.5 0.831414
+66  1.3001 0.168586 0.5
+67  1.3001 0.831414 0.5
+68  0.168586 1.3001 0.5
+69  0.831414 1.3001 0.5
+70  0.5 1.3001 0.168586
+71  0.5 1.3001 0.831414
+72  0.5 0.168586 1.3001
+73  0.5 0.831414 1.3001
+74  0.168586 0.5 1.3001
+75  0.831414 0.5 1.3001
+76  0.5 0.5 0.933013
+77  0.5 0.5 0.0669873
+78  0.933013 0.5 0.5
+79  0.0669873 0.5 0.5
+80  0.5 0.933013 0.5
+81  0.5 0.0669873 0.5
+82  0.183268 -0.241189 0.183268
+83  0.183268 -0.241189 0.816732
+84  0.816732 -0.241189 0.183268
+85  0.816732 -0.241189 0.816732
+86  0.183268 0.183268 -0.241189
+87  0.816732 0.183268 -0.241189
+88  0.183268 0.816732 -0.241189
+89  0.816732 0.816732 -0.241189
+90  -0.241189 0.183268 0.183268
+91  -0.241189 0.816732 0.183268
+92  -0.241189 0.183268 0.816732
+93  -0.241189 0.816732 0.816732
+94  1.24119 0.183268 0.183268
+95  1.24119 0.816732 0.183268
+96  1.24119 0.183268 0.816732
+97  1.24119 0.816732 0.816732
+98  0.183268 1.24119 0.183268
+99  0.183268 1.24119 0.816732
+100  0.816732 1.24119 0.183268
+101  0.816732 1.24119 0.816732
+102  0.183268 0.183268 1.24119
+103  0.816732 0.183268 1.24119
+104  0.183268 0.816732 1.24119
+105  0.816732 0.816732 1.24119
+106  0.880366 0.880366 0.5
+107  0.119634 0.880366 0.5
+108  0.880366 0.119634 0.5
+109  0.119634 0.119634 0.5
+110  0.880366 0.5 0.880366
+111  0.880366 0.5 0.119634
+112  0.119634 0.5 0.880366
+113  0.119634 0.5 0.119634
+114  0.5 0.880366 0.880366
+115  0.5 0.119634 0.880366
+116  0.5 0.880366 0.119634
+117  0.5 0.119634 0.119634
+118  0.151399 0.151399 0.151399
+119  0.848601 0.151399 0.151399
+120  0.151399 0.848601 0.151399
+121  0.848601 0.848601 0.151399
+122  0.151399 0.151399 0.848601
+123  0.848601 0.151399 0.848601
+124  0.151399 0.848601 0.848601
+125  0.848601 0.848601 0.848601
+$ENDNOD
+$ELM
+64
+1 5 0 0 8 1 28 82 32 30 86 118 90 
+2 5 0 0 8 28 9 54 82 86 56 117 118 
+3 5 0 0 8 30 86 118 90 10 58 113 60 
+4 5 0 0 8 86 56 117 118 58 22 77 113 
+5 5 0 0 8 32 82 52 11 90 118 109 62 
+6 5 0 0 8 82 54 21 52 118 117 81 109 
+7 5 0 0 8 90 118 109 62 60 113 79 23 
+8 5 0 0 8 118 117 81 109 113 77 27 79 
+9 5 0 0 8 9 29 84 54 56 87 119 117 
+10 5 0 0 8 29 2 36 84 87 34 94 119 
+11 5 0 0 8 56 87 119 117 22 59 111 77 
+12 5 0 0 8 87 34 94 119 59 12 64 111 
+13 5 0 0 8 54 84 53 21 117 119 108 81 
+14 5 0 0 8 84 36 13 53 119 94 66 108 
+15 5 0 0 8 117 119 108 81 77 111 78 27 
+16 5 0 0 8 119 94 66 108 111 64 24 78 
+17 5 0 0 8 10 58 113 60 31 88 120 91 
+18 5 0 0 8 58 22 77 113 88 57 116 120 
+19 5 0 0 8 31 88 120 91 3 38 98 40 
+20 5 0 0 8 88 57 116 120 38 14 70 98 
+21 5 0 0 8 60 113 79 23 91 120 107 63 
+22 5 0 0 8 113 77 27 79 120 116 80 107 
+23 5 0 0 8 91 120 107 63 40 98 68 15 
+24 5 0 0 8 120 116 80 107 98 70 25 68 
+25 5 0 0 8 22 59 111 77 57 89 121 116 
+26 5 0 0 8 59 12 64 111 89 35 95 121 
+27 5 0 0 8 57 89 121 116 14 39 100 70 
+28 5 0 0 8 89 35 95 121 39 4 42 100 
+29 5 0 0 8 77 111 78 27 116 121 106 80 
+30 5 0 0 8 111 64 24 78 121 95 67 106 
+31 5 0 0 8 116 121 106 80 70 100 69 25 
+32 5 0 0 8 121 95 67 106 100 42 16 69 
+33 5 0 0 8 11 52 83 33 62 109 122 92 
+34 5 0 0 8 52 21 55 83 109 81 115 122 
+35 5 0 0 8 62 109 122 92 23 79 112 61 
+36 5 0 0 8 109 81 115 122 79 27 76 112 
+37 5 0 0 8 33 83 44 5 92 122 102 46 
+38 5 0 0 8 83 55 17 44 122 115 72 102 
+39 5 0 0 8 92 122 102 46 61 112 74 18 
+40 5 0 0 8 122 115 72 102 112 76 26 74 
+41 5 0 0 8 21 53 85 55 81 108 123 115 
+42 5 0 0 8 53 13 37 85 108 66 96 123 
+43 5 0 0 8 81 108 123 115 27 78 110 76 
+44 5 0 0 8 108 66 96 123 78 24 65 110 
+45 5 0 0 8 55 85 45 17 115 123 103 72 
+46 5 0 0 8 85 37 6 45 123 96 48 103 
+47 5 0 0 8 115 123 103 72 76 110 75 26 
+48 5 0 0 8 123 96 48 103 110 65 19 75 
+49 5 0 0 8 23 79 112 61 63 107 124 93 
+50 5 0 0 8 79 27 76 112 107 80 114 124 
+51 5 0 0 8 63 107 124 93 15 68 99 41 
+52 5 0 0 8 107 80 114 124 68 25 71 99 
+53 5 0 0 8 61 112 74 18 93 124 104 47 
+54 5 0 0 8 112 76 26 74 124 114 73 104 
+55 5 0 0 8 93 124 104 47 41 99 50 7 
+56 5 0 0 8 124 114 73 104 99 71 20 50 
+57 5 0 0 8 27 78 110 76 80 106 125 114 
+58 5 0 0 8 78 24 65 110 106 67 97 125 
+59 5 0 0 8 80 106 125 114 25 69 101 71 
+60 5 0 0 8 106 67 97 125 69 16 43 101 
+61 5 0 0 8 76 110 75 26 114 125 105 73 
+62 5 0 0 8 110 65 19 75 125 97 49 105 
+63 5 0 0 8 114 125 105 73 71 101 51 20 
+64 5 0 0 8 125 97 49 105 101 43 8 51 
+$ENDELM
diff --git a/tests/opencascade/plane_intersection_01.cc b/tests/opencascade/plane_intersection_01.cc
new file mode 100644 (file)
index 0000000..726a38e
--- /dev/null
@@ -0,0 +1,86 @@
+//-----------------------------------------------------------
+//
+//    Copyright (C) 2014 by the deal.II authors 
+//
+//    This file is subject to LGPL and may not be distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//-----------------------------------------------------------
+
+// Create a sphere, intersect it with a plane, and use the resulting
+// object as a boundary.
+
+#include "../tests.h"
+
+#include <deal.II/grid/occ_utilities.h>
+#include <deal.II/grid/occ_boundary_lib.h>
+
+#include <fstream>
+#include <base/logstream.h>
+#include <grid/tria.h>
+#include <grid/tria_accessor.h>
+#include <grid/grid_out.h>
+#include <grid/tria_iterator.h>
+#include <grid/grid_generator.h>
+#include <grid/tria_boundary_lib.h>
+
+#include <gp_Pnt.hxx>
+#include <gp_Dir.hxx>
+#include <gp_Ax2.hxx>
+#include <GC_MakeCircle.hxx>
+#include <BRepBuilderAPI_MakeEdge.hxx>
+#include <BRepBuilderAPI_MakeFace.hxx>
+#include <BRepBuilderAPI_MakeWire.hxx>
+#include <BRepPrimAPI_MakeSphere.hxx>
+#include <TopExp_Explorer.hxx>
+#include <TopoDS_Wire.hxx>
+#include <TopoDS_Edge.hxx>
+#include <TopoDS_Face.hxx>
+#include <TopoDS.hxx>
+using namespace OpenCASCADE;
+
+int main()
+{
+  std::ofstream logfile("output");
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+
+  // Create a sphere
+  gp_Pnt center(.5,.5,.5);
+  Standard_Real radius(Point<3>().distance(Pnt(center)));
+
+  TopoDS_Face face = BRepPrimAPI_MakeSphere(center, radius);
+  
+  // Intersect the sphere with the plane z = 0
+  TopoDS_Shape intersect = intersect_plane(face, 0, 0, 1, 0);
+  // This shape should be the standard circle. Only one edge should be
+  // in intersect: our edge.
+  TopExp_Explorer exp(intersect, TopAbs_EDGE);
+  TopoDS_Edge edge = TopoDS::Edge(exp.Current());
+  
+  ArclengthProjectionLineManifold<2,3> boundary_line(edge);
+
+  // The unit square.
+  Triangulation<2,3> tria;
+  GridGenerator::hyper_cube(tria);
+
+  // Set the exterior boundary
+  tria.set_manifold(1, boundary_line);
+
+  // This is here to ignore the points created in the interior of the
+  // face.
+  tria.begin()->set_all_manifold_ids(1);
+  tria.begin()->set_manifold_id(0);
+
+  // We refine twice, and expect the outer points to end up on a
+  // smooth curve interpolating the square vertices.
+  tria.refine_global(2);
+
+
+  // You can open the generated file with gmsh.
+  GridOut gridout;
+  gridout.write_msh (tria, logfile);
+  return 0;
+}
diff --git a/tests/opencascade/plane_intersection_01.output b/tests/opencascade/plane_intersection_01.output
new file mode 100644 (file)
index 0000000..6203ac0
--- /dev/null
@@ -0,0 +1,48 @@
+
+$NOD
+25
+1  0.00000 0.00000 0.00000
+2  1.00000 0.00000 0.00000
+3  0.00000 1.00000 0.00000
+4  1.00000 1.00000 0.00000
+5  0.500000 -0.207107 0.00000
+6  -0.207107 0.500000 0.00000
+7  1.20711 0.500000 0.00000
+8  0.500000 1.20711 0.00000
+9  0.500000 0.500000 0.00000
+10  0.229402 -0.153281 0.00000
+11  0.770598 -0.153281 0.00000
+12  -0.153281 0.229402 0.00000
+13  -0.153281 0.770598 0.00000
+14  1.15328 0.229402 0.00000
+15  1.15328 0.770598 0.00000
+16  0.229402 1.15328 0.00000
+17  0.770598 1.15328 0.00000
+18  0.500000 0.146447 0.00000
+19  0.500000 0.853553 0.00000
+20  0.146447 0.500000 0.00000
+21  0.853553 0.500000 0.00000
+22  0.189433 0.189433 0.00000
+23  0.810567 0.189433 0.00000
+24  0.189433 0.810567 0.00000
+25  0.810567 0.810567 0.00000
+$ENDNOD
+$ELM
+16
+1 3 0 0 4 1 10 22 12 
+2 3 0 0 4 10 5 18 22 
+3 3 0 0 4 12 22 20 6 
+4 3 0 0 4 22 18 9 20 
+5 3 0 0 4 5 11 23 18 
+6 3 0 0 4 11 2 14 23 
+7 3 0 0 4 18 23 21 9 
+8 3 0 0 4 23 14 7 21 
+9 3 0 0 4 6 20 24 13 
+10 3 0 0 4 20 9 19 24 
+11 3 0 0 4 13 24 16 3 
+12 3 0 0 4 24 19 8 16 
+13 3 0 0 4 9 21 25 19 
+14 3 0 0 4 21 7 15 25 
+15 3 0 0 4 19 25 17 8 
+16 3 0 0 4 25 15 4 17 
+$ENDELM

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.