From 82f0cbca98c36061cfaea3e1d15c2f2ab414d9b6 Mon Sep 17 00:00:00 2001 From: Daniel Arndt Date: Tue, 12 May 2020 12:47:18 -0400 Subject: [PATCH] Instantiate OpenCASCADE --- source/opencascade/boundary_lib.inst.in | 2 + source/opencascade/manifold_lib.cc | 348 +++++++++++++----------- source/opencascade/manifold_lib.inst.in | 2 + 3 files changed, 191 insertions(+), 161 deletions(-) diff --git a/source/opencascade/boundary_lib.inst.in b/source/opencascade/boundary_lib.inst.in index 10c4a204dc..a6b81548d5 100644 --- a/source/opencascade/boundary_lib.inst.in +++ b/source/opencascade/boundary_lib.inst.in @@ -21,6 +21,8 @@ for (deal_II_dimension : DIMENSIONS) template class DirectionalProjectionBoundary; template class NormalToMeshProjectionBoundary; #if deal_II_dimension <= 2 + template class NormalProjectionBoundary; template class DirectionalProjectionBoundary; + template class NormalToMeshProjectionBoundary; #endif } diff --git a/source/opencascade/manifold_lib.cc b/source/opencascade/manifold_lib.cc index 701eaada99..819d9fb8de 100644 --- a/source/opencascade/manifold_lib.cc +++ b/source/opencascade/manifold_lib.cc @@ -75,8 +75,7 @@ namespace OpenCASCADE } } // namespace - /*============================== NormalProjectionManifold - * ==============================*/ + /*======================= NormalProjectionManifold =========================*/ template NormalProjectionManifold::NormalProjectionManifold( const TopoDS_Shape &sh, @@ -117,8 +116,7 @@ namespace OpenCASCADE } - /*============================== DirectionalProjectionManifold - * ==============================*/ + /*===================== DirectionalProjectionManifold ======================*/ template DirectionalProjectionManifold::DirectionalProjectionManifold( const TopoDS_Shape & sh, @@ -162,8 +160,7 @@ namespace OpenCASCADE - /*============================== NormalToMeshProjectionManifold - * ==============================*/ + /*===================== NormalToMeshProjectionManifold =====================*/ template NormalToMeshProjectionManifold::NormalToMeshProjectionManifold( const TopoDS_Shape &sh, @@ -187,170 +184,199 @@ namespace OpenCASCADE } + namespace + { + template + Point + internal_project_to_manifold(const TopoDS_Shape &, + const double, + const ArrayView> &, + const Point &) + { + Assert(false, ExcNotImplemented()); + return {}; + } + + template <> + Point<3> + internal_project_to_manifold( + const TopoDS_Shape & sh, + const double tolerance, + const ArrayView> &surrounding_points, + const Point<3> & candidate) + { + constexpr int spacedim = 3; + TopoDS_Shape out_shape; + Tensor<1, spacedim> average_normal; +# ifdef DEBUG + for (unsigned int i = 0; i < surrounding_points.size(); ++i) + { + Assert(closest_point(sh, surrounding_points[i], tolerance) + .distance(surrounding_points[i]) < + std::max(tolerance * surrounding_points[i].norm(), + tolerance), + ExcPointNotOnManifold(surrounding_points[i])); + } +# endif + + switch (surrounding_points.size()) + { + case 2: + { + for (unsigned int i = 0; i < surrounding_points.size(); ++i) + { + std::tuple, Tensor<1, 3>, double, double> + p_and_diff_forms = closest_point_and_differential_forms( + sh, surrounding_points[i], tolerance); + average_normal += std::get<1>(p_and_diff_forms); + } + + average_normal /= 2.0; + + Assert( + average_normal.norm() > 1e-4, + ExcMessage( + "Failed to refine cell: the average of the surface normals at the surrounding edge turns out to be a null vector, making the projection direction undetermined.")); + + Tensor<1, 3> T = surrounding_points[0] - surrounding_points[1]; + T /= T.norm(); + average_normal = average_normal - (average_normal * T) * T; + 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]; + 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(); + u = surrounding_points[4] - surrounding_points[7]; + v = surrounding_points[6] - surrounding_points[7]; + const double n3_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> n3(n3_coords); + n3 = n3 / n3.norm(); + u = surrounding_points[6] - surrounding_points[7]; + v = surrounding_points[5] - surrounding_points[7]; + 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(); + + average_normal = (n1 + n2 + n3 + n4) / 4.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; + } + default: + { + // 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; + } + } + + return line_intersection(sh, candidate, average_normal, tolerance); + } + } // namespace + + template Point NormalToMeshProjectionManifold::project_to_manifold( const ArrayView> &surrounding_points, const Point & candidate) const { - TopoDS_Shape out_shape; - Tensor<1, 3> average_normal; -# ifdef DEBUG - for (unsigned int i = 0; i < surrounding_points.size(); ++i) - { - Assert(closest_point(sh, surrounding_points[i], tolerance) - .distance(surrounding_points[i]) < - std::max(tolerance * surrounding_points[i].norm(), tolerance), - ExcPointNotOnManifold(surrounding_points[i])); - } -# endif - - switch (surrounding_points.size()) - { - case 2: - { - for (unsigned int i = 0; i < surrounding_points.size(); ++i) - { - std::tuple, Tensor<1, 3>, double, double> - p_and_diff_forms = - closest_point_and_differential_forms(sh, - surrounding_points[i], - tolerance); - average_normal += std::get<1>(p_and_diff_forms); - } - - average_normal /= 2.0; - - Assert( - average_normal.norm() > 1e-4, - ExcMessage( - "Failed to refine cell: the average of the surface normals at the surrounding edge turns out to be a null vector, making the projection direction undetermined.")); - - Tensor<1, 3> T = surrounding_points[0] - surrounding_points[1]; - T /= T.norm(); - average_normal = average_normal - (average_normal * T) * T; - 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]; - 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(); - u = surrounding_points[4] - surrounding_points[7]; - v = surrounding_points[6] - surrounding_points[7]; - const double n3_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> n3(n3_coords); - n3 = n3 / n3.norm(); - u = surrounding_points[6] - surrounding_points[7]; - v = surrounding_points[5] - surrounding_points[7]; - 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(); - - average_normal = (n1 + n2 + n3 + n4) / 4.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; - } - default: - { - // 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; - } - } - - return line_intersection(sh, candidate, average_normal, tolerance); + return internal_project_to_manifold(sh, + tolerance, + surrounding_points, + candidate); } - /*============================== ArclengthProjectionLineManifold - * ==============================*/ + /*==================== ArclengthProjectionLineManifold =====================*/ template ArclengthProjectionLineManifold:: ArclengthProjectionLineManifold(const TopoDS_Shape &sh, diff --git a/source/opencascade/manifold_lib.inst.in b/source/opencascade/manifold_lib.inst.in index 937221b2ab..d4e529f5db 100644 --- a/source/opencascade/manifold_lib.inst.in +++ b/source/opencascade/manifold_lib.inst.in @@ -23,7 +23,9 @@ for (deal_II_dimension : DIMENSIONS) template class ArclengthProjectionLineManifold; template class NURBSPatchManifold; #if deal_II_dimension <= 2 + template class NormalProjectionManifold; template class DirectionalProjectionManifold; + template class NormalToMeshProjectionManifold; template class ArclengthProjectionLineManifold; template class NURBSPatchManifold; #endif -- 2.39.5