From 73d02cb86ea2a677d330796b99de43ac2c0a4721 Mon Sep 17 00:00:00 2001 From: Luca Heltai Date: Sat, 9 Dec 2017 11:59:19 +0100 Subject: [PATCH] OpenCASCADE classes and utilities now work also in spacedim=2 --- doc/news/changes/minor/20171208LucaHeltai | 3 + include/deal.II/opencascade/utilities.h | 82 ++++++---- source/opencascade/CMakeLists.txt | 1 + source/opencascade/boundary_lib.cc | 55 ++++--- source/opencascade/boundary_lib.inst.in | 6 + source/opencascade/utilities.cc | 144 +++++++++-------- source/opencascade/utilities.inst.in | 62 +++++++ tests/opencascade/arclength_boundary_02.cc | 56 +++++++ .../opencascade/arclength_boundary_02.output | 151 ++++++++++++++++++ 9 files changed, 437 insertions(+), 123 deletions(-) create mode 100644 doc/news/changes/minor/20171208LucaHeltai create mode 100644 source/opencascade/utilities.inst.in create mode 100644 tests/opencascade/arclength_boundary_02.cc create mode 100644 tests/opencascade/arclength_boundary_02.output diff --git a/doc/news/changes/minor/20171208LucaHeltai b/doc/news/changes/minor/20171208LucaHeltai new file mode 100644 index 0000000000..be924f86ac --- /dev/null +++ b/doc/news/changes/minor/20171208LucaHeltai @@ -0,0 +1,3 @@ +New: Now opencascade also works for spacedim == 2. +
+(Luca Heltai, 2017/12/08) diff --git a/include/deal.II/opencascade/utilities.h b/include/deal.II/opencascade/utilities.h index 8d16be67e2..2fbe26d691 100644 --- a/include/deal.II/opencascade/utilities.h +++ b/include/deal.II/opencascade/utilities.h @@ -94,7 +94,11 @@ DEAL_II_NAMESPACE_OPEN * 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. + * If you whish to use these tools when the dimension of the space is two, then + * make sure your cad files are actually flat and that all z coordinates are equal + * to zero, as otherwise you will get many exceptions. + * + * @author Luca Heltai, Andrea Mola, 2011--2017. */ namespace OpenCASCADE { @@ -179,9 +183,9 @@ namespace OpenCASCADE 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 (which is of - * type TopoDS_Edge). The points are reordered internally according to their + * Creates a smooth BSpline curve passing through the points in the assigned + * vector, and store it in the returned TopoDS_Shape (which is of type + * TopoDS_Edge). The points are reordered internally according to their * scalar product with the direction, if direction is different from zero, * otherwise they are used as passed. Notice that this function changes the * input points if required by the algorithm. @@ -195,8 +199,9 @@ namespace OpenCASCADE * points. If the algorithm fails in generating such a curve, an exception * is thrown. */ - TopoDS_Edge interpolation_curve(std::vector > &curve_points, - const Tensor<1,3> &direction=Tensor<1,3>(), + template + TopoDS_Edge interpolation_curve(std::vector > &curve_points, + const Tensor<1, dim> &direction=Tensor<1,dim>(), const bool closed=false, const double tolerance=1e-7); @@ -211,16 +216,20 @@ namespace OpenCASCADE std::vector &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 + * Create a triangulation from a single face. This class extracts the first u + * and v parameter of the parametric surface making up this face, and creates + * a Triangulation<2,spacedim> 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. + * + * If you call this function with a Triangulation<2,2>, make sure that the + * input face has all z coordinates set to zero, or you'll get an exception. */ + template void create_triangulation(const TopoDS_Face &face, - Triangulation<2,3> &tria); + Triangulation<2,spacedim> &tria); /** @@ -274,9 +283,10 @@ namespace OpenCASCADE * the u coordinate and the v coordinate (which is different from zero only * if the resulting shape is a face). */ - std::tuple, TopoDS_Shape, double, double> + template + std::tuple, TopoDS_Shape, double, double> project_point_and_pull_back(const TopoDS_Shape &in_shape, - const Point<3> &origin, + const Point &origin, const double tolerance=1e-7); /** @@ -285,9 +295,10 @@ namespace OpenCASCADE * are iterated, faces first, then edges, and the returned point is the * closest one to the @p in_shape, regardless of its type. */ - Point<3> closest_point(const TopoDS_Shape &in_shape, - const Point<3> &origin, - const double tolerance=1e-7); + template + Point closest_point(const TopoDS_Shape &in_shape, + const Point &origin, + const double tolerance=1e-7); /** * Given an elementary shape @p in_shape and the reference coordinates @@ -297,9 +308,10 @@ namespace OpenCASCADE * used as input to this function. If this is not the case, an Exception is * thrown. */ - Point<3> push_forward(const TopoDS_Shape &in_shape, - const double u, - const double v); + template + Point push_forward(const TopoDS_Shape &in_shape, + const double u, + const double v); /** @@ -334,22 +346,29 @@ namespace OpenCASCADE * * The optional @p tolerance parameter is used to compute distances. */ - Point<3> line_intersection(const TopoDS_Shape &in_shape, - const Point<3> &origin, - const Tensor<1,3> &direction, - const double tolerance=1e-7); + template + Point line_intersection(const TopoDS_Shape &in_shape, + const Point &origin, + const Tensor<1,dim> &direction, + const double tolerance=1e-7); /** - * Convert OpenCASCADE point into a Point<3>. + * Convert OpenCASCADE point into a Point. + * + * The tolerance argument is used to check if the non used components of the + * OpenCASCADE point are close to zero. If this is not the case, an assertion + * is thrown in debug mode. */ - Point<3> point(const gp_Pnt &p); + template + Point point(const gp_Pnt &p, const double &tolerance=1e-10); /** * Convert Point<3> into OpenCASCADE point. */ - gp_Pnt point(const Point<3> &p); + template + gp_Pnt point(const Point &p); /** @@ -358,9 +377,10 @@ namespace OpenCASCADE * optional parameter is used as a relative tolerance when comparing * objects. */ - bool point_compare(const Point<3> &p1, - const Point<3> &p2, - const Tensor<1,3> &direction = Tensor<1,3>(), + template + bool point_compare(const Point &p1, + const Point &p2, + const Tensor<1,dim> &direction = Tensor<1,dim>(), const double tolerance = 1e-10); @@ -368,16 +388,18 @@ namespace OpenCASCADE * Exception thrown when the point specified as argument does not lie * between @p tolerance from the given TopoDS_Shape. */ + template DeclException1 (ExcPointNotOnManifold, - Point<3>, + Point, <<"The point [ "< DeclException1 (ExcProjectionFailed, - Point<3>, + Point, <<"Projection of point [ "<< arg1 << " ] failed."); diff --git a/source/opencascade/CMakeLists.txt b/source/opencascade/CMakeLists.txt index 9a13042353..0d9905db34 100644 --- a/source/opencascade/CMakeLists.txt +++ b/source/opencascade/CMakeLists.txt @@ -22,6 +22,7 @@ SET(_src SET(_inst boundary_lib.inst.in + utilities.inst.in ) FILE(GLOB _header diff --git a/source/opencascade/boundary_lib.cc b/source/opencascade/boundary_lib.cc index 8125c5bfd2..1545c0f4c3 100644 --- a/source/opencascade/boundary_lib.cc +++ b/source/opencascade/boundary_lib.cc @@ -99,7 +99,7 @@ namespace OpenCASCADE Assert(closest_point(sh, surrounding_points[i], tolerance) .distance(surrounding_points[i]) < std::max(tolerance*surrounding_points[i].norm(), tolerance), - ExcPointNotOnManifold(surrounding_points[i])); + ExcPointNotOnManifold(surrounding_points[i])); #endif return closest_point(sh, candidate,tolerance); } @@ -129,7 +129,7 @@ namespace OpenCASCADE Assert(closest_point(sh, surrounding_points[i],tolerance) .distance(surrounding_points[i]) < std::max(tolerance*surrounding_points[i].norm(), tolerance), - ExcPointNotOnManifold(surrounding_points[i])); + ExcPointNotOnManifold(surrounding_points[i])); #endif return line_intersection(sh, candidate, direction, tolerance); } @@ -162,7 +162,7 @@ namespace OpenCASCADE Assert(closest_point(sh, surrounding_points[i], tolerance) .distance(surrounding_points[i]) < std::max(tolerance*surrounding_points[i].norm(), tolerance), - ExcPointNotOnManifold(surrounding_points[i])); + ExcPointNotOnManifold(surrounding_points[i])); } #endif @@ -248,9 +248,7 @@ namespace OpenCASCADE curve(curve_adaptor(sh)), tolerance(tolerance), length(shape_length(sh)) - { - Assert(spacedim == 3, ExcNotImplemented()); - } + {} template @@ -261,7 +259,7 @@ namespace OpenCASCADE ShapeAnalysis_Curve curve_analysis; gp_Pnt proj; const double dist = curve_analysis.Project(curve->GetCurve(), point(space_point), tolerance, proj, t, true); - Assert(dist < tolerance*length, ExcPointNotOnManifold(space_point)); + Assert(dist < tolerance*length, ExcPointNotOnManifold(space_point)); (void)dist; // Silence compiler warning in Release mode. return Point<1>(GCPnts_AbscissaPoint::Length(curve->GetCurve(),curve->GetCurve().FirstParameter(),t)); } @@ -274,7 +272,7 @@ namespace OpenCASCADE { GCPnts_AbscissaPoint AP(curve->GetCurve(), chart_point[0], curve->GetCurve().FirstParameter()); gp_Pnt P = curve->GetCurve().Value(AP.Parameter()); - return point(P); + return point(P); } template @@ -286,10 +284,9 @@ namespace OpenCASCADE tolerance(tolerance) {} - template <> - Point<2> - NURBSPatchManifold<2, 3>:: - pull_back(const Point<3> &space_point) const + template Point<2> + NURBSPatchManifold:: + pull_back(const Point &space_point) const { Handle(Geom_Surface) SurfToProj = BRep_Tool::Surface(face); @@ -302,20 +299,20 @@ namespace OpenCASCADE return Point<2>(u,v); } - template <> - Point<3> - NURBSPatchManifold<2, 3>:: + template + Point + NURBSPatchManifold:: push_forward(const Point<2> &chart_point) const { - return ::dealii::OpenCASCADE::push_forward(face, chart_point[0], chart_point[1]); + return ::dealii::OpenCASCADE::push_forward(face, chart_point[0], chart_point[1]); } - template <> - DerivativeForm<1,2,3> - NURBSPatchManifold<2, 3>:: + template + DerivativeForm<1,2,spacedim> + NURBSPatchManifold:: push_forward_gradient(const Point<2> &chart_point) const { - DerivativeForm<1,2,3> DX; + DerivativeForm<1,2,spacedim> DX; Handle(Geom_Surface) surf = BRep_Tool::Surface(face); gp_Pnt q; @@ -324,17 +321,24 @@ namespace OpenCASCADE DX[0][0] = Du.X(); DX[1][0] = Du.Y(); - DX[2][0] = Du.Z(); + if (spacedim>2) + DX[2][0] = Du.Z(); + else + Assert(std::abs(Du.Z()) < tolerance, + ExcMessage("Expecting derivative along Z to be zero! Bailing out.")); DX[0][1] = Dv.X(); DX[1][1] = Dv.Y(); - DX[2][1] = Dv.Z(); - + if (spacedim>2) + DX[2][1] = Dv.Z(); + else + Assert(std::abs(Dv.Z()) < tolerance, + ExcMessage("Expecting derivative along Z to be zero! Bailing out.")); return DX; } - template <> + template std::tuple - NURBSPatchManifold<2, 3>:: + NURBSPatchManifold:: get_uv_bounds() const { Standard_Real umin, umax, vmin, vmax; @@ -343,7 +347,6 @@ namespace OpenCASCADE } // Explicit instantiations - template class NURBSPatchManifold<2, 3 >; #include "boundary_lib.inst" } // end namespace OpenCASCADE diff --git a/source/opencascade/boundary_lib.inst.in b/source/opencascade/boundary_lib.inst.in index a43f51b688..6514422184 100644 --- a/source/opencascade/boundary_lib.inst.in +++ b/source/opencascade/boundary_lib.inst.in @@ -21,6 +21,12 @@ for (deal_II_dimension : DIMENSIONS) template class DirectionalProjectionBoundary; template class NormalToMeshProjectionBoundary; template class ArclengthProjectionLineManifold; + template class NURBSPatchManifold; +#if deal_II_dimension <= 2 + template class DirectionalProjectionBoundary; + template class ArclengthProjectionLineManifold; + template class NURBSPatchManifold; +#endif } diff --git a/source/opencascade/utilities.cc b/source/opencascade/utilities.cc index eefb3b44f3..bd83e5f571 100644 --- a/source/opencascade/utilities.cc +++ b/source/opencascade/utilities.cc @@ -165,27 +165,52 @@ namespace OpenCASCADE } } - gp_Pnt point(const Point<3> &p) + template + gp_Pnt point(const Point &p) { - return gp_Pnt(p(0), p(1), p(2)); + switch (spacedim) + { + case 1: + return gp_Pnt(p[0], 0, 0); + case 2: + return gp_Pnt(p[0], p[1], 0); + case 3: + return gp_Pnt(p[0], p[1], p[2]); + } } - - Point<3> point(const gp_Pnt &p) + template + Point point(const gp_Pnt &p, const double &tolerance) { - return Point<3>(p.X(), p.Y(), p.Z()); + (void) tolerance; + switch (spacedim) + { + case 1: + Assert(std::abs(p.Y()) <= tolerance, + ExcMessage("Cannot convert OpenCASCADE point to 1d if p.Y() != 0.")); + Assert(std::abs(p.Z()) <= tolerance, + ExcMessage("Cannot convert OpenCASCADE point to 1d if p.Z() != 0.")); + return Point(p.X()); + case 2: + Assert(std::abs(p.Z()) <= tolerance, + ExcMessage("Cannot convert OpenCASCADE point to 2d if p.Z() != 0.")); + return Point(p.X(), p.Y()); + case 3: + return Point(p.X(), p.Y(), p.Z()); + } } - bool point_compare(const Point<3> &p1, - const Point<3> &p2, - const Tensor<1,3> &direction, + template + bool point_compare(const Point &p1, + const Point &p2, + const Tensor<1,dim> &direction, const double tolerance) { const double rel_tol=std::max(tolerance, std::max(p1.norm(), p2.norm())*tolerance); if (direction.norm() > 0.0) return (p1*direction < p2*direction-rel_tol); else - for (int d=2; d>=0; --d) + for (int d=dim; d>=0; --d) if (p1[d] < p2[d]-rel_tol) return true; else if (p2[d] < p1[d]-rel_tol) @@ -387,16 +412,16 @@ namespace OpenCASCADE return out_shape; } - - Point<3> line_intersection(const TopoDS_Shape &in_shape, - const Point<3> &origin, - const Tensor<1,3> &direction, - const double tolerance) + template + Point line_intersection(const TopoDS_Shape &in_shape, + const Point &origin, + const Tensor<1,dim> &direction, + const double tolerance) { // translating original Point to gp point gp_Pnt P0 = point(origin); - gp_Ax1 gpaxis(P0, gp_Dir(direction[0], direction[1], direction[2])); + gp_Ax1 gpaxis(P0, gp_Dir(direction[0], dim > 1 ? direction[1] : 0, dim>2 ? direction[2] : 0)); gp_Lin line(gpaxis); // destination point @@ -413,7 +438,7 @@ namespace OpenCASCADE Assert(Inters.IsDone(), ExcMessage("Could not project point.")); double minDistance = 1e7; - Point<3> result; + Point result; for (int i=0; i(Inters.Pnt(i+1)); } } return result; } - TopoDS_Edge interpolation_curve(std::vector > &curve_points, - const Tensor<1,3> &direction, + template + TopoDS_Edge interpolation_curve(std::vector > &curve_points, + const Tensor<1,dim> &direction, const bool closed, const double tolerance) { @@ -439,7 +465,7 @@ namespace OpenCASCADE if (direction*direction > 0) { std::sort(curve_points.begin(), curve_points.end(), - [&](const Point<3> &p1, const Point<3> &p2) + [&](const Point &p1, const Point &p2) { return OpenCASCADE::point_compare(p1, p2, direction, tolerance); }); @@ -527,15 +553,12 @@ namespace OpenCASCADE unsigned int point_index = start_point_index; // point_index and face_index always run together - std::vector > pointlist; + 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])); + pointlist.push_back(current_point); // Get next point if (face_to_verts[face_index].first != point_index ) @@ -551,7 +574,7 @@ namespace OpenCASCADE } while (point_index != start_point_index); - interpolation_curves.push_back(interpolation_curve(pointlist, Tensor<1,3>(), true)); + interpolation_curves.push_back(interpolation_curve(pointlist, Tensor<1,spacedim>(), true)); finished = true; for (const auto &f : visited_faces) @@ -566,11 +589,10 @@ namespace OpenCASCADE } - - std::tuple, TopoDS_Shape, double, double> - project_point_and_pull_back(const TopoDS_Shape &in_shape, - const Point<3> &origin, - const double tolerance) + template + std::tuple, TopoDS_Shape, double, double> project_point_and_pull_back(const TopoDS_Shape &in_shape, + const Point &origin, + const double tolerance) { TopExp_Explorer exp; gp_Pnt Pproj = point(origin); @@ -598,7 +620,7 @@ namespace OpenCASCADE SurfToProj->D0(proj_params.X(),proj_params.Y(),tmp_proj); - double distance = point(tmp_proj).distance(origin); + double distance = point(tmp_proj).distance(origin); if (distance < minDistance) { minDistance = distance; @@ -643,16 +665,17 @@ namespace OpenCASCADE } Assert(counter > 0, ExcMessage("Could not find projection points.")); - return std::tuple, TopoDS_Shape, double, double> - (point(Pproj),out_shape, u, v); + return std::tuple, TopoDS_Shape, double, double> + (point(Pproj),out_shape, u, v); } - Point<3> closest_point(const TopoDS_Shape &in_shape, - const Point<3> &origin, - const double tolerance) + template + Point closest_point(const TopoDS_Shape &in_shape, + const Point &origin, + const double tolerance) { - std::tuple, TopoDS_Shape, double, double> + std::tuple, TopoDS_Shape, double, double> ref = project_point_and_pull_back(in_shape, origin, tolerance); return std::get<0>(ref); } @@ -690,26 +713,27 @@ namespace OpenCASCADE return push_forward_and_differential_forms(face, u, v, tolerance); } - Point<3> push_forward(const TopoDS_Shape &in_shape, - const double u, - const double v) + template + Point push_forward(const TopoDS_Shape &in_shape, + const double u, + const double v) { switch (in_shape.ShapeType()) { case TopAbs_FACE: { BRepAdaptor_Surface surf(TopoDS::Face(in_shape)); - return point(surf.Value(u,v)); + return point(surf.Value(u,v)); } case TopAbs_EDGE: { BRepAdaptor_Curve curve(TopoDS::Edge(in_shape)); - return point(curve.Value(u)); + return point(curve.Value(u)); } default: Assert(false, ExcUnsupportedShape()); } - return Point<3>(); + return Point(); } std::tuple, Tensor<1,3>, double, double> @@ -739,13 +763,14 @@ namespace OpenCASCADE Max_Curvature *= -1; } - return std::tuple, Tensor<1,3>, double, double>(point(Value), normal, Min_Curvature, Max_Curvature); + return std::tuple, Tensor<1,3>, double, double>(point<3>(Value), normal, Min_Curvature, Max_Curvature); } + template void create_triangulation(const TopoDS_Face &face, - Triangulation<2,3> &tria) + Triangulation<2,spacedim> &tria) { BRepAdaptor_Surface surf(face); const double u0 = surf.FirstUParameter(); @@ -754,13 +779,13 @@ namespace OpenCASCADE const double v1 = surf.LastVParameter(); std::vector > cells; - std::vector > vertices; + std::vector > vertices; SubCellData t; - vertices.push_back(point(surf.Value(u0,v0))); - vertices.push_back(point(surf.Value(u1,v0))); - vertices.push_back(point(surf.Value(u0,v1))); - vertices.push_back(point(surf.Value(u1,v1))); + vertices.push_back(point(surf.Value(u0,v0))); + vertices.push_back(point(surf.Value(u1,v0))); + vertices.push_back(point(surf.Value(u0,v1))); + vertices.push_back(point(surf.Value(u1,v1))); CellData<2> cell; for (unsigned int i=0; i<4; ++i) @@ -770,22 +795,7 @@ namespace OpenCASCADE tria.create_triangulation(vertices, cells, t); } - - - // Explicit instantiations - template - std::vector create_curves_from_triangulation_boundary - (const Triangulation<2, 2> &triangulation, - const Mapping<2, 2> &mapping); - - - - template - std::vector create_curves_from_triangulation_boundary - (const Triangulation<2, 3> &triangulation, - const Mapping<2, 3> &mapping); - - +#include "utilities.inst" } // end namespace diff --git a/source/opencascade/utilities.inst.in b/source/opencascade/utilities.inst.in new file mode 100644 index 0000000000..7fae0b9abf --- /dev/null +++ b/source/opencascade/utilities.inst.in @@ -0,0 +1,62 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 1998 - 2017 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + +for (deal_II_dimension : DIMENSIONS) +{ +#if deal_II_dimension > 1 + // Explicit instantiations for dim = 2 and 3 + template + std::vector create_curves_from_triangulation_boundary + (const Triangulation<2, deal_II_dimension> &triangulation, + const Mapping<2, deal_II_dimension> &mapping); + + template bool point_compare(const Point& p1, + const Point& p2, + const Tensor<1,deal_II_dimension>& direction, + const double tolerance ); + + template Point point(const gp_Pnt &p, const double &tolerance); + + template gp_Pnt point(const Point &p); + + template + TopoDS_Edge interpolation_curve(std::vector >& curve_points, + const Tensor<1, deal_II_dimension>& direction=Tensor<1,deal_II_dimension>(), + const bool closed, + const double tolerance); + + template Point push_forward(const TopoDS_Shape &in_shape, + const double u, + const double v); + + template Point line_intersection(const TopoDS_Shape &in_shape, + const Point &origin, + const Tensor<1,deal_II_dimension> &direction, + const double tolerance); + + template void create_triangulation(const TopoDS_Face &face, + Triangulation<2,deal_II_dimension> &tria); + + template std::tuple, TopoDS_Shape, double, double> + project_point_and_pull_back(const TopoDS_Shape &in_shape, + const Point& origin, + const double tolerance); + + template Point closest_point(const TopoDS_Shape &in_shape, + const Point &origin, + const double tolerance); + +#endif +} diff --git a/tests/opencascade/arclength_boundary_02.cc b/tests/opencascade/arclength_boundary_02.cc new file mode 100644 index 0000000000..47ffb03e95 --- /dev/null +++ b/tests/opencascade/arclength_boundary_02.cc @@ -0,0 +1,56 @@ +//----------------------------------------------------------- +// +// 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. +// +//----------------------------------------------------------- + +#include "../tests.h" + +#include +#include +#include +#include +#include + +#include +#include +#include + +// 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 > pts; + pts.push_back(Point<2>(0,0)); + pts.push_back(Point<2>(0,1)); + pts.push_back(Point<2>(1,1)); + pts.push_back(Point<2>(1,0)); + + TopoDS_Edge edge = interpolation_curve(pts, Tensor<1,2>(), true); + ArclengthProjectionLineManifold<2,2> manifold(edge); + + Triangulation<2> tria; + GridGenerator::hyper_cube(tria); + + tria.set_all_manifold_ids_on_boundary(1); + tria.begin_active()->face(2)->set_manifold_id(numbers::invalid_manifold_id); + tria.set_manifold(1, manifold); + + tria.refine_global(3); + 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 index 0000000000..fea54f20a8 --- /dev/null +++ b/tests/opencascade/arclength_boundary_02.output @@ -0,0 +1,151 @@ +$NOD +81 +1 0 0 0 +2 1 0 0 +3 0 1 0 +4 1 1 0 +5 0.5 0 0 +6 -0.251337 0.515857 0 +7 1.16406 0.505545 0 +8 0.504074 1.14672 0 +9 0.4896 0.521016 0 +10 0.25 0 0 +11 0.75 0 0 +12 -0.158882 0.242372 0 +13 -0.20277 0.796991 0 +14 1.11495 0.242077 0 +15 1.13535 0.7715 0 +16 0.242844 1.10455 0 +17 0.767092 1.12113 0 +18 0.4948 0.260508 0 +19 0.496837 0.833869 0 +20 0.119131 0.518436 0 +21 0.826831 0.51328 0 +22 0.180414 0.257274 0 +23 0.79253 0.255303 0 +24 0.174797 0.80468 0 +25 0.79798 0.801633 0 +26 0.125 0 0 +27 0.375 0 0 +28 0.625 0 0 +29 0.875 0 0 +30 -0.08483 0.11767 0 +31 -0.217128 0.375124 0 +32 -0.249712 0.660431 0 +33 -0.114709 0.911629 0 +34 1.06551 0.117217 0 +35 1.1481 0.372213 0 +36 1.16089 0.639761 0 +37 1.08298 0.894921 0 +38 0.117841 1.06055 0 +39 0.372227 1.13328 0 +40 0.636548 1.14359 0 +41 0.891132 1.07494 0 +42 0.4974 0.130254 0 +43 0.4922 0.390762 0 +44 0.493218 0.677443 0 +45 0.500456 0.990296 0 +46 -0.0661032 0.517147 0 +47 0.304365 0.519726 0 +48 0.658215 0.517148 0 +49 0.995446 0.509412 0 +50 0.215207 0.128637 0 +51 0.149773 0.387855 0 +52 0.0107661 0.249823 0 +53 0.337607 0.258891 0 +54 0.771265 0.127652 0 +55 0.80968 0.384292 0 +56 0.643665 0.257905 0 +57 0.953738 0.24869 0 +58 0.146964 0.661558 0 +59 0.208821 0.954615 0 +60 -0.0139862 0.800836 0 +61 0.335817 0.819275 0 +62 0.812405 0.657456 0 +63 0.782536 0.961381 0 +64 0.647409 0.817751 0 +65 0.966663 0.786567 0 +66 0.0672094 0.124472 0 +67 0.356303 0.129445 0 +68 -0.033678 0.381489 0 +69 0.320986 0.389308 0 +70 0.634332 0.128953 0 +71 0.915373 0.123867 0 +72 0.65094 0.387527 0 +73 0.978892 0.378252 0 +74 -0.0513737 0.660995 0 +75 0.320091 0.6695 0 +76 0.0516048 0.929231 0 +77 0.354022 0.976278 0 +78 0.652812 0.667449 0 +79 0.986646 0.648609 0 +80 0.641978 0.980671 0 +81 0.927966 0.926509 0 +$ENDNOD +$ELM +64 +1 3 0 0 4 1 26 66 30 +2 3 0 0 4 26 10 50 66 +3 3 0 0 4 30 66 52 12 +4 3 0 0 4 66 50 22 52 +5 3 0 0 4 10 27 67 50 +6 3 0 0 4 27 5 42 67 +7 3 0 0 4 50 67 53 22 +8 3 0 0 4 67 42 18 53 +9 3 0 0 4 12 52 68 31 +10 3 0 0 4 52 22 51 68 +11 3 0 0 4 31 68 46 6 +12 3 0 0 4 68 51 20 46 +13 3 0 0 4 22 53 69 51 +14 3 0 0 4 53 18 43 69 +15 3 0 0 4 51 69 47 20 +16 3 0 0 4 69 43 9 47 +17 3 0 0 4 5 28 70 42 +18 3 0 0 4 28 11 54 70 +19 3 0 0 4 42 70 56 18 +20 3 0 0 4 70 54 23 56 +21 3 0 0 4 11 29 71 54 +22 3 0 0 4 29 2 34 71 +23 3 0 0 4 54 71 57 23 +24 3 0 0 4 71 34 14 57 +25 3 0 0 4 18 56 72 43 +26 3 0 0 4 56 23 55 72 +27 3 0 0 4 43 72 48 9 +28 3 0 0 4 72 55 21 48 +29 3 0 0 4 23 57 73 55 +30 3 0 0 4 57 14 35 73 +31 3 0 0 4 55 73 49 21 +32 3 0 0 4 73 35 7 49 +33 3 0 0 4 6 46 74 32 +34 3 0 0 4 46 20 58 74 +35 3 0 0 4 32 74 60 13 +36 3 0 0 4 74 58 24 60 +37 3 0 0 4 20 47 75 58 +38 3 0 0 4 47 9 44 75 +39 3 0 0 4 58 75 61 24 +40 3 0 0 4 75 44 19 61 +41 3 0 0 4 13 60 76 33 +42 3 0 0 4 60 24 59 76 +43 3 0 0 4 33 76 38 3 +44 3 0 0 4 76 59 16 38 +45 3 0 0 4 24 61 77 59 +46 3 0 0 4 61 19 45 77 +47 3 0 0 4 59 77 39 16 +48 3 0 0 4 77 45 8 39 +49 3 0 0 4 9 48 78 44 +50 3 0 0 4 48 21 62 78 +51 3 0 0 4 44 78 64 19 +52 3 0 0 4 78 62 25 64 +53 3 0 0 4 21 49 79 62 +54 3 0 0 4 49 7 36 79 +55 3 0 0 4 62 79 65 25 +56 3 0 0 4 79 36 15 65 +57 3 0 0 4 19 64 80 45 +58 3 0 0 4 64 25 63 80 +59 3 0 0 4 45 80 40 8 +60 3 0 0 4 80 63 17 40 +61 3 0 0 4 25 65 81 63 +62 3 0 0 4 65 15 37 81 +63 3 0 0 4 63 81 41 17 +64 3 0 0 4 81 37 4 41 +$ENDELM -- 2.39.5