From: Luca Heltai Date: Tue, 5 Dec 2017 12:07:52 +0000 (+0100) Subject: OpenCASCADE::create_curves_from_triangulation_boundary X-Git-Tag: v9.0.0-rc1~687^2~2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=e6b0506f13ac3a376968bafc672fd6c2d1ea835e;p=dealii.git OpenCASCADE::create_curves_from_triangulation_boundary --- diff --git a/doc/news/changes/minor/20171205LucaHeltai b/doc/news/changes/minor/20171205LucaHeltai new file mode 100644 index 0000000000..39468e6536 --- /dev/null +++ b/doc/news/changes/minor/20171205LucaHeltai @@ -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. +
+(Dirk Peschka, Luca Heltai, 2017/12/05) diff --git a/include/deal.II/opencascade/utilities.h b/include/deal.II/opencascade/utilities.h index ba788d1cf5..e51525d215 100644 --- a/include/deal.II/opencascade/utilities.h +++ b/include/deal.II/opencascade/utilities.h @@ -23,6 +23,7 @@ #include #include +#include #include @@ -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 + std::vector 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 diff --git a/source/opencascade/utilities.cc b/source/opencascade/utilities.cc index 821a7ab564..77429dbd42 100644 --- a/source/opencascade/utilities.cc +++ b/source/opencascade/utilities.cc @@ -462,6 +462,112 @@ namespace OpenCASCADE return out_shape; } + + + template + std::vector 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 > vert_to_faces; + std::map > face_to_verts; + std::map visited_faces; + std::map > vert_to_point; + + unsigned int face_index; + + for (auto cell: triangulation.active_cell_iterators()) + for (unsigned int f=0; f::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,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 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 > 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, 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 index 0000000000..6e6ad0b9f8 --- /dev/null +++ b/tests/opencascade/interpolation_curve_boundary_01.cc @@ -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 + +#include +#include +#include +#include +#include +#include +#include + +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 +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 index 0000000000..fca148838a --- /dev/null +++ b/tests/opencascade/interpolation_curve_boundary_01.output @@ -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 +