From dfb466fadf126e14df4158fd2085cf90578daed6 Mon Sep 17 00:00:00 2001 From: Luca Heltai Date: Wed, 21 Mar 2018 18:09:11 +0100 Subject: [PATCH] Fix bug in normal to mesh projection due to Boundary -> Manifold --- source/opencascade/boundary_lib.cc | 25 +++- .../normal_to_mesh_projection_03.cc | 107 ++++++++++++++++++ .../normal_to_mesh_projection_03.output | 52 +++++++++ 3 files changed, 180 insertions(+), 4 deletions(-) create mode 100644 tests/opencascade/normal_to_mesh_projection_03.cc create mode 100644 tests/opencascade/normal_to_mesh_projection_03.output diff --git a/source/opencascade/boundary_lib.cc b/source/opencascade/boundary_lib.cc index 1a3df42f5c..773a30b7c0 100644 --- a/source/opencascade/boundary_lib.cc +++ b/source/opencascade/boundary_lib.cc @@ -189,6 +189,27 @@ namespace OpenCASCADE average_normal /= average_normal.norm(); break; } + case 4: + { + Tensor<1,3> u = surrounding_points[1]-surrounding_points[0]; + Tensor<1,3> v = surrounding_points[2]-surrounding_points[0]; + const double n1_coords[3] = {u[1] *v[2]-u[2] *v[1],u[2] *v[0]-u[0] *v[2],u[0] *v[1]-u[1] *v[0]}; + Tensor<1,3> n1(n1_coords); + n1 = n1/n1.norm(); + u = surrounding_points[2]-surrounding_points[3]; + v = surrounding_points[1]-surrounding_points[3]; + const double n2_coords[3] = {u[1] *v[2]-u[2] *v[1],u[2] *v[0]-u[0] *v[2],u[0] *v[1]-u[1] *v[0]}; + Tensor<1,3> n2(n2_coords); + n2 = n2/n2.norm(); + + average_normal = (n1+n2)/2.0; + + Assert(average_normal.norm() > tolerance, + ExcMessage("Failed to refine cell: the normal estimated via the surrounding points turns out to be a null vector, making the projection direction undetermined.")); + + average_normal /= average_normal.norm(); + break; + } case 8: { Tensor<1,3> u = surrounding_points[1]-surrounding_points[0]; @@ -211,10 +232,6 @@ namespace OpenCASCADE const double n4_coords[3] = {u[1] *v[2]-u[2] *v[1],u[2] *v[0]-u[0] *v[2],u[0] *v[1]-u[1] *v[0]}; Tensor<1,3> n4(n4_coords); n4 = n4/n4.norm(); - //for (unsigned int i=0; i +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include + +using namespace OpenCASCADE; + +int main () +{ + initlog(); + + // Create a bspline passing through the points + std::vector > pts; + pts.push_back(Point<3>(0,0,0)); + pts.push_back(Point<3>(1,0,0)); + pts.push_back(Point<3>(1,0,-1)); + pts.push_back(Point<3>(0,0,-1)); + + TopoDS_Edge edge1 = interpolation_curve(pts); + for (unsigned int i=0; i(0,1,0); + TopoDS_Edge edge2 = interpolation_curve(pts); + + TopoDS_Face face = BRepFill::Face (edge1, edge2); + + + NormalToMeshProjectionBoundary<2,3> manifold(face); + + Triangulation<2,3> tria; + GridGenerator::hyper_cube(tria); + + tria.begin_active()->set_all_manifold_ids(1); + tria.set_manifold(1, manifold); + + tria.refine_global(1); + + Tensor<1,3> direction({.1, .01, .001}); + double tolerance = 1e-10; + + { + FE_Q<2,3> fe(2); + DoFHandler<2,3> dh(tria); + dh.distribute_dofs(fe); + MappingQGeneric<2,3> mapping2(2); + std::vector> spoints(dh.n_dofs()); + DoFTools::map_dofs_to_support_points(mapping2, dh, spoints); + + std::sort(spoints.begin(), spoints.end(), + [&](const Point<3> &p1, const Point<3> &p2) + { + return OpenCASCADE::point_compare(p1, p2, direction, tolerance); + }); + + for (auto p: spoints) + deallog << p << std::endl; + deallog << "============" << std::endl; + } + + tria.refine_global(1); + { + FE_Q<2,3> fe(1); + DoFHandler<2,3> dh(tria); + dh.distribute_dofs(fe); + std::vector> spoints(dh.n_dofs()); + DoFTools::map_dofs_to_support_points(StaticMappingQ1<2,3>::mapping, + dh, spoints); + + std::sort(spoints.begin(), spoints.end(), + [&](const Point<3> &p1, const Point<3> &p2) + { + return OpenCASCADE::point_compare(p1, p2, direction, tolerance); + }); + + for (auto p: spoints) + deallog << p << std::endl; + } +} diff --git a/tests/opencascade/normal_to_mesh_projection_03.output b/tests/opencascade/normal_to_mesh_projection_03.output new file mode 100644 index 0000000000..87f13a9267 --- /dev/null +++ b/tests/opencascade/normal_to_mesh_projection_03.output @@ -0,0 +1,52 @@ + +DEAL::0.00000 0.00000 0.00000 +DEAL::0.00000 0.250000 0.00000 +DEAL::0.00000 0.500000 0.00000 +DEAL::0.00000 0.750000 0.00000 +DEAL::0.00000 1.00000 0.00000 +DEAL::0.234806 0.00000 0.153640 +DEAL::0.234806 0.250000 0.153640 +DEAL::0.234806 0.500000 0.153640 +DEAL::0.234806 0.750000 0.153640 +DEAL::0.234806 1.00000 0.153640 +DEAL::0.500000 0.00000 0.245356 +DEAL::0.500000 0.250000 0.245356 +DEAL::0.500000 0.500000 0.245356 +DEAL::0.500000 0.750000 0.245356 +DEAL::0.500000 1.00000 0.245356 +DEAL::0.789901 0.00000 0.203990 +DEAL::0.789901 0.250000 0.203990 +DEAL::0.789901 0.500000 0.203990 +DEAL::0.789901 0.750000 0.203990 +DEAL::0.789901 1.00000 0.203990 +DEAL::1.00000 0.00000 0.00000 +DEAL::1.00000 0.250000 -2.79105e-17 +DEAL::1.00000 0.500000 -4.62817e-17 +DEAL::1.00000 0.750000 -3.17705e-17 +DEAL::1.00000 1.00000 0.00000 +DEAL::============ +DEAL::0.00000 0.00000 0.00000 +DEAL::0.00000 0.250000 0.00000 +DEAL::0.00000 0.500000 0.00000 +DEAL::0.00000 0.750000 0.00000 +DEAL::0.00000 1.00000 0.00000 +DEAL::0.234806 0.00000 0.153640 +DEAL::0.234806 0.250000 0.153640 +DEAL::0.234806 0.500000 0.153640 +DEAL::0.234806 0.750000 0.153640 +DEAL::0.234806 1.00000 0.153640 +DEAL::0.500000 0.00000 0.245356 +DEAL::0.500000 0.250000 0.245356 +DEAL::0.500000 0.500000 0.245356 +DEAL::0.500000 0.750000 0.245356 +DEAL::0.500000 1.00000 0.245356 +DEAL::0.789901 0.00000 0.203990 +DEAL::0.789901 0.250000 0.203990 +DEAL::0.789901 0.500000 0.203990 +DEAL::0.789901 0.750000 0.203990 +DEAL::0.789901 1.00000 0.203990 +DEAL::1.00000 0.00000 0.00000 +DEAL::1.00000 0.250000 -2.79105e-17 +DEAL::1.00000 0.500000 -4.62817e-17 +DEAL::1.00000 0.750000 -3.17705e-17 +DEAL::1.00000 1.00000 0.00000 -- 2.39.5