]> https://gitweb.dealii.org/ - dealii.git/commitdiff
OpenCASCADE::create_curves_from_triangulation_boundary
authorLuca Heltai <luca.heltai@sissa.it>
Tue, 5 Dec 2017 12:07:52 +0000 (13:07 +0100)
committerLuca Heltai <luca.heltai@sissa.it>
Tue, 5 Dec 2017 12:07:52 +0000 (13:07 +0100)
doc/news/changes/minor/20171205LucaHeltai [new file with mode: 0644]
include/deal.II/opencascade/utilities.h
source/opencascade/utilities.cc
tests/opencascade/interpolation_curve_boundary_01.cc [new file with mode: 0644]
tests/opencascade/interpolation_curve_boundary_01.output [new file with mode: 0644]

diff --git a/doc/news/changes/minor/20171205LucaHeltai b/doc/news/changes/minor/20171205LucaHeltai
new file mode 100644 (file)
index 0000000..39468e6
--- /dev/null
@@ -0,0 +1,5 @@
+New: Created a new method OpenCASCADE::Utilities::create_curves_from_triangulation_boundary
+that smoothly interpolates the boundary of two dimensional triangulations into a vector of 
+OpenCASCADE TopoDS_Edge closed curves, representing the connected components of the boundary. 
+<br> 
+(Dirk Peschka, Luca Heltai, 2017/12/05)
index ba788d1cf591d0f03c34e22d0b2a43fc021ccd08..e51525d215356b0b08c96ca0aced0d93bc326504 100644 (file)
@@ -23,6 +23,7 @@
 
 #include <deal.II/grid/tria.h>
 #include <deal.II/base/point.h>
+#include <deal.II/fe/mapping_q1.h>
 
 #include <string>
 
@@ -88,8 +89,10 @@ DEAL_II_NAMESPACE_OPEN
  * TopoDS_Shapes of one type or another, and provide interfaces to common
  * deal.II objects, like Triangulation, Manifold, and so on.
  *
- * Notice that these tools are only useful when spacedim is equal to three,
- * since OpenCASCADE only operates in three-dimensional mode.
+ * Notice that most of these tools are only useful when spacedim is equal to
+ * three, since OpenCASCADE only operates in three-dimensional mode. In some
+ * cases they can be used in two dimensions as well, and the third dimension
+ * will be set to zero.
  *
  * @author Luca Heltai, Andrea Mola, 2011--2014.
  */
@@ -219,6 +222,28 @@ namespace OpenCASCADE
   void create_triangulation(const TopoDS_Face &face,
                             Triangulation<2,3> &tria);
 
+
+  /**
+   * Given a Triangulation and an optional Mapping, create a vector of smooth
+   * curves that interpolate the connected parts of the boundary vertices of
+   * the Triangulation and return them as a vector of TopoDS_Edge.
+   *
+   * The returned curves are ordered with respect to the indices of the faces
+   * that make up the triangulation boundary, i.e., the first curve is the one
+   * extracted starting from the face with the lowest index, and so on.
+   *
+   * @param[input] triangulation Input triangulation
+   * @param[input] mapping Optional input mapping
+   * @return An std::vector of TopoDS_Edge objects representing the smooth
+   *  interpolation of the boundary of the `triangulation`
+   *
+   * @author Dirk Peschka, Luca Heltai, 2017.
+   */
+  template <int spacedim>
+  std::vector<TopoDS_Edge> create_curves_from_triangulation_boundary
+  (const Triangulation<2,spacedim> &triangulation,
+   const Mapping<2,spacedim> &mapping=StaticMappingQ1<2,spacedim>::mapping);
+
   /**
    * 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
index 821a7ab564be10d4b3fa56a7aaff72006b9263f9..77429dbd427b96509592d3737bb57a6e1bc64c7a 100644 (file)
@@ -462,6 +462,112 @@ namespace OpenCASCADE
     return out_shape;
   }
 
+
+
+  template<int spacedim>
+  std::vector<TopoDS_Edge> create_curves_from_triangulation_boundary
+  (const Triangulation<2, spacedim> &triangulation,
+   const Mapping<2, spacedim> &mapping)
+
+  {
+    // store maps from global vertex index to pairs of global face   indices
+    //        and from global face   index to pairs of global vertex indices
+    std::map<unsigned int, std::pair<unsigned int,unsigned int> > vert_to_faces;
+    std::map<unsigned int, std::pair<unsigned int,unsigned int> > face_to_verts;
+    std::map<unsigned int, bool>                                  visited_faces;
+    std::map<unsigned int, Point<spacedim> >                      vert_to_point;
+
+    unsigned int face_index;
+
+    for (auto cell: triangulation.active_cell_iterators())
+      for (unsigned int f=0; f<GeometryInfo<2>::faces_per_cell; ++f)
+        if (cell->face(f)->at_boundary())
+          {
+            // get global face and vertex indices
+            face_index = cell->face(f)->index();
+            unsigned int v0 = cell->face(f)->vertex_index(0);
+            unsigned int v1 = cell->face(f)->vertex_index(1);
+            face_to_verts[face_index].first  = v0;
+            face_to_verts[face_index].second = v1;
+            visited_faces[face_index] = false;
+
+            // extract mapped vertex locations
+            std::array<Point<spacedim>,GeometryInfo<2>::vertices_per_cell> verts = mapping.get_vertices(cell);
+            vert_to_point[v0]=verts[GeometryInfo<2>::face_to_cell_vertices(f,0,true,false,false)];
+            vert_to_point[v1]=verts[GeometryInfo<2>::face_to_cell_vertices(f,1,true,false,false)];
+
+            // distribute indices into maps
+            if (vert_to_faces.find(v0) == vert_to_faces.end())
+              {
+                vert_to_faces[v0].first  = face_index;
+              }
+            else
+              {
+                vert_to_faces[v0].second = face_index;
+              }
+            if (vert_to_faces.find(v1) == vert_to_faces.end())
+              {
+                vert_to_faces[v1].first  = face_index;
+              }
+            else
+              {
+                vert_to_faces[v1].second = face_index;
+              }
+          }
+
+    // run through maps in an orderly fashion, i.e., through the
+    // boundary in one cycle and add points to pointlist.
+    std::vector<TopoDS_Edge> interpolation_curves;
+    bool finished = (face_to_verts.size() == 0);
+    face_index = finished ? 0 : face_to_verts.begin()->first;
+
+    while (finished == false)
+      {
+        const unsigned int start_point_index = face_to_verts[face_index].first;
+        unsigned int point_index             = start_point_index;
+
+        // point_index and face_index always run together
+        std::vector<Point<3> > pointlist;
+        do
+          {
+            visited_faces[face_index]            = true;
+            auto current_point = vert_to_point[point_index];
+            if (spacedim==2)
+              pointlist.push_back(Point<3>(current_point[0],current_point[1],0));
+            else
+              pointlist.push_back(Point<3>(current_point[0],current_point[1],current_point[2]));
+
+            // Get next point
+            if (face_to_verts[face_index].first != point_index )
+              point_index = face_to_verts[face_index].first;
+            else
+              point_index = face_to_verts[face_index].second;
+
+            // Get next face
+            if (vert_to_faces[point_index].first != face_index)
+              face_index = vert_to_faces[point_index].first;
+            else
+              face_index = vert_to_faces[point_index].second;
+          }
+        while (point_index != start_point_index);
+
+        interpolation_curves.push_back(interpolation_curve(pointlist, Tensor<1,3>(), true));
+        pointlist.clear();
+
+        finished = true;
+        for (const auto &f : visited_faces)
+          if (f.second == false)
+            {
+              face_index = f.first;
+              finished = false;
+              break;
+            }
+      }
+    return interpolation_curves;
+  }
+
+
+
   std::tuple<Point<3>, TopoDS_Shape, double, double>
   project_point_and_pull_back(const TopoDS_Shape &in_shape,
                               const Point<3> &origin,
diff --git a/tests/opencascade/interpolation_curve_boundary_01.cc b/tests/opencascade/interpolation_curve_boundary_01.cc
new file mode 100644 (file)
index 0000000..6e6ad0b
--- /dev/null
@@ -0,0 +1,106 @@
+//-----------------------------------------------------------
+//
+//    Copyright (C) 2014 - 2015 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/opencascade/utilities.h>
+
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/grid_out.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria_boundary_lib.h>
+#include <deal.II/grid/grid_tools.h>
+
+using namespace OpenCASCADE;
+
+void remove_iges_header(const std::string &in_filename,
+                        const std::string &out_filename)
+{
+  std::ifstream in(in_filename);
+  std::ofstream out(out_filename);
+  std::string line;
+  unsigned int counter = 5;
+  while (counter--) std::getline(in,line);
+  while (std::getline(in,line))
+    out << line << std::endl;
+  in.close();
+  out.close();
+}
+
+template<int spacedim>
+void test()
+{
+  deallog << "Testing <2," << spacedim << ">" << std::endl;
+
+  Triangulation<2,spacedim> tria1;
+  Triangulation<2,spacedim> tria2;
+  Triangulation<2,spacedim> tria3;
+
+  GridGenerator::hyper_cube(tria1, 0, 1);
+  GridGenerator::hyper_cube(tria2, 2, 3);
+
+  GridGenerator::merge_triangulations(tria1, tria2, tria3);
+
+  tria3.refine_global(1);
+
+  auto v1 = create_curves_from_triangulation_boundary(tria1);
+  auto v3 = create_curves_from_triangulation_boundary(tria3);
+
+  if (spacedim == 3)
+    {
+      Triangulation<2,3> tria4;
+      GridGenerator::hyper_sphere(tria4);
+      auto v4 = create_curves_from_triangulation_boundary(tria4);
+
+      if (v4.size() != 0)
+        deallog << "Not OK!: size is " << v4.size() << std::endl;
+      else
+        deallog << "OK" << std::endl;
+    }
+
+  if (v1.size() != 1)
+    deallog << "Not OK: size is " << v1.size() << std::endl;
+  else
+    deallog << "OK" << std::endl;
+
+  if (v3.size() != 2)
+    deallog << "Not OK!: size is " << v3.size() << std::endl;
+  else
+    deallog << "OK" << std::endl;
+
+  write_IGES(v1[0], "edge1.iges");
+  write_IGES(v3[0], "edge30.iges");
+  write_IGES(v3[1], "edge31.iges");
+
+  remove_iges_header("edge1.iges", "edge1_noheader.iges");
+  remove_iges_header("edge30.iges", "edge30_noheader.iges");
+  remove_iges_header("edge31.iges", "edge31_noheader.iges");
+
+  cat_file("edge1_noheader.iges");
+  cat_file("edge30_noheader.iges");
+  cat_file("edge31_noheader.iges");
+}
+
+int main ()
+{
+  initlog();
+
+  test<2>();
+  test<3>();
+
+  return 0;
+}
+
diff --git a/tests/opencascade/interpolation_curve_boundary_01.output b/tests/opencascade/interpolation_curve_boundary_01.output
new file mode 100644 (file)
index 0000000..fca1488
--- /dev/null
@@ -0,0 +1,72 @@
+
+DEAL::Testing <2,2>
+DEAL::OK
+DEAL::OK
+     126       1       0       0       0       0       0       000000000D0000001
+     126       0       0       4       0                               0D0000002
+126,6,3,1,1,1,1,-2.,-1.,0.,0.,1.,2.,3.,4.,4.,5.,6.,1.,1.,1.,1.,  0000001P0000001
+1.,1.,1.,-0.5,0.388888889,0.,0.5,-0.388888889,0.,1.142857143,    0000001P0000002
+-0.19047619,0.,1.25,1.25,0.,-0.142857143,1.19047619,0.,-0.5,     0000001P0000003
+0.388888889,0.,0.5,-0.388888889,0.,0.,4.,0.,0.,1.;               0000001P0000004
+S      1G      4D      2P      4                                        T0000001
+
+     126       1       0       0       0       0       0       000000000D0000001
+     126       0       0       8       0                               0D0000002
+126,10,3,1,1,1,1,-1.,-0.5,0.,0.,0.5,1.,1.5,2.,2.5,3.,3.5,4.,4.,  0000001P0000001
+4.5,5.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,-0.111111111,            0000001P0000002
+-5.555555556E-02,0.,0.111111111,5.555555556E-02,0.,0.481590574,  0000001P0000003
+1.251840943E-02,0.,1.147766323,-0.127147766,0.,0.927344134,      0000001P0000004
+0.496072656,0.,1.142857143,1.142857143,0.,0.501227295,           0000001P0000005
+0.932498773,0.,-0.147766323,1.127147766,0.,8.983799705E-02,      0000001P0000006
+0.558910162,0.,-0.111111111,-5.555555556E-02,0.,0.111111111,     0000001P0000007
+5.555555556E-02,0.,0.,4.,0.,0.,1.;                               0000001P0000008
+S      1G      4D      2P      8                                        T0000001
+
+     126       1       0       0       0       0       0       000000000D0000001
+     126       0       0       8       0                               0D0000002
+126,10,3,1,1,1,1,-1.,-0.5,0.,0.,0.5,1.,1.5,2.,2.5,3.,3.5,4.,4.,  0000001P0000001
+4.5,5.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.888888889,1.944444444, 0000001P0000002
+0.,2.111111111,2.055555556,0.,2.481590574,2.012518409,0.,        0000001P0000003
+3.147766323,1.872852234,0.,2.927344134,2.496072656,0.,           0000001P0000004
+3.142857143,3.142857143,0.,2.501227295,2.932498773,0.,           0000001P0000005
+1.852233677,3.127147766,0.,2.089837997,2.558910162,0.,           0000001P0000006
+1.888888889,1.944444444,0.,2.111111111,2.055555556,0.,0.,4.,0.,  0000001P0000007
+0.,1.;                                                           0000001P0000008
+S      1G      4D      2P      8                                        T0000001
+
+DEAL::Testing <2,3>
+DEAL::OK
+DEAL::OK
+DEAL::OK
+     126       1       0       0       0       0       0       000000000D0000001
+     126       0       0       4       0                               0D0000002
+126,6,3,1,1,1,1,-2.,-1.,0.,0.,1.,2.,3.,4.,4.,5.,6.,1.,1.,1.,1.,  0000001P0000001
+1.,1.,1.,-0.5,0.388888889,0.,0.5,-0.388888889,0.,1.142857143,    0000001P0000002
+-0.19047619,0.,1.25,1.25,0.,-0.142857143,1.19047619,0.,-0.5,     0000001P0000003
+0.388888889,0.,0.5,-0.388888889,0.,0.,4.,0.,0.,1.;               0000001P0000004
+S      1G      4D      2P      4                                        T0000001
+
+     126       1       0       0       0       0       0       000000000D0000001
+     126       0       0       8       0                               0D0000002
+126,10,3,1,1,1,1,-1.,-0.5,0.,0.,0.5,1.,1.5,2.,2.5,3.,3.5,4.,4.,  0000001P0000001
+4.5,5.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,-0.111111111,            0000001P0000002
+-5.555555556E-02,0.,0.111111111,5.555555556E-02,0.,0.481590574,  0000001P0000003
+1.251840943E-02,0.,1.147766323,-0.127147766,0.,0.927344134,      0000001P0000004
+0.496072656,0.,1.142857143,1.142857143,0.,0.501227295,           0000001P0000005
+0.932498773,0.,-0.147766323,1.127147766,0.,8.983799705E-02,      0000001P0000006
+0.558910162,0.,-0.111111111,-5.555555556E-02,0.,0.111111111,     0000001P0000007
+5.555555556E-02,0.,0.,4.,0.,0.,1.;                               0000001P0000008
+S      1G      4D      2P      8                                        T0000001
+
+     126       1       0       0       0       0       0       000000000D0000001
+     126       0       0       8       0                               0D0000002
+126,10,3,1,1,1,1,-1.,-0.5,0.,0.,0.5,1.,1.5,2.,2.5,3.,3.5,4.,4.,  0000001P0000001
+4.5,5.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.888888889,1.944444444, 0000001P0000002
+0.,2.111111111,2.055555556,0.,2.481590574,2.012518409,0.,        0000001P0000003
+3.147766323,1.872852234,0.,2.927344134,2.496072656,0.,           0000001P0000004
+3.142857143,3.142857143,0.,2.501227295,2.932498773,0.,           0000001P0000005
+1.852233677,3.127147766,0.,2.089837997,2.558910162,0.,           0000001P0000006
+1.888888889,1.944444444,0.,2.111111111,2.055555556,0.,0.,4.,0.,  0000001P0000007
+0.,1.;                                                           0000001P0000008
+S      1G      4D      2P      8                                        T0000001
+

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.