From: Nicola Giuliani Date: Mon, 16 Dec 2019 08:58:10 +0000 (+0100) Subject: added default case for normal to mesh projection manifold X-Git-Tag: v9.2.0-rc1~776^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F9174%2Fhead;p=dealii.git added default case for normal to mesh projection manifold and indent --- diff --git a/source/opencascade/manifold_lib.cc b/source/opencascade/manifold_lib.cc index dd873d0d15..701eaada99 100644 --- a/source/opencascade/manifold_lib.cc +++ b/source/opencascade/manifold_lib.cc @@ -302,7 +302,45 @@ namespace OpenCASCADE } default: { - AssertThrow(false, ExcNotImplemented()); + // Given an arbitrary number of points we compute all the possible + // normal vectors + for (unsigned int i = 0; i < surrounding_points.size(); ++i) + for (unsigned int j = 0; j < surrounding_points.size(); ++j) + if (j != i) + for (unsigned int k = 0; k < surrounding_points.size(); ++k) + if (k != j && k != i) + { + Tensor<1, 3> u = + surrounding_points[i] - surrounding_points[j]; + Tensor<1, 3> v = + surrounding_points[i] - surrounding_points[k]; + const double n_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(n_coords); + if (n1.norm() > tolerance) + { + n1 = n1 / n1.norm(); + if (average_normal.norm() < tolerance) + average_normal = n1; + else + { + auto dot_prod = n1 * average_normal; + // We check that the direction of the normal + // vector w.r.t the current average, and make + // sure we flip it if it is opposite + if (dot_prod > 0) + average_normal += n1; + else + average_normal -= n1; + } + } + } + Assert( + average_normal.norm() > tolerance, + ExcMessage( + "Failed to compute a normal: the normal estimated via the surrounding points turns out to be a null vector, making the projection direction undetermined.")); + average_normal = average_normal / average_normal.norm(); break; } } diff --git a/tests/opencascade/normal_to_mesh_projection_04.cc b/tests/opencascade/normal_to_mesh_projection_04.cc new file mode 100644 index 0000000000..9f303abfd7 --- /dev/null +++ b/tests/opencascade/normal_to_mesh_projection_04.cc @@ -0,0 +1,101 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2019 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.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + + +// Create a BSpline surface, and test NormalToMeshProjectionManifold projection +// using 5 surrounding points. + +#include +#include +#include + +#include +#include +#include + +#include +#include +#include + +#include +#include + +#include +#include +#include +#include +#include +#include + +#include "../tests.h" + +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>(0.5, 2, -0.5)); + 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 < pts.size(); ++i) + pts[i] += Point<3>(0, 1, 0); + TopoDS_Edge edge2 = interpolation_curve(pts); + + TopoDS_Face face = BRepFill::Face(edge1, edge2); + + NormalToMeshProjectionManifold<2, 3> manifold(face); + FlatManifold<2, 3> flat_manifold; + std::vector> surrounding_points; + surrounding_points.push_back(Point<3>(0., 0., 0.)); + surrounding_points.push_back(Point<3>(1., 0., 0.)); + Vector weights(surrounding_points.size()); + weights = 1. / surrounding_points.size(); + auto intermediate_point = + manifold.get_new_point(ArrayView>(&surrounding_points[0], + surrounding_points.size()), + ArrayView(&weights[0], weights.size())); + surrounding_points.push_back(intermediate_point); + surrounding_points.push_back(Point<3>(0., 1., 0.)); + surrounding_points.push_back(Point<3>(1., 1., 0.)); + weights.reinit(surrounding_points.size()); + weights = 1. / surrounding_points.size(); + + + auto new_point = + manifold.get_new_point(ArrayView>(&surrounding_points[0], + surrounding_points.size()), + ArrayView(&weights[0], weights.size())); + auto new_point_flat = + flat_manifold.get_new_point(ArrayView>(&surrounding_points[0], + surrounding_points.size()), + ArrayView(&weights[0], weights.size())); + + for (unsigned int i = 0; i < surrounding_points.size(); ++i) + { + deallog << "Surrunding point " << i << " " << surrounding_points[i] + << " weight " << weights[i] << std::endl; + } + deallog << "Mean point " << new_point_flat << std::endl; + deallog << "Projected point " + << " " << new_point << std::endl; +} diff --git a/tests/opencascade/normal_to_mesh_projection_04.output b/tests/opencascade/normal_to_mesh_projection_04.output new file mode 100644 index 0000000000..0a8857d5a7 --- /dev/null +++ b/tests/opencascade/normal_to_mesh_projection_04.output @@ -0,0 +1,8 @@ + +DEAL::Surrunding point 0 0.00000 0.00000 0.00000 weight 0.200000 +DEAL::Surrunding point 1 1.00000 0.00000 0.00000 weight 0.200000 +DEAL::Surrunding point 2 0.500000 0.00000 0.245356 weight 0.200000 +DEAL::Surrunding point 3 0.00000 1.00000 0.00000 weight 0.200000 +DEAL::Surrunding point 4 1.00000 1.00000 0.00000 weight 0.200000 +DEAL::Mean point 0.500000 0.400000 0.0490712 +DEAL::Projected point 0.500000 0.447721 0.245356