diff --git a/applications/IgaApplication/custom_utilities/create_breps_sbm_utilities.h b/applications/IgaApplication/custom_utilities/create_breps_sbm_utilities.h new file mode 100644 index 000000000000..8a8170d1713a --- /dev/null +++ b/applications/IgaApplication/custom_utilities/create_breps_sbm_utilities.h @@ -0,0 +1,438 @@ +// | / | +// ' / __| _` | __| _ \ __| +// . \ | ( | | ( |\__ ` +// _|\_\_| \__,_|\__|\___/ ____/ +// Multi-Physics +// +// License: BSD License +// Kratos default license: kratos/license.txt +// +// Main authors: Nicolo' Antonelli +// Andrea Gorgi +// + +# pragma once + +// System includes + +// External includes + +// Project includes +#include "includes/kratos_parameters.h" +#include "includes/model_part.h" + +// Geometries +#include "geometries/nurbs_curve_geometry.h" +#include "geometries/brep_surface.h" +#include "geometries/brep_curve.h" +#include "geometries/brep_curve_on_surface.h" + + +namespace Kratos +{ + +///@name Kratos Classes +///@{ +template +class CreateBrepsSbmUtilities : public IO +{ + public: + + ///@} + ///@name Type Definitions + ///@{ + + /// Pointer definition of CreateBrepsSbmUtilities + KRATOS_CLASS_POINTER_DEFINITION(CreateBrepsSbmUtilities); + + using SizeType = std::size_t; + using IndexType = std::size_t; + + using GeometryType = Geometry; + using GeometryPointerType = typename GeometryType::Pointer; + using GeometrySurrogateArrayType = DenseVector; + + using ContainerNodeType = PointerVector; + using ContainerEmbeddedNodeType = PointerVector; + + using NurbsSurfaceGeometryType = NurbsSurfaceGeometry<3, PointerVector>; + using NurbsSurfaceGeometryPointerType = typename NurbsSurfaceGeometryType::Pointer; + + using BrepSurfaceType = BrepSurface; + using BrepCurveOnSurfaceType = BrepCurveOnSurface; + + using BrepCurveOnSurfaceLoopType = DenseVector; + using BrepCurveOnSurfaceLoopArrayType = DenseVector>; + + ///@} + ///@name Life Cycle + ///@{ + + /// Constructor with path to input file. + CreateBrepsSbmUtilities( + SizeType EchoLevel = 0) + : mEchoLevel(EchoLevel) + { + } + + /// Destructor. + ~CreateBrepsSbmUtilities() = default; + + ///@} + ///@name Python exposed Functions + ///@{ + + + /** + * @brief Create a Surrogate Boundary object + * + * @param pSurface + * @param rSurrogateModelPartInner + * @param rSurrogateModelPartOuter + * @param rCoordsA + * @param rCoordsB + * @param rModelPart + */ + void CreateSurrogateBoundary(NurbsSurfaceGeometryPointerType& pSurface, + const ModelPart& rSurrogateModelPartInner, + const ModelPart& rSurrogateModelPartOuter, + const Point& rCoordsA, + const Point& rCoordsB, + ModelPart& rModelPart) + { + CreateBrepSurface(pSurface, rSurrogateModelPartInner, rSurrogateModelPartOuter, rModelPart, mEchoLevel); + CreateBrepCurveOnSurfaces(pSurface, rSurrogateModelPartInner, rSurrogateModelPartOuter, rCoordsA, rCoordsB, rModelPart, mEchoLevel); + } + + /** + * @brief Create a Surrogate Boundary object + * + * @param pSurface + * @param rCoordsA + * @param rCoordsB + * @param rModelPart + */ + void CreateSurrogateBoundary(NurbsSurfaceGeometryPointerType& pSurface, + const Point& rCoordsA, + const Point& rCoordsB, + ModelPart& rModelPart) + { + CreateBrepSurface(pSurface, rModelPart, mEchoLevel); + IndexType id_brep_curve_on_surface = 2; // because id 1 is the brep surface + CreateBrepCurvesOnRectangle(pSurface, rCoordsA, rCoordsB, id_brep_curve_on_surface, rModelPart); + } + + +private: + + /** + * @brief Create a Brep Surface object + * + * @param pSurface + * @param rModelPart + * @param EchoLevel + */ + static void CreateBrepSurface( + NurbsSurfaceGeometryPointerType pSurface, + ModelPart& rModelPart, + const SizeType EchoLevel = 0) + { + KRATOS_INFO_IF("ReadBrepSurface", (EchoLevel > 3)) + << "Creating BrepSurface \""<< std::endl; + + BrepCurveOnSurfaceLoopArrayType outer_loops, inner_loops; + + auto p_brep_surface = + Kratos::make_shared( + pSurface, + outer_loops, + inner_loops, + false); + + // Sets the brep as geometry parent of the nurbs surface. + pSurface->SetGeometryParent(p_brep_surface.get()); + p_brep_surface->SetId(1); + rModelPart.AddGeometry(p_brep_surface); + } + + /** + * @brief Create a Brep Surface object + * + * @param pSurface + * @param rSurrogateModelPartInner + * @param rSurrogateModelPartOuter + * @param rModelPart + * @param EchoLevel + */ + static void CreateBrepSurface( + NurbsSurfaceGeometryPointerType pSurface, + const ModelPart& rSurrogateModelPartInner, + const ModelPart& rSurrogateModelPartOuter, + ModelPart& rModelPart, + const SizeType EchoLevel = 0) + { + KRATOS_INFO_IF("ReadBrepSurface", (EchoLevel > 3)) + << "Creating BrepSurface \""<< std::endl; + + BrepCurveOnSurfaceLoopArrayType outer_loops, inner_loops; + + GeometrySurrogateArrayType surrogate_outer_loop_geometries(rSurrogateModelPartOuter.NumberOfConditions()); + GeometrySurrogateArrayType surrogate_inner_loop_geometries(rSurrogateModelPartInner.NumberOfConditions()); + int count = 0; + for (auto i_cond : rSurrogateModelPartOuter.Conditions()) + { + surrogate_outer_loop_geometries[count] = i_cond.pGetGeometry(); + count++; + } + + count = 0; + for (auto i_cond : rSurrogateModelPartInner.Conditions()) + { + surrogate_inner_loop_geometries[count] = i_cond.pGetGeometry(); + count++; + } + + auto p_brep_surface = + Kratos::make_shared( + pSurface, + outer_loops, + inner_loops); + + p_brep_surface->SetSurrogateInnerLoopGeometries(surrogate_inner_loop_geometries); + p_brep_surface->SetSurrogateOuterLoopGeometries(surrogate_outer_loop_geometries); + + // Sets the brep as geometry parent of the nurbs surface. + pSurface->SetGeometryParent(p_brep_surface.get()); + p_brep_surface->SetId(1); + rModelPart.AddGeometry(p_brep_surface); + } + + /** + * @brief Create a Brep Curve On Surfaces object + * + * @param pSurface + * @param rSurrogateModelPartInner + * @param rSurrogateModelPartOuter + * @param rCoordsA + * @param rCoordsB + * @param rModelPart + * @param EchoLevel + */ + static void CreateBrepCurveOnSurfaces( + const NurbsSurfaceGeometryPointerType pSurface, + const ModelPart& rSurrogateModelPartInner, + const ModelPart& rSurrogateModelPartOuter, + const Point& rCoordsA, + const Point& rCoordsB, + ModelPart& rModelPart, + const SizeType EchoLevel = 0) { + // OUTER + IndexType id_brep_curve_on_surface = 2; // because id 1 is the brep surface + + if (rSurrogateModelPartOuter.NumberOfConditions() > 0) + { + for (auto &i_cond : rSurrogateModelPartOuter.Conditions()) { + Geometry::PointsArrayType points; + Point first_point = i_cond.GetGeometry()[0]; + Point second_point = i_cond.GetGeometry()[1]; + Vector knot_vector = ZeroVector(2); + + // define active range: the first point and second point are in order + Vector active_range_knot_vector = ZeroVector(2); + // Compute the knot vector needed + if (first_point[0] == second_point[0]) { + // the brep curve is vertical + active_range_knot_vector[0] = first_point[1]; + active_range_knot_vector[1] = second_point[1]; + } else { + // the brep curve is horizontal + active_range_knot_vector[0] = first_point[0]; + active_range_knot_vector[1] = second_point[0]; + } + + // check if the brep is entering or exiting + if (i_cond.Is(BOUNDARY)) + { + // if true is exiting -> reverse the order of the + Point temp = first_point; + first_point = second_point; + second_point = temp; + } + + Point::Pointer first_brep_point = Kratos::make_shared(first_point[0], first_point[1], 0.0); + Point::Pointer p_second_brep_point = Kratos::make_shared(second_point[0], second_point[1], 0.0); + + NurbsCurveGeometry<2, PointerVector>::Pointer p_trimming_curve = CreateBrepCurve(first_brep_point, p_second_brep_point, active_range_knot_vector); + + NurbsInterval brep_active_range(active_range_knot_vector[0], active_range_knot_vector[1]); + + bool curve_direction = true; + auto p_brep_curve_on_surface = Kratos::make_shared( + pSurface, p_trimming_curve, brep_active_range, curve_direction); + + p_brep_curve_on_surface->SetId(id_brep_curve_on_surface); + + rModelPart.AddGeometry(p_brep_curve_on_surface); + id_brep_curve_on_surface++; + } + + } else { + CreateBrepCurvesOnRectangle(pSurface, rCoordsA, rCoordsB, id_brep_curve_on_surface, rModelPart); + } + + + // INNER + for (IndexType i_element = 1; i_element < rSurrogateModelPartInner.Elements().size()+1; i_element++) { + IndexType first_surrogate_cond_id = rSurrogateModelPartInner.pGetElement(i_element)->GetGeometry()[0].Id(); // Element 1 because is the only surrogate loop + IndexType last_surrogate_cond_id = rSurrogateModelPartInner.pGetElement(i_element)->GetGeometry()[1].Id(); // Element 1 because is the only surrogate loop + + for (IndexType cond_id = first_surrogate_cond_id; cond_id <= last_surrogate_cond_id; ++cond_id) { + + auto p_cond = rSurrogateModelPartInner.pGetCondition(cond_id); + Geometry::PointsArrayType points; + Point first_point = p_cond->GetGeometry()[0]; + Point second_point = p_cond->GetGeometry()[1]; + Vector knot_vector = ZeroVector(2); + + // define active range: the first point and second point are in order + Vector active_range_knot_vector = ZeroVector(2); + // Compute the knot vector needed + if (first_point[0] == second_point[0]) { + // the brep curve is vertical + active_range_knot_vector[0] = first_point[1]; + active_range_knot_vector[1] = second_point[1]; + } else { + // the brep curve is horizontal + active_range_knot_vector[0] = first_point[0]; + active_range_knot_vector[1] = second_point[0]; + } + + // check if the brep is entering or exiting + if (p_cond->Is(BOUNDARY)) + { + // if true is exiting -> reverse the order of the + Point temp = first_point; + first_point = second_point; + second_point = temp; + } + + Point::Pointer p_first_brep_point = Kratos::make_shared(first_point[0], first_point[1], 0.0); + Point::Pointer p_second_brep_point = Kratos::make_shared(second_point[0], second_point[1], 0.0); + NurbsCurveGeometry<2, PointerVector>::Pointer p_trimming_curve = CreateBrepCurve(p_first_brep_point, p_second_brep_point, active_range_knot_vector); + NurbsInterval brep_active_range(active_range_knot_vector[0], active_range_knot_vector[1]); + + bool curve_direction = true; + auto p_brep_curve_on_surface = Kratos::make_shared( + pSurface, p_trimming_curve, brep_active_range, curve_direction); + + p_brep_curve_on_surface->SetId(id_brep_curve_on_surface); + rModelPart.AddGeometry(p_brep_curve_on_surface); + id_brep_curve_on_surface++; + } + } + } + + /** + * @brief Create a Single Brep object + * + * @param pFirstBrepPoint + * @param pSecondBrepPoint + * @param rActiveRangeKnotVector + * @return NurbsCurveGeometry<2, PointerVector>::Pointer + */ + static typename NurbsCurveGeometry<2, PointerVector>::Pointer CreateBrepCurve( + const Point::Pointer pFirstBrepPoint, + const Point::Pointer pSecondBrepPoint, + const Vector& rActiveRangeKnotVector) + { + // Create the data for the trimming curves + PointerVector control_points; + control_points.push_back(pFirstBrepPoint); + control_points.push_back(pSecondBrepPoint); + const int polynomial_degree = 1; + Vector knot_vector = ZeroVector(4) ; + knot_vector[0] = rActiveRangeKnotVector[0] ; + knot_vector[1] = rActiveRangeKnotVector[0] ; + knot_vector[2] = rActiveRangeKnotVector[1] ; + knot_vector[3] = rActiveRangeKnotVector[1] ; + // Create the trimming curves + typename NurbsCurveGeometry<2, PointerVector>::Pointer p_trimming_curve( + new NurbsCurveGeometry<2, PointerVector>( + control_points, + polynomial_degree, + knot_vector)); + return p_trimming_curve; + } + + /** + * @brief Create a Brep Curves On Rectangle object + * + * @param pSurfaceGeometry + * @param rCoordsA + * @param rCoordsB + * @param rLastGeometryId + * @param rModelPart + */ + static void CreateBrepCurvesOnRectangle(const NurbsSurfaceGeometryPointerType pSurfaceGeometry, + const Point& rCoordsA, + const Point& rCoordsB, + IndexType& rLastGeometryId, + ModelPart& rModelPart) { + Vector knot_vector = ZeroVector(2); + knot_vector[0] = 0.0; + knot_vector[1] = std::abs(rCoordsB[0] - rCoordsA[0]); + const int p = 1; + + NurbsCurveGeometry<2, PointerVector>::PointsArrayType segment1; + segment1.push_back(Point::Pointer(new Point(rCoordsA[0], rCoordsA[1]))); + segment1.push_back(Point::Pointer(new Point(rCoordsB[0], rCoordsA[1]))); + + NurbsCurveGeometry<2, PointerVector>::PointsArrayType segment2; + segment2.push_back(Point::Pointer(new Point(rCoordsB[0], rCoordsA[1]))); + segment2.push_back(Point::Pointer(new Point(rCoordsB[0], rCoordsB[1]))); + + NurbsCurveGeometry<2, PointerVector>::PointsArrayType segment3; + segment3.push_back(Point::Pointer(new Point(rCoordsB[0], rCoordsB[1]))); + segment3.push_back(Point::Pointer(new Point(rCoordsA[0], rCoordsB[1]))); + + NurbsCurveGeometry<2, PointerVector>::PointsArrayType segment4; + segment4.push_back(Point::Pointer(new Point(rCoordsA[0], rCoordsB[1]))); + segment4.push_back(Point::Pointer(new Point(rCoordsA[0], rCoordsA[1]))); + + auto p_curve_1 = Kratos::make_shared>>(segment1, p, knot_vector); + auto p_curve_2 = Kratos::make_shared>>(segment2, p, knot_vector); + auto p_curve_3 = Kratos::make_shared>>(segment3, p, knot_vector); + auto p_curve_4 = Kratos::make_shared>>(segment4, p, knot_vector); + + auto brep_curve_on_surface = BrepCurveOnSurface< PointerVector, true, PointerVector>(pSurfaceGeometry, p_curve_1); + auto p_brep_curve_on_surface = Kratos::make_shared(pSurfaceGeometry, p_curve_1); + p_brep_curve_on_surface->SetId(rLastGeometryId); + rModelPart.AddGeometry(p_brep_curve_on_surface); + + brep_curve_on_surface = BrepCurveOnSurface< PointerVector, true, PointerVector>(pSurfaceGeometry, p_curve_2); + p_brep_curve_on_surface = Kratos::make_shared(pSurfaceGeometry, p_curve_2); + p_brep_curve_on_surface->SetId(++rLastGeometryId); + rModelPart.AddGeometry(p_brep_curve_on_surface); + + brep_curve_on_surface = BrepCurveOnSurface< PointerVector, true, PointerVector>(pSurfaceGeometry, p_curve_3); + p_brep_curve_on_surface = Kratos::make_shared(pSurfaceGeometry, p_curve_3); + p_brep_curve_on_surface->SetId(++rLastGeometryId); + rModelPart.AddGeometry(p_brep_curve_on_surface); + + brep_curve_on_surface = BrepCurveOnSurface< PointerVector, true, PointerVector>(pSurfaceGeometry, p_curve_4); + p_brep_curve_on_surface = Kratos::make_shared(pSurfaceGeometry, p_curve_4); + p_brep_curve_on_surface->SetId(++rLastGeometryId); + rModelPart.AddGeometry(p_brep_curve_on_surface); + + rLastGeometryId++; + } + + ///@} + ///@name Utility functions + ///@{ + + int mEchoLevel; + + ///@} +}; // Class CreateBrepsSbmUtilities +} // namespace Kratos. diff --git a/applications/IgaApplication/tests/cpp_tests/test_sbm_create_brep_utilities.cpp b/applications/IgaApplication/tests/cpp_tests/test_sbm_create_brep_utilities.cpp new file mode 100644 index 000000000000..e4515a4eca71 --- /dev/null +++ b/applications/IgaApplication/tests/cpp_tests/test_sbm_create_brep_utilities.cpp @@ -0,0 +1,225 @@ +// | / | +// ' / __| _` | __| _ \ __| +// . \ | ( | | ( |\__ ` +// _|\_\_| \__,_|\__|\___/ ____/ +// Multi-Physics +// +// License: BSD License +// Kratos default license: kratos/license.txt +// +// Main authors: Nicolo' Antonelli +// Main authors: Andrea Gorgi + +// Project includes +#include "containers/model.h" +#include "testing/testing.h" +#include "custom_utilities/create_breps_sbm_utilities.h" +#include "includes/kratos_parameters.h" +#include "custom_modelers/nurbs_geometry_modeler.h" + +namespace Kratos::Testing +{ + using NurbsSurfaceGeometryType = NurbsSurfaceGeometry<3, PointerVector>; + using NurbsSurfaceGeometryPointerType = typename NurbsSurfaceGeometryType::Pointer; + +// Tests the SnakeSbmUtilities with a square outer geometry +KRATOS_TEST_CASE_IN_SUITE(TestCreateBrepsSbmUtilitiesOuter, KratosIgaFastSuite) +{ + using IndexType = std::size_t; + using GeometryType = Geometry; + using ContainerNodeType = PointerVector; + using ContainerEmbeddedNodeType = PointerVector; + using BrepCurveOnSurfaceType = BrepCurveOnSurface; + using CoordinatesArrayType = GeometryType::CoordinatesArrayType; + + // Initialization + std::size_t EchoLevel = 0; + Model model; + ModelPart& iga_model_part = model.CreateModelPart("IgaModelPart"); + ModelPart& surrogate_sub_model_part_outer = iga_model_part.CreateSubModelPart("surrogate_outer"); + ModelPart& surrogate_sub_model_part_inner = iga_model_part.CreateSubModelPart("surrogate_inner"); + + surrogate_sub_model_part_outer.CreateNewProperties(0); + + surrogate_sub_model_part_outer.CreateNewNode(1, -1.0, -1.0, 0.0); + surrogate_sub_model_part_outer.CreateNewNode(2, -1.0, -0.9, 0.0); + surrogate_sub_model_part_outer.CreateNewNode(3, -1.0, -0.8, 0.0); + + surrogate_sub_model_part_outer.CreateNewNode(4, -0.9, -0.8, 0.0); + surrogate_sub_model_part_outer.CreateNewNode(5, -0.8, -0.8, 0.0); + + surrogate_sub_model_part_outer.CreateNewNode(6, -0.8, -0.9, 0.0); + surrogate_sub_model_part_outer.CreateNewNode(7, -0.8, -1.0, 0.0); + + surrogate_sub_model_part_outer.CreateNewNode(8, -0.9, -1.0, 0.0); + + + Properties::Pointer p_prop = surrogate_sub_model_part_outer.pGetProperties(0); + surrogate_sub_model_part_outer.CreateNewCondition("LineCondition2D2N", 1, {{1, 2}}, p_prop); + surrogate_sub_model_part_outer.CreateNewCondition("LineCondition2D2N", 2, {{2, 3}}, p_prop); + surrogate_sub_model_part_outer.CreateNewCondition("LineCondition2D2N", 3, {{3, 4}}, p_prop); + surrogate_sub_model_part_outer.CreateNewCondition("LineCondition2D2N", 4, {{4, 5}}, p_prop); + surrogate_sub_model_part_outer.CreateNewCondition("LineCondition2D2N", 5, {{5, 6}}, p_prop); + surrogate_sub_model_part_outer.CreateNewCondition("LineCondition2D2N", 6, {{6, 7}}, p_prop); + surrogate_sub_model_part_outer.CreateNewCondition("LineCondition2D2N", 7, {{7, 8}}, p_prop); + surrogate_sub_model_part_outer.CreateNewCondition("LineCondition2D2N", 8, {{8, 1}}, p_prop); + + for (auto &cond : surrogate_sub_model_part_outer.Conditions()) + { + cond.Set(BOUNDARY, true); + } + + Kratos::Parameters parameters(R"( + { + "echo_level": 0, + "model_part_name": "IgaModelPart", + "geometry_name" : "nurbs_surface", + "lower_point_xyz": [-1, -1, 0], + "upper_point_xyz": [1, 1, 0], + "lower_point_uvw": [-1, -1, 0], + "upper_point_uvw": [1, 1, 0], + "polynomial_order": [2, 2], + "number_of_knot_spans": [20, 20] + } + )"); + + const Point& A_uvw = Point(-1, -1, 0); + const Point& B_uvw = Point(1, 1, 0); + + NurbsGeometryModeler nurbs_modeler(model, parameters); + nurbs_modeler.SetupGeometryModel(); + + NurbsSurfaceGeometryPointerType surface = std::dynamic_pointer_cast(iga_model_part.pGetGeometry("nurbs_surface")); + + // Create the breps for the outer sbm boundary + CreateBrepsSbmUtilities CreateBrepsSbmUtilities(EchoLevel); + CreateBrepsSbmUtilities.CreateSurrogateBoundary(surface, surrogate_sub_model_part_inner, surrogate_sub_model_part_outer, A_uvw, B_uvw, iga_model_part); + + const double tolerance = 1e-12; + + for (IndexType current_id = 2; current_id < iga_model_part.NumberOfGeometries(); current_id++) { + + BrepCurveOnSurfaceType::Pointer brep_curve = std::dynamic_pointer_cast(iga_model_part.pGetGeometry(current_id)); + auto ppp = brep_curve->DomainInterval(); + CoordinatesArrayType vertex_1_on_physical; + CoordinatesArrayType vertex_2_on_physical; + CoordinatesArrayType local_coord_vertex_1 = ZeroVector(3); + CoordinatesArrayType local_coord_vertex_2 = ZeroVector(3); + + local_coord_vertex_1[0] = ppp.GetT0(); local_coord_vertex_2[0] = ppp.GetT1(); + + brep_curve->GlobalCoordinates(vertex_1_on_physical, local_coord_vertex_1); + brep_curve->GlobalCoordinates(vertex_2_on_physical, local_coord_vertex_2); + + KRATOS_EXPECT_NEAR(vertex_2_on_physical[0], surrogate_sub_model_part_outer.GetNode(current_id-1).X() , tolerance); + KRATOS_EXPECT_NEAR(vertex_2_on_physical[1], surrogate_sub_model_part_outer.GetNode(current_id-1).Y() , tolerance); + } + + // Ensure the number of nodes matches expectation + KRATOS_EXPECT_EQ(iga_model_part.Geometries().size(), 10); + +} + + +KRATOS_TEST_CASE_IN_SUITE(TestCreateBrepsSbmUtilitiesInner, KratosIgaFastSuite) +{ + using IndexType = std::size_t; + using GeometryType = Geometry; + using ContainerNodeType = PointerVector; + using ContainerEmbeddedNodeType = PointerVector; + using BrepCurveOnSurfaceType = BrepCurveOnSurface; + using CoordinatesArrayType = GeometryType::CoordinatesArrayType; + + // Initialization + std::size_t EchoLevel = 0; + Model model; + ModelPart& iga_model_part = model.CreateModelPart("IgaModelPart"); + ModelPart& surrogate_sub_model_part_outer = iga_model_part.CreateSubModelPart("surrogate_outer"); + ModelPart& surrogate_sub_model_part_inner = iga_model_part.CreateSubModelPart("surrogate_inner"); + + surrogate_sub_model_part_inner.CreateNewProperties(0); + + surrogate_sub_model_part_inner.CreateNewNode(1, -1.0, -1.0, 0.0); + surrogate_sub_model_part_inner.CreateNewNode(2, -1.0, -0.9, 0.0); + surrogate_sub_model_part_inner.CreateNewNode(3, -1.0, -0.8, 0.0); + + surrogate_sub_model_part_inner.CreateNewNode(4, -0.9, -0.8, 0.0); + surrogate_sub_model_part_inner.CreateNewNode(5, -0.8, -0.8, 0.0); + + surrogate_sub_model_part_inner.CreateNewNode(6, -0.8, -0.9, 0.0); + surrogate_sub_model_part_inner.CreateNewNode(7, -0.8, -1.0, 0.0); + + surrogate_sub_model_part_inner.CreateNewNode(8, -0.9, -1.0, 0.0); + + + Properties::Pointer p_prop = surrogate_sub_model_part_inner.pGetProperties(0); + surrogate_sub_model_part_inner.CreateNewCondition("LineCondition2D2N", 1, {{1, 2}}, p_prop); + surrogate_sub_model_part_inner.CreateNewCondition("LineCondition2D2N", 2, {{2, 3}}, p_prop); + surrogate_sub_model_part_inner.CreateNewCondition("LineCondition2D2N", 3, {{3, 4}}, p_prop); + surrogate_sub_model_part_inner.CreateNewCondition("LineCondition2D2N", 4, {{4, 5}}, p_prop); + surrogate_sub_model_part_inner.CreateNewCondition("LineCondition2D2N", 5, {{5, 6}}, p_prop); + surrogate_sub_model_part_inner.CreateNewCondition("LineCondition2D2N", 6, {{6, 7}}, p_prop); + surrogate_sub_model_part_inner.CreateNewCondition("LineCondition2D2N", 7, {{7, 8}}, p_prop); + surrogate_sub_model_part_inner.CreateNewCondition("LineCondition2D2N", 8, {{8,1}}, p_prop); + + for (auto &cond : surrogate_sub_model_part_inner.Conditions()) + { + cond.Set(BOUNDARY, false); + } + + // Create "fictituos element" to memorize starting and ending node id for each surrogate boundary loop + std::vector elem_nodes{1, 8}; + surrogate_sub_model_part_inner.CreateNewElement("Element2D2N", 1, elem_nodes, p_prop); + + Kratos::Parameters parameters(R"( + { + "echo_level": 0, + "model_part_name": "IgaModelPart", + "geometry_name" : "nurbs_surface", + "lower_point_xyz": [-2, -2, 0], + "upper_point_xyz": [0, 0, 0], + "lower_point_uvw": [-2, -2, 0], + "upper_point_uvw": [0, 0, 0], + "polynomial_order": [2, 2], + "number_of_knot_spans": [20, 20] + } + )"); + + const Point& A_uvw = Point(-2, -2, 0); + const Point& B_uvw = Point(0, 0, 0); + + NurbsGeometryModeler nurbs_modeler(model, parameters); + nurbs_modeler.SetupGeometryModel(); + + NurbsSurfaceGeometryPointerType surface = std::dynamic_pointer_cast(iga_model_part.pGetGeometry("nurbs_surface")); + + // Create the breps for the outer sbm boundary + CreateBrepsSbmUtilities CreateBrepsSbmUtilities(EchoLevel); + CreateBrepsSbmUtilities.CreateSurrogateBoundary(surface, surrogate_sub_model_part_inner, surrogate_sub_model_part_outer, A_uvw, B_uvw, iga_model_part); + + const double tolerance = 1e-12; + + for (IndexType current_id = 6; current_id < iga_model_part.NumberOfGeometries(); current_id++) { + + BrepCurveOnSurfaceType::Pointer brep_curve = std::dynamic_pointer_cast(iga_model_part.pGetGeometry(current_id)); + auto ppp = brep_curve->DomainInterval(); + CoordinatesArrayType vertex_1_on_physical; + CoordinatesArrayType vertex_2_on_physical; + CoordinatesArrayType local_coord_vertex_1 = ZeroVector(3); + CoordinatesArrayType local_coord_vertex_2 = ZeroVector(3); + + local_coord_vertex_1[0] = ppp.GetT0(); local_coord_vertex_2[0] = ppp.GetT1(); + + brep_curve->GlobalCoordinates(vertex_1_on_physical, local_coord_vertex_1); + brep_curve->GlobalCoordinates(vertex_2_on_physical, local_coord_vertex_2); + + KRATOS_EXPECT_NEAR(vertex_1_on_physical[0], surrogate_sub_model_part_inner.GetNode(current_id-5).X() , tolerance); + KRATOS_EXPECT_NEAR(vertex_1_on_physical[1], surrogate_sub_model_part_inner.GetNode(current_id-5).Y() , tolerance); + } + + // Ensure the number of nodes matches expectation + KRATOS_EXPECT_EQ(iga_model_part.Geometries().size(), 14); + +} + +} \ No newline at end of file