diff --git a/kratos/geometries/nurbs_curve_on_surface_geometry.h b/kratos/geometries/nurbs_curve_on_surface_geometry.h index 54a2b6026b3f..b49cfbcf87ab 100644 --- a/kratos/geometries/nurbs_curve_on_surface_geometry.h +++ b/kratos/geometries/nurbs_curve_on_surface_geometry.h @@ -260,11 +260,127 @@ class NurbsCurveOnSurfaceGeometry : public GeometrySpansLocalSpace(surface_spans_u, 0); mpNurbsSurface->SpansLocalSpace(surface_spans_v, 1); - CurveAxisIntersection::ComputeAxisIntersection( - rSpans, - *(mpNurbsCurve.get()), Start, End, - surface_spans_u, surface_spans_v, - 1e-6); + // In the SBM case the NurbsSurface is marked as IS_SBM = true + const bool is_sbm = mpNurbsSurface->GetValue(IS_SBM); + + if (is_sbm) { + // SBM case: compute only knot spans intersections + ComputeAxisIntersectionSBM( + rSpans, + Start, End, + surface_spans_u, surface_spans_v); + } + else { + // trimming case: compute also inner knot spans intersections + CurveAxisIntersection::ComputeAxisIntersection( + rSpans, + *(mpNurbsCurve.get()), Start, End, + surface_spans_u, surface_spans_v, + 1e-6); + } + } + + /** + * @brief + * SBM case: compute only knot spans intersections + * @param rSpans + * @param Start + * @param End + * @param surface_spans_u + * @param surface_spans_v + * @param Tolerance + * @return int + */ + void ComputeAxisIntersectionSBM( + std::vector& rSpans, + double Start, + double End, + const std::vector& surface_spans_u, + const std::vector& surface_spans_v) const + { + CoordinatesArrayType physical_coord_1 = ZeroVector(3); + CoordinatesArrayType local_coord_1 = ZeroVector(3); + local_coord_1[0] = Start; + mpNurbsCurve->GlobalCoordinates(physical_coord_1, local_coord_1); + + CoordinatesArrayType physical_coord_2 = ZeroVector(3); + CoordinatesArrayType local_coord_2 = ZeroVector(3); + local_coord_2[0] = End; + mpNurbsCurve->GlobalCoordinates(physical_coord_2, local_coord_2); + + // 2D need to re-order + bool reverse_order = false; + if ((physical_coord_2[0] - physical_coord_1[0] < 0 ) || + (physical_coord_2[1] - physical_coord_1[1] < 0 )) { + reverse_order = true; + } + std::vector tempSpans; + const double tolerance_orientation = 1e-10; + const double tolerance_intersection = 1e-15; + + // Scale factor between volume coordinate and surface one + double physical_length_segment = norm_2(physical_coord_1-physical_coord_2); + double parameter_length_segment = norm_2(local_coord_1-local_coord_2); + double scale_factor = parameter_length_segment/physical_length_segment; + + std::vector knot_interval(2); + + // Compute constant coordinate -> understand segment orientation (vertical/horizontal) + if (std::abs(physical_coord_1[0]-physical_coord_2[0]) > tolerance_orientation) { + // horizontal case + knot_interval[0] = physical_coord_1[0]; + knot_interval[1] = physical_coord_2[0]; + std::sort(knot_interval.begin(), knot_interval.end()); + + knot_interval[0] -= tolerance_intersection; + knot_interval[1] += tolerance_intersection; + // Compare with volume_spans_u + for (IndexType i = 0; i < surface_spans_u.size(); i++) { + double curr_knot_value = surface_spans_u[i]; + if (curr_knot_value < knot_interval[0]) {continue;} + if (std::abs(curr_knot_value - knot_interval[0]) < tolerance_intersection*10) knot_interval[0] = curr_knot_value; + if (curr_knot_value > knot_interval[1]) {break;} + if (std::abs(curr_knot_value - knot_interval[1]) < tolerance_intersection*10) knot_interval[1] = curr_knot_value; + double knot_value_in_curve_parameter = Start + (curr_knot_value-knot_interval[0]) * scale_factor; + + tempSpans.push_back(knot_value_in_curve_parameter); + } + + } else if (std::abs(physical_coord_1[1]-physical_coord_2[1]) > tolerance_orientation) { + // vertical case + knot_interval[0] = physical_coord_1[1]; + knot_interval[1] = physical_coord_2[1]; + std::sort(knot_interval.begin(), knot_interval.end()); + + knot_interval[0] -= tolerance_intersection; + knot_interval[1] += tolerance_intersection; + // Compare with volume_spans_v + for (IndexType i = 0; i < surface_spans_v.size(); i++) { + double curr_knot_value = surface_spans_v[i]; + if (curr_knot_value < knot_interval[0]) {continue;} + if (std::abs(curr_knot_value - knot_interval[0]) < tolerance_intersection*10) knot_interval[0] = curr_knot_value; + if (curr_knot_value > knot_interval[1]) {break;} + if (std::abs(curr_knot_value - knot_interval[1]) < tolerance_intersection*10) knot_interval[1] = curr_knot_value; + double knot_value_in_curve_parameter = Start + (curr_knot_value-knot_interval[0]) * scale_factor; + + tempSpans.push_back(knot_value_in_curve_parameter); + } + } else { + KRATOS_ERROR << "ComputeAxisIntersectionSBM :: brep not parallel to any surface knot vector"; + } + + rSpans.resize(tempSpans.size()); + if (tempSpans.size() > 2){ + if (reverse_order) { + for (IndexType i = 0; i < tempSpans.size(); i++) { + rSpans[i] = End - tempSpans[tempSpans.size()-1-i]; + } + } else { + rSpans = tempSpans; + } + } else rSpans = tempSpans; + + KRATOS_ERROR_IF(rSpans.size()<2) << "ComputeAxisIntersectionSBM :: Wrong number of intersection found (<2)" << std::endl; } /* @brief Provides the nurbs boundaries of the NURBS/B-Spline curve. @@ -449,7 +565,58 @@ class NurbsCurveOnSurfaceGeometry : public GeometryGetValue(IS_SBM) ; + bool is_brep_internal = false; + IndexType InternalBrepSpanU; + IndexType InternalBrepSpanV; + + if (is_sbm) { + // If IS_SBM -> check to which axis the brep is alligned + std::vector first_integration_point(2); // first integration point of the brep in the parameter space + std::vector last_integration_point(2); // last integration point of the brep in the parameter space + mpNurbsCurve->GlobalSpaceDerivatives(first_integration_point,rIntegrationPoints[0],1); + mpNurbsCurve->GlobalSpaceDerivatives(last_integration_point,rIntegrationPoints[rIntegrationPoints.size()-1],1); + + // check if the brep is internal (external breps do not need to be computed in a different knot span) + is_brep_internal = true; + const double tolerance = 1e-14; + // At this point all the true are boundary GPs + if (std::abs(first_integration_point[0][0]-mpNurbsSurface->KnotsU()[0]) < tolerance) + is_brep_internal = false; + if (std::abs(first_integration_point[0][0]-mpNurbsSurface->KnotsU()[mpNurbsSurface->KnotsU().size()-1]) < tolerance) + is_brep_internal = false; + if (std::abs(first_integration_point[0][1]-mpNurbsSurface->KnotsV()[0]) < tolerance) + is_brep_internal = false; + if (std::abs(first_integration_point[0][1]-mpNurbsSurface->KnotsV()[mpNurbsSurface->KnotsV().size()-1]) < tolerance) + is_brep_internal = false; + + // If is_brep_internal -> the current GP is a surrogate GP otherwise it is an external GP + if (is_brep_internal) { + // brep = knot spans edge -> same SpanU and SpanV for all its integration points + InternalBrepSpanU = NurbsUtilities::GetLowerSpan(mpNurbsSurface->PolynomialDegreeU(), mpNurbsSurface->KnotsU(), first_integration_point[0][0]); + InternalBrepSpanV = NurbsUtilities::GetLowerSpan(mpNurbsSurface->PolynomialDegreeV(), mpNurbsSurface->KnotsV(), first_integration_point[0][1]); + + // Understand if the knot-boundary is along U or V + bool gp_is_along_V; + if (std::abs(first_integration_point[0][0] - last_integration_point[0][0]) < tolerance) { + gp_is_along_V = true; + } + else {gp_is_along_V = false;} + + if (gp_is_along_V && first_integration_point[0][1] > last_integration_point[0][1]) { + // Is decreasing along y -> Add 1 at InternalBrepSpanU + InternalBrepSpanU++; + } + if (!gp_is_along_V && first_integration_point[0][0] < last_integration_point[0][0]) { + // Is increasing along x -> Add 1 at InternalBrepSpanV + InternalBrepSpanV++; + } + } + } + + // Loop over the integration points of the brep + for (IndexType i = 0; i < rIntegrationPoints.size(); ++i) { std::vector global_space_derivatives(2); mpNurbsCurve->GlobalSpaceDerivatives( @@ -457,15 +624,28 @@ class NurbsCurveOnSurfaceGeometry : public GeometryIsRational()) { + if (mpNurbsSurface->IsRational()) { shape_function_container.ComputeNurbsShapeFunctionValues( mpNurbsSurface->KnotsU(), mpNurbsSurface->KnotsV(), mpNurbsSurface->Weights(), global_space_derivatives[0][0], global_space_derivatives[0][1]); } else { - shape_function_container.ComputeBSplineShapeFunctionValues( - mpNurbsSurface->KnotsU(), mpNurbsSurface->KnotsV(), - global_space_derivatives[0][0], global_space_derivatives[0][1]); + if (is_sbm && is_brep_internal) { + // SBM case on internal brep + shape_function_container.ComputeBSplineShapeFunctionValuesAtSpan( + mpNurbsSurface->KnotsU(), + mpNurbsSurface->KnotsV(), + InternalBrepSpanU, + InternalBrepSpanV, + global_space_derivatives[0][0], + global_space_derivatives[0][1]); + } + else { + // Trimming case or external brep + shape_function_container.ComputeBSplineShapeFunctionValues( + mpNurbsSurface->KnotsU(), mpNurbsSurface->KnotsV(), + global_space_derivatives[0][0], global_space_derivatives[0][1]); + } } /// Get List of Control Points diff --git a/kratos/includes/variables.h b/kratos/includes/variables.h index d601268fafd6..fd9981490a20 100644 --- a/kratos/includes/variables.h +++ b/kratos/includes/variables.h @@ -485,6 +485,8 @@ namespace Kratos KRATOS_DEFINE_VARIABLE(double, VARIATIONAL_REDISTANCE_COEFFICIENT_FIRST) KRATOS_DEFINE_VARIABLE(double, VARIATIONAL_REDISTANCE_COEFFICIENT_SECOND) + // SBM variables + KRATOS_DEFINE_VARIABLE(bool, IS_SBM) } // namespace Kratos. #undef KRATOS_EXPORT_MACRO diff --git a/kratos/sources/variables.cpp b/kratos/sources/variables.cpp index 729b9bec5c73..74dc59edd22a 100644 --- a/kratos/sources/variables.cpp +++ b/kratos/sources/variables.cpp @@ -477,6 +477,9 @@ KRATOS_CREATE_3D_VARIABLE_WITH_COMPONENTS(PARAMETER_2D_COORDINATES) KRATOS_CREATE_VARIABLE(double, VARIATIONAL_REDISTANCE_COEFFICIENT_FIRST) KRATOS_CREATE_VARIABLE(double, VARIATIONAL_REDISTANCE_COEFFICIENT_SECOND) +// SBM variables +KRATOS_CREATE_VARIABLE(bool, IS_SBM) + //------------------------------------------------------------------------------// //------------------------------------------------------------------------------// //------------------------------------------------------------------------------// @@ -954,5 +957,8 @@ void KratosApplication::RegisterVariables() { // Variational redistance KRATOS_REGISTER_VARIABLE(VARIATIONAL_REDISTANCE_COEFFICIENT_FIRST) KRATOS_REGISTER_VARIABLE(VARIATIONAL_REDISTANCE_COEFFICIENT_SECOND) + + // SBM variables + KRATOS_REGISTER_VARIABLE(IS_SBM) } } // namespace Kratos. diff --git a/kratos/tests/cpp_tests/geometries/test_nurbs_curve_on_surface_geometry.cpp b/kratos/tests/cpp_tests/geometries/test_nurbs_curve_on_surface_geometry.cpp index 08485ecf79d0..b3a2fcd819bd 100644 --- a/kratos/tests/cpp_tests/geometries/test_nurbs_curve_on_surface_geometry.cpp +++ b/kratos/tests/cpp_tests/geometries/test_nurbs_curve_on_surface_geometry.cpp @@ -313,6 +313,55 @@ typedef Node NodeType; return NurbsCurveOnSurfaceGeometry<3, PointerVector, PointerVector>(surface, curve); } + NurbsCurveOnSurfaceGeometry<3, PointerVector, PointerVector> GenerateReferenceNurbs2dforKnotIntersections(std::vector& brep_coordinates) + { + Vector knot_vector_u = ZeroVector(4); + Vector knot_vector_v = ZeroVector(4); + const SizeType p = 1; + const SizeType q = 1; + + // lower + PointerVector points_surface; + + knot_vector_u[0] = 0.0; knot_vector_u[1] = 1.0; knot_vector_u[2] = 2.0; knot_vector_u[3] = 3.0; + knot_vector_v[0] = 0.0; knot_vector_v[1] = 1.0; knot_vector_v[2] = 2.0; knot_vector_v[3] = 3.0; + + points_surface.push_back(Point::Pointer(new Point(0.0, 0.0, 0.0))); + points_surface.push_back(Point::Pointer(new Point(1.0, 0.0, 0.0))); + points_surface.push_back(Point::Pointer(new Point(2.0, 0.0, 0.0))); + points_surface.push_back(Point::Pointer(new Point(3.0, 0.0, 0.0))); + // + points_surface.push_back(Point::Pointer(new Point(0.0, 0.0, 0.0))); + points_surface.push_back(Point::Pointer(new Point(1.0, 1.0, 0.0))); + points_surface.push_back(Point::Pointer(new Point(2.0, 1.0, 0.0))); + points_surface.push_back(Point::Pointer(new Point(3.0, 1.0, 0.0))); + // + points_surface.push_back(Point::Pointer(new Point(0.0, 2.0, 0.0))); + points_surface.push_back(Point::Pointer(new Point(1.0, 2.0, 0.0))); + points_surface.push_back(Point::Pointer(new Point(2.0, 2.0, 0.0))); + points_surface.push_back(Point::Pointer(new Point(3.0, 2.0, 0.0))); + // + points_surface.push_back(Point::Pointer(new Point(0.0, 3.0, 0.0))); + points_surface.push_back(Point::Pointer(new Point(1.0, 3.0, 0.0))); + points_surface.push_back(Point::Pointer(new Point(2.0, 3.0, 0.0))); + points_surface.push_back(Point::Pointer(new Point(3.0, 3.0, 0.0))); + + auto p_surface = Kratos::make_shared>>(points_surface, p, q, knot_vector_u, knot_vector_v); + + Vector knot_vector_curve = ZeroVector(2); + knot_vector_curve[0] = 0.0; + knot_vector_curve[1] = 3.0; + + NurbsCurveGeometry<2, PointerVector>::PointsArrayType segment; + segment.push_back(Point::Pointer(new Point(brep_coordinates[0][0], brep_coordinates[0][1]))); + segment.push_back(Point::Pointer(new Point(brep_coordinates[1][0], brep_coordinates[1][1]))); + + auto p_curve = Kratos::make_shared>>(segment, p, knot_vector_curve); + + // Create and return a curve on surface geometry + return NurbsCurveOnSurfaceGeometry<3, PointerVector, PointerVector>(p_surface, p_curve); + } + ///// Tests // Ported from the ANurbs library (https://github.com/oberbichler/ANurbs) KRATOS_TEST_CASE_IN_SUITE(BSplineCurveOnSurfaceBSpline, KratosCoreNurbsGeometriesFastSuite) @@ -537,5 +586,244 @@ typedef Node NodeType; KRATOS_EXPECT_EQ(quadrature_points[i].GetGeometryType(), geometry_type); } } + + // test intersection with background surface in the SBM scenario + KRATOS_TEST_CASE_IN_SUITE(NurbsCurveOnSurfaceIntersectionSpansExternalBodyFitted, KratosCoreNurbsGeometriesFastSuite) + { + // Create a Nurbs curve on a Nurbs surface + std::vector brep_coordinates(2); + brep_coordinates[0].resize(2); + brep_coordinates[1].resize(2); + + brep_coordinates[0][0] = 0.0; + brep_coordinates[0][1] = 0.0; + brep_coordinates[1][0] = 0.0; + brep_coordinates[1][1] = 3.0; + + auto curve_on_surface = GenerateReferenceNurbs2dforKnotIntersections(brep_coordinates); + + auto p_surface = curve_on_surface.pGetGeometryPart(GeometryType::BACKGROUND_GEOMETRY_INDEX); + + p_surface->SetValue(IS_SBM, true); + + std::vector spans; + + curve_on_surface.SpansLocalSpace(spans); + + // Test size + KRATOS_EXPECT_EQ(spans.size(), 4); + + // Compare each value + KRATOS_EXPECT_NEAR(spans[0], 0.0, TOLERANCE); + KRATOS_EXPECT_NEAR(spans[1], 1.0, TOLERANCE); + KRATOS_EXPECT_NEAR(spans[2], 2.0, TOLERANCE); + KRATOS_EXPECT_NEAR(spans[3], 3.0, TOLERANCE); + + const auto geometry_family = GeometryData::KratosGeometryFamily::Kratos_Nurbs; + const auto geometry_type = GeometryData::KratosGeometryType::Kratos_Nurbs_Curve_On_Surface; + KRATOS_EXPECT_EQ(curve_on_surface.GetGeometryFamily(), geometry_family); + KRATOS_EXPECT_EQ(curve_on_surface.GetGeometryType(), geometry_type); + } + + // test intersection with background surface in the SBM scenario + KRATOS_TEST_CASE_IN_SUITE(NurbsCurveOnSurfaceIntersectionSpansExternalSBM, KratosCoreNurbsGeometriesFastSuite) + { + // Create a Nurbs curve on a Nurbs surface + std::vector brep_coordinates(2); + brep_coordinates[0].resize(2); + brep_coordinates[1].resize(2); + + brep_coordinates[0][0] = 1.0; + brep_coordinates[0][1] = 0.0; + brep_coordinates[1][0] = 2.0; + brep_coordinates[1][1] = 0.0; + + auto curve_on_surface = GenerateReferenceNurbs2dforKnotIntersections(brep_coordinates); + + auto p_surface = curve_on_surface.pGetGeometryPart(GeometryType::BACKGROUND_GEOMETRY_INDEX); + + p_surface->SetValue(IS_SBM, true); + + std::vector spans; + + curve_on_surface.SpansLocalSpace(spans); + + // Test size + KRATOS_EXPECT_EQ(spans.size(), 2); + + // Compare each value + KRATOS_EXPECT_NEAR(spans[0], 0.0, TOLERANCE); + KRATOS_EXPECT_NEAR(spans[1], 3.0, TOLERANCE); + + const auto geometry_family = GeometryData::KratosGeometryFamily::Kratos_Nurbs; + const auto geometry_type = GeometryData::KratosGeometryType::Kratos_Nurbs_Curve_On_Surface; + KRATOS_EXPECT_EQ(curve_on_surface.GetGeometryFamily(), geometry_family); + KRATOS_EXPECT_EQ(curve_on_surface.GetGeometryType(), geometry_type); + } + + // test intersection with background surface in the SBM scenario + KRATOS_TEST_CASE_IN_SUITE(NurbsCurveOnSurfaceIntersectionSpansInternalSBM, KratosCoreNurbsGeometriesFastSuite) + { + // Create a Nurbs curve on a Nurbs surface + std::vector brep_coordinates(2); + brep_coordinates[0].resize(2); + brep_coordinates[1].resize(2); + + brep_coordinates[0][0] = 1.0; + brep_coordinates[0][1] = 1.0; + brep_coordinates[1][0] = 2.0; + brep_coordinates[1][1] = 1.0; + + auto curve_on_surface = GenerateReferenceNurbs2dforKnotIntersections(brep_coordinates); + auto p_surface = curve_on_surface.pGetGeometryPart(GeometryType::BACKGROUND_GEOMETRY_INDEX); + p_surface->SetValue(IS_SBM, true); + + std::vector spans; + + curve_on_surface.SpansLocalSpace(spans); + + // Test size + KRATOS_EXPECT_EQ(spans.size(), 2); + + // Compare each value + KRATOS_EXPECT_NEAR(spans[0], 0.0, TOLERANCE); + KRATOS_EXPECT_NEAR(spans[1], 3.0, TOLERANCE); + + const auto geometry_family = GeometryData::KratosGeometryFamily::Kratos_Nurbs; + const auto geometry_type = GeometryData::KratosGeometryType::Kratos_Nurbs_Curve_On_Surface; + KRATOS_EXPECT_EQ(curve_on_surface.GetGeometryFamily(), geometry_family); + KRATOS_EXPECT_EQ(curve_on_surface.GetGeometryType(), geometry_type); + } + + // test quadrature points of curve on surface + KRATOS_TEST_CASE_IN_SUITE(NurbsCurveOnSurfaceCreateQuadraturePointsSBMInternal, KratosCoreNurbsGeometriesFastSuite) + { + // Create a Nurbs curve on a Nurbs surface + std::vector brep_coordinates(2); + brep_coordinates[0].resize(2); + brep_coordinates[1].resize(2); + + brep_coordinates[0][0] = 1.0; + brep_coordinates[0][1] = 1.0; + brep_coordinates[1][0] = 2.0; + brep_coordinates[1][1] = 1.0; + + auto curve_on_surface = GenerateReferenceNurbs2dforKnotIntersections(brep_coordinates); + auto p_surface = curve_on_surface.pGetGeometryPart(GeometryType::BACKGROUND_GEOMETRY_INDEX); + p_surface->SetValue(IS_SBM, true); + + // Check general information, input to ouput + typename Geometry::IntegrationPointsArrayType integration_points; + IntegrationInfo integration_info = curve_on_surface.GetDefaultIntegrationInfo(); + curve_on_surface.CreateIntegrationPoints(integration_points, integration_info); + + typename Geometry::GeometriesArrayType quadrature_points; + curve_on_surface.CreateQuadraturePointGeometries(quadrature_points, 3, integration_points, integration_info); + + KRATOS_EXPECT_EQ(quadrature_points.size(), 3); + + array_1d global_coords; + array_1d local_coords; + local_coords[0] = integration_points[2][0]; + local_coords[1] = integration_points[2][1]; + curve_on_surface.GlobalCoordinates(global_coords, local_coords); + + KRATOS_EXPECT_VECTOR_NEAR(quadrature_points[2].Center(), global_coords, TOLERANCE); + + std::vector expected_cps(4); + + expected_cps[0].resize(2); expected_cps[1].resize(2); expected_cps[2].resize(2); expected_cps[3].resize(2); + + expected_cps[0][0] = 1.0; expected_cps[0][1] = 1.0; + expected_cps[1][0] = 2.0; expected_cps[1][1] = 1.0; + expected_cps[2][0] = 1.0; expected_cps[2][1] = 2.0; + expected_cps[3][0] = 2.0; expected_cps[3][1] = 2.0; + + // Check the correct SpanU++ + KRATOS_EXPECT_EQ(quadrature_points[0].size(), 4); + for (IndexType i = 0; i < quadrature_points[0].size(); i++){ + + Vector cp_coordinates(2); + cp_coordinates[0] = quadrature_points[0][i].X(); + cp_coordinates[1] = quadrature_points[0][i].Y(); + + KRATOS_EXPECT_NEAR(norm_2(cp_coordinates-expected_cps[i]), 0.0, TOLERANCE); + + } + const auto geometry_family = GeometryData::KratosGeometryFamily::Kratos_Quadrature_Geometry; + const auto geometry_type = GeometryData::KratosGeometryType::Kratos_Quadrature_Point_Curve_On_Surface_Geometry; + + for (IndexType i = 0; i < quadrature_points.size(); ++i) { + KRATOS_EXPECT_EQ(quadrature_points[i].GetGeometryFamily(), geometry_family); + KRATOS_EXPECT_EQ(quadrature_points[i].GetGeometryType(), geometry_type); + } + } + + // test quadrature points of curve on surface + KRATOS_TEST_CASE_IN_SUITE(NurbsCurveOnSurfaceCreateQuadraturePointsSBMExternal, KratosCoreNurbsGeometriesFastSuite) + { + // Create a Nurbs curve on a Nurbs surface + std::vector brep_coordinates(2); + brep_coordinates[0].resize(2); + brep_coordinates[1].resize(2); + + brep_coordinates[0][0] = 1.0; + brep_coordinates[0][1] = 0.0; + brep_coordinates[1][0] = 2.0; + brep_coordinates[1][1] = 0.0; + + auto curve_on_surface = GenerateReferenceNurbs2dforKnotIntersections(brep_coordinates); + auto p_surface = curve_on_surface.pGetGeometryPart(GeometryType::BACKGROUND_GEOMETRY_INDEX); + p_surface->SetValue(IS_SBM, true); + + // Check general information, input to ouput + typename Geometry::IntegrationPointsArrayType integration_points; + IntegrationInfo integration_info = curve_on_surface.GetDefaultIntegrationInfo(); + curve_on_surface.CreateIntegrationPoints(integration_points, integration_info); + + typename Geometry::GeometriesArrayType quadrature_points; + curve_on_surface.CreateQuadraturePointGeometries(quadrature_points, 3, integration_points, integration_info); + + KRATOS_EXPECT_EQ(quadrature_points.size(), 3); + + array_1d global_coords; + array_1d local_coords; + local_coords[0] = integration_points[2][0]; + local_coords[1] = integration_points[2][1]; + curve_on_surface.GlobalCoordinates(global_coords, local_coords); + + KRATOS_EXPECT_VECTOR_NEAR(quadrature_points[2].Center(), global_coords, TOLERANCE); + + std::vector expected_cps(4); + + expected_cps[0].resize(2); expected_cps[1].resize(2); expected_cps[2].resize(2); expected_cps[3].resize(2); + + expected_cps[0][0] = 1.0; expected_cps[0][1] = 0.0; + expected_cps[1][0] = 2.0; expected_cps[1][1] = 0.0; + expected_cps[2][0] = 1.0; expected_cps[2][1] = 1.0; + expected_cps[3][0] = 2.0; expected_cps[3][1] = 1.0; + + // Check the correct SpanU++ + KRATOS_EXPECT_EQ(quadrature_points[0].size(), 4); + for (IndexType i = 0; i < quadrature_points[0].size(); i++){ + + Vector cp_coordinates(2); + cp_coordinates[0] = quadrature_points[0][i].X(); + cp_coordinates[1] = quadrature_points[0][i].Y(); + + KRATOS_EXPECT_NEAR(norm_2(cp_coordinates-expected_cps[i]), 0.0, TOLERANCE); + + } + const auto geometry_family = GeometryData::KratosGeometryFamily::Kratos_Quadrature_Geometry; + const auto geometry_type = GeometryData::KratosGeometryType::Kratos_Quadrature_Point_Curve_On_Surface_Geometry; + + for (IndexType i = 0; i < quadrature_points.size(); ++i) { + KRATOS_EXPECT_EQ(quadrature_points[i].GetGeometryFamily(), geometry_family); + KRATOS_EXPECT_EQ(quadrature_points[i].GetGeometryType(), geometry_type); + } + } + + + } // namespace Testing. } // namespace Kratos.