From db53ed6c223ecf21c52c4646070d132a073d1477 Mon Sep 17 00:00:00 2001 From: Luca Heltai Date: Tue, 12 Aug 2014 15:54:12 +0200 Subject: [PATCH] Added OCC::NormalProjectionBoundary Added interpolation curve boundary test. Updated documentation. --- include/deal.II/grid/occ_boundary_lib.h | 110 ++++++++++++++++++ include/deal.II/grid/occ_utilities.h | 35 ++++-- source/grid/CMakeLists.txt | 2 + source/grid/occ_boundary_lib.cc | 42 +++++++ source/grid/occ_boundary_lib.inst.in | 25 ++++ source/grid/occ_utilities.cc | 14 +-- tests/opencascade/circle_normal_projection.cc | 102 ++++++++++++++++ .../circle_normal_projection.output | 43 +++++++ .../interpolation_curve_boundary.cc | 87 ++++++++++++++ .../interpolation_curve_boundary.output | 48 ++++++++ 10 files changed, 490 insertions(+), 18 deletions(-) create mode 100644 include/deal.II/grid/occ_boundary_lib.h create mode 100644 source/grid/occ_boundary_lib.cc create mode 100644 source/grid/occ_boundary_lib.inst.in create mode 100644 tests/opencascade/circle_normal_projection.cc create mode 100644 tests/opencascade/circle_normal_projection.output create mode 100644 tests/opencascade/interpolation_curve_boundary.cc create mode 100644 tests/opencascade/interpolation_curve_boundary.output diff --git a/include/deal.II/grid/occ_boundary_lib.h b/include/deal.II/grid/occ_boundary_lib.h new file mode 100644 index 0000000000..e46be218fe --- /dev/null +++ b/include/deal.II/grid/occ_boundary_lib.h @@ -0,0 +1,110 @@ +// --------------------------------------------------------------------- +// $Id$ +// +// Copyright (C) 2014 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + +#ifndef __deal2__occ_boundary_lib_h +#define __deal2__occ_boundary_lib_h + +#include + +#ifdef DEAL_II_WITH_OPENCASCADE + +#include +#include + +DEAL_II_NAMESPACE_OPEN + +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. + * + * This class makes no assumptions on the shape you pass to it, and + * the topological dimension of the Manifold is inferred from the + * TopoDS_Shape itself. In debug mode there is a sanity check to + * make sure that the surrounding points (the ones used in + * project_to_manifold()) actually live on the Manifold, i.e., + * calling OpenCASCADE::closest_point() on those points leaves them + * untouched. If this is not the case, an ExcPointNotOnManifold is + * thrown. + * + * This could happen, for example, if you are trying to use a shape + * of type TopoDS_Edge when projecting on a face. In this case, the + * vertices of the face would be collapsed to the edge, and your + * surrounding points would not be on lying on the given shape, + * raising an exception. + * + * @author Luca Heltai, Andrea Mola, 2011--2014. + */ + template + class NormalProjectionBoundary : public Boundary { + public: + 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 + * within tolerance from the given TopoDS_Shape. If this is not + * the case, an exception is thrown. + * + * The projected point is computed using OpenCASCADE normal + * projection algorithms. + */ + virtual Point + project_to_manifold (const std::vector > &surrounding_points, + const Point &candidate) const; + + /** + * Exception thrown when the point specified as argument does not + * lie between #tolerance from the given TopoDS_Shape. + */ + DeclException1 (ExcPointNotOnManifold, + Point, + <<"The point [ "< + #ifdef DEAL_II_WITH_OPENCASCADE #include @@ -27,7 +29,7 @@ #include #include -#include +#include DEAL_II_NAMESPACE_OPEN @@ -123,18 +125,28 @@ namespace OpenCASCADE /** * Creates a 3D smooth BSpline curve passing through the points in - * the assigned vector. The points are reordered internally - * according to their scalar product with the direction, if - * direction is different from zero, otherwise they are used as + * 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. * * This class is used to interpolate a BsplineCurve passing through - * an array of points, with a C2 Continuity. + * an array of points, with a C2 Continuity. If the optional + * parameter #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 + * input points. If the algorithm fails in generating such a curve, + * an exception is thrown. */ - TopoDS_Shape interpolation_curve(std::vector > &curve_points, - const dealii::Point<3> direction=dealii::Point<3>(), - const double tolerance=1e-7); + TopoDS_Edge interpolation_curve(std::vector > &curve_points, + const dealii::Point<3> direction=dealii::Point<3>(), + const bool closed=false, + const double tolerance=1e-7); + /** * Get the closest point to the given topological shape. If the @@ -155,19 +167,20 @@ namespace OpenCASCADE /** * Convert OpenCASCADE point into a Point<3>. */ - Point<3> Pnt(const gp_Pnt &p); + inline Point<3> Pnt(const gp_Pnt &p); /** * Convert Point<3> into OpenCASCADE point. */ - gp_Pnt Pnt(const Point<3> &p); + inline 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. + * 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>(), diff --git a/source/grid/CMakeLists.txt b/source/grid/CMakeLists.txt index 738fc8c43b..8fa6c75d02 100644 --- a/source/grid/CMakeLists.txt +++ b/source/grid/CMakeLists.txt @@ -26,6 +26,7 @@ SET(_src manifold.cc manifold_lib.cc occ_utilities.cc + occ_boundary_lib.cc persistent_tria.cc tria_accessor.cc tria_boundary.cc @@ -45,6 +46,7 @@ SET(_inst intergrid_map.inst.in manifold.inst.in manifold_lib.inst.in + occ_boundary_lib.inst.in tria_accessor.inst.in tria_boundary.inst.in tria_boundary_lib.inst.in diff --git a/source/grid/occ_boundary_lib.cc b/source/grid/occ_boundary_lib.cc new file mode 100644 index 0000000000..2ef10faf2e --- /dev/null +++ b/source/grid/occ_boundary_lib.cc @@ -0,0 +1,42 @@ +#include + +#ifdef DEAL_II_WITH_OPENCASCADE + +DEAL_II_NAMESPACE_OPEN + +namespace OpenCASCADE +{ + template + NormalProjectionBoundary::NormalProjectionBoundary(const TopoDS_Shape sh, + const double tolerance) : + sh(sh), + tolerance(tolerance) + { + Assert(spacedim == 3, ExcImpossibleInDim(spacedim)); + } + + + template + Point NormalProjectionBoundary:: + project_to_manifold (const std::vector > &surrounding_points, + const Point &candidate) const { + TopoDS_Shape out_shape; + double u=0, v=0; + for(unsigned int i=0; i0 ? + tolerance*surrounding_points[i].norm() : + tolerance), + ExcPointNotOnManifold(surrounding_points[i])); + + return closest_point(sh, candidate, out_shape, u, v); + } + + // Explicit instantiations +#include "occ_boundary_lib.inst" + +} // end namespace OpenCASCADE + +DEAL_II_NAMESPACE_CLOSE + +#endif diff --git a/source/grid/occ_boundary_lib.inst.in b/source/grid/occ_boundary_lib.inst.in new file mode 100644 index 0000000000..6e97bc7391 --- /dev/null +++ b/source/grid/occ_boundary_lib.inst.in @@ -0,0 +1,25 @@ +// --------------------------------------------------------------------- +// $Id$ +// +// Copyright (C) 1998 - 2013 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + + +for (deal_II_dimension : DIMENSIONS) + { + template class NormalProjectionBoundary; + } + + + diff --git a/source/grid/occ_utilities.cc b/source/grid/occ_utilities.cc index 61c62445c9..73b04d3ed1 100644 --- a/source/grid/occ_utilities.cc +++ b/source/grid/occ_utilities.cc @@ -1,5 +1,4 @@ #include -#include #ifdef DEAL_II_WITH_OPENCASCADE @@ -106,10 +105,10 @@ namespace OpenCASCADE } 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) + const dealii::Point<3> direction, + const double tolerance) { - const double rel_tol=std::max(p1.norm(), p2.norm())*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 @@ -245,8 +244,9 @@ namespace OpenCASCADE - TopoDS_Shape interpolation_curve(std::vector > &curve_points, + TopoDS_Edge interpolation_curve(std::vector > &curve_points, const Point<3> direction, + const bool closed, const double tolerance) { @@ -266,12 +266,12 @@ namespace OpenCASCADE } - GeomAPI_Interpolate bspline_generator(vertices, false, tolerance); + 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_Shape out_shape = BRepBuilderAPI_MakeEdge(bspline); + TopoDS_Edge out_shape = BRepBuilderAPI_MakeEdge(bspline); return out_shape; } diff --git a/tests/opencascade/circle_normal_projection.cc b/tests/opencascade/circle_normal_projection.cc new file mode 100644 index 0000000000..f40c3edeb7 --- /dev/null +++ b/tests/opencascade/circle_normal_projection.cc @@ -0,0 +1,102 @@ + +//---------------------------- 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 +// +// 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. +// +//---------------------------- occ_circle_projection.cc --------------------------- + + +// Create a circle, a Triangulation, and try to project normally on +// it. + +#include "../tests.h" + +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +using namespace OpenCASCADE; + +int main () +{ + std::ofstream logfile("output"); + deallog.attach(logfile); + deallog.depth_console(0); + + // The circle passing through the + // vertices of the unit square + gp_Dir z_axis(0.,0.,1.); + gp_Pnt center(.5,.5,0.); + gp_Ax2 axis(center, z_axis); + Standard_Real radius(std::sqrt(2.)/2.); + + Handle(Geom_Curve) circle = GC_MakeCircle(axis, radius); + TopoDS_Edge edge = BRepBuilderAPI_MakeEdge(circle); + + + // Create a boundary projector. + NormalProjectionBoundary<2,3> boundary_line(edge); + + // This one is for checking: This + // is what deal.II would do for a + // circle. + HyperBallBoundary<2,3> boundary_line_deal (Point<3>(.5,.5,0), + std::sqrt(2.)/2.); + + + + // The unit square. + Triangulation<2,3> tria; + GridGenerator::hyper_cube(tria); + + // 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 + // of the face. + tria.begin()->set_material_id(1); + + // We refine twice, and expect the + // outer points to end up on the + // circle. + tria.refine_global(2); + + + // You can open the generated file + // with paraview. + GridOut gridout; + gridout.write_ucd (tria, logfile); + + return 0; +} + diff --git a/tests/opencascade/circle_normal_projection.output b/tests/opencascade/circle_normal_projection.output new file mode 100644 index 0000000000..3291307f5a --- /dev/null +++ b/tests/opencascade/circle_normal_projection.output @@ -0,0 +1,43 @@ + +25 16 0 0 0 +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 +1 1 quad 1 10 22 12 +2 1 quad 10 5 18 22 +3 1 quad 12 22 20 6 +4 1 quad 22 18 9 20 +5 1 quad 5 11 23 18 +6 1 quad 11 2 14 23 +7 1 quad 18 23 21 9 +8 1 quad 23 14 7 21 +9 1 quad 6 20 24 13 +10 1 quad 20 9 19 24 +11 1 quad 13 24 16 3 +12 1 quad 24 19 8 16 +13 1 quad 9 21 25 19 +14 1 quad 21 7 15 25 +15 1 quad 19 25 17 8 +16 1 quad 25 15 4 17 diff --git a/tests/opencascade/interpolation_curve_boundary.cc b/tests/opencascade/interpolation_curve_boundary.cc new file mode 100644 index 0000000000..9d43fa31cf --- /dev/null +++ b/tests/opencascade/interpolation_curve_boundary.cc @@ -0,0 +1,87 @@ + +//---------------------------- 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 +// +// 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. +// +//---------------------------- occ_circle_projection.cc --------------------------- + + +// 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 +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +using namespace OpenCASCADE; + +int main () +{ + std::ofstream logfile("output"); + deallog.attach(logfile); + deallog.depth_console(0); + + // The curve passing through the vertices of the unit square. + std::vector > vertices; + vertices.push_back(Point<3>(0,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)); + TopoDS_Shape shape = interpolation_curve(vertices, Point<3>(), true); + + // Create a boundary projector. + NormalProjectionBoundary<2,3> boundary_line(shape); + + // The unit square. + Triangulation<2,3> tria; + GridGenerator::hyper_cube(tria); + + // Set the exterior boundary + tria.set_boundary(0, boundary_line); + + // This is here to ignore the points created in the interior of the + // face. + tria.begin()->set_material_id(1); + + // 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 paraview. + 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/interpolation_curve_boundary.output b/tests/opencascade/interpolation_curve_boundary.output new file mode 100644 index 0000000000..894bc26860 --- /dev/null +++ b/tests/opencascade/interpolation_curve_boundary.output @@ -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.439917 -0.236739 0.00000 +6 -0.310022 0.473832 0.00000 +7 1.14617 0.492348 0.00000 +8 0.488859 1.16307 0.00000 +9 0.470616 0.486564 0.00000 +10 0.207735 -0.140071 0.00000 +11 0.764326 -0.218873 0.00000 +12 -0.202060 0.202837 0.00000 +13 -0.224611 0.783020 0.00000 +14 1.10302 0.237503 0.00000 +15 1.12328 0.758622 0.00000 +16 0.234185 1.11236 0.00000 +17 0.760208 1.13858 0.00000 +18 0.455267 0.124913 0.00000 +19 0.479738 0.824816 0.00000 +20 0.0802967 0.480198 0.00000 +21 0.808395 0.489456 0.00000 +22 0.142719 0.173942 0.00000 +23 0.773464 0.171896 0.00000 +24 0.152383 0.790482 0.00000 +25 0.784659 0.794181 0.00000 +$ENDNOD +$ELM +16 +1 3 1 0 4 1 10 22 12 +2 3 1 0 4 10 5 18 22 +3 3 1 0 4 12 22 20 6 +4 3 1 0 4 22 18 9 20 +5 3 1 0 4 5 11 23 18 +6 3 1 0 4 11 2 14 23 +7 3 1 0 4 18 23 21 9 +8 3 1 0 4 23 14 7 21 +9 3 1 0 4 6 20 24 13 +10 3 1 0 4 20 9 19 24 +11 3 1 0 4 13 24 16 3 +12 3 1 0 4 24 19 8 16 +13 3 1 0 4 9 21 25 19 +14 3 1 0 4 21 7 15 25 +15 3 1 0 4 19 25 17 8 +16 3 1 0 4 25 15 4 17 +$ENDELM -- 2.39.5