From 92fd157c38ada95c2f24aa310ebccf6a160e6157 Mon Sep 17 00:00:00 2001 From: andrewgorgi Date: Tue, 14 Jan 2025 11:29:15 +0100 Subject: [PATCH 1/8] Ready for PR --- .../nurbs_geometry_modeler.cpp | 38 +- .../custom_modelers/nurbs_geometry_modeler.h | 33 +- .../nurbs_geometry_modeler_sbm.cpp | 122 +++ .../nurbs_geometry_modeler_sbm.h | 129 +++ .../create_breps_sbm_utilities.h | 406 +++++++++ .../IgaApplication/iga_application.cpp | 1 + applications/IgaApplication/iga_application.h | 2 + .../tests/test_IgaApplication.py | 2 + .../IgaApplication/tests/test_modelers_sbm.py | 285 ++++++ kratos/geometries/brep_surface.h | 22 +- .../brep_trimming_utilities.cpp | 10 +- .../brep_trimming_utilities.h | 7 +- .../nurbs_utilities/snake_sbm_utilities.cpp | 839 ++++++++++++++++++ .../nurbs_utilities/snake_sbm_utilities.h | 179 ++++ 14 files changed, 2045 insertions(+), 30 deletions(-) create mode 100644 applications/IgaApplication/custom_modelers/nurbs_geometry_modeler_sbm.cpp create mode 100644 applications/IgaApplication/custom_modelers/nurbs_geometry_modeler_sbm.h create mode 100644 applications/IgaApplication/custom_utilities/create_breps_sbm_utilities.h create mode 100644 applications/IgaApplication/tests/test_modelers_sbm.py create mode 100644 kratos/utilities/nurbs_utilities/snake_sbm_utilities.cpp create mode 100644 kratos/utilities/nurbs_utilities/snake_sbm_utilities.h diff --git a/applications/IgaApplication/custom_modelers/nurbs_geometry_modeler.cpp b/applications/IgaApplication/custom_modelers/nurbs_geometry_modeler.cpp index 9d51adbf1986..7d6489162a8d 100644 --- a/applications/IgaApplication/custom_modelers/nurbs_geometry_modeler.cpp +++ b/applications/IgaApplication/custom_modelers/nurbs_geometry_modeler.cpp @@ -106,7 +106,7 @@ namespace Kratos ///@name Private Operations ///@{ void NurbsGeometryModeler::CreateAndAddRegularGrid2D( ModelPart& r_model_part, const Point& A_xyz, const Point& B_xyz, - const Point& A_uvw, const Point& B_uvw, SizeType OrderU, SizeType OrderV,SizeType NumKnotSpansU, SizeType NumKnotSpansV) + const Point& A_uvw, const Point& B_uvw, SizeType OrderU, SizeType OrderV,SizeType NumKnotSpansU, SizeType NumKnotSpansV, bool add_surface_to_model_part) { KRATOS_ERROR_IF( B_xyz.X() <= A_xyz.X() || B_xyz.Y() <= A_xyz.Y() ) << "NurbsGeometryModeler: " << "The two Points A_xyz and B_xyz must meet the following requirement: (B_xyz-A_xyz) > (0,0,0). However, (B_xyz-A_xyz)=" << B_xyz-A_xyz << std::endl; @@ -175,22 +175,25 @@ namespace Kratos // Set up weights. Required in case no refinement is performed. Vector WeightsRefined(points.size(),1.0); - - // Add geometry to model part - if( mParameters.Has("geometry_name") ){ - p_surface_geometry->SetId(mParameters["geometry_name"].GetString()); - } else { - const SizeType number_of_geometries = r_model_part.NumberOfGeometries(); - SizeType last_geometry_id = 0; - if( number_of_geometries > 0 ){ - for( auto it = r_model_part.GeometriesBegin(); it!= r_model_part.GeometriesEnd(); ++it){ - last_geometry_id = it->Id(); + + // In some cases there is no need of add the surface geometry, add_surface_to_model_part is true by default + if (add_surface_to_model_part) { + // Add geometry to model part + if( mParameters.Has("geometry_name") ){ + p_surface_geometry->SetId(mParameters["geometry_name"].GetString()); + } else { + const SizeType number_of_geometries = r_model_part.NumberOfGeometries(); + SizeType last_geometry_id = 0; + if( number_of_geometries > 0 ){ + for( auto it = r_model_part.GeometriesBegin(); it!= r_model_part.GeometriesEnd(); ++it){ + last_geometry_id = it->Id(); + } } + p_surface_geometry->SetId(last_geometry_id+1); } - p_surface_geometry->SetId(last_geometry_id+1); + r_model_part.AddGeometry(p_surface_geometry); } - r_model_part.AddGeometry(p_surface_geometry); - + // Perform knot refinement. PointerVector PointsRefined = p_surface_geometry->Points(); if( NumKnotSpansU > 1) { @@ -237,6 +240,13 @@ namespace Kratos p_surface_geometry->SetInternals(PointsRefined, p_surface_geometry->PolynomialDegreeU(), p_surface_geometry->PolynomialDegreeV(), p_surface_geometry->KnotsU(), p_surface_geometry->KnotsV(), WeightsRefined); + + // assign the value p_surface_geometry to the class member mpSurface for derived classes + mpSurface = p_surface_geometry; + mInsertKnotsU = insert_knots_u; + mInsertKnotsV = insert_knots_v; + mKnotVectorU = knot_vector_u; + mKnotVectorV = knot_vector_v; } diff --git a/applications/IgaApplication/custom_modelers/nurbs_geometry_modeler.h b/applications/IgaApplication/custom_modelers/nurbs_geometry_modeler.h index bd97bf17f76c..6496b90a0fbb 100644 --- a/applications/IgaApplication/custom_modelers/nurbs_geometry_modeler.h +++ b/applications/IgaApplication/custom_modelers/nurbs_geometry_modeler.h @@ -82,15 +82,8 @@ class KRATOS_API(IGA_APPLICATION) NurbsGeometryModeler ///@} -private: - ///@name Private Member Variables - ///@{ - - Model* mpModel; - - ///@} - ///@name Private Operations - ///@{ +protected: +///@{ /** * @brief Creates a regular grid composed out of bivariant B-splines. @@ -100,8 +93,26 @@ class KRATOS_API(IGA_APPLICATION) NurbsGeometryModeler * @param NumKnotSpans Number of equidistant elements/knot spans in each direction u,v. * @note The CP'S are defined as nodes and added to the rModelPart. **/ - void CreateAndAddRegularGrid2D( ModelPart& r_model_part, const Point& A_xyz, const Point& B_xyz, const Point& A_uvw, const Point& B_uvw, - SizeType OrderU, SizeType OrderV, SizeType NumKnotSpansU, SizeType NumKnotSpansV ); + virtual void CreateAndAddRegularGrid2D( ModelPart& r_model_part, const Point& A_xyz, const Point& B_xyz, const Point& A_uvw, const Point& B_uvw, + SizeType OrderU, SizeType OrderV, SizeType NumKnotSpansU, SizeType NumKnotSpansV, bool add_surface_to_model_part = true); + + NurbsSurfaceGeometryPointerType mpSurface; + + Vector mKnotVectorU; + Vector mKnotVectorV; + std::vector mInsertKnotsU; + std::vector mInsertKnotsV; + + +private: + ///@name Private Member Variables + ///@{ + + Model* mpModel; + + ///@} + ///@name Private Operations + /** * @brief Creates a cartesian grid composed out of trivariant B-spline cubes. diff --git a/applications/IgaApplication/custom_modelers/nurbs_geometry_modeler_sbm.cpp b/applications/IgaApplication/custom_modelers/nurbs_geometry_modeler_sbm.cpp new file mode 100644 index 000000000000..322f8795c83e --- /dev/null +++ b/applications/IgaApplication/custom_modelers/nurbs_geometry_modeler_sbm.cpp @@ -0,0 +1,122 @@ +// | / | +// ' / __| _` | __| _ \ __| +// . \ | ( | | ( |\__ ` +// _|\_\_| \__,_|\__|\___/ ____/ +// Multi-Physics +// +// License: BSD License +// Kratos default license: kratos/license.txt +// +// Main authors: Nicolo' Antonelli +// Andrea Gorgi +// + +// Project includes +#include "includes/define.h" +#include "nurbs_geometry_modeler_sbm.h" +#include "geometries/nurbs_shape_function_utilities/nurbs_volume_refinement_utilities.h" +#include "custom_utilities/create_breps_sbm_utilities.h" + +namespace Kratos +{ + + ///@name Stages + ///@{ + + void NurbsGeometryModelerSbm::SetupGeometryModel(){ + + //---------------------------------------------------------------------------------------------------------------- + KRATOS_INFO_IF("NurbsGeometryModelerSbm", mEchoLevel > 1) << "[NURBS MODELER SBM]:: calling NurbsGeometryModeler" << std::endl; + + // Call the SetupGeometryModel method of the base class NurbsGeometryModeler + NurbsGeometryModeler::SetupGeometryModel(); + } + + ///@} + ///@name Private Operations + ///@{ + void NurbsGeometryModelerSbm::CreateAndAddRegularGrid2D( ModelPart& r_model_part, const Point& A_xyz, const Point& B_xyz, + const Point& A_uvw, const Point& B_uvw, SizeType OrderU, SizeType OrderV,SizeType NumKnotSpansU, SizeType NumKnotSpansV, bool add_surface_to_model_part) + { + + // Call the CreateAndAddRegularGrid2D method of the base class NurbsGeometryModeler + NurbsGeometryModeler::CreateAndAddRegularGrid2D(r_model_part, A_xyz, B_xyz, + A_uvw, B_uvw, OrderU, OrderV, NumKnotSpansU, NumKnotSpansV, false); + + // Create the Domain/Iga Model Part + const std::string iga_model_part_name = mParameters["model_part_name"].GetString(); + ModelPart& iga_model_part = mpModel->HasModelPart(iga_model_part_name) + ? mpModel->GetModelPart(iga_model_part_name) + : mpModel->CreateModelPart(iga_model_part_name); + + // Create the True Model Part -> contains all the true boundary features + std::string skin_model_part_inner_initial_name = "skin_model_part_inner_initial_name"; + std::string skin_model_part_outer_initial_name = "skin_model_part_outer_initial_name"; + std::string skin_model_part_name; + if (mParameters.Has("skin_model_part_inner_initial_name")) { + skin_model_part_inner_initial_name = mParameters["skin_model_part_inner_initial_name"].GetString(); + + if (!mpModel->HasModelPart(skin_model_part_inner_initial_name)) + KRATOS_ERROR << "The skin_model_part '" << skin_model_part_inner_initial_name << "' was not created in the model.\n" + << "Check the reading of the mdpa file in the import mdpa modeler."<< std::endl; + } + if (mParameters.Has("skin_model_part_outer_initial_name")) { + skin_model_part_outer_initial_name = mParameters["skin_model_part_outer_initial_name"].GetString(); + + if (!mpModel->HasModelPart(skin_model_part_outer_initial_name)) + KRATOS_ERROR << "The skin_model_part '" << skin_model_part_outer_initial_name << "' was not created in the model.\n" + << "Check the reading of the mdpa file in the import mdpa modeler."<< std::endl; + } + // If there is not neither skin_inner nor skin_outer throw an error since you are using the sbm modeler + if (!(mParameters.Has("skin_model_part_inner_initial_name") || mParameters.Has("skin_model_part_outer_initial_name"))){ + KRATOS_ERROR << "None of the 'skin_model_part_name' have not been defined " << + "in the nurbs_geometry_modeler_sbm, define it " << + "in the project paramer json" << std::endl; + } + + if (mParameters.Has("skin_model_part_name")) + skin_model_part_name = mParameters["skin_model_part_name"].GetString(); + else + KRATOS_ERROR << "The skin_model_part name '" << skin_model_part_name << "' was not defined in the project parameters.\n" << std::endl; + + // inner + ModelPart& skin_model_part_inner_initial = mpModel->HasModelPart(skin_model_part_inner_initial_name) + ? mpModel->GetModelPart(skin_model_part_inner_initial_name) + : mpModel->CreateModelPart(skin_model_part_inner_initial_name); + // outer + ModelPart& skin_model_part_outer_initial = mpModel->HasModelPart(skin_model_part_outer_initial_name) + ? mpModel->GetModelPart(skin_model_part_outer_initial_name) + : mpModel->CreateModelPart(skin_model_part_outer_initial_name); + + // Create the surrogate sub model parts inner and outer + ModelPart& surrogate_sub_model_part_inner = iga_model_part.CreateSubModelPart("surrogate_inner"); + ModelPart& surrogate_sub_model_part_outer = iga_model_part.CreateSubModelPart("surrogate_outer"); + + // Skin model part refined after Snake Process + ModelPart& skin_model_part = mpModel->CreateModelPart(skin_model_part_name); + ModelPart& skin_sub_model_part_in = skin_model_part.CreateSubModelPart("inner"); + ModelPart& skin_sub_model_part_out = skin_model_part.CreateSubModelPart("outer"); + + + // compute unique_knot_vector_u + Vector unique_knot_vector_u(2+(NumKnotSpansU-1)); + unique_knot_vector_u[0] = mKnotVectorU[0]; unique_knot_vector_u[NumKnotSpansU] = mKnotVectorU[mKnotVectorU.size()-1]; + for (SizeType i_knot_insertion = 0; i_knot_insertion < NumKnotSpansU-1; i_knot_insertion++) { + unique_knot_vector_u[i_knot_insertion+1] = mInsertKnotsU[i_knot_insertion]; + } + + // compute unique_knot_vector_v + Vector unique_knot_vector_v(2+(NumKnotSpansV-1)); + unique_knot_vector_v[0] = mKnotVectorV[0]; unique_knot_vector_v[NumKnotSpansV] = mKnotVectorV[mKnotVectorV.size()-1]; + + for (SizeType i_knot_insertion = 0; i_knot_insertion < NumKnotSpansV-1; i_knot_insertion++) { + unique_knot_vector_v[i_knot_insertion+1] = mInsertKnotsV[i_knot_insertion]; + } + SnakeSBMUtilities::CreateTheSnakeCoordinates(iga_model_part, skin_model_part_inner_initial, skin_model_part_outer_initial, skin_model_part, mEchoLevel, + unique_knot_vector_u, unique_knot_vector_v, mParameters) ; + // Create the breps for the outer sbm boundary + CreateBrepsSBMUtilities CreateBrepsSBMUtilities(mEchoLevel); + CreateBrepsSBMUtilities.CreateSurrogateBoundary(mpSurface, r_model_part, surrogate_sub_model_part_inner, surrogate_sub_model_part_outer, A_uvw, B_uvw); + } + +} // end namespace kratos diff --git a/applications/IgaApplication/custom_modelers/nurbs_geometry_modeler_sbm.h b/applications/IgaApplication/custom_modelers/nurbs_geometry_modeler_sbm.h new file mode 100644 index 000000000000..92525dc60e04 --- /dev/null +++ b/applications/IgaApplication/custom_modelers/nurbs_geometry_modeler_sbm.h @@ -0,0 +1,129 @@ +// | / | +// ' / __| _` | __| _ \ __| +// . \ | ( | | ( |\__ ` +// _|\_\_| \__,_|\__|\___/ ____/ +// Multi-Physics +// +// License: BSD License +// Kratos default license: kratos/license.txt +// +// Main authors: Nicolo' Antonelli +// Andrea Gorgi +// + +#if !defined(KRATOS_NURBS_GEOMETRY_MODELER_SBM_H_INCLUDED ) +#define KRATOS_NURBS_GEOMETRY_MODELER_SBM_H_INCLUDED + +// System includes + +// External includes + +// Project includes +#include "includes/model_part.h" +#include "nurbs_geometry_modeler.h" +#include "geometries/nurbs_volume_geometry.h" +#include "geometries/nurbs_surface_geometry.h" +#include "geometries/nurbs_shape_function_utilities/nurbs_surface_refinement_utilities.h" +#include "geometries/brep_curve_on_surface.h" +#include "utilities/nurbs_utilities/snake_sbm_utilities.h" + +namespace Kratos { + +class KRATOS_API(IGA_APPLICATION) NurbsGeometryModelerSbm + : public NurbsGeometryModeler +{ +public: + ///@name Type Definitions + ///@{ + KRATOS_CLASS_POINTER_DEFINITION( NurbsGeometryModelerSbm ); + + typedef std::size_t IndexType; + typedef std::size_t SizeType; + typedef Node NodeType; + + typedef Geometry GeometryType; + typedef typename GeometryType::Pointer GeometryPointerType; + + typedef NurbsSurfaceGeometry<3, PointerVector> NurbsSurfaceGeometryType; + typedef typename NurbsSurfaceGeometryType::Pointer NurbsSurfaceGeometryPointerType; + + typedef NurbsVolumeGeometry> NurbsVolumeGeometryType; + typedef typename NurbsVolumeGeometryType::Pointer NurbsVolumeGeometryPointerType; + + typedef PointerVector ContainerNodeType; + typedef PointerVector ContainerEmbeddedNodeType; + + ///@} + ///@name Life Cycle + ///@{ + + /// Default constructor. + NurbsGeometryModelerSbm() + : NurbsGeometryModeler() {} + + /// Constructor. + NurbsGeometryModelerSbm( + Model & rModel, + const Parameters ModelerParameters = Parameters()) + : NurbsGeometryModeler(rModel, ModelerParameters) + , mpModel(&rModel) + { + } + + /// Destructor. + virtual ~NurbsGeometryModelerSbm() = default; + + /// Creates the Modeler Pointer + Modeler::Pointer Create(Model& rModel, const Parameters ModelParameters) const override + { + return Kratos::make_shared(rModel, ModelParameters); + } + + ///@} + ///@name Stages + ///@{ + + void SetupGeometryModel() override; + + ///@} + +protected: + + /** + * @brief Creates a regular grid composed out of bivariant B-splines. + * @param PointA Lower point of bounding box. + * @param PointB Upper point of bounding box. + * @param Order Polynomial degree in each direction u,v. + * @param NumKnotSpans Number of equidistant elements/knot spans in each direction u,v. + * @note The CP'S are defined as nodes and added to the rModelPart. + **/ + void CreateAndAddRegularGrid2D( ModelPart& r_model_part, const Point& A_xyz, const Point& B_xyz, const Point& A_uvw, const Point& B_uvw, + SizeType OrderU, SizeType OrderV, SizeType NumKnotSpansU, SizeType NumKnotSpansV, bool add_surface_to_model_part ) override; + +private: + + ///@name Private Member Variables + ///@{ + + Model* mpModel; + + ///@} + ///@name Private Operations + ///@{ + + /** + * @brief Creates a cartesian grid composed out of trivariant B-spline cubes. + * @param PointA Lower point of bounding box. + * @param PointB Upper point of bounding box. + * @param Order Polynomial degree in each direction u,v,w. + * @param NumKnotSpans Number of equidistant elements/knot spans in each direction u,v,w. + * @note The CP'S are defined as nodes and added to the rModelPart. + **/ + void CreateAndAddRegularGrid3D( ModelPart& r_model_part, const Point& A_xyz, const Point& B_xyz, const Point& A_uvw, const Point& B_uvw, + SizeType OrderU, SizeType OrderV, SizeType OrderW, SizeType NumKnotSpansU, SizeType NumKnotSpansV, SizeType NumKnotSpansW ); + + Parameters ReadParamatersFile(const std::string& rDataFileName) const; +}; + +} // End namesapce Kratos +#endif // KRATOS_NURBS_GEOMETRY_MODELER_H_INCLUDED \ No newline at end of file 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..8e508eea7d30 --- /dev/null +++ b/applications/IgaApplication/custom_utilities/create_breps_sbm_utilities.h @@ -0,0 +1,406 @@ +// | / | +// ' / __| _` | __| _ \ __| +// . \ | ( | | ( |\__ ` +// _|\_\_| \__,_|\__|\___/ ____/ +// Multi-Physics +// +// License: BSD License +// Kratos default license: kratos/license.txt +// +// Main authors: Nicolo' Antonelli +// Andrea Gorgi +// + +#if !defined(KRATOS_CREATE_BREPS_SBM_UTILITIES_INCLUDED ) +#define KRATOS_CREATE_BREPS_SBM_UTILITIES_INCLUDED + + +// System includes + +// External includes + +// Project includes +#include "includes/io.h" +#include "includes/kratos_parameters.h" +#include "includes/model_part.h" + +// Geometries +#include "geometries/coupling_geometry.h" +#include "geometries/point_on_geometry.h" + +#include "geometries/nurbs_curve_geometry.h" +#include "geometries/nurbs_surface_geometry.h" + +#include "geometries/brep_surface.h" +#include "geometries/brep_curve.h" + +#include "geometries/nurbs_curve_on_surface_geometry.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); + + typedef std::size_t SizeType; + typedef std::size_t IndexType; + + typedef Geometry GeometryType; + typedef typename GeometryType::Pointer GeometryPointerType; + + typedef PointerVector ContainerNodeType; + typedef PointerVector ContainerEmbeddedNodeType; + + typedef CouplingGeometry CouplingGeometryType; + + typedef NurbsSurfaceGeometry<3, ContainerNodeType> NurbsSurfaceType; + typedef NurbsCurveGeometry<2, ContainerEmbeddedNodeType> NurbsTrimmingCurveType; + + typedef typename NurbsSurfaceType::Pointer NurbsSurfacePointerType; + typedef typename NurbsTrimmingCurveType::Pointer NurbsTrimmingCurvePointerType; + + typedef NurbsSurfaceGeometry<3, PointerVector> NurbsSurfaceGeometryType; + typedef typename NurbsSurfaceGeometryType::Pointer NurbsSurfaceGeometryPointerType; + + typedef BrepSurface BrepSurfaceType; + typedef BrepCurveOnSurface BrepCurveOnSurfaceType; + + typedef BrepCurve BrepCurveType; + typedef PointOnGeometry PointOnGeometryOnSurfaceType; + typedef PointOnGeometry PointOnGeometryOnCurveType; + + typedef DenseVector BrepCurveOnSurfaceArrayType; + typedef DenseVector BrepCurveOnSurfaceLoopType; + typedef DenseVector> BrepCurveOnSurfaceLoopArrayType; + + ///@} + ///@name Life Cycle + ///@{ + + /// Constructor with path to input file. + CreateBrepsSBMUtilities( + SizeType EchoLevel = 0) + : mEchoLevel(EchoLevel) + { + } + + /// Destructor. + ~CreateBrepsSBMUtilities() = default; + + ///@} + ///@name Python exposed Functions + ///@{ + + /// Adds all CAD geometries to the herin provided model_part. + void CreateSurrogateBoundary(NurbsSurfaceGeometryPointerType& p_surface, ModelPart& rModelPart, ModelPart& rSurrogateModelPart_inner, ModelPart& rSurrogateModelPart_outer, const Point& A_uvw, const Point& B_uvw) + { + CreateBrepSurface(p_surface, rModelPart, rSurrogateModelPart_inner, rSurrogateModelPart_outer, mEchoLevel); + CreateBrepCurveOnSurfaces(p_surface, rModelPart, rSurrogateModelPart_inner, rSurrogateModelPart_outer, A_uvw, B_uvw, mEchoLevel); + } + + +private: + + // 2D + static void CreateBrepSurface( + NurbsSurfaceGeometryPointerType p_surface, + ModelPart& rModelPart, + ModelPart& rSurrogateModelPart_inner, + ModelPart& rSurrogateModelPart_outer, + SizeType EchoLevel = 0) + { + KRATOS_INFO_IF("ReadBrepSurface", (EchoLevel > 3)) + << "Creating BrepSurface \""<< std::endl; + + bool case1 = rSurrogateModelPart_outer.Nodes().size() > 0; + bool case2 = rSurrogateModelPart_inner.Nodes().size() > 0; + BrepCurveOnSurfaceLoopArrayType outer_loops, inner_loops; + + auto p_brep_surface = + Kratos::make_shared( + p_surface, + outer_loops, + inner_loops, + rSurrogateModelPart_inner, + rSurrogateModelPart_outer); + + // Sets the brep as geometry parent of the nurbs surface. + p_surface->SetGeometryParent(p_brep_surface.get()); + + SizeType last_geometry_id = rModelPart.GetParentModelPart().Geometries().size(); + p_brep_surface->SetId(1); + rModelPart.AddGeometry(p_brep_surface); + + } + + + static void CreateBrepCurveOnSurfaces( + NurbsSurfaceGeometryPointerType p_surface, + ModelPart& rModelPart, + ModelPart& rSurrogateModelPart_inner, + ModelPart& rSurrogateModelPart_outer, + const Point& A_uvw, const Point& B_uvw, + SizeType EchoLevel = 0) { + + // OUTER : + // Each element in the surrogate_model_part represents a surrogate boundary loop. First "node" is the initial ID of the first surrogate node and + // the second "node" is the last surrogate node of that loop. (We have done this in the case we have multiple surrogate boundaries and 1 model part) + + int Id_brep_curve_on_surface = 2; + + if (rSurrogateModelPart_outer.Nodes().size() > 0) { + int sizeSurrogateLoop_outer = rSurrogateModelPart_outer.Nodes().size(); + std::vector surrogatecoord_x_outer(sizeSurrogateLoop_outer); + std::vector surrogatecoord_y_outer(sizeSurrogateLoop_outer); + int countSurrogateLoop_outer = 0; + for (auto i_node = rSurrogateModelPart_outer.NodesEnd()-1; i_node != rSurrogateModelPart_outer.NodesBegin()-1; i_node--) { + surrogatecoord_x_outer[countSurrogateLoop_outer] = i_node->X(); + surrogatecoord_y_outer[countSurrogateLoop_outer] = i_node->Y(); + countSurrogateLoop_outer++; + } + // OUTER + std::vector>::Pointer> trimming_curves_GPT_outer; + for (std::size_t i = 0; i < surrogatecoord_x_outer.size(); ++i) { + Vector active_range_knot_vector = ZeroVector(2); + + Point::Pointer point1 = Kratos::make_shared(surrogatecoord_x_outer[i], surrogatecoord_y_outer[i], 0.0); + Point::Pointer point2 = Kratos::make_shared( + surrogatecoord_x_outer[(i + 1) % surrogatecoord_x_outer.size()], // Wrap around for the last point + surrogatecoord_y_outer[(i + 1) % surrogatecoord_y_outer.size()], // Wrap around for the last point + 0.0); + // Compute the knot vector needed + if (surrogatecoord_x_outer[(i + 1) % surrogatecoord_x_outer.size()]==surrogatecoord_x_outer[i]) { + active_range_knot_vector[0] = surrogatecoord_y_outer[i]; + active_range_knot_vector[1] = surrogatecoord_y_outer[(i + 1) % surrogatecoord_y_outer.size()]; + } + else { + active_range_knot_vector[0] = surrogatecoord_x_outer[i]; + active_range_knot_vector[1] = surrogatecoord_x_outer[(i + 1) % surrogatecoord_x_outer.size()]; + } + //// Order the active_range_knot_vector + if (active_range_knot_vector[0] > active_range_knot_vector[1]) { + double temp = active_range_knot_vector[1]; + active_range_knot_vector[1] = active_range_knot_vector[0] ; + active_range_knot_vector[0] = temp ; + } + NurbsCurveGeometry<2, PointerVector>::Pointer p_trimming_curve = CreateSingleBrep(point1, point2, active_range_knot_vector); + trimming_curves_GPT_outer.push_back(p_trimming_curve); + } + + BrepCurveOnSurfaceLoopType trimming_brep_curve_vector_outer(surrogatecoord_x_outer.size()); + + for (std::size_t i = 0; i < trimming_curves_GPT_outer.size(); ++i) { + Vector active_range_vector = ZeroVector(2); + if (surrogatecoord_x_outer[(i + 1) % surrogatecoord_x_outer.size()]==surrogatecoord_x_outer[i]) { + active_range_vector[0] = surrogatecoord_y_outer[i]; + active_range_vector[1] = surrogatecoord_y_outer[(i + 1) % surrogatecoord_x_outer.size()]; + } + else { + active_range_vector[0] = surrogatecoord_x_outer[i]; + active_range_vector[1] = surrogatecoord_x_outer[(i + 1) % surrogatecoord_x_outer.size()]; + } + bool curve_direction = true; + + // Re-order + if (active_range_vector[0] > active_range_vector[1]) { + double temp = active_range_vector[1]; + active_range_vector[1] = active_range_vector[0] ; + active_range_vector[0] = temp ; + } + + NurbsInterval brep_active_range(active_range_vector[0], active_range_vector[1]); + + auto p_brep_curve_on_surface = Kratos::make_shared( + p_surface, trimming_curves_GPT_outer[i], 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(rModelPart, p_surface, A_uvw, B_uvw, Id_brep_curve_on_surface); + } + + // INNER + for (int iel = 1; iel < rSurrogateModelPart_inner.Elements().size()+1; iel++) { + int firstSurrogateNodeId = rSurrogateModelPart_inner.pGetElement(iel)->GetGeometry()[0].Id(); // Element 1 because is the only surrogate loop + int lastSurrogateNodeId = rSurrogateModelPart_inner.pGetElement(iel)->GetGeometry()[1].Id(); // Element 1 because is the only surrogate loop + int sizeSurrogateLoop = lastSurrogateNodeId - firstSurrogateNodeId + 1; + + int countSurrogateLoop = 0; + std::vector surrogatecoord_x(sizeSurrogateLoop); + std::vector surrogatecoord_y(sizeSurrogateLoop); + for (int id_node = firstSurrogateNodeId; id_node < lastSurrogateNodeId+1; id_node++) { + surrogatecoord_x[countSurrogateLoop] = rSurrogateModelPart_inner.GetNode(id_node).X(); + surrogatecoord_y[countSurrogateLoop] = rSurrogateModelPart_inner.GetNode(id_node).Y(); + countSurrogateLoop++; + } + std::vector>::Pointer> trimming_curves_GPT; + + for (std::size_t i = 0; i < surrogatecoord_x.size(); ++i) { + Vector active_range_knot_vector = ZeroVector(2); + + Point::Pointer point1 = Kratos::make_shared(surrogatecoord_x[i], surrogatecoord_y[i], 0.0); + Point::Pointer point2 = Kratos::make_shared( + surrogatecoord_x[(i + 1) % surrogatecoord_x.size()], // Wrap around for the last point + surrogatecoord_y[(i + 1) % surrogatecoord_y.size()], // Wrap around for the last point + 0.0); + // Compute the knot vector needed + if (surrogatecoord_x[(i + 1) % surrogatecoord_x.size()]==surrogatecoord_x[i]) { + active_range_knot_vector[0] = surrogatecoord_y[i]; + active_range_knot_vector[1] = surrogatecoord_y[(i + 1) % surrogatecoord_y.size()]; + } + else { + active_range_knot_vector[0] = surrogatecoord_x[i]; + active_range_knot_vector[1] = surrogatecoord_x[(i + 1) % surrogatecoord_x.size()]; + } + // Order the active_range_knot_vector + if (active_range_knot_vector[0] > active_range_knot_vector[1]) { + double temp = active_range_knot_vector[1]; + active_range_knot_vector[1] = active_range_knot_vector[0] ; + active_range_knot_vector[0] = temp ; + } + NurbsCurveGeometry<2, PointerVector>::Pointer p_trimming_curve = CreateSingleBrep(point1, point2, active_range_knot_vector); + trimming_curves_GPT.push_back(p_trimming_curve); + } + + BrepCurveOnSurfaceLoopType trimming_brep_curve_vector(surrogatecoord_x.size()); + + + for (std::size_t i = 0; i < trimming_curves_GPT.size(); ++i) { + Vector active_range_vector = ZeroVector(2); + if (surrogatecoord_x[(i + 1) % surrogatecoord_x.size()]==surrogatecoord_x[i]) { + active_range_vector[0] = surrogatecoord_y[i]; + active_range_vector[1] = surrogatecoord_y[(i + 1) % surrogatecoord_x.size()]; + } + else { + active_range_vector[0] = surrogatecoord_x[i]; + active_range_vector[1] = surrogatecoord_x[(i + 1) % surrogatecoord_x.size()]; + } + bool curve_direction = true; + + // re-order + if (active_range_vector[0] > active_range_vector[1]) { + double temp = active_range_vector[1]; + active_range_vector[1] = active_range_vector[0] ; + active_range_vector[0] = temp ; + } + + NurbsInterval brep_active_range(active_range_vector[0], active_range_vector[1]); + + auto p_brep_curve_on_surface = Kratos::make_shared( + p_surface, trimming_curves_GPT[i], 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++; + trimming_brep_curve_vector[i] = p_brep_curve_on_surface ; + + } + + } + } + + + + ///@} + ///@name Utility functions + ///@{ + + int mEchoLevel; + + static typename NurbsCurveGeometry<2, PointerVector>::Pointer CreateSingleBrep( + Point::Pointer point1, Point::Pointer point2, Vector active_range_knot_vector) + { + // Create the data for the trimming curves + PointerVector control_points; + control_points.push_back(point1); + control_points.push_back(point2); + int polynomial_degree = 1; + Vector knot_vector = ZeroVector(4) ; + knot_vector[0] = active_range_knot_vector[0] ; + knot_vector[1] = active_range_knot_vector[0] ; + knot_vector[2] = active_range_knot_vector[1] ; + knot_vector[3] = active_range_knot_vector[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; + } + + + + static void CreateBrepCurvesOnRectangle(ModelPart& r_model_part, NurbsSurfaceGeometryPointerType p_surface_geometry, const Point& A_uvw, const Point& B_uvw, int &last_geometry_id) { + Vector knot_vector = ZeroVector(2); + knot_vector[0] = 0.0; + knot_vector[1] = std::abs(B_uvw[0] - A_uvw[0]); + int p = 1; + + NurbsCurveGeometry<2, PointerVector>::PointsArrayType segment1; + segment1.push_back(Point::Pointer(new Point(A_uvw[0], A_uvw[1]))); + segment1.push_back(Point::Pointer(new Point(B_uvw[0], A_uvw[1]))); + + NurbsCurveGeometry<2, PointerVector>::PointsArrayType segment2; + segment2.push_back(Point::Pointer(new Point(B_uvw[0], A_uvw[1]))); + segment2.push_back(Point::Pointer(new Point(B_uvw[0], B_uvw[1]))); + + NurbsCurveGeometry<2, PointerVector>::PointsArrayType segment3; + segment3.push_back(Point::Pointer(new Point(B_uvw[0], B_uvw[1]))); + segment3.push_back(Point::Pointer(new Point(A_uvw[0], B_uvw[1]))); + + NurbsCurveGeometry<2, PointerVector>::PointsArrayType segment4; + segment4.push_back(Point::Pointer(new Point(A_uvw[0], B_uvw[1]))); + segment4.push_back(Point::Pointer(new Point(A_uvw[0], A_uvw[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>(p_surface_geometry, p_curve_1); + auto p_brep_curve_on_surface = Kratos::make_shared(p_surface_geometry, p_curve_1); + p_brep_curve_on_surface->SetId(last_geometry_id); + r_model_part.AddGeometry(p_brep_curve_on_surface); + + brep_curve_on_surface = BrepCurveOnSurface< PointerVector, true, PointerVector>(p_surface_geometry, p_curve_2); + p_brep_curve_on_surface = Kratos::make_shared(p_surface_geometry, p_curve_2); + p_brep_curve_on_surface->SetId(++last_geometry_id); + r_model_part.AddGeometry(p_brep_curve_on_surface); + + brep_curve_on_surface = BrepCurveOnSurface< PointerVector, true, PointerVector>(p_surface_geometry, p_curve_3); + p_brep_curve_on_surface = Kratos::make_shared(p_surface_geometry, p_curve_3); + p_brep_curve_on_surface->SetId(++last_geometry_id); + r_model_part.AddGeometry(p_brep_curve_on_surface); + + brep_curve_on_surface = BrepCurveOnSurface< PointerVector, true, PointerVector>(p_surface_geometry, p_curve_4); + p_brep_curve_on_surface = Kratos::make_shared(p_surface_geometry, p_curve_4); + p_brep_curve_on_surface->SetId(++last_geometry_id); + r_model_part.AddGeometry(p_brep_curve_on_surface); + + last_geometry_id++; + } + + ///@} +}; // Class CreateBrepsSBMUtilities +} // namespace Kratos. + +#endif diff --git a/applications/IgaApplication/iga_application.cpp b/applications/IgaApplication/iga_application.cpp index e7465c6e28cb..2c3fdf636892 100644 --- a/applications/IgaApplication/iga_application.cpp +++ b/applications/IgaApplication/iga_application.cpp @@ -90,6 +90,7 @@ KRATOS_INFO("") << " KRATOS _____ _____\n" KRATOS_REGISTER_MODELER("IgaModeler", mIgaModeler); KRATOS_REGISTER_MODELER("RefinementModeler", mRefinementModeler); KRATOS_REGISTER_MODELER("NurbsGeometryModeler", mNurbsGeometryModeler); + KRATOS_REGISTER_MODELER("NurbsGeometryModelerSbm", mNurbsGeometryModelerSbm); // VARIABLES KRATOS_REGISTER_VARIABLE(CROSS_AREA) diff --git a/applications/IgaApplication/iga_application.h b/applications/IgaApplication/iga_application.h index 59fa1ad32e55..f543fd51d992 100644 --- a/applications/IgaApplication/iga_application.h +++ b/applications/IgaApplication/iga_application.h @@ -43,6 +43,7 @@ #include "custom_modelers/iga_modeler.h" #include "custom_modelers/refinement_modeler.h" #include "custom_modelers/nurbs_geometry_modeler.h" +#include "custom_modelers/nurbs_geometry_modeler_sbm.h" namespace Kratos { @@ -141,6 +142,7 @@ class KRATOS_API(IGA_APPLICATION) KratosIgaApplication : public KratosApplicatio const IgaModeler mIgaModeler; const RefinementModeler mRefinementModeler; const NurbsGeometryModeler mNurbsGeometryModeler; + const NurbsGeometryModelerSbm mNurbsGeometryModelerSbm; ///@} ///@name Private methods diff --git a/applications/IgaApplication/tests/test_IgaApplication.py b/applications/IgaApplication/tests/test_IgaApplication.py index 6df962b14fc2..82efb0db1550 100644 --- a/applications/IgaApplication/tests/test_IgaApplication.py +++ b/applications/IgaApplication/tests/test_IgaApplication.py @@ -45,6 +45,7 @@ from test_nurbs_volume_element import TestNurbsVolumeElement as TTestNurbsVolumeElements # Modelers tests from test_modelers import TestModelers as TTestModelers +from test_modelers_sbm import TestModelers as TTestModelersSbm # Processes tests from test_map_nurbs_volume_results_to_embedded_geometry_process import TestMapNurbsVolumeResultsToEmbeddedGeometryProcess as TTestMapNurbsVolumeResultsToEmbeddedGeometryProcess @@ -93,6 +94,7 @@ def AssembleTestSuites(): TTestNurbsVolumeElements, # Modelers TTestModelers, + TTestModelersSbm, TTestMapNurbsVolumeResultsToEmbeddedGeometryProcess ])) diff --git a/applications/IgaApplication/tests/test_modelers_sbm.py b/applications/IgaApplication/tests/test_modelers_sbm.py new file mode 100644 index 000000000000..c37efcd6ce7c --- /dev/null +++ b/applications/IgaApplication/tests/test_modelers_sbm.py @@ -0,0 +1,285 @@ +import KratosMultiphysics +import KratosMultiphysics.IgaApplication +import KratosMultiphysics.KratosUnittest as KratosUnittest + +def run_modelers(current_model, modelers_list): + from KratosMultiphysics.modeler_factory import KratosModelerFactory + factory = KratosModelerFactory() + list_of_modelers = factory.ConstructListOfModelers(current_model, modelers_list) + + for modeler in list_of_modelers: + modeler.SetupGeometryModel() + + for modeler in list_of_modelers: + modeler.PrepareGeometryModel() + + for modeler in list_of_modelers: + modeler.SetupModelPart() + +class TestModelersSbm(KratosUnittest.TestCase): + + # test for SBM + def test_nurbs_geometry_2d_modeler_sbm_outer_skin_boundary(self): + current_model = KratosMultiphysics.Model() + modeler_settings = KratosMultiphysics.Parameters(""" + [{ + "modeler_name": "NurbsGeometryModelerSbm", + "Parameters": { + "model_part_name" : "IgaModelPart", + "lower_point_xyz": [0.0,0.0,0.0], + "upper_point_xyz": [2.0,2.0,0.0], + "lower_point_uvw": [0.0,0.0,0.0], + "upper_point_uvw": [2.0,2.0,0.0], + "polynomial_order" : [1, 1], + "number_of_knot_spans" : [5,5], + "sbm_parameters": { + "lambda_outer": 0.5, + "number_of_inner_loops": 0 + }, + "skin_model_part_outer_initial_name": "skinModelPart_outer_initial", + "skin_model_part_name": "skinModelPart" + } + }] + """) + # New model + current_model = KratosMultiphysics.Model() + skin_model_part_outer_initial = current_model.CreateModelPart("skinModelPart_outer_initial") + + skin_model_part_outer_initial.CreateNewProperties(1) + + skin_model_part_outer_initial.CreateNewNode(1, 0.0, 0.0, 0.0) + skin_model_part_outer_initial.CreateNewNode(2, 2.0, 0.0, 0.0) + skin_model_part_outer_initial.CreateNewNode(3, 2.0, 2.0, 0.0) + skin_model_part_outer_initial.CreateNewNode(4, 0.0, 2.0, 0.0) + + skin_model_part_outer_initial.CreateNewCondition("LineCondition2D2N", 1, [1, 2], skin_model_part_outer_initial.GetProperties()[1]) + skin_model_part_outer_initial.CreateNewCondition("LineCondition2D2N", 2, [2, 3], skin_model_part_outer_initial.GetProperties()[1]) + skin_model_part_outer_initial.CreateNewCondition("LineCondition2D2N", 3, [3, 4], skin_model_part_outer_initial.GetProperties()[1]) + skin_model_part_outer_initial.CreateNewCondition("LineCondition2D2N", 4, [4, 1], skin_model_part_outer_initial.GetProperties()[1]) + + run_modelers(current_model, modeler_settings) + + # test surrogate nodes + self.assertEqual(current_model["IgaModelPart.surrogate_outer"].GetNodes()[43].X, 0.4) + self.assertEqual(current_model["IgaModelPart.surrogate_outer"].GetNodes()[43].Y, 2.0) + self.assertEqual(current_model["IgaModelPart.surrogate_outer"].GetNodes()[43].Z, 0.0) + + self.assertEqual(current_model["IgaModelPart.surrogate_outer"].GetNodes()[51].X, 2.0) + self.assertEqual(current_model["IgaModelPart.surrogate_outer"].GetNodes()[51].Y, 0.4) + self.assertEqual(current_model["IgaModelPart.surrogate_outer"].GetNodes()[51].Z, 0.0) + + # test the creation of the breps + self.assertEqual(current_model["IgaModelPart"].NumberOfGeometries(), 21) + self.assertAlmostEqual(current_model["IgaModelPart"].GetGeometry(1).Center().X, 1.0) + self.assertAlmostEqual(current_model["IgaModelPart"].GetGeometry(1).Center().Y, 1.0) + self.assertAlmostEqual(current_model["IgaModelPart"].GetGeometry(1).Center().Z, 0.0) + + + # inner boundary single inner loop + def test_nurbs_geometry_2d_modeler_sbm_inner_skin_boundary(self): + current_model = KratosMultiphysics.Model() + modeler_settings = KratosMultiphysics.Parameters(""" + [{ + "modeler_name": "NurbsGeometryModelerSbm", + "Parameters": { + "model_part_name" : "IgaModelPart", + "lower_point_xyz": [0.0,0.0,0.0], + "upper_point_xyz": [4.0,6.0,0.0], + "lower_point_uvw": [0.0,0.0,0.0], + "upper_point_uvw": [4.0,6.0,0.0], + "polynomial_order" : [2, 2], + "number_of_knot_spans" : [6,9], + "sbm_parameters": { + "lambda_inner": 0.5, + "number_of_inner_loops": 1 + }, + "skin_model_part_inner_initial_name": "skinModelPart_inner_initial", + "skin_model_part_name": "skinModelPart" + } + }] + """) + # New model + current_model = KratosMultiphysics.Model() + skin_model_part_inner_initial = current_model.CreateModelPart("skinModelPart_inner_initial") + + skin_model_part_inner_initial.CreateNewProperties(1) + + skin_model_part_inner_initial.CreateNewNode(1, 0.5, 0.5, 0.0) + skin_model_part_inner_initial.CreateNewNode(2, 2.0, 0.5, 0.0) + skin_model_part_inner_initial.CreateNewNode(3, 2.0, 3.5, 0.0) + skin_model_part_inner_initial.CreateNewNode(4, 0.5, 3.5, 0.0) + + skin_model_part_inner_initial.CreateNewCondition("LineCondition2D2N", 1, [1, 2], skin_model_part_inner_initial.GetProperties()[1]) + skin_model_part_inner_initial.CreateNewCondition("LineCondition2D2N", 2, [2, 3], skin_model_part_inner_initial.GetProperties()[1]) + skin_model_part_inner_initial.CreateNewCondition("LineCondition2D2N", 3, [3, 4], skin_model_part_inner_initial.GetProperties()[1]) + skin_model_part_inner_initial.CreateNewCondition("LineCondition2D2N", 4, [4, 1], skin_model_part_inner_initial.GetProperties()[1]) + + run_modelers(current_model, modeler_settings) + + # test surrogate nodes + self.assertEqual(current_model["IgaModelPart.surrogate_inner"].GetNodes()[92].X, 2/3) + self.assertEqual(current_model["IgaModelPart.surrogate_inner"].GetNodes()[92].Y, 8/3) + self.assertEqual(current_model["IgaModelPart.surrogate_inner"].GetNodes()[92].Z, 0.0) + + self.assertEqual(current_model["IgaModelPart.surrogate_inner"].GetNodes()[98].X, 2.0) + self.assertEqual(current_model["IgaModelPart.surrogate_inner"].GetNodes()[98].Y, 4/3) + self.assertEqual(current_model["IgaModelPart.surrogate_inner"].GetNodes()[98].Z, 0.0) + + # test the creation of the breps + self.assertEqual(current_model["IgaModelPart"].NumberOfGeometries(), 17) + self.assertAlmostEqual(current_model["IgaModelPart"].GetGeometry(1).Center().X, 2.0) + self.assertAlmostEqual(current_model["IgaModelPart"].GetGeometry(1).Center().Y, 3.0) + self.assertAlmostEqual(current_model["IgaModelPart"].GetGeometry(1).Center().Z, 0.0) + + + # inner boundary double inner loops + def test_nurbs_geometry_2d_modeler_sbm_two_inner_skin_boundary(self): + current_model = KratosMultiphysics.Model() + modeler_settings = KratosMultiphysics.Parameters(""" + [{ + "modeler_name": "NurbsGeometryModelerSbm", + "Parameters": { + "model_part_name" : "IgaModelPart", + "lower_point_xyz": [0.0,0.0,0.0], + "upper_point_xyz": [4.0,6.0,0.0], + "lower_point_uvw": [0.0,0.0,0.0], + "upper_point_uvw": [4.0,6.0,0.0], + "polynomial_order" : [2, 2], + "number_of_knot_spans" : [10,15], + "sbm_parameters": { + "lambda_inner": 0.5, + "number_of_inner_loops": 2 + }, + "skin_model_part_inner_initial_name": "skinModelPart_inner_initial", + "skin_model_part_name": "skinModelPart" + } + }] + """) + + current_model = KratosMultiphysics.Model() + + skin_model_part_inner_initial = current_model.CreateModelPart("skinModelPart_inner_initial") + skin_model_part_inner_initial.CreateNewProperties(1) + + # object 1 + skin_model_part_inner_initial.CreateNewNode(1, 0.5, 0.5, 0.0) + skin_model_part_inner_initial.CreateNewNode(2, 2.0, 0.5, 0.0) + skin_model_part_inner_initial.CreateNewNode(3, 2.0, 3.5, 0.0) + skin_model_part_inner_initial.CreateNewNode(4, 0.5, 3.5, 0.0) + + skin_model_part_inner_initial.CreateNewCondition("LineCondition2D2N", 1, [1, 2], skin_model_part_inner_initial.GetProperties()[1]) + skin_model_part_inner_initial.CreateNewCondition("LineCondition2D2N", 2, [2, 3], skin_model_part_inner_initial.GetProperties()[1]) + skin_model_part_inner_initial.CreateNewCondition("LineCondition2D2N", 3, [3, 4], skin_model_part_inner_initial.GetProperties()[1]) + skin_model_part_inner_initial.CreateNewCondition("LineCondition2D2N", 4, [4, 1], skin_model_part_inner_initial.GetProperties()[1]) + + # object 2 + skin_model_part_inner_initial.CreateNewNode(5, 1.5, 4.6, 0.0) + skin_model_part_inner_initial.CreateNewNode(6, 3.7, 4.7, 0.0) + skin_model_part_inner_initial.CreateNewNode(7, 3.0, 5.5, 0.0) + + skin_model_part_inner_initial.CreateNewCondition("LineCondition2D2N", 5, [5, 6], skin_model_part_inner_initial.GetProperties()[1]) + skin_model_part_inner_initial.CreateNewCondition("LineCondition2D2N", 6, [6, 7], skin_model_part_inner_initial.GetProperties()[1]) + skin_model_part_inner_initial.CreateNewCondition("LineCondition2D2N", 7, [7, 5], skin_model_part_inner_initial.GetProperties()[1]) + + run_modelers(current_model, modeler_settings) + + # test surrogate nodes + self.assertAlmostEqual(current_model["IgaModelPart.surrogate_inner"].GetNodes()[226].X, 1.6) + self.assertAlmostEqual(current_model["IgaModelPart.surrogate_inner"].GetNodes()[226].Y, 0.4) + self.assertAlmostEqual(current_model["IgaModelPart.surrogate_inner"].GetNodes()[226].Z, 0.0) + + self.assertAlmostEqual(current_model["IgaModelPart.surrogate_inner"].GetNodes()[235].X, 3.2) + self.assertAlmostEqual(current_model["IgaModelPart.surrogate_inner"].GetNodes()[235].Y, 5.2) + self.assertAlmostEqual(current_model["IgaModelPart.surrogate_inner"].GetNodes()[235].Z, 0.0) + + self.assertAlmostEqual(current_model["IgaModelPart.surrogate_inner"].GetNodes()[213].X, 0.4) + self.assertAlmostEqual(current_model["IgaModelPart.surrogate_inner"].GetNodes()[213].Y, 3.6) + self.assertAlmostEqual(current_model["IgaModelPart.surrogate_inner"].GetNodes()[213].Z, 0.0) + + # test the creation of the breps + self.assertEqual(current_model["IgaModelPart"].NumberOfGeometries(), 41) + self.assertAlmostEqual(current_model["IgaModelPart"].GetGeometry(1).Center().X, 2.0) + self.assertAlmostEqual(current_model["IgaModelPart"].GetGeometry(1).Center().Y, 3.0) + self.assertAlmostEqual(current_model["IgaModelPart"].GetGeometry(1).Center().Z, 0.0) + + + # inner boundary single inner loop + outer + def test_nurbs_geometry_2d_modeler_sbm_inner_outer_skin_boundary(self): + current_model = KratosMultiphysics.Model() + modeler_settings = KratosMultiphysics.Parameters(""" + [{ + "modeler_name": "NurbsGeometryModelerSbm", + "Parameters": { + "model_part_name" : "IgaModelPart", + "lower_point_xyz": [0.0,0.0,0.0], + "upper_point_xyz": [4.0,6.0,0.0], + "lower_point_uvw": [0.0,0.0,0.0], + "upper_point_uvw": [4.0,6.0,0.0], + "polynomial_order" : [2, 2], + "number_of_knot_spans" : [20,10], + "sbm_parameters": { + "lambda_inner": 1.0, + "lambda_outer": 0.5, + "number_of_inner_loops": 1 + }, + "skin_model_part_inner_initial_name": "skinModelPart_inner_initial", + "skin_model_part_outer_initial_name": "skinModelPart_outer_initial", + "skin_model_part_name": "skinModelPart" + } + }] + """) + current_model = KratosMultiphysics.Model() + + # Create Inner skin boundary + + skin_model_part_inner_initial = current_model.CreateModelPart("skinModelPart_inner_initial") + skin_model_part_inner_initial.CreateNewProperties(1) + skin_model_part_inner_initial.CreateNewNode(1, 1.5, 0.5, 0.0) + skin_model_part_inner_initial.CreateNewNode(2, 3.0, 0.5, 0.0) + skin_model_part_inner_initial.CreateNewNode(3, 3.0, 3.5, 0.0) + skin_model_part_inner_initial.CreateNewNode(4, 1.5, 3.5, 0.0) + skin_model_part_inner_initial.CreateNewCondition("LineCondition2D2N", 1, [1, 2], skin_model_part_inner_initial.GetProperties()[1]) + skin_model_part_inner_initial.CreateNewCondition("LineCondition2D2N", 2, [2, 3], skin_model_part_inner_initial.GetProperties()[1]) + skin_model_part_inner_initial.CreateNewCondition("LineCondition2D2N", 3, [3, 4], skin_model_part_inner_initial.GetProperties()[1]) + skin_model_part_inner_initial.CreateNewCondition("LineCondition2D2N", 4, [4, 1], skin_model_part_inner_initial.GetProperties()[1]) + + # Create Outer skin boundary + skin_model_part_outer_initial = current_model.CreateModelPart("skinModelPart_outer_initial") + skin_model_part_outer_initial.CreateNewProperties(1) + skin_model_part_outer_initial.CreateNewNode(1, 0.4, 0.0, 0.0) + skin_model_part_outer_initial.CreateNewNode(2, 3.5, 0.0, 0.0) + skin_model_part_outer_initial.CreateNewNode(3, 3.2, 5.0, 0.0) + skin_model_part_outer_initial.CreateNewNode(4, 1.0, 5.5, 0.0) + skin_model_part_outer_initial.CreateNewCondition("LineCondition2D2N", 1, [1, 2], skin_model_part_outer_initial.GetProperties()[1]) + skin_model_part_outer_initial.CreateNewCondition("LineCondition2D2N", 2, [2, 3], skin_model_part_outer_initial.GetProperties()[1]) + skin_model_part_outer_initial.CreateNewCondition("LineCondition2D2N", 3, [3, 4], skin_model_part_outer_initial.GetProperties()[1]) + skin_model_part_outer_initial.CreateNewCondition("LineCondition2D2N", 4, [4, 1], skin_model_part_outer_initial.GetProperties()[1]) + + run_modelers(current_model, modeler_settings) + + # test surrogate nodes + self.assertAlmostEqual(current_model["IgaModelPart.surrogate_inner"].GetNodes()[281].X, 2.8) + self.assertAlmostEqual(current_model["IgaModelPart.surrogate_inner"].GetNodes()[281].Y, 0.6) + self.assertAlmostEqual(current_model["IgaModelPart.surrogate_inner"].GetNodes()[281].Z, 0.0) + + self.assertAlmostEqual(current_model["IgaModelPart.surrogate_inner"].GetNodes()[269].X, 1.6) + self.assertAlmostEqual(current_model["IgaModelPart.surrogate_inner"].GetNodes()[269].Y, 3.0) + self.assertAlmostEqual(current_model["IgaModelPart.surrogate_inner"].GetNodes()[269].Z, 0.0) + + self.assertAlmostEqual(current_model["IgaModelPart.surrogate_outer"].GetNodes()[295].X, 0.8) + self.assertAlmostEqual(current_model["IgaModelPart.surrogate_outer"].GetNodes()[295].Y, 3.6) + self.assertAlmostEqual(current_model["IgaModelPart.surrogate_outer"].GetNodes()[295].Z, 0.0) + + self.assertAlmostEqual(current_model["IgaModelPart.surrogate_outer"].GetNodes()[320].X, 3.6) + self.assertAlmostEqual(current_model["IgaModelPart.surrogate_outer"].GetNodes()[320].Y, 0.6) + self.assertAlmostEqual(current_model["IgaModelPart.surrogate_outer"].GetNodes()[320].Z, 0.0) + + # test the creation of the breps + self.assertEqual(current_model["IgaModelPart"].NumberOfGeometries(), 73) + self.assertAlmostEqual(current_model["IgaModelPart"].GetGeometry(1).Center().X, 2.0) + self.assertAlmostEqual(current_model["IgaModelPart"].GetGeometry(1).Center().Y, 3.0) + self.assertAlmostEqual(current_model["IgaModelPart"].GetGeometry(1).Center().Z, 0.0) + + +if __name__ == '__main__': + KratosUnittest.main() diff --git a/kratos/geometries/brep_surface.h b/kratos/geometries/brep_surface.h index af74c470e587..ee48fe8e41a0 100644 --- a/kratos/geometries/brep_surface.h +++ b/kratos/geometries/brep_surface.h @@ -63,6 +63,7 @@ class BrepSurface typedef NurbsSurfaceGeometry<3, TContainerPointType> NurbsSurfaceType; typedef BrepCurveOnSurface BrepCurveOnSurfaceType; + typedef BrepTrimmingUtilities BrepTrimmingUtilitiesType; typedef DenseVector BrepCurveOnSurfaceArrayType; typedef DenseVector BrepCurveOnSurfaceLoopType; @@ -117,6 +118,22 @@ class BrepSurface { } + // Constructor for SBM + BrepSurface( + typename NurbsSurfaceType::Pointer pSurface, + BrepCurveOnSurfaceLoopArrayType& BrepOuterLoopArray, + BrepCurveOnSurfaceLoopArrayType& BrepInnerLoopArray, + ModelPart& rSurrogateModelPart_inner, ModelPart& rSurrogateModelPart_outer) + : BaseType(PointsArrayType(), &msGeometryData) + , mpNurbsSurface(pSurface) + , mOuterLoopArray(BrepOuterLoopArray) + , mInnerLoopArray(BrepInnerLoopArray) + , mpSurrogateModelPart_inner(&rSurrogateModelPart_inner) + , mpSurrogateModelPart_outer(&rSurrogateModelPart_outer) + { + mIsTrimmed = false; + } + explicit BrepSurface(const PointsArrayType& ThisPoints) : BaseType(ThisPoints, &msGeometryData) { @@ -467,7 +484,7 @@ class BrepSurface mpNurbsSurface->SpansLocalSpace(spans_u, 0); mpNurbsSurface->SpansLocalSpace(spans_v, 1); - BrepTrimmingUtilities::CreateBrepSurfaceTrimmingIntegrationPoints( + BrepTrimmingUtilitiesType::CreateBrepSurfaceTrimmingIntegrationPoints( rIntegrationPoints, mOuterLoopArray, mInnerLoopArray, spans_u, spans_v, @@ -594,6 +611,9 @@ class BrepSurface BrepCurveOnSurfaceArrayType mEmbeddedEdgesArray; + ModelPart* mpSurrogateModelPart_inner = nullptr; + ModelPart* mpSurrogateModelPart_outer = nullptr; + /** IsTrimmed is used to optimize processes as * e.g. creation of integration domain. */ diff --git a/kratos/utilities/geometry_utilities/brep_trimming_utilities.cpp b/kratos/utilities/geometry_utilities/brep_trimming_utilities.cpp index 0d682b17d244..7f98d2b5424a 100644 --- a/kratos/utilities/geometry_utilities/brep_trimming_utilities.cpp +++ b/kratos/utilities/geometry_utilities/brep_trimming_utilities.cpp @@ -6,6 +6,9 @@ // // License: BSD License // Kratos default license: kratos/license.txt +// +// Main authors: Nicolo' Antonelli +// Andrea Gorgi // header includes #include "brep_trimming_utilities.h" @@ -14,8 +17,8 @@ namespace Kratos { ///@name Kratos Classes ///@{ - - void BrepTrimmingUtilities::CreateBrepSurfaceTrimmingIntegrationPoints( + template + void BrepTrimmingUtilities::CreateBrepSurfaceTrimmingIntegrationPoints( IntegrationPointsArrayType& rIntegrationPoints, const DenseVector>& rOuterLoops, const DenseVector>& rInnerLoops, @@ -152,4 +155,7 @@ namespace Kratos // const std::vector& rSpansV, // IntegrationInfo& rIntegrationInfo); + template class KRATOS_API(KRATOS_CORE) BrepTrimmingUtilities; + template class KRATOS_API(KRATOS_CORE) BrepTrimmingUtilities; + } // namespace Kratos. diff --git a/kratos/utilities/geometry_utilities/brep_trimming_utilities.h b/kratos/utilities/geometry_utilities/brep_trimming_utilities.h index 2ee856f1b8b8..5162ac6842e6 100644 --- a/kratos/utilities/geometry_utilities/brep_trimming_utilities.h +++ b/kratos/utilities/geometry_utilities/brep_trimming_utilities.h @@ -6,6 +6,9 @@ // // License: BSD License // Kratos default license: kratos/license.txt +// +// Main authors: Nicolo' Antonelli +// Andrea Gorgi #if !defined(KRATOS_BREP_TRIMMING_UTILITIES_H_INCLUDED) #define KRATOS_BREP_TRIMMING_UTILITIES_H_INCLUDED @@ -32,7 +35,7 @@ namespace Kratos ///@{ //using namespace ClipperLib; - + template class KRATOS_API(KRATOS_CORE) BrepTrimmingUtilities { public: @@ -52,7 +55,7 @@ namespace Kratos typedef signed long long cInt; - using BrepCurveOnSurfacePointerType = typename BrepCurveOnSurface, false, PointerVector>::Pointer; + using BrepCurveOnSurfacePointerType = typename BrepCurveOnSurface, TShiftedBoundary, PointerVector>::Pointer; //template static void CreateBrepSurfaceTrimmingIntegrationPoints( diff --git a/kratos/utilities/nurbs_utilities/snake_sbm_utilities.cpp b/kratos/utilities/nurbs_utilities/snake_sbm_utilities.cpp new file mode 100644 index 000000000000..7c993bff971a --- /dev/null +++ b/kratos/utilities/nurbs_utilities/snake_sbm_utilities.cpp @@ -0,0 +1,839 @@ +// | / | +// ' / __| _` | __| _ \ __| +// . \ | ( | | ( |\__ ` +// _|\_\_| \__,_|\__|\___/ ____/ +// Multi-Physics +// +// License: BSD License +// Kratos default license: kratos/license.txt +// + +// System includes + +// External includes + +// Project includes +#include "snake_sbm_utilities.h" + +namespace Kratos +{ + + void SnakeSBMUtilities::CreateTheSnakeCoordinates(ModelPart& iga_model_part, + ModelPart& skin_model_part_inner_initial, + ModelPart& skin_model_part_outer_initial, + ModelPart& skin_model_part, + int rEchoLevel, + Vector& knot_vector_u, + Vector& knot_vector_v, + const Parameters mParameters) { + + // Check + KRATOS_ERROR_IF_NOT(mParameters.Has("sbm_parameters")) << "sbm_parameters has not been defined in the nurbs modeler" << std::endl; + + // Initilize the property of skin_model_part_in and out + if (skin_model_part_inner_initial.Nodes().size()>0) { + + CreateTheSnakeCoordinates(iga_model_part, skin_model_part_inner_initial, skin_model_part, rEchoLevel, knot_vector_u, knot_vector_v, mParameters, true); + + } + if (skin_model_part_outer_initial.Nodes().size()>0) { + + CreateTheSnakeCoordinates(iga_model_part, skin_model_part_outer_initial, skin_model_part, rEchoLevel, knot_vector_u, knot_vector_v, mParameters, false); + } + } + + void SnakeSBMUtilities::CreateTheSnakeCoordinates(ModelPart& iga_model_part, + ModelPart& skin_model_part_initial, + ModelPart& skin_model_part, + int rEchoLevel, + Vector& knot_vector_u, + Vector& knot_vector_v, + const Parameters mParameters, + bool is_inner) + { + + std::string surrogate_sub_model_part_name; + std::string skin_sub_model_part_name; + // ModelPart skin_sub_model_part; + if (is_inner) { + surrogate_sub_model_part_name = "surrogate_inner"; + skin_sub_model_part_name = "inner"; + } + else { + surrogate_sub_model_part_name = "surrogate_outer"; + skin_sub_model_part_name = "outer"; + } + + ModelPart& skin_sub_model_part = skin_model_part.GetSubModelPart(skin_sub_model_part_name); + ModelPart& surrogate_sub_model_part = iga_model_part.GetSubModelPart(surrogate_sub_model_part_name); + + if (!skin_model_part_initial.HasProperties(0)) skin_model_part_initial.CreateNewProperties(0); + if (!skin_model_part.HasProperties(0)) skin_model_part.CreateNewProperties(0); + Properties::Pointer p_cond_prop_in = skin_model_part_initial.pGetProperties(0); + skin_model_part_initial.AddProperties(p_cond_prop_in); + + Vector knot_step_uv(2); + knot_step_uv[0] = std::abs(knot_vector_u[int(knot_vector_u.size()/2) +1] - knot_vector_u[int(knot_vector_u.size()/2)] ) ; + knot_step_uv[1] = std::abs(knot_vector_v[int(knot_vector_v.size()/2) +1] - knot_vector_v[int(knot_vector_v.size()/2)] ) ; + + Vector meshSizes_uv(2); + meshSizes_uv[0] = knot_step_uv[0]; + meshSizes_uv[1] = knot_step_uv[1]; + ModelPart& surrogate_model_part = iga_model_part.GetSubModelPart(surrogate_sub_model_part_name); + surrogate_model_part.GetProcessInfo().SetValue(MARKER_MESHES, meshSizes_uv); + + Vector starting_pos_uv(2); + starting_pos_uv[0] = knot_vector_u[0]; + starting_pos_uv[1] = knot_vector_v[0]; + + Vector parameterExternalCoordinates(4); + parameterExternalCoordinates[0] = knot_vector_u[0]; + parameterExternalCoordinates[1] = knot_vector_v[0]; + parameterExternalCoordinates[2] = knot_vector_u[knot_vector_u.size()-1]; + parameterExternalCoordinates[3] = knot_vector_v[knot_vector_v.size()-1]; + + surrogate_model_part.GetProcessInfo().SetValue(LOAD_MESHES, parameterExternalCoordinates); + + // Create the matrix of active/inactive knot spans, one for inner and one for outer loop + unsigned int numberOfLoops; + if (is_inner) + numberOfLoops = mParameters["sbm_parameters"]["number_of_inner_loops"].GetInt(); + else + numberOfLoops = 1; + + std::vector n_knot_spans_uv(2); + n_knot_spans_uv[0] = knot_vector_u.size()-1; + n_knot_spans_uv[1] = knot_vector_v.size()-1; + + std::vector>> knot_spans_available; + knot_spans_available.reserve(numberOfLoops); + + for (int i = 0; i < numberOfLoops; ++i) { + std::vector> matrix; + matrix.reserve(n_knot_spans_uv[1]); + for (int j = 0; j <= n_knot_spans_uv[1]; ++j) { + std::vector row(n_knot_spans_uv[0]); + matrix.push_back(row); + } + knot_spans_available.push_back(matrix); + } + + // Optimized Snake -> for inner loops + int idMatrixKnotSpansAvailable = 0; + int idFirstNode; + bool newInnerLoop = true; + + + KRATOS_INFO_IF("::[SnakeSBMUtilities]::", rEchoLevel > 0 && is_inner) << "Inner :: Starting SnakeStep" << std::endl; + KRATOS_INFO_IF("::[SnakeSBMUtilities]::", rEchoLevel > 0 && !is_inner) << "Outer :: Starting SnakeStep" << std::endl; + + if (skin_model_part_initial.Conditions().size() > 0) { + + // CREATE FIRST NODE FOR SKIN SUB MODEL PART + auto initial_condition = skin_model_part_initial.GetCondition(1); + double x_true_boundary0 = initial_condition.GetGeometry()[0].X(); + double y_true_boundary0 = initial_condition.GetGeometry()[0].Y(); + + const int id_first_node = skin_model_part.GetRootModelPart().Nodes().size()+1; + skin_sub_model_part.CreateNewNode(id_first_node, x_true_boundary0, y_true_boundary0, 0.0); + + for (auto &i_cond : skin_model_part_initial.Conditions()) { + if (newInnerLoop) { + idFirstNode = i_cond.GetGeometry()[0].Id(); + newInnerLoop = false; + } + // Collect the coordinates of the points of the i_cond + double x_true_boundary1 = i_cond.GetGeometry()[0].X(); + double y_true_boundary1 = i_cond.GetGeometry()[0].Y(); + double x_true_boundary2 = i_cond.GetGeometry()[1].X(); + double y_true_boundary2 = i_cond.GetGeometry()[1].Y(); + + std::vector> xy_coord_i_cond(2); + xy_coord_i_cond[0].resize(2); xy_coord_i_cond[1].resize(2); + + xy_coord_i_cond[0][0] = i_cond.GetGeometry()[0].X(); // x_true_boundary1 + xy_coord_i_cond[1][0] = i_cond.GetGeometry()[0].Y(); // y_true_boundary1 + xy_coord_i_cond[0][1] = i_cond.GetGeometry()[1].X(); // x_true_boundary2 + xy_coord_i_cond[1][1] = i_cond.GetGeometry()[1].Y(); // y_true_boundary2 + + // Collect the intersections of the skin boundary with the knot values + std::vector> knot_span_uv(2); + knot_span_uv[0].resize(2); knot_span_uv[1].resize(2); + + knot_span_uv[0][0] = (x_true_boundary1-starting_pos_uv[0]) / knot_step_uv[0]; // knot_span_u_1st_point + knot_span_uv[1][0] = (y_true_boundary1-starting_pos_uv[1]) / knot_step_uv[1]; // knot_span_v_1st_point + knot_span_uv[0][1] = (x_true_boundary2-starting_pos_uv[0]) / knot_step_uv[0]; // knot_span_u_2nd_point + knot_span_uv[1][1] = (y_true_boundary2-starting_pos_uv[1]) / knot_step_uv[1]; // knot_span_v_2nd_point + + // In the inner case, check is the immersed object is inside the rectangular domain + if (is_inner && + (knot_span_uv[0][0] < 0 || knot_span_uv[0][0] >= n_knot_spans_uv[0] || + knot_span_uv[1][0] < 0 || knot_span_uv[1][0] >= n_knot_spans_uv[1] || + knot_span_uv[0][1] < 0 || knot_span_uv[0][1] >= n_knot_spans_uv[0] || + knot_span_uv[1][1] < 0 || knot_span_uv[1][1] >= n_knot_spans_uv[1]) ) + KRATOS_ERROR << "[SnakeSbmUtilities]:: The skin boundary provided is bigger than the background geometry in the parameter space." << std::endl; + + // additional check knot_span_uv computation on the domain border [especially for outer boundary] + if (knot_span_uv[0][0] == n_knot_spans_uv[0]) knot_span_uv[0][0]--; + if (knot_span_uv[1][0] == n_knot_spans_uv[1]) knot_span_uv[1][0]--; + if (knot_span_uv[0][1] == n_knot_spans_uv[0]) knot_span_uv[0][1]--; + if (knot_span_uv[1][1] == n_knot_spans_uv[1]) knot_span_uv[1][1]--; + + SnakeStep(skin_sub_model_part, knot_spans_available, idMatrixKnotSpansAvailable, + knot_span_uv, xy_coord_i_cond, knot_step_uv, starting_pos_uv); + + if (i_cond.GetGeometry()[1].Id() == idFirstNode) { + idMatrixKnotSpansAvailable++; + newInnerLoop = true; + } + } + } + + PointVector points; + for (auto &i_cond : skin_sub_model_part.Conditions()) { + points.push_back(PointTypePointer(new PointType(i_cond.Id(), i_cond.GetGeometry()[0].X(), i_cond.GetGeometry()[0].Y(), i_cond.GetGeometry()[0].Z()))); + } + DynamicBins testBins(points.begin(), points.end()); + + KRATOS_INFO_IF("::[SnakeSBMUtilities]::", rEchoLevel > 0 && is_inner) << "Inner :: Ending SnakeStep" << std::endl; + KRATOS_INFO_IF("::[SnakeSBMUtilities]::", rEchoLevel > 0 && !is_inner) << "Outer :: Ending SnakeStep" << std::endl; + + // Read lambda parameters: 0.0 -> External, 0.5 -> Optimal + double lambda; + if (is_inner) + if (mParameters["sbm_parameters"].Has("lambda_inner")) + lambda = mParameters["sbm_parameters"]["lambda_inner"].GetDouble(); + else + lambda = 0.5; + else + if (mParameters["sbm_parameters"].Has("lambda_outer")) + lambda = mParameters["sbm_parameters"]["lambda_outer"].GetDouble(); + else + lambda = 0.5; + + KRATOS_INFO_IF("::[SnakeSBMUtilities]::", rEchoLevel > 0 && is_inner) << "Inner :: MarkKnotSpansAvailable" << std::endl; + KRATOS_INFO_IF("::[SnakeSBMUtilities]::", rEchoLevel > 0 && !is_inner) << "Outer :: MarkKnotSpansAvailable" << std::endl; + + for (int i = 0; i < numberOfLoops; i++) { + int idInnerLoop = i; + // Mark the knot_spans_available's for inner and outer loops + MarkKnotSpansAvailable(knot_spans_available, idInnerLoop, testBins, skin_sub_model_part, lambda, + n_knot_spans_uv, knot_step_uv, starting_pos_uv); + + if (is_inner) { + CreateSurrogateBuondaryFromSnake_inner (knot_spans_available, idInnerLoop, surrogate_sub_model_part, + n_knot_spans_uv, knot_vector_u, knot_vector_v, starting_pos_uv ); + KRATOS_INFO_IF("::[SnakeSBMUtilities]::", rEchoLevel > 0) << "Inner :: Snake process has finished" << std::endl; + } + else { + CreateSurrogateBuondaryFromSnake_outer (testBins, skin_sub_model_part, knot_spans_available, idInnerLoop, surrogate_sub_model_part, + n_knot_spans_uv, knot_vector_u, knot_vector_v, starting_pos_uv); + KRATOS_INFO_IF("::[SnakeSBMUtilities]::", rEchoLevel > 0) << "Outer :: Snake process has finished" << std::endl; + } + } + + KRATOS_INFO_IF("::[SnakeSBMUtilities]::", rEchoLevel > 0 && is_inner) << "Inner :: Loop finished" << std::endl; + KRATOS_INFO_IF("::[SnakeSBMUtilities]::", rEchoLevel > 0 && !is_inner) << "Outer :: Loop finished" << std::endl; + } + + + void SnakeSBMUtilities::SnakeStep(ModelPart& skin_model_part, + std::vector>> &knot_spans_available, + int idMatrix, + std::vector> knot_spans_uv, + std::vector> xy_coord_i_cond, + Vector knot_step_uv, + Vector starting_pos_uv) + { + bool isSplitted = false; + + if (knot_spans_uv[0][0] != knot_spans_uv[0][1] || knot_spans_uv[1][0] != knot_spans_uv[1][1]) { // INTERSECTION BETWEEN TRUE AND SURROGATE BOUNDARY + // Check if we are jumping some cut knot spans. If yes we split the true segment + if (std::abs(knot_spans_uv[1][0]-knot_spans_uv[1][1]) > 1 || std::abs(knot_spans_uv[0][0]-knot_spans_uv[0][1]) > 1 || + (knot_spans_uv[0][0] != knot_spans_uv[0][1] && knot_spans_uv[1][0] != knot_spans_uv[1][1]) ) { + isSplitted = true; + + // Split the segment and do it recursively + double x_true_boundary_split = (xy_coord_i_cond[0][0]+xy_coord_i_cond[0][1]) / 2; + double y_true_boundary_split = (xy_coord_i_cond[1][0]+xy_coord_i_cond[1][1]) / 2; + int knot_span_u_point_split = (x_true_boundary_split-starting_pos_uv[0]) / knot_step_uv[0] ; + int knot_span_v_point_split = (y_true_boundary_split-starting_pos_uv[1]) / knot_step_uv[1] ; + + // update xy_coord for the first split segment + std::vector> xy_coord_i_cond_split(2); + xy_coord_i_cond_split[0].resize(2); xy_coord_i_cond_split[1].resize(2); + xy_coord_i_cond_split[0][0] = xy_coord_i_cond[0][0]; // x_true_boundary1 + xy_coord_i_cond_split[1][0] = xy_coord_i_cond[1][0]; // y_true_boundary1 + xy_coord_i_cond_split[0][1] = x_true_boundary_split; // x_true_boundary_split + xy_coord_i_cond_split[1][1] = y_true_boundary_split; // y_true_boundary_split + // update knot_span_uv for the first split segment + std::vector> knot_span_uv_split(2); + knot_span_uv_split[0].resize(2); knot_span_uv_split[1].resize(2); + knot_span_uv_split[0][0] = knot_spans_uv[0][0]; // knot_span_u_1st_point + knot_span_uv_split[1][0] = knot_spans_uv[1][0]; // knot_span_v_1st_point + knot_span_uv_split[0][1] = knot_span_u_point_split; // knot_span_u_point_split + knot_span_uv_split[1][1] = knot_span_v_point_split; // knot_span_v_point_split + + // __We do it recursively first split__ + SnakeStep(skin_model_part, knot_spans_available, idMatrix, knot_span_uv_split, + xy_coord_i_cond_split, knot_step_uv, starting_pos_uv ); + + // update xy_coord for the second split segment + xy_coord_i_cond_split[0][0] = x_true_boundary_split; // x_true_boundary_split + xy_coord_i_cond_split[1][0] = y_true_boundary_split; // y_true_boundary_split + xy_coord_i_cond_split[0][1] = xy_coord_i_cond[0][1]; // x_true_boundary2 + xy_coord_i_cond_split[1][1] = xy_coord_i_cond[1][1]; // y_true_boundary2 + // update knot_span_uv for the first split segment + knot_span_uv_split[0][0] = knot_span_u_point_split; // knot_span_u_point_split + knot_span_uv_split[1][0] = knot_span_v_point_split; // knot_span_v_point_split + knot_span_uv_split[0][1] = knot_spans_uv[0][1]; // knot_span_u_2nd_point + knot_span_uv_split[1][1] = knot_spans_uv[1][1]; // knot_span_v_2nd_point + + // __We do it recursively second split__ + SnakeStep(skin_model_part, knot_spans_available, idMatrix, knot_span_uv_split, + xy_coord_i_cond_split, knot_step_uv, starting_pos_uv); + } + // Check if the true boundary crosses an u or a v knot value + else if (knot_spans_uv[0][0] != knot_spans_uv[0][1]) { // u knot value is crossed + // Find the "knot_spans_available" using the intersection + knot_spans_available[idMatrix][knot_spans_uv[1][0]][knot_spans_uv[0][0]] = 2; + knot_spans_available[idMatrix][knot_spans_uv[1][0]][knot_spans_uv[0][1]] = 2; + + } + else if (knot_spans_uv[1][0] != knot_spans_uv[1][1]) { // v knot value is crossed + // Find the "knot_spans_available" using the intersection (Snake_coordinate classic -> External Boundary) + knot_spans_available[idMatrix][knot_spans_uv[1][0]][knot_spans_uv[0][0]] = 2; + knot_spans_available[idMatrix][knot_spans_uv[1][1]][knot_spans_uv[0][0]] = 2; + } + } + if (!isSplitted) { + // Call the root model part for the Ids of the node + auto idNode1 = skin_model_part.GetRootModelPart().Nodes().size(); + auto idNode2 = idNode1+1; + // Create two nodes and two conditions for each skin condition + skin_model_part.CreateNewNode(idNode2, (xy_coord_i_cond[0][0]+xy_coord_i_cond[0][1] ) / 2, (xy_coord_i_cond[1][0]+xy_coord_i_cond[1][1] ) / 2, 0.0); + skin_model_part.CreateNewNode(idNode2+1, xy_coord_i_cond[0][1], xy_coord_i_cond[1][1], 0.0); + Properties::Pointer p_cond_prop = skin_model_part.pGetProperties(0); + Condition::Pointer p_cond1 = skin_model_part.CreateNewCondition("LineCondition2D2N", idNode1, {{idNode1, idNode2}}, p_cond_prop ); + Condition::Pointer p_cond2 = skin_model_part.CreateNewCondition("LineCondition2D2N", idNode2, {{idNode2, idNode2+1}}, p_cond_prop ); + skin_model_part.AddCondition(p_cond1); + skin_model_part.AddCondition(p_cond2); + } + } + + + bool SnakeSBMUtilities::isPointInsideSkinBoundary(Point& point1, DynamicBins& testBins, ModelPart& skin_model_part) + { + // Get the nearest point of the true boundary + PointerType pointToSearch = PointerType(new PointType(1000000, point1.X(), point1.Y(), 0.0)); + PointerType nearestPoint = testBins.SearchNearestPoint(*pointToSearch); + + // Get the closest Condition the initial_skin_model_part_in.Conditions + int id1 = nearestPoint->Id(); + auto nearestCondition1 = skin_model_part.GetCondition(id1); + // Check if the condition is the first one and therefore the previous one does not exist + int id2 = id1 - 1; + if (id1 == skin_model_part.ConditionsBegin()->Id()) { + int nConditions = skin_model_part.Conditions().size(); + id2 = id1 + nConditions - 1; + } + auto nearestCondition2 = skin_model_part.GetCondition(id2); + // The two candidates nodes + Point candidatePoint1 = nearestCondition1.GetGeometry()[1]; + Point candidatePoint2 = nearestCondition2.GetGeometry()[0]; + + Point point2 = nearestCondition1.GetGeometry()[0]; // FIRST POINT IN TRUE GEOM + Point point3 = candidatePoint1;// SECOND POINT IN TRUE GEOM + if (MathUtils::Norm(candidatePoint1-point1) > MathUtils::Norm(candidatePoint2-point1)){ + // Need to invert the order to preserve the positivity of the area + point3 = point2; + point2 = candidatePoint2; + } + + array_1d v_1 = point2 - point1; + array_1d v_2 = point3 - point1; + array_1d crossProduct; + MathUtils::CrossProduct(crossProduct, v_1, v_2); + + bool isInside = false; + if (crossProduct[2] > 0) {isInside = true;} + + return isInside; + } + + + + void SnakeSBMUtilities::MarkKnotSpansAvailable(std::vector>> & knot_spans_available, int idMatrix,DynamicBins& testBin, + ModelPart& skin_model_part, double lambda, std::vector& n_knot_spans_uv, + Vector& knot_step_uv, const Vector& starting_pos_uv) { + for (int i = 0; i < n_knot_spans_uv[1]; i++) { + for (int j = 0; j < n_knot_spans_uv[0]; j++) { + if (knot_spans_available[idMatrix][i][j] == 2) { + // Check the 8 neighbor knot spans -> Is there any completely inside? Note that we can just check 1 point. + + // right node + if (i != n_knot_spans_uv[1]-1) + if (knot_spans_available[idMatrix][i+1][j] == 0) { + Point gaussPoint = Point((j+0.5) * knot_step_uv[0] + starting_pos_uv[0], (i+1+0.5) * knot_step_uv[1] +starting_pos_uv[1], 0); + if (isPointInsideSkinBoundary(gaussPoint, testBin, skin_model_part)) {knot_spans_available[idMatrix][i+1][j] = 1;} + } + // left node + if (i != 0) + if (knot_spans_available[idMatrix][i-1][j] == 0) { + Point gaussPoint = Point((j+0.5) * knot_step_uv[0]+starting_pos_uv[0], (i-1+0.5) * knot_step_uv[1] + starting_pos_uv[1], 0); + if (isPointInsideSkinBoundary(gaussPoint, testBin, skin_model_part)) {knot_spans_available[idMatrix][i-1][j] = 1;} + } + // up node + if (j != n_knot_spans_uv[0]-1) + if (knot_spans_available[idMatrix][i][j+1] == 0) { + Point gaussPoint = Point((j+1+0.5) * knot_step_uv[0]+starting_pos_uv[0], (i+0.5) * knot_step_uv[1]+starting_pos_uv[1], 0); + if (isPointInsideSkinBoundary(gaussPoint, testBin, skin_model_part)) {knot_spans_available[idMatrix][i][j+1] = 1;} + } + //down node + if (j != 0) + if (knot_spans_available[idMatrix][i][j-1] == 0) { + Point gaussPoint = Point((j-1+0.5) * knot_step_uv[0]+starting_pos_uv[0], (i+0.5) * knot_step_uv[1]+starting_pos_uv[1], 0); + if (isPointInsideSkinBoundary(gaussPoint, testBin, skin_model_part)) {knot_spans_available[idMatrix][i][j-1] = 1;} + } + + // corner right-down node + if (j != 0 && i != n_knot_spans_uv[1]-1) + if (knot_spans_available[idMatrix][i+1][j-1] == 0) { + Point gaussPoint = Point((j-1+0.5) * knot_step_uv[0]+starting_pos_uv[0], (i+1+0.5) * knot_step_uv[1]+starting_pos_uv[1], 0); + if (isPointInsideSkinBoundary(gaussPoint, testBin, skin_model_part)) {knot_spans_available[idMatrix][i+1][j-1] = 1;} + } + // corner left-down node + if (j != 0 && i != 0) + if (knot_spans_available[idMatrix][i-1][j-1] == 0) { + Point gaussPoint = Point((j-1+0.5) * knot_step_uv[0]+starting_pos_uv[0], (i-1+0.5) * knot_step_uv[1]+starting_pos_uv[1], 0); + if (isPointInsideSkinBoundary(gaussPoint, testBin, skin_model_part)) {knot_spans_available[idMatrix][i-1][j-1] = 1;} + } + // corner right-up node + if (j != n_knot_spans_uv[0]-1 && i != n_knot_spans_uv[1]-1) + if (knot_spans_available[idMatrix][i+1][j+1] == 0) { + Point gaussPoint = Point((j+1+0.5) * knot_step_uv[0]+starting_pos_uv[0], (i+1+0.5) * knot_step_uv[1]+starting_pos_uv[1], 0); + if (isPointInsideSkinBoundary(gaussPoint, testBin, skin_model_part)) {knot_spans_available[idMatrix][i+1][j+1] = 1;} + } + // corner left-up node + if (j != n_knot_spans_uv[0]-1 && i != 0) + if (knot_spans_available[idMatrix][i-1][j+1] == 0) { + Point gaussPoint = Point((j+1+0.5) * knot_step_uv[0]+starting_pos_uv[0], (i-1+0.5) * knot_step_uv[1]+starting_pos_uv[1], 0); + if (isPointInsideSkinBoundary(gaussPoint, testBin, skin_model_part)) {knot_spans_available[idMatrix][i-1][j+1] = 1;} + } + + // Create 25 "fake" GaussPoints to check if the majority are inside or outside + const int numFakeGaussPoints = 5; + int numberOfInsideGaussianPoints = 0; + for (int i_GPx = 0; i_GPx < numFakeGaussPoints; i_GPx++){ + double x_coord = j*knot_step_uv[0] + knot_step_uv[0]/(numFakeGaussPoints+1)*(i_GPx+1) + starting_pos_uv[0]; + + // NOTE:: The v-knot spans are upside down in the matrix!! + for (int i_GPy = 0; i_GPy < numFakeGaussPoints; i_GPy++) + { + double y_coord = i*knot_step_uv[1] + knot_step_uv[1]/(numFakeGaussPoints+1)*(i_GPy+1) + starting_pos_uv[1]; + Point gaussPoint = Point(x_coord, y_coord, 0); // GAUSSIAN POINT + if (isPointInsideSkinBoundary(gaussPoint, testBin, skin_model_part)) { + // Sum over the number of numFakeGaussPoints per knot span + numberOfInsideGaussianPoints++; + } + } + + } + + // Mark the knot span as available or not depending on the number of Gauss Points Inside/Outside + if (numberOfInsideGaussianPoints < lambda*numFakeGaussPoints*numFakeGaussPoints) { + knot_spans_available[idMatrix][i][j] = -1; // Cut knot spans that have been checked + } + else{ + knot_spans_available[idMatrix][i][j] = 1; // The knot span is considered DEACTIVE + } + + // exit(0); + + } + } + } + } + + + void SnakeSBMUtilities::CreateSurrogateBuondaryFromSnake_inner (std::vector>> & knot_spans_available, int idMatrix, + ModelPart& surrogate_model_part_inner, std::vector& n_knot_spans_uv, + Vector& knot_vector_u, Vector& knot_vector_v, + const Vector& starting_pos_uv){ + // Find the start of the Snake + int start_i = 0; + int start_j = 0; + for (int i = 0; i < (n_knot_spans_uv[0]); i++) { + for (int j = 0; j < (n_knot_spans_uv[1]); j++ ) { + if (knot_spans_available[idMatrix][j][i] == 1 ) { + // Check if one of the neighouring knot span is also == 1, + // otherwise the knot span is isolated I cannot be taken. + if ((knot_spans_available[idMatrix][j+1][i] == 1) || (knot_spans_available[idMatrix][j][i+1] == 1)) + { + start_i = i; + start_j = j; + break; + } + } + } + if (knot_spans_available[idMatrix][start_j][start_i] == 1 ) { break; } + } + + if (!surrogate_model_part_inner.HasProperties(0)) {surrogate_model_part_inner.CreateNewProperties(0);} + + Properties::Pointer p_cond_prop = surrogate_model_part_inner.pGetProperties(0); + Node snakeNode(1 , knot_vector_u[start_i], knot_vector_v[start_j], 0.0); + + const int id_first_node = surrogate_model_part_inner.GetRootModelPart().Nodes().size() + 1; + surrogate_model_part_inner.CreateNewNode(id_first_node, snakeNode); + IndexType idSnakeNode = id_first_node+1; + + // Follow the clockwise loop + int end = 0; + // We are going orizontally + int direction = 0; + // 0 = up_vertical, 1 = right_orizontal, 2 = down_vertical, 3 = left_orizontal + int i = start_i; int j = start_j; + int I = start_i; int J = start_j; + int steps = 0; + + const int max_number_of_steps = 1e5; + while (end == 0 && steps < max_number_of_steps) { + steps ++ ; + int is_special_case = 1; // Variable to check if we are in a super special case and we shouldn't go out of the while loop + // Try to go left w.r.t. the current direction + if (direction == 0) {I-- ; } + if (direction == 1) {J++; } + if (direction == 2) {I++ ; } + if (direction == 3) {J-- ; } + if (knot_spans_available[idMatrix][J][I] == 1) { + direction = direction - 1; + if (direction == -1) {direction = 3;} + } + else { + // Need to try to go straight or right w.r.t. the current direction; risetto e muovo + if (direction == 0) {I++ ; J++ ;} + if (direction == 1) {J-- ; I++ ;} + if (direction == 2) {I-- ; J-- ;} + if (direction == 3) {J++ ; I-- ;} + if (knot_spans_available[idMatrix][J][I] == 1) { + // going straight + if (direction == 0) {j++ ; } + if (direction == 1) {i++ ; } + if (direction == 2) {j-- ; } + if (direction == 3) {i-- ; } + + Node snakeNode(idSnakeNode , knot_vector_u[i], knot_vector_v[j], 0.0); + surrogate_model_part_inner.CreateNewNode(idSnakeNode, snakeNode); + Condition::Pointer pcond1 = surrogate_model_part_inner.CreateNewCondition("LineCondition2D2N", idSnakeNode-1, {{idSnakeNode-1, idSnakeNode}}, p_cond_prop ); + idSnakeNode++; + } + else { + // Need to search to hte right; rese e muovo + if (direction == 0) {J-- ; I++ ;} + if (direction == 1) {I-- ; J-- ;} + if (direction == 2) {J++ ; I-- ;} + if (direction == 3) {I++ ; J++ ;} + if (knot_spans_available[idMatrix][J][I] == 1) { + // We are moving to the right -> First move straight, store, the move to the right (i,j), store again + if (direction == 0) {j++ ; } + if (direction == 1) {i++ ; } + if (direction == 2) {j-- ; } + if (direction == 3) {i-- ; } + Node snakeNode1(idSnakeNode , knot_vector_u[i], knot_vector_v[j], 0.0); + surrogate_model_part_inner.CreateNewNode(idSnakeNode, snakeNode1); + Condition::Pointer pcond1 = surrogate_model_part_inner.CreateNewCondition("LineCondition2D2N", idSnakeNode-1, {{idSnakeNode-1, idSnakeNode}}, p_cond_prop ); + idSnakeNode++; + if (direction == 0) {i++ ; } + if (direction == 1) {j-- ; } + if (direction == 2) {i-- ; } + if (direction == 3) {j++ ; } + Node snakeNode2(idSnakeNode , knot_vector_u[i], knot_vector_v[j], 0.0); + surrogate_model_part_inner.CreateNewNode(idSnakeNode, snakeNode2); + Condition::Pointer pcond2 = surrogate_model_part_inner.CreateNewCondition("LineCondition2D2N", idSnakeNode-1, {{idSnakeNode-1, idSnakeNode}}, p_cond_prop ); + idSnakeNode++; + direction = direction + 1; + if (direction == 4) {direction = 0;} + } + else { // Special case of "isolated" knot span to be circumnavigated + is_special_case = 0; + // reset and move (I,J) backward + if (direction == 0) {I-- ; J-- ;} + if (direction == 1) {J++ ; I-- ;} + if (direction == 2) {I++ ; J++ ;} + if (direction == 3) {J-- ; I++ ;} + // Check + if (knot_spans_available[idMatrix][J][I] == 1) { + // First straight, store, then move to the right, store again, then move to the right and store again + if (direction == 0) {j++ ; } + if (direction == 1) {i++ ; } + if (direction == 2) {j-- ; } + if (direction == 3) {i-- ; } + + Node snakeNode(idSnakeNode , knot_vector_u[i], knot_vector_v[j], 0.0); + surrogate_model_part_inner.CreateNewNode(idSnakeNode, snakeNode); + Condition::Pointer pcond1 = surrogate_model_part_inner.CreateNewCondition("LineCondition2D2N", idSnakeNode-1, {{idSnakeNode-1, idSnakeNode}}, p_cond_prop ); + idSnakeNode++; + + if (direction == 0) {i++ ; } + if (direction == 1) {j-- ; } + if (direction == 2) {i-- ; } + if (direction == 3) {j++ ; } + + Node snakeNode2(idSnakeNode , knot_vector_u[i], knot_vector_v[j], 0.0); + surrogate_model_part_inner.CreateNewNode(idSnakeNode, snakeNode2); + Condition::Pointer pcond2 = surrogate_model_part_inner.CreateNewCondition("LineCondition2D2N", idSnakeNode-1, {{idSnakeNode-1, idSnakeNode}}, p_cond_prop ); + idSnakeNode++; + + if (direction == 0) {j-- ; } + if (direction == 1) {i-- ; } + if (direction == 2) {j++ ; } + if (direction == 3) {i++ ; } + + Node snakeNode3(idSnakeNode , knot_vector_u[i], knot_vector_v[j], 0.0); + surrogate_model_part_inner.CreateNewNode(idSnakeNode, snakeNode3); + Condition::Pointer pcond3 = surrogate_model_part_inner.CreateNewCondition("LineCondition2D2N", idSnakeNode-1, {{idSnakeNode-1, idSnakeNode}}, p_cond_prop ); + idSnakeNode++; + + // Finally rotate the direction of 180 degrees + direction = direction + 2; + if (direction == 4) {direction = 0;} + if (direction == 5) {direction = 1;} + } + else + KRATOS_ERROR << "SnakeSBMUtilities:: " << "error in the Snakes Coordinates" << std::endl; + } + } + } + // Check if we have close the snake + if (I == start_i && J == start_j && is_special_case == 1 ) { + // End of the while loop + end = 1; + KRATOS_INFO("number of steps in the snake step") << steps << std::endl; + } + } + KRATOS_ERROR_IF(steps >= max_number_of_steps-1) << "SnakeSBMUtilities:: " << "reached maximum number of steps" << std::endl; + + // Create "fictituos element" to memorize starting and ending node id for each surrogate boundary loop + IndexType initialId = id_first_node; + if (surrogate_model_part_inner.Elements().size()>0) { + // Check if it is not the first inner loop + initialId = surrogate_model_part_inner.GetElement(surrogate_model_part_inner.GetRootModelPart().Elements().size()).GetGeometry()[1].Id()+1; + } + std::vector elem_nodes{initialId, idSnakeNode-1}; + surrogate_model_part_inner.CreateNewElement("Element2D2N", surrogate_model_part_inner.GetRootModelPart().Elements().size()+1, elem_nodes, p_cond_prop); + } + + + void SnakeSBMUtilities::CreateSurrogateBuondaryFromSnake_outer (DynamicBins& testBin_out, ModelPart& initial_skin_model_part_out, + std::vector>> & knot_spans_available, int idMatrix, + ModelPart& surrogate_model_part_outer, std::vector& n_knot_spans_uv, + Vector& knot_vector_u, Vector& knot_vector_v, + const Vector& starting_pos_uv){ + + // CHECK ALL THE EXTERNAL KNOT SPANS + + // LEFT BOUNDARY + double knot_step_u = knot_vector_u[1]-knot_vector_u[0]; + double knot_step_v = knot_vector_v[1]-knot_vector_v[0]; + + for (int i = 0; i<2; i++) { + for (int j = 0; j < (n_knot_spans_uv[1]); j++ ) { + Point centroidKnotSpan = Point((j+0.5)*knot_step_u+starting_pos_uv[0], (i+0.5)*knot_step_v+starting_pos_uv[1], 0); + if (isPointInsideSkinBoundary(centroidKnotSpan, testBin_out, initial_skin_model_part_out) && knot_spans_available[idMatrix][i][j] != -1) { + knot_spans_available[idMatrix][i][j] = 1; + } + } + } + // TOP BOUNDARY + for (int j = knot_spans_available[idMatrix][0].size()-1; j > knot_spans_available[idMatrix][0].size()-3; j--) { + for (int i = 0; i < (n_knot_spans_uv[0]); i++) { + Point centroidKnotSpan = Point((j+0.5)*knot_step_u+starting_pos_uv[0], (i+0.5)*knot_step_v+starting_pos_uv[1], 0); + if (isPointInsideSkinBoundary(centroidKnotSpan, testBin_out, initial_skin_model_part_out) && knot_spans_available[idMatrix][i][j] != -1) { + knot_spans_available[idMatrix][i][j] = 1; + } + } + } + // RIGHT BOUNDARY + for (int i = knot_spans_available[idMatrix][0].size()-1; i > knot_spans_available[idMatrix][0].size()-3; i--) { + for (int j = n_knot_spans_uv[1]-1; j > -1; j-- ) { + Point centroidKnotSpan = Point((j+0.5)*knot_step_u+starting_pos_uv[0], (i+0.5)*knot_step_v+starting_pos_uv[1], 0); + if (isPointInsideSkinBoundary(centroidKnotSpan, testBin_out, initial_skin_model_part_out) && knot_spans_available[idMatrix][i][j] != -1) { + knot_spans_available[idMatrix][i][j] = 1; + } + } + } + // BOTTOM BOUNDARY + for (int j = 0; j<2; j++) { + for (int i = n_knot_spans_uv[0]-1; i > -1 ; i--) { + Point centroidKnotSpan = Point((j+0.5)*knot_step_u+starting_pos_uv[0], (i+0.5)*knot_step_v+starting_pos_uv[1], 0); + if (isPointInsideSkinBoundary(centroidKnotSpan, testBin_out, initial_skin_model_part_out) && knot_spans_available[idMatrix][i][j] != -1) { + knot_spans_available[idMatrix][i][j] = 1; + } + } + } + + // Find the start of the Snake + int start_i = 0; + int start_j = 0; + for (int i = 0; i < (n_knot_spans_uv[0]); i++) { + for (int j = 0; j < (n_knot_spans_uv[1]); j++ ) { + if (knot_spans_available[idMatrix][j][i] == 1 ) { + // Check if one of the neighouring knot span is also == 1, + // otherwise the knot span is isolated I cannot be taken. + if ((knot_spans_available[idMatrix][j+1][i] == 1) || (knot_spans_available[idMatrix][j][i+1] == 1)) + { + start_i = i; + start_j = j; + break; + } + } + } + if (knot_spans_available[idMatrix][start_j][start_i] == 1 ) { break; } + } + + // EXTEND THE MATRIX + //------------------- + std::vector> knot_spans_available_extended(n_knot_spans_uv[1]+2, std::vector(n_knot_spans_uv[0]+2)); + + for (int i = 0; i < knot_spans_available[idMatrix].size(); i++){ + for (int j = 0; j < knot_spans_available[idMatrix][0].size(); j++) { + knot_spans_available_extended[i+1][j+1] = knot_spans_available[idMatrix][i][j]; + } + } + + Properties::Pointer p_cond_prop = surrogate_model_part_outer.CreateNewProperties(1001); + Node snakeNode(1 , knot_vector_u[start_i], knot_vector_v[start_j], 0.0); + + const int id_first_node = surrogate_model_part_outer.GetRootModelPart().Nodes().size() + 1; + + surrogate_model_part_outer.CreateNewNode(id_first_node, snakeNode); + IndexType idSnakeNode = id_first_node+1; + + // Follow the clockwise loop + int end = 0; + // We are going orizontally + int direction = 0 ; // 0 = up_vertical, 1 = right_orizontal, 2 = down_vertical, 3 = left_orizontal + int I = start_i+1; + int J = start_j+1; + int i = start_i; + int j = start_j; + int steps = 0; + const int max_number_of_steps = 1e5; + while (end == 0 && steps < max_number_of_steps) { + steps ++ ; + int is_special_case = 1; // Variable to check if we are in a super special case and we shouldn't go out of the while loop + // Try to go left w.r.t. the current direction + if (direction == 0) {I-- ; } + if (direction == 1) {J++; } + if (direction == 2) {I++ ; } + if (direction == 3) {J-- ; } + if (knot_spans_available_extended[J][I] == 1) { + direction = direction - 1; + if (direction == -1) {direction = 3;} + } + else { + // Need to try to go straight or right w.r.t. the current direction; reset and move + if (direction == 0) {I++ ; J++ ;} + if (direction == 1) {J-- ; I++ ;} + if (direction == 2) {I-- ; J-- ;} + if (direction == 3) {J++ ; I-- ;} + + if (knot_spans_available_extended[J][I] == 1) { + + // Going straight! -> Don't store and move (i,j) + if (direction == 0) {j++ ; } + if (direction == 1) {i++ ; } + if (direction == 2) {j-- ; } + if (direction == 3) {i-- ; } + + Node snakeNode(idSnakeNode , knot_vector_u[i], knot_vector_v[j], 0.0); + surrogate_model_part_outer.CreateNewNode(idSnakeNode, snakeNode); + Condition::Pointer pcond1 = surrogate_model_part_outer.CreateNewCondition("LineCondition2D2N", idSnakeNode-1, {{idSnakeNode-1, idSnakeNode}}, p_cond_prop ); + idSnakeNode++; + } + else { + // Need to search to the right; Reset and move + if (direction == 0) {J-- ; I++ ;} + if (direction == 1) {I-- ; J-- ;} + if (direction == 2) {J++ ; I-- ;} + if (direction == 3) {I++ ; J++ ;} + + if (knot_spans_available_extended[J][I] == 1) { + + // Going to the right! -> first go forward, store, more right (i,j), store again + if (direction == 0) {j++ ; } + if (direction == 1) {i++ ; } + if (direction == 2) {j-- ; } + if (direction == 3) {i-- ; } + Node snakeNode1(idSnakeNode , knot_vector_u[i], knot_vector_v[j], 0.0); + surrogate_model_part_outer.CreateNewNode(idSnakeNode, snakeNode1); + Condition::Pointer pcond1 = surrogate_model_part_outer.CreateNewCondition("LineCondition2D2N", idSnakeNode-1, {{idSnakeNode-1, idSnakeNode}}, p_cond_prop ); + idSnakeNode++; + if (direction == 0) {i++ ; } + if (direction == 1) {j-- ; } + if (direction == 2) {i-- ; } + if (direction == 3) {j++ ; } + Node snakeNode2(idSnakeNode , knot_vector_u[i], knot_vector_v[j], 0.0); + surrogate_model_part_outer.CreateNewNode(idSnakeNode, snakeNode2); + Condition::Pointer pcond2 = surrogate_model_part_outer.CreateNewCondition("LineCondition2D2N", idSnakeNode-1, {{idSnakeNode-1, idSnakeNode}}, p_cond_prop ); + idSnakeNode++; + direction = direction + 1; + if (direction == 4) {direction = 0;} + } + else { // Special case of "isolated" knot span to be circumnavigated + is_special_case = 0; + // Reset and move (I,J) "backward" + if (direction == 0) {I-- ; J-- ;} + if (direction == 1) {J++ ; I-- ;} + if (direction == 2) {I++ ; J++ ;} + if (direction == 3) {J-- ; I++ ;} + // Check + if (knot_spans_available_extended[J][I] == 1) { + // First go forward, store, then move to right, then store again, then move to the right and store again + if (direction == 0) {j++ ; } + if (direction == 1) {i++ ; } + if (direction == 2) {j-- ; } + if (direction == 3) {i-- ; } + Node snakeNode(idSnakeNode , knot_vector_u[i], knot_vector_v[j], 0.0); + surrogate_model_part_outer.CreateNewNode(idSnakeNode, snakeNode); + Condition::Pointer pcond1 = surrogate_model_part_outer.CreateNewCondition("LineCondition2D2N", idSnakeNode-1, {{idSnakeNode-1, idSnakeNode}}, p_cond_prop ); + idSnakeNode++; + if (direction == 0) {i++ ; } + if (direction == 1) {j-- ; } + if (direction == 2) {i-- ; } + if (direction == 3) {j++ ; } + Node snakeNode2(idSnakeNode , knot_vector_u[i], knot_vector_v[j], 0.0); + surrogate_model_part_outer.CreateNewNode(idSnakeNode, snakeNode2); + Condition::Pointer pcond2 = surrogate_model_part_outer.CreateNewCondition("LineCondition2D2N", idSnakeNode-1, {{idSnakeNode-1, idSnakeNode}}, p_cond_prop ); + idSnakeNode++; + if (direction == 0) {j-- ; } + if (direction == 1) {i-- ; } + if (direction == 2) {j++ ; } + if (direction == 3) {i++ ; } + Node snakeNode3(idSnakeNode , knot_vector_u[i], knot_vector_v[j], 0.0); + surrogate_model_part_outer.CreateNewNode(idSnakeNode, snakeNode3); + Condition::Pointer pcond3 = surrogate_model_part_outer.CreateNewCondition("LineCondition2D2N", idSnakeNode-1, {{idSnakeNode-1, idSnakeNode}}, p_cond_prop ); + idSnakeNode++; + // Finally rotate the direction of 180 degrees + direction = direction + 2; + if (direction == 4) {direction = 0;} + if (direction == 5) {direction = 1;} + } + else + KRATOS_ERROR << "SnakeSBMUtilities:: " << "error in the Snakes Coordinates" << std::endl; + } + } + } + // Check if we have close the snake + if (I == start_i+1 && J == start_j+1 && is_special_case == 1 ) { + // End of the while loop + end = 1; + } + } + // Create "fictituos element" to memorize starting and ending node id for each surrogate boundary loop + std::vector elem_nodes{1, idSnakeNode-1}; + int elem_id = surrogate_model_part_outer.GetRootModelPart().Elements().size()+1; + surrogate_model_part_outer.CreateNewElement("Element2D2N", elem_id, elem_nodes, p_cond_prop); + } +} // namespace Kratos. diff --git a/kratos/utilities/nurbs_utilities/snake_sbm_utilities.h b/kratos/utilities/nurbs_utilities/snake_sbm_utilities.h new file mode 100644 index 000000000000..ae456ed545a5 --- /dev/null +++ b/kratos/utilities/nurbs_utilities/snake_sbm_utilities.h @@ -0,0 +1,179 @@ +// | / | +// ' / __| _` | __| _ \ __| +// . \ | ( | | ( |\__ ` +// _|\_\_| \__,_|\__|\___/ ____/ +// Multi-Physics +// +// License: BSD License +// Kratos default license: kratos/license.txt +// +// Main authors: Nicolo' Antonelli +// Andrea Gorgi +// + +#if !defined(KRATOS_SNAKE_SBM_UTILITIES_H_INCLUDED ) +#define KRATOS_SNAKE_SBM_UTILITIES_H_INCLUDED + +// Project includes +#include "includes/model_part.h" +#include "spatial_containers/bins_dynamic.h" + +namespace Kratos +{ + ///@name Kratos Classes + ///@{ + + /** + * @class SnakeSBMUtilities + * @brief Utility class for implementing Snake-based Surrogate Boundary Method (SBM). + * This class provides various functions to create and manipulate snake-like surrogate boundaries + * for finite element models in Kratos. + */ + class KRATOS_API(IGA_APPLICATION) SnakeSBMUtilities + { + public: + /** + * @brief Creates the initial snake coordinates for 2D or 3D skin. + * @param iga_model_part Model part containing the IGA geometry. + * @param skin_model_part_inner_initial Model part for the initial inner skin boundary. + * @param skin_model_part_outer_initial Model part for the initial outer skin boundary. + * @param skin_model_part Model part where the snake coordinates will be stored. + * @param rEchoLevel Verbosity level for logging. + * @param knot_vector_u Knot vector in the U-direction. + * @param knot_vector_v Knot vector in the V-direction. + * @param mParameters Input parameters for the process. + */ + static void CreateTheSnakeCoordinates(ModelPart& iga_model_part, + ModelPart& skin_model_part_inner_initial, + ModelPart& skin_model_part_outer_initial, + ModelPart& skin_model_part, + int rEchoLevel, + Vector& knot_vector_u, + Vector& knot_vector_v, + const Parameters mParameters); + + /** + * @brief Creates the snake coordinates for a specific loop (inner or outer). + * @param iga_model_part Model part containing the IGA geometry. + * @param skin_model_part_initial Model part for the initial skin boundary. + * @param skin_model_part Model part where the snake coordinates will be stored. + * @param rEchoLevel Verbosity level for logging. + * @param knot_vector_u Knot vector in the U-direction. + * @param knot_vector_v Knot vector in the V-direction. + * @param mParameters Input parameters for the process. + * @param is_inner_loop Boolean flag to indicate if the loop is inner or outer. + */ + static void CreateTheSnakeCoordinates(ModelPart& iga_model_part, + ModelPart& skin_model_part_initial, + ModelPart& skin_model_part, + int rEchoLevel, + Vector& knot_vector_u, + Vector& knot_vector_v, + const Parameters mParameters, + bool is_inner_loop); + + /** + * @brief Performs a single step in the snake algorithm for 2D models. + * @param skin_model_part The skin model part to be updated. + * @param knot_spans_available Knot spans available for the snake. + * @param idMatrixKnotSpansAvailable ID of the matrix tracking available knot spans. + * @param knot_spans_uv Knot spans in UV coordinates. + * @param xy_coord_i_cond XY coordinates of control points. + * @param knot_step_uv Step size in UV space. + * @param starting_pos_uv Starting position in UV space. + */ + static void SnakeStep(ModelPart& skin_model_part, + std::vector>> &knot_spans_available, + int idMatrixKnotSpansAvailable, + std::vector> knot_spans_uv, + std::vector> xy_coord_i_cond, + Vector knot_step_uv, + Vector starting_pos_uv); + + private: + ///@name IGA functionalities + ///@{ + using PointType = Node; + using PointTypePointer = Node::Pointer; + using PointVector = std::vector; + using PointIterator = std::vector::iterator; + using DistanceVector = std::vector; + using DistanceIterator = std::vector::iterator; + using DynamicBins = BinsDynamic<3, PointType, PointVector, PointTypePointer, PointIterator, DistanceIterator>; + using PointerType = DynamicBins::PointerType; + + /** + * @brief Checks if a point is inside the given skin boundary. + * @param pointToSearch The point to be checked. + * @param testBins Spatial bins to search for nearby points. + * @param skin_model_part_in Input skin model part. + * @return True if the point is inside the boundary, false otherwise. + */ + static bool isPointInsideSkinBoundary(Point& pointToSearch, + DynamicBins& testBins, + ModelPart& skin_model_part_in); + + /** + * @brief Marks the knot spans as available for the snake algorithm. + * @param knot_spans_available Reference to the knot spans availability matrix. + * @param idMatrix ID of the matrix being updated. + * @param testBin Spatial bins to search for nearby points. + * @param skin_model_part Skin model part being analyzed. + * @param lambda Relaxation parameter. + * @param n_knot_spans_uv Number of knot spans in UV directions. + * @param knot_step_uv Step size in UV space. + * @param starting_pos_uv Starting position in UV space. + */ + static void MarkKnotSpansAvailable(std::vector>> & knot_spans_available, + int idMatrix,DynamicBins& testBin, + ModelPart& skin_model_part, + double lambda, + std::vector& n_knot_spans_uv, + Vector& knot_step_uv, + const Vector& starting_pos_uv); + + /** + * @brief Creates the inner surrogate boundary from the snake algorithm. + * @param knot_spans_available Reference to the knot spans availability matrix. + * @param idMatrix ID of the matrix being updated. + * @param surrogate_model_part_inner Model part representing the inner boundary. + * @param n_knot_spans_uv Number of knot spans in UV directions. + * @param knot_vector_u Knot vector in the U-direction. + * @param knot_vector_v Knot vector in the V-direction. + * @param starting_pos_uv Starting position in UV space. + */ + static void CreateSurrogateBuondaryFromSnake_inner (std::vector>> & knot_spans_available, + int idMatrix, + ModelPart& surrogate_model_part_inner, + std::vector& n_knot_spans_uv, + Vector& knot_vector_u, + Vector& knot_vector_v, + const Vector& starting_pos_uv); + + /** + * @brief Creates the outer surrogate boundary from the snake algorithm. + * @param testBin_out Spatial bins for testing. + * @param initial_skin_model_part_out Initial outer skin model part. + * @param knot_spans_available Reference to the knot spans availability matrix. + * @param idMatrix ID of the matrix being updated. + * @param surrogate_model_part_outer Model part representing the outer boundary. + * @param n_knot_spans_uv Number of knot spans in UV directions. + * @param knot_vector_u Knot vector in the U-direction. + * @param knot_vector_v Knot vector in the V-direction. + * @param starting_pos_uv Starting position in UV space. + */ + static void CreateSurrogateBuondaryFromSnake_outer (DynamicBins& testBin_out, + ModelPart& initial_skin_model_part_out, + std::vector>> & knot_spans_available, + int idMatrix, + ModelPart& surrogate_model_part_outer, + std::vector& n_knot_spans_uv, + Vector& knot_vector_u, + Vector& knot_vector_v, + const Vector& starting_pos_uv); + + }; // Class SnakeSBMUtilities + +} // namespace Kratos. + +#endif // KRATOS_DIRECTOR_UTILITIES_H_INCLUDED defined From 859abbde8cf5e3940fea312db6ea74f19523ae0e Mon Sep 17 00:00:00 2001 From: andrewgorgi Date: Tue, 14 Jan 2025 16:00:47 +0100 Subject: [PATCH 2/8] Modified int->IndexType + Body fitted case in the NurbsGeometryModelerSbm MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-authored-by: Nicolò Antonelli --- .../nurbs_geometry_modeler_sbm.cpp | 15 ++-- .../create_breps_sbm_utilities.h | 69 +++++++++++++++---- .../IgaApplication/tests/test_modelers_sbm.py | 57 +++++++++++++++ .../nurbs_utilities/snake_sbm_utilities.cpp | 26 +++---- 4 files changed, 135 insertions(+), 32 deletions(-) diff --git a/applications/IgaApplication/custom_modelers/nurbs_geometry_modeler_sbm.cpp b/applications/IgaApplication/custom_modelers/nurbs_geometry_modeler_sbm.cpp index 322f8795c83e..1a76a866a453 100644 --- a/applications/IgaApplication/custom_modelers/nurbs_geometry_modeler_sbm.cpp +++ b/applications/IgaApplication/custom_modelers/nurbs_geometry_modeler_sbm.cpp @@ -69,9 +69,14 @@ namespace Kratos } // If there is not neither skin_inner nor skin_outer throw an error since you are using the sbm modeler if (!(mParameters.Has("skin_model_part_inner_initial_name") || mParameters.Has("skin_model_part_outer_initial_name"))){ - KRATOS_ERROR << "None of the 'skin_model_part_name' have not been defined " << - "in the nurbs_geometry_modeler_sbm, define it " << - "in the project paramer json" << std::endl; + + // Create the breps for the outer sbm boundary + CreateBrepsSBMUtilities CreateBrepsSBMUtilities(mEchoLevel); + CreateBrepsSBMUtilities.CreateSurrogateBoundary(mpSurface, r_model_part, A_uvw, B_uvw); + + KRATOS_WARNING("None of the 'skin_model_part_name' have not been defined ") << + "in the nurbs_geometry_modeler_sbm in the project paramer json" << std::endl; + return; } if (mParameters.Has("skin_model_part_name")) @@ -94,8 +99,8 @@ namespace Kratos // Skin model part refined after Snake Process ModelPart& skin_model_part = mpModel->CreateModelPart(skin_model_part_name); - ModelPart& skin_sub_model_part_in = skin_model_part.CreateSubModelPart("inner"); - ModelPart& skin_sub_model_part_out = skin_model_part.CreateSubModelPart("outer"); + skin_model_part.CreateSubModelPart("inner"); + skin_model_part.CreateSubModelPart("outer"); // compute unique_knot_vector_u diff --git a/applications/IgaApplication/custom_utilities/create_breps_sbm_utilities.h b/applications/IgaApplication/custom_utilities/create_breps_sbm_utilities.h index 8e508eea7d30..0505b711ad5f 100644 --- a/applications/IgaApplication/custom_utilities/create_breps_sbm_utilities.h +++ b/applications/IgaApplication/custom_utilities/create_breps_sbm_utilities.h @@ -104,17 +104,63 @@ class CreateBrepsSBMUtilities : public IO ///@name Python exposed Functions ///@{ - /// Adds all CAD geometries to the herin provided model_part. + /// Adds the surface geometry to the herin provided model_part and create the boundary breps for the SBM case. void CreateSurrogateBoundary(NurbsSurfaceGeometryPointerType& p_surface, ModelPart& rModelPart, ModelPart& rSurrogateModelPart_inner, ModelPart& rSurrogateModelPart_outer, const Point& A_uvw, const Point& B_uvw) { CreateBrepSurface(p_surface, rModelPart, rSurrogateModelPart_inner, rSurrogateModelPart_outer, mEchoLevel); CreateBrepCurveOnSurfaces(p_surface, rModelPart, rSurrogateModelPart_inner, rSurrogateModelPart_outer, A_uvw, B_uvw, mEchoLevel); } + /// Adds the surface geometry to the herin provided model_part and create the boundary breps when SBM is not needed. + void CreateSurrogateBoundary(NurbsSurfaceGeometryPointerType& p_surface, ModelPart& rModelPart, const Point& A_uvw, const Point& B_uvw) + { + CreateBrepSurface(p_surface, rModelPart, mEchoLevel); + int id_brep_curve_on_surface = 2; + CreateBrepCurvesOnRectangle(rModelPart, p_surface, A_uvw, B_uvw, id_brep_curve_on_surface); + } + private: - // 2D + /** + * @brief Create a Brep Surface object + * + * @param p_surface + * @param rModelPart + * @param EchoLevel + */ + static void CreateBrepSurface( + NurbsSurfaceGeometryPointerType p_surface, + ModelPart& rModelPart, + 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( + p_surface, + outer_loops, + inner_loops, + false); + + // Sets the brep as geometry parent of the nurbs surface. + p_surface->SetGeometryParent(p_brep_surface.get()); + p_brep_surface->SetId(1); + rModelPart.AddGeometry(p_brep_surface); + } + + /** + * @brief Create a Brep Surface object for the SBM case + * + * @param p_surface + * @param rModelPart + * @param rSurrogateModelPart_inner + * @param rSurrogateModelPart_outer + * @param EchoLevel + */ static void CreateBrepSurface( NurbsSurfaceGeometryPointerType p_surface, ModelPart& rModelPart, @@ -125,8 +171,6 @@ class CreateBrepsSBMUtilities : public IO KRATOS_INFO_IF("ReadBrepSurface", (EchoLevel > 3)) << "Creating BrepSurface \""<< std::endl; - bool case1 = rSurrogateModelPart_outer.Nodes().size() > 0; - bool case2 = rSurrogateModelPart_inner.Nodes().size() > 0; BrepCurveOnSurfaceLoopArrayType outer_loops, inner_loops; auto p_brep_surface = @@ -139,11 +183,8 @@ class CreateBrepsSBMUtilities : public IO // Sets the brep as geometry parent of the nurbs surface. p_surface->SetGeometryParent(p_brep_surface.get()); - - SizeType last_geometry_id = rModelPart.GetParentModelPart().Geometries().size(); p_brep_surface->SetId(1); rModelPart.AddGeometry(p_brep_surface); - } @@ -159,7 +200,7 @@ class CreateBrepsSBMUtilities : public IO // Each element in the surrogate_model_part represents a surrogate boundary loop. First "node" is the initial ID of the first surrogate node and // the second "node" is the last surrogate node of that loop. (We have done this in the case we have multiple surrogate boundaries and 1 model part) - int Id_brep_curve_on_surface = 2; + int id_brep_curve_on_surface = 2; if (rSurrogateModelPart_outer.Nodes().size() > 0) { int sizeSurrogateLoop_outer = rSurrogateModelPart_outer.Nodes().size(); @@ -225,19 +266,19 @@ class CreateBrepsSBMUtilities : public IO auto p_brep_curve_on_surface = Kratos::make_shared( p_surface, trimming_curves_GPT_outer[i], brep_active_range, curve_direction); - p_brep_curve_on_surface->SetId(Id_brep_curve_on_surface); + p_brep_curve_on_surface->SetId(id_brep_curve_on_surface); rModelPart.AddGeometry(p_brep_curve_on_surface); - Id_brep_curve_on_surface++; + id_brep_curve_on_surface++; } } else { - CreateBrepCurvesOnRectangle(rModelPart, p_surface, A_uvw, B_uvw, Id_brep_curve_on_surface); + CreateBrepCurvesOnRectangle(rModelPart, p_surface, A_uvw, B_uvw, id_brep_curve_on_surface); } // INNER - for (int iel = 1; iel < rSurrogateModelPart_inner.Elements().size()+1; iel++) { + for (IndexType iel = 1; iel < rSurrogateModelPart_inner.Elements().size()+1; iel++) { int firstSurrogateNodeId = rSurrogateModelPart_inner.pGetElement(iel)->GetGeometry()[0].Id(); // Element 1 because is the only surrogate loop int lastSurrogateNodeId = rSurrogateModelPart_inner.pGetElement(iel)->GetGeometry()[1].Id(); // Element 1 because is the only surrogate loop int sizeSurrogateLoop = lastSurrogateNodeId - firstSurrogateNodeId + 1; @@ -305,9 +346,9 @@ class CreateBrepsSBMUtilities : public IO auto p_brep_curve_on_surface = Kratos::make_shared( p_surface, trimming_curves_GPT[i], brep_active_range, curve_direction); - p_brep_curve_on_surface->SetId(Id_brep_curve_on_surface); + p_brep_curve_on_surface->SetId(id_brep_curve_on_surface); rModelPart.AddGeometry(p_brep_curve_on_surface); - Id_brep_curve_on_surface++; + id_brep_curve_on_surface++; trimming_brep_curve_vector[i] = p_brep_curve_on_surface ; } diff --git a/applications/IgaApplication/tests/test_modelers_sbm.py b/applications/IgaApplication/tests/test_modelers_sbm.py index c37efcd6ce7c..a72a1d429a89 100644 --- a/applications/IgaApplication/tests/test_modelers_sbm.py +++ b/applications/IgaApplication/tests/test_modelers_sbm.py @@ -18,6 +18,63 @@ def run_modelers(current_model, modelers_list): class TestModelersSbm(KratosUnittest.TestCase): + # test the call to the nurbs_geometry_modeler_sbm to create a rectangle + the breps + def test_nurbs_geometry_2d_modeler_control_points(self): + current_model = KratosMultiphysics.Model() + modeler_settings = KratosMultiphysics.Parameters(""" + [{ + "modeler_name": "NurbsGeometryModelerSbm", + "Parameters": { + "model_part_name" : "IgaModelPart", + "lower_point_xyz": [0.0,0.0,0.0], + "upper_point_xyz": [1.0,1.0,0.0], + "lower_point_uvw": [0.0,0.0,0.0], + "upper_point_uvw": [1.0,1.0,0.0], + "polynomial_order" : [4, 1], + "number_of_knot_spans" : [3,2] + } + }] + """) + run_modelers(current_model, modeler_settings) + + model_part = current_model.GetModelPart("IgaModelPart") + # Check control points + nodes_mp = model_part.NodesArray(0) + + for i in range(3): + self.assertAlmostEqual(nodes_mp[i*7+0].X,0.0) + self.assertAlmostEqual(nodes_mp[i*7+1].X,0.083333333334) + self.assertAlmostEqual(nodes_mp[i*7+2].X,0.25) + self.assertAlmostEqual(nodes_mp[i*7+3].X,0.5) + self.assertAlmostEqual(nodes_mp[i*7+4].X,0.75) + self.assertAlmostEqual(nodes_mp[i*7+5].X,0.916666666667) + self.assertAlmostEqual(nodes_mp[i*7+6].X,1.0) + for i in range(7): + self.assertAlmostEqual(nodes_mp[i].Y,0.0) + self.assertAlmostEqual(nodes_mp[i].Z,0.0) + for i in range(7,14): + self.assertAlmostEqual(nodes_mp[i].Y,0.5) + self.assertAlmostEqual(nodes_mp[i].Z,0.0) + for i in range(14,21): + self.assertAlmostEqual(nodes_mp[i].Y,1.0) + self.assertAlmostEqual(nodes_mp[i].Z,0.0) + + geometry = model_part.GetGeometry(1) + + # Check if geometry holds same nodes as model part + for i, node_geo in enumerate(geometry): + self.assertEqual(node_geo.Id, nodes_mp[i].Id ) + self.assertEqual(node_geo.X, nodes_mp[i].X ) + self.assertEqual(node_geo.Y, nodes_mp[i].Y ) + self.assertEqual(node_geo.Z, nodes_mp[i].Z ) + + # test the creation of the breps + self.assertEqual(current_model["IgaModelPart"].NumberOfGeometries(), 5) + self.assertAlmostEqual(current_model["IgaModelPart"].GetGeometry(1).Center().X, 0.5) + self.assertAlmostEqual(current_model["IgaModelPart"].GetGeometry(1).Center().Y, 0.5) + self.assertAlmostEqual(current_model["IgaModelPart"].GetGeometry(1).Center().Z, 0.0) + + # test for SBM def test_nurbs_geometry_2d_modeler_sbm_outer_skin_boundary(self): current_model = KratosMultiphysics.Model() diff --git a/kratos/utilities/nurbs_utilities/snake_sbm_utilities.cpp b/kratos/utilities/nurbs_utilities/snake_sbm_utilities.cpp index 7c993bff971a..4281c2fba5f2 100644 --- a/kratos/utilities/nurbs_utilities/snake_sbm_utilities.cpp +++ b/kratos/utilities/nurbs_utilities/snake_sbm_utilities.cpp @@ -108,7 +108,7 @@ namespace Kratos std::vector>> knot_spans_available; knot_spans_available.reserve(numberOfLoops); - for (int i = 0; i < numberOfLoops; ++i) { + for (IndexType i = 0; i < numberOfLoops; ++i) { std::vector> matrix; matrix.reserve(n_knot_spans_uv[1]); for (int j = 0; j <= n_knot_spans_uv[1]; ++j) { @@ -120,7 +120,7 @@ namespace Kratos // Optimized Snake -> for inner loops int idMatrixKnotSpansAvailable = 0; - int idFirstNode; + IndexType idFirstNode; bool newInnerLoop = true; @@ -214,8 +214,8 @@ namespace Kratos KRATOS_INFO_IF("::[SnakeSBMUtilities]::", rEchoLevel > 0 && is_inner) << "Inner :: MarkKnotSpansAvailable" << std::endl; KRATOS_INFO_IF("::[SnakeSBMUtilities]::", rEchoLevel > 0 && !is_inner) << "Outer :: MarkKnotSpansAvailable" << std::endl; - for (int i = 0; i < numberOfLoops; i++) { - int idInnerLoop = i; + for (IndexType i = 0; i < numberOfLoops; i++) { + IndexType idInnerLoop = i; // Mark the knot_spans_available's for inner and outer loops MarkKnotSpansAvailable(knot_spans_available, idInnerLoop, testBins, skin_sub_model_part, lambda, n_knot_spans_uv, knot_step_uv, starting_pos_uv); @@ -329,10 +329,10 @@ namespace Kratos PointerType nearestPoint = testBins.SearchNearestPoint(*pointToSearch); // Get the closest Condition the initial_skin_model_part_in.Conditions - int id1 = nearestPoint->Id(); + IndexType id1 = nearestPoint->Id(); auto nearestCondition1 = skin_model_part.GetCondition(id1); // Check if the condition is the first one and therefore the previous one does not exist - int id2 = id1 - 1; + IndexType id2 = id1 - 1; if (id1 == skin_model_part.ConditionsBegin()->Id()) { int nConditions = skin_model_part.Conditions().size(); id2 = id1 + nConditions - 1; @@ -424,11 +424,11 @@ namespace Kratos // Create 25 "fake" GaussPoints to check if the majority are inside or outside const int numFakeGaussPoints = 5; int numberOfInsideGaussianPoints = 0; - for (int i_GPx = 0; i_GPx < numFakeGaussPoints; i_GPx++){ + for (IndexType i_GPx = 0; i_GPx < numFakeGaussPoints; i_GPx++){ double x_coord = j*knot_step_uv[0] + knot_step_uv[0]/(numFakeGaussPoints+1)*(i_GPx+1) + starting_pos_uv[0]; // NOTE:: The v-knot spans are upside down in the matrix!! - for (int i_GPy = 0; i_GPy < numFakeGaussPoints; i_GPy++) + for (IndexType i_GPy = 0; i_GPy < numFakeGaussPoints; i_GPy++) { double y_coord = i*knot_step_uv[1] + knot_step_uv[1]/(numFakeGaussPoints+1)*(i_GPy+1) + starting_pos_uv[1]; Point gaussPoint = Point(x_coord, y_coord, 0); // GAUSSIAN POINT @@ -646,7 +646,7 @@ namespace Kratos } } // TOP BOUNDARY - for (int j = knot_spans_available[idMatrix][0].size()-1; j > knot_spans_available[idMatrix][0].size()-3; j--) { + for (int j = int(knot_spans_available[idMatrix][0].size())-1; j > int(knot_spans_available[idMatrix][0].size())-3; j--) { for (int i = 0; i < (n_knot_spans_uv[0]); i++) { Point centroidKnotSpan = Point((j+0.5)*knot_step_u+starting_pos_uv[0], (i+0.5)*knot_step_v+starting_pos_uv[1], 0); if (isPointInsideSkinBoundary(centroidKnotSpan, testBin_out, initial_skin_model_part_out) && knot_spans_available[idMatrix][i][j] != -1) { @@ -655,7 +655,7 @@ namespace Kratos } } // RIGHT BOUNDARY - for (int i = knot_spans_available[idMatrix][0].size()-1; i > knot_spans_available[idMatrix][0].size()-3; i--) { + for (int i = int(knot_spans_available[idMatrix][0].size())-1; i > int(knot_spans_available[idMatrix][0].size())-3; i--) { for (int j = n_knot_spans_uv[1]-1; j > -1; j-- ) { Point centroidKnotSpan = Point((j+0.5)*knot_step_u+starting_pos_uv[0], (i+0.5)*knot_step_v+starting_pos_uv[1], 0); if (isPointInsideSkinBoundary(centroidKnotSpan, testBin_out, initial_skin_model_part_out) && knot_spans_available[idMatrix][i][j] != -1) { @@ -696,8 +696,8 @@ namespace Kratos //------------------- std::vector> knot_spans_available_extended(n_knot_spans_uv[1]+2, std::vector(n_knot_spans_uv[0]+2)); - for (int i = 0; i < knot_spans_available[idMatrix].size(); i++){ - for (int j = 0; j < knot_spans_available[idMatrix][0].size(); j++) { + for (IndexType i = 0; i < knot_spans_available[idMatrix].size(); i++){ + for (IndexType j = 0; j < knot_spans_available[idMatrix][0].size(); j++) { knot_spans_available_extended[i+1][j+1] = knot_spans_available[idMatrix][i][j]; } } @@ -833,7 +833,7 @@ namespace Kratos } // Create "fictituos element" to memorize starting and ending node id for each surrogate boundary loop std::vector elem_nodes{1, idSnakeNode-1}; - int elem_id = surrogate_model_part_outer.GetRootModelPart().Elements().size()+1; + IndexType elem_id = surrogate_model_part_outer.GetRootModelPart().Elements().size()+1; surrogate_model_part_outer.CreateNewElement("Element2D2N", elem_id, elem_nodes, p_cond_prop); } } // namespace Kratos. From 7a518dd9d1e0b2ffa14fd9aa4d8c246b9844e365 Mon Sep 17 00:00:00 2001 From: andrewgorgi Date: Wed, 15 Jan 2025 11:20:10 +0100 Subject: [PATCH 3/8] Correct the typo in test_IgaApplication.py --- applications/IgaApplication/tests/test_IgaApplication.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/applications/IgaApplication/tests/test_IgaApplication.py b/applications/IgaApplication/tests/test_IgaApplication.py index 82efb0db1550..18e8969600c7 100644 --- a/applications/IgaApplication/tests/test_IgaApplication.py +++ b/applications/IgaApplication/tests/test_IgaApplication.py @@ -45,7 +45,7 @@ from test_nurbs_volume_element import TestNurbsVolumeElement as TTestNurbsVolumeElements # Modelers tests from test_modelers import TestModelers as TTestModelers -from test_modelers_sbm import TestModelers as TTestModelersSbm +from test_modelers_sbm import TestModelersSbm as TTestModelersSbm # Processes tests from test_map_nurbs_volume_results_to_embedded_geometry_process import TestMapNurbsVolumeResultsToEmbeddedGeometryProcess as TTestMapNurbsVolumeResultsToEmbeddedGeometryProcess From 9f49410b93c46624bc8eec4e21b43e92d70a0e27 Mon Sep 17 00:00:00 2001 From: andrewgorgi Date: Thu, 16 Jan 2025 18:05:57 +0100 Subject: [PATCH 4/8] Suggestions by Alejandro MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-authored-by: Nicolò Antonelli --- .../nurbs_geometry_modeler_sbm.cpp | 11 ++-- .../nurbs_geometry_modeler_sbm.h | 34 +++++-------- .../create_breps_sbm_utilities.h | 50 +++++++++---------- .../brep_trimming_utilities.cpp | 3 -- .../brep_trimming_utilities.h | 7 +-- .../nurbs_utilities/snake_sbm_utilities.cpp | 47 ++++++++--------- .../nurbs_utilities/snake_sbm_utilities.h | 7 +-- 7 files changed, 65 insertions(+), 94 deletions(-) diff --git a/applications/IgaApplication/custom_modelers/nurbs_geometry_modeler_sbm.cpp b/applications/IgaApplication/custom_modelers/nurbs_geometry_modeler_sbm.cpp index 1a76a866a453..ef2738013065 100644 --- a/applications/IgaApplication/custom_modelers/nurbs_geometry_modeler_sbm.cpp +++ b/applications/IgaApplication/custom_modelers/nurbs_geometry_modeler_sbm.cpp @@ -12,10 +12,9 @@ // // Project includes -#include "includes/define.h" #include "nurbs_geometry_modeler_sbm.h" -#include "geometries/nurbs_shape_function_utilities/nurbs_volume_refinement_utilities.h" #include "custom_utilities/create_breps_sbm_utilities.h" +#include "utilities/nurbs_utilities/snake_sbm_utilities.h" namespace Kratos { @@ -56,15 +55,15 @@ namespace Kratos if (mParameters.Has("skin_model_part_inner_initial_name")) { skin_model_part_inner_initial_name = mParameters["skin_model_part_inner_initial_name"].GetString(); - if (!mpModel->HasModelPart(skin_model_part_inner_initial_name)) - KRATOS_ERROR << "The skin_model_part '" << skin_model_part_inner_initial_name << "' was not created in the model.\n" + KRATOS_ERROR_IF_NOT(mpModel->HasModelPart(skin_model_part_inner_initial_name)) + << "The skin_model_part '" << skin_model_part_inner_initial_name << "' was not created in the model.\n" << "Check the reading of the mdpa file in the import mdpa modeler."<< std::endl; } if (mParameters.Has("skin_model_part_outer_initial_name")) { skin_model_part_outer_initial_name = mParameters["skin_model_part_outer_initial_name"].GetString(); - if (!mpModel->HasModelPart(skin_model_part_outer_initial_name)) - KRATOS_ERROR << "The skin_model_part '" << skin_model_part_outer_initial_name << "' was not created in the model.\n" + KRATOS_ERROR_IF_NOT(mpModel->HasModelPart(skin_model_part_outer_initial_name)) + << "The skin_model_part '" << skin_model_part_outer_initial_name << "' was not created in the model.\n" << "Check the reading of the mdpa file in the import mdpa modeler."<< std::endl; } // If there is not neither skin_inner nor skin_outer throw an error since you are using the sbm modeler diff --git a/applications/IgaApplication/custom_modelers/nurbs_geometry_modeler_sbm.h b/applications/IgaApplication/custom_modelers/nurbs_geometry_modeler_sbm.h index 92525dc60e04..3d3e2a1a643d 100644 --- a/applications/IgaApplication/custom_modelers/nurbs_geometry_modeler_sbm.h +++ b/applications/IgaApplication/custom_modelers/nurbs_geometry_modeler_sbm.h @@ -11,21 +11,14 @@ // Andrea Gorgi // -#if !defined(KRATOS_NURBS_GEOMETRY_MODELER_SBM_H_INCLUDED ) -#define KRATOS_NURBS_GEOMETRY_MODELER_SBM_H_INCLUDED +# pragma once // System includes // External includes // Project includes -#include "includes/model_part.h" #include "nurbs_geometry_modeler.h" -#include "geometries/nurbs_volume_geometry.h" -#include "geometries/nurbs_surface_geometry.h" -#include "geometries/nurbs_shape_function_utilities/nurbs_surface_refinement_utilities.h" -#include "geometries/brep_curve_on_surface.h" -#include "utilities/nurbs_utilities/snake_sbm_utilities.h" namespace Kratos { @@ -37,21 +30,21 @@ class KRATOS_API(IGA_APPLICATION) NurbsGeometryModelerSbm ///@{ KRATOS_CLASS_POINTER_DEFINITION( NurbsGeometryModelerSbm ); - typedef std::size_t IndexType; - typedef std::size_t SizeType; - typedef Node NodeType; + using IndexType = std::size_t; + using SizeType = std::size_t; + using NodeType = Node; - typedef Geometry GeometryType; - typedef typename GeometryType::Pointer GeometryPointerType; + using GeometryType = Geometry; + using GeometryPointerType = GeometryType::Pointer; - typedef NurbsSurfaceGeometry<3, PointerVector> NurbsSurfaceGeometryType; - typedef typename NurbsSurfaceGeometryType::Pointer NurbsSurfaceGeometryPointerType; + using NurbsSurfaceGeometryType = NurbsSurfaceGeometry<3, PointerVector>; + using NurbsSurfaceGeometryPointerType = NurbsSurfaceGeometryType::Pointer; - typedef NurbsVolumeGeometry> NurbsVolumeGeometryType; - typedef typename NurbsVolumeGeometryType::Pointer NurbsVolumeGeometryPointerType; + using NurbsVolumeGeometryType = NurbsVolumeGeometry>; + using NurbsVolumeGeometryPointerType = NurbsVolumeGeometryType::Pointer; - typedef PointerVector ContainerNodeType; - typedef PointerVector ContainerEmbeddedNodeType; + using ContainerNodeType = PointerVector; + using ContainerEmbeddedNodeType = PointerVector; ///@} ///@name Life Cycle @@ -125,5 +118,4 @@ class KRATOS_API(IGA_APPLICATION) NurbsGeometryModelerSbm Parameters ReadParamatersFile(const std::string& rDataFileName) const; }; -} // End namesapce Kratos -#endif // KRATOS_NURBS_GEOMETRY_MODELER_H_INCLUDED \ No newline at end of file +} // End namesapce Kratos \ No newline at end of file diff --git a/applications/IgaApplication/custom_utilities/create_breps_sbm_utilities.h b/applications/IgaApplication/custom_utilities/create_breps_sbm_utilities.h index 0505b711ad5f..5c41ab4c6128 100644 --- a/applications/IgaApplication/custom_utilities/create_breps_sbm_utilities.h +++ b/applications/IgaApplication/custom_utilities/create_breps_sbm_utilities.h @@ -11,9 +11,7 @@ // Andrea Gorgi // -#if !defined(KRATOS_CREATE_BREPS_SBM_UTILITIES_INCLUDED ) -#define KRATOS_CREATE_BREPS_SBM_UTILITIES_INCLUDED - +# pragma once // System includes @@ -55,36 +53,36 @@ class CreateBrepsSBMUtilities : public IO /// Pointer definition of CreateBrepsSBMUtilities KRATOS_CLASS_POINTER_DEFINITION(CreateBrepsSBMUtilities); - typedef std::size_t SizeType; - typedef std::size_t IndexType; + using SizeType = std::size_t; + using IndexType = std::size_t; - typedef Geometry GeometryType; - typedef typename GeometryType::Pointer GeometryPointerType; + using GeometryType = Geometry; + using GeometryPointerType = typename GeometryType::Pointer; - typedef PointerVector ContainerNodeType; - typedef PointerVector ContainerEmbeddedNodeType; + using ContainerNodeType = PointerVector; + using ContainerEmbeddedNodeType = PointerVector; - typedef CouplingGeometry CouplingGeometryType; + using CouplingGeometryType = CouplingGeometry; - typedef NurbsSurfaceGeometry<3, ContainerNodeType> NurbsSurfaceType; - typedef NurbsCurveGeometry<2, ContainerEmbeddedNodeType> NurbsTrimmingCurveType; + using NurbsSurfaceType = NurbsSurfaceGeometry<3, ContainerNodeType>; + using NurbsTrimmingCurveType = NurbsCurveGeometry<2, ContainerEmbeddedNodeType>; - typedef typename NurbsSurfaceType::Pointer NurbsSurfacePointerType; - typedef typename NurbsTrimmingCurveType::Pointer NurbsTrimmingCurvePointerType; + using NurbsSurfacePointerType = typename NurbsSurfaceType::Pointer; + using NurbsTrimmingCurvePointerType = typename NurbsTrimmingCurveType::Pointer; - typedef NurbsSurfaceGeometry<3, PointerVector> NurbsSurfaceGeometryType; - typedef typename NurbsSurfaceGeometryType::Pointer NurbsSurfaceGeometryPointerType; + using NurbsSurfaceGeometryType = NurbsSurfaceGeometry<3, PointerVector>; + using NurbsSurfaceGeometryPointerType = typename NurbsSurfaceGeometryType::Pointer; - typedef BrepSurface BrepSurfaceType; - typedef BrepCurveOnSurface BrepCurveOnSurfaceType; - - typedef BrepCurve BrepCurveType; - typedef PointOnGeometry PointOnGeometryOnSurfaceType; - typedef PointOnGeometry PointOnGeometryOnCurveType; + using BrepSurfaceType = BrepSurface; + using BrepCurveOnSurfaceType = BrepCurveOnSurface; - typedef DenseVector BrepCurveOnSurfaceArrayType; - typedef DenseVector BrepCurveOnSurfaceLoopType; - typedef DenseVector> BrepCurveOnSurfaceLoopArrayType; + using BrepCurveType = BrepCurve; + using PointOnGeometryOnSurfaceType = PointOnGeometry; + using PointOnGeometryOnCurveType = PointOnGeometry; + + using BrepCurveOnSurfaceArrayType = DenseVector; + using BrepCurveOnSurfaceLoopType = DenseVector; + using BrepCurveOnSurfaceLoopArrayType = DenseVector>; ///@} ///@name Life Cycle @@ -443,5 +441,3 @@ class CreateBrepsSBMUtilities : public IO ///@} }; // Class CreateBrepsSBMUtilities } // namespace Kratos. - -#endif diff --git a/kratos/utilities/geometry_utilities/brep_trimming_utilities.cpp b/kratos/utilities/geometry_utilities/brep_trimming_utilities.cpp index 7f98d2b5424a..d77b2368a967 100644 --- a/kratos/utilities/geometry_utilities/brep_trimming_utilities.cpp +++ b/kratos/utilities/geometry_utilities/brep_trimming_utilities.cpp @@ -6,9 +6,6 @@ // // License: BSD License // Kratos default license: kratos/license.txt -// -// Main authors: Nicolo' Antonelli -// Andrea Gorgi // header includes #include "brep_trimming_utilities.h" diff --git a/kratos/utilities/geometry_utilities/brep_trimming_utilities.h b/kratos/utilities/geometry_utilities/brep_trimming_utilities.h index 5162ac6842e6..dcf4fc6086a4 100644 --- a/kratos/utilities/geometry_utilities/brep_trimming_utilities.h +++ b/kratos/utilities/geometry_utilities/brep_trimming_utilities.h @@ -6,12 +6,8 @@ // // License: BSD License // Kratos default license: kratos/license.txt -// -// Main authors: Nicolo' Antonelli -// Andrea Gorgi -#if !defined(KRATOS_BREP_TRIMMING_UTILITIES_H_INCLUDED) -#define KRATOS_BREP_TRIMMING_UTILITIES_H_INCLUDED +#pragma once // Std includes #include @@ -329,4 +325,3 @@ namespace Kratos ///@} // Kratos Classes } // namespace Kratos. -#endif // KRATOS_BREP_TRIMMING_UTILITIES_H_INCLUDED defined diff --git a/kratos/utilities/nurbs_utilities/snake_sbm_utilities.cpp b/kratos/utilities/nurbs_utilities/snake_sbm_utilities.cpp index 4281c2fba5f2..5f3b02e80d2f 100644 --- a/kratos/utilities/nurbs_utilities/snake_sbm_utilities.cpp +++ b/kratos/utilities/nurbs_utilities/snake_sbm_utilities.cpp @@ -72,7 +72,7 @@ namespace Kratos Properties::Pointer p_cond_prop_in = skin_model_part_initial.pGetProperties(0); skin_model_part_initial.AddProperties(p_cond_prop_in); - Vector knot_step_uv(2); + array_1d knot_step_uv(2); knot_step_uv[0] = std::abs(knot_vector_u[int(knot_vector_u.size()/2) +1] - knot_vector_u[int(knot_vector_u.size()/2)] ) ; knot_step_uv[1] = std::abs(knot_vector_v[int(knot_vector_v.size()/2) +1] - knot_vector_v[int(knot_vector_v.size()/2)] ) ; @@ -82,7 +82,7 @@ namespace Kratos ModelPart& surrogate_model_part = iga_model_part.GetSubModelPart(surrogate_sub_model_part_name); surrogate_model_part.GetProcessInfo().SetValue(MARKER_MESHES, meshSizes_uv); - Vector starting_pos_uv(2); + array_1d starting_pos_uv; starting_pos_uv[0] = knot_vector_u[0]; starting_pos_uv[1] = knot_vector_v[0]; @@ -143,27 +143,25 @@ namespace Kratos newInnerLoop = false; } // Collect the coordinates of the points of the i_cond - double x_true_boundary1 = i_cond.GetGeometry()[0].X(); - double y_true_boundary1 = i_cond.GetGeometry()[0].Y(); - double x_true_boundary2 = i_cond.GetGeometry()[1].X(); - double y_true_boundary2 = i_cond.GetGeometry()[1].Y(); + const auto& r_coords_true_boundary1 = i_cond.GetGeometry()[0].Coordinates(); + const auto& r_coords_true_boundary2 = i_cond.GetGeometry()[1].Coordinates(); std::vector> xy_coord_i_cond(2); xy_coord_i_cond[0].resize(2); xy_coord_i_cond[1].resize(2); - xy_coord_i_cond[0][0] = i_cond.GetGeometry()[0].X(); // x_true_boundary1 - xy_coord_i_cond[1][0] = i_cond.GetGeometry()[0].Y(); // y_true_boundary1 - xy_coord_i_cond[0][1] = i_cond.GetGeometry()[1].X(); // x_true_boundary2 - xy_coord_i_cond[1][1] = i_cond.GetGeometry()[1].Y(); // y_true_boundary2 + xy_coord_i_cond[0][0] = r_coords_true_boundary1[0]; + xy_coord_i_cond[1][0] = r_coords_true_boundary1[1]; + xy_coord_i_cond[0][1] = r_coords_true_boundary2[0]; + xy_coord_i_cond[1][1] = r_coords_true_boundary2[1]; // Collect the intersections of the skin boundary with the knot values std::vector> knot_span_uv(2); knot_span_uv[0].resize(2); knot_span_uv[1].resize(2); - knot_span_uv[0][0] = (x_true_boundary1-starting_pos_uv[0]) / knot_step_uv[0]; // knot_span_u_1st_point - knot_span_uv[1][0] = (y_true_boundary1-starting_pos_uv[1]) / knot_step_uv[1]; // knot_span_v_1st_point - knot_span_uv[0][1] = (x_true_boundary2-starting_pos_uv[0]) / knot_step_uv[0]; // knot_span_u_2nd_point - knot_span_uv[1][1] = (y_true_boundary2-starting_pos_uv[1]) / knot_step_uv[1]; // knot_span_v_2nd_point + knot_span_uv[0][0] = (r_coords_true_boundary1[0]-starting_pos_uv[0]) / knot_step_uv[0]; // knot_span_u_1st_point + knot_span_uv[1][0] = (r_coords_true_boundary1[1]-starting_pos_uv[1]) / knot_step_uv[1]; // knot_span_v_1st_point + knot_span_uv[0][1] = (r_coords_true_boundary2[0]-starting_pos_uv[0]) / knot_step_uv[0]; // knot_span_u_2nd_point + knot_span_uv[1][1] = (r_coords_true_boundary2[1]-starting_pos_uv[1]) / knot_step_uv[1]; // knot_span_v_2nd_point // In the inner case, check is the immersed object is inside the rectangular domain if (is_inner && @@ -365,7 +363,7 @@ namespace Kratos void SnakeSBMUtilities::MarkKnotSpansAvailable(std::vector>> & knot_spans_available, int idMatrix,DynamicBins& testBin, ModelPart& skin_model_part, double lambda, std::vector& n_knot_spans_uv, - Vector& knot_step_uv, const Vector& starting_pos_uv) { + array_1d& knot_step_uv, const Vector& starting_pos_uv) { for (int i = 0; i < n_knot_spans_uv[1]; i++) { for (int j = 0; j < n_knot_spans_uv[0]; j++) { if (knot_spans_available[idMatrix][i][j] == 2) { @@ -447,9 +445,6 @@ namespace Kratos else{ knot_spans_available[idMatrix][i][j] = 1; // The knot span is considered DEACTIVE } - - // exit(0); - } } } @@ -489,16 +484,16 @@ namespace Kratos IndexType idSnakeNode = id_first_node+1; // Follow the clockwise loop - int end = 0; + bool end = false; // We are going orizontally - int direction = 0; + std::size_t direction = 0; // 0 = up_vertical, 1 = right_orizontal, 2 = down_vertical, 3 = left_orizontal int i = start_i; int j = start_j; int I = start_i; int J = start_j; int steps = 0; const int max_number_of_steps = 1e5; - while (end == 0 && steps < max_number_of_steps) { + while (!end && steps < max_number_of_steps) { steps ++ ; int is_special_case = 1; // Variable to check if we are in a super special case and we shouldn't go out of the while loop // Try to go left w.r.t. the current direction @@ -608,7 +603,7 @@ namespace Kratos // Check if we have close the snake if (I == start_i && J == start_j && is_special_case == 1 ) { // End of the while loop - end = 1; + end = true; KRATOS_INFO("number of steps in the snake step") << steps << std::endl; } } @@ -711,16 +706,16 @@ namespace Kratos IndexType idSnakeNode = id_first_node+1; // Follow the clockwise loop - int end = 0; + bool end = false; // We are going orizontally - int direction = 0 ; // 0 = up_vertical, 1 = right_orizontal, 2 = down_vertical, 3 = left_orizontal + std::size_t direction = 0 ; // 0 = up_vertical, 1 = right_orizontal, 2 = down_vertical, 3 = left_orizontal int I = start_i+1; int J = start_j+1; int i = start_i; int j = start_j; int steps = 0; const int max_number_of_steps = 1e5; - while (end == 0 && steps < max_number_of_steps) { + while (!end && steps < max_number_of_steps) { steps ++ ; int is_special_case = 1; // Variable to check if we are in a super special case and we shouldn't go out of the while loop // Try to go left w.r.t. the current direction @@ -828,7 +823,7 @@ namespace Kratos // Check if we have close the snake if (I == start_i+1 && J == start_j+1 && is_special_case == 1 ) { // End of the while loop - end = 1; + end = true; } } // Create "fictituos element" to memorize starting and ending node id for each surrogate boundary loop diff --git a/kratos/utilities/nurbs_utilities/snake_sbm_utilities.h b/kratos/utilities/nurbs_utilities/snake_sbm_utilities.h index ae456ed545a5..008c643a76c6 100644 --- a/kratos/utilities/nurbs_utilities/snake_sbm_utilities.h +++ b/kratos/utilities/nurbs_utilities/snake_sbm_utilities.h @@ -11,8 +11,7 @@ // Andrea Gorgi // -#if !defined(KRATOS_SNAKE_SBM_UTILITIES_H_INCLUDED ) -#define KRATOS_SNAKE_SBM_UTILITIES_H_INCLUDED +#pragma once // Project includes #include "includes/model_part.h" @@ -129,7 +128,7 @@ namespace Kratos ModelPart& skin_model_part, double lambda, std::vector& n_knot_spans_uv, - Vector& knot_step_uv, + array_1d& knot_step_uv, const Vector& starting_pos_uv); /** @@ -175,5 +174,3 @@ namespace Kratos }; // Class SnakeSBMUtilities } // namespace Kratos. - -#endif // KRATOS_DIRECTOR_UTILITIES_H_INCLUDED defined From 625d3af8a885106d7f255fd8cc51a52dcea1d17a Mon Sep 17 00:00:00 2001 From: andrewgorgi Date: Thu, 16 Jan 2025 20:15:26 +0100 Subject: [PATCH 5/8] direction variable converted to integer (not SizeType) --- kratos/utilities/nurbs_utilities/snake_sbm_utilities.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/kratos/utilities/nurbs_utilities/snake_sbm_utilities.cpp b/kratos/utilities/nurbs_utilities/snake_sbm_utilities.cpp index 5f3b02e80d2f..c774bcf3c1b1 100644 --- a/kratos/utilities/nurbs_utilities/snake_sbm_utilities.cpp +++ b/kratos/utilities/nurbs_utilities/snake_sbm_utilities.cpp @@ -486,7 +486,7 @@ namespace Kratos // Follow the clockwise loop bool end = false; // We are going orizontally - std::size_t direction = 0; + int direction = 0; // 0 = up_vertical, 1 = right_orizontal, 2 = down_vertical, 3 = left_orizontal int i = start_i; int j = start_j; int I = start_i; int J = start_j; @@ -708,7 +708,7 @@ namespace Kratos // Follow the clockwise loop bool end = false; // We are going orizontally - std::size_t direction = 0 ; // 0 = up_vertical, 1 = right_orizontal, 2 = down_vertical, 3 = left_orizontal + int direction = 0 ; // 0 = up_vertical, 1 = right_orizontal, 2 = down_vertical, 3 = left_orizontal int I = start_i+1; int J = start_j+1; int i = start_i; From 05c029036524580229e28b3a11664979a6419582 Mon Sep 17 00:00:00 2001 From: andrewgorgi Date: Wed, 29 Jan 2025 16:05:14 +0100 Subject: [PATCH 6/8] Correction to SnakeSbmUtilities + PR suggestions by Ricky --- .../custom_modelers/nurbs_geometry_modeler.h | 2 +- .../nurbs_geometry_modeler_sbm.h | 2 +- .../create_breps_sbm_utilities.h | 38 ++++++------ kratos/geometries/brep_surface.h | 10 ++-- .../nurbs_utilities.h | 11 +++- .../nurbs_utilities/snake_sbm_utilities.cpp | 60 ++++++++++--------- .../nurbs_utilities/snake_sbm_utilities.h | 6 +- 7 files changed, 70 insertions(+), 59 deletions(-) diff --git a/applications/IgaApplication/custom_modelers/nurbs_geometry_modeler.h b/applications/IgaApplication/custom_modelers/nurbs_geometry_modeler.h index 6496b90a0fbb..735f34957007 100644 --- a/applications/IgaApplication/custom_modelers/nurbs_geometry_modeler.h +++ b/applications/IgaApplication/custom_modelers/nurbs_geometry_modeler.h @@ -66,7 +66,7 @@ class KRATOS_API(IGA_APPLICATION) NurbsGeometryModeler } /// Destructor. - virtual ~NurbsGeometryModeler() = default; + ~NurbsGeometryModeler() = default; /// Creates the Modeler Pointer Modeler::Pointer Create(Model& rModel, const Parameters ModelParameters) const override diff --git a/applications/IgaApplication/custom_modelers/nurbs_geometry_modeler_sbm.h b/applications/IgaApplication/custom_modelers/nurbs_geometry_modeler_sbm.h index 3d3e2a1a643d..9593506c9dda 100644 --- a/applications/IgaApplication/custom_modelers/nurbs_geometry_modeler_sbm.h +++ b/applications/IgaApplication/custom_modelers/nurbs_geometry_modeler_sbm.h @@ -64,7 +64,7 @@ class KRATOS_API(IGA_APPLICATION) NurbsGeometryModelerSbm } /// Destructor. - virtual ~NurbsGeometryModelerSbm() = default; + ~NurbsGeometryModelerSbm() = default; /// Creates the Modeler Pointer Modeler::Pointer Create(Model& rModel, const Parameters ModelParameters) const override diff --git a/applications/IgaApplication/custom_utilities/create_breps_sbm_utilities.h b/applications/IgaApplication/custom_utilities/create_breps_sbm_utilities.h index 5c41ab4c6128..a1ba88435ce6 100644 --- a/applications/IgaApplication/custom_utilities/create_breps_sbm_utilities.h +++ b/applications/IgaApplication/custom_utilities/create_breps_sbm_utilities.h @@ -103,10 +103,10 @@ class CreateBrepsSBMUtilities : public IO ///@{ /// Adds the surface geometry to the herin provided model_part and create the boundary breps for the SBM case. - void CreateSurrogateBoundary(NurbsSurfaceGeometryPointerType& p_surface, ModelPart& rModelPart, ModelPart& rSurrogateModelPart_inner, ModelPart& rSurrogateModelPart_outer, const Point& A_uvw, const Point& B_uvw) + void CreateSurrogateBoundary(NurbsSurfaceGeometryPointerType& p_surface, ModelPart& rModelPart, ModelPart& rSurrogateModelPartInner, ModelPart& rSurrogateModelPartOuter, const Point& A_uvw, const Point& B_uvw) { - CreateBrepSurface(p_surface, rModelPart, rSurrogateModelPart_inner, rSurrogateModelPart_outer, mEchoLevel); - CreateBrepCurveOnSurfaces(p_surface, rModelPart, rSurrogateModelPart_inner, rSurrogateModelPart_outer, A_uvw, B_uvw, mEchoLevel); + CreateBrepSurface(p_surface, rModelPart, rSurrogateModelPartInner, rSurrogateModelPartOuter, mEchoLevel); + CreateBrepCurveOnSurfaces(p_surface, rModelPart, rSurrogateModelPartInner, rSurrogateModelPartOuter, A_uvw, B_uvw, mEchoLevel); } /// Adds the surface geometry to the herin provided model_part and create the boundary breps when SBM is not needed. @@ -155,15 +155,15 @@ class CreateBrepsSBMUtilities : public IO * * @param p_surface * @param rModelPart - * @param rSurrogateModelPart_inner - * @param rSurrogateModelPart_outer + * @param rSurrogateModelPartInner + * @param rSurrogateModelPartOuter * @param EchoLevel */ static void CreateBrepSurface( NurbsSurfaceGeometryPointerType p_surface, ModelPart& rModelPart, - ModelPart& rSurrogateModelPart_inner, - ModelPart& rSurrogateModelPart_outer, + ModelPart& rSurrogateModelPartInner, + ModelPart& rSurrogateModelPartOuter, SizeType EchoLevel = 0) { KRATOS_INFO_IF("ReadBrepSurface", (EchoLevel > 3)) @@ -176,8 +176,8 @@ class CreateBrepsSBMUtilities : public IO p_surface, outer_loops, inner_loops, - rSurrogateModelPart_inner, - rSurrogateModelPart_outer); + rSurrogateModelPartInner, + rSurrogateModelPartOuter); // Sets the brep as geometry parent of the nurbs surface. p_surface->SetGeometryParent(p_brep_surface.get()); @@ -189,8 +189,8 @@ class CreateBrepsSBMUtilities : public IO static void CreateBrepCurveOnSurfaces( NurbsSurfaceGeometryPointerType p_surface, ModelPart& rModelPart, - ModelPart& rSurrogateModelPart_inner, - ModelPart& rSurrogateModelPart_outer, + ModelPart& rSurrogateModelPartInner, + ModelPart& rSurrogateModelPartOuter, const Point& A_uvw, const Point& B_uvw, SizeType EchoLevel = 0) { @@ -200,12 +200,12 @@ class CreateBrepsSBMUtilities : public IO int id_brep_curve_on_surface = 2; - if (rSurrogateModelPart_outer.Nodes().size() > 0) { - int sizeSurrogateLoop_outer = rSurrogateModelPart_outer.Nodes().size(); + if (rSurrogateModelPartOuter.Nodes().size() > 0) { + int sizeSurrogateLoop_outer = rSurrogateModelPartOuter.Nodes().size(); std::vector surrogatecoord_x_outer(sizeSurrogateLoop_outer); std::vector surrogatecoord_y_outer(sizeSurrogateLoop_outer); int countSurrogateLoop_outer = 0; - for (auto i_node = rSurrogateModelPart_outer.NodesEnd()-1; i_node != rSurrogateModelPart_outer.NodesBegin()-1; i_node--) { + for (auto i_node = rSurrogateModelPartOuter.NodesEnd()-1; i_node != rSurrogateModelPartOuter.NodesBegin()-1; i_node--) { surrogatecoord_x_outer[countSurrogateLoop_outer] = i_node->X(); surrogatecoord_y_outer[countSurrogateLoop_outer] = i_node->Y(); countSurrogateLoop_outer++; @@ -276,17 +276,17 @@ class CreateBrepsSBMUtilities : public IO } // INNER - for (IndexType iel = 1; iel < rSurrogateModelPart_inner.Elements().size()+1; iel++) { - int firstSurrogateNodeId = rSurrogateModelPart_inner.pGetElement(iel)->GetGeometry()[0].Id(); // Element 1 because is the only surrogate loop - int lastSurrogateNodeId = rSurrogateModelPart_inner.pGetElement(iel)->GetGeometry()[1].Id(); // Element 1 because is the only surrogate loop + for (IndexType iel = 1; iel < rSurrogateModelPartInner.Elements().size()+1; iel++) { + int firstSurrogateNodeId = rSurrogateModelPartInner.pGetElement(iel)->GetGeometry()[0].Id(); // Element 1 because is the only surrogate loop + int lastSurrogateNodeId = rSurrogateModelPartInner.pGetElement(iel)->GetGeometry()[1].Id(); // Element 1 because is the only surrogate loop int sizeSurrogateLoop = lastSurrogateNodeId - firstSurrogateNodeId + 1; int countSurrogateLoop = 0; std::vector surrogatecoord_x(sizeSurrogateLoop); std::vector surrogatecoord_y(sizeSurrogateLoop); for (int id_node = firstSurrogateNodeId; id_node < lastSurrogateNodeId+1; id_node++) { - surrogatecoord_x[countSurrogateLoop] = rSurrogateModelPart_inner.GetNode(id_node).X(); - surrogatecoord_y[countSurrogateLoop] = rSurrogateModelPart_inner.GetNode(id_node).Y(); + surrogatecoord_x[countSurrogateLoop] = rSurrogateModelPartInner.GetNode(id_node).X(); + surrogatecoord_y[countSurrogateLoop] = rSurrogateModelPartInner.GetNode(id_node).Y(); countSurrogateLoop++; } std::vector>::Pointer> trimming_curves_GPT; diff --git a/kratos/geometries/brep_surface.h b/kratos/geometries/brep_surface.h index ee48fe8e41a0..38106fa2aa3c 100644 --- a/kratos/geometries/brep_surface.h +++ b/kratos/geometries/brep_surface.h @@ -123,13 +123,13 @@ class BrepSurface typename NurbsSurfaceType::Pointer pSurface, BrepCurveOnSurfaceLoopArrayType& BrepOuterLoopArray, BrepCurveOnSurfaceLoopArrayType& BrepInnerLoopArray, - ModelPart& rSurrogateModelPart_inner, ModelPart& rSurrogateModelPart_outer) + ModelPart& rSurrogateModelPartInner, ModelPart& rSurrogateModelPartOuter) : BaseType(PointsArrayType(), &msGeometryData) , mpNurbsSurface(pSurface) , mOuterLoopArray(BrepOuterLoopArray) , mInnerLoopArray(BrepInnerLoopArray) - , mpSurrogateModelPart_inner(&rSurrogateModelPart_inner) - , mpSurrogateModelPart_outer(&rSurrogateModelPart_outer) + , mpSurrogateModelPartInner(&rSurrogateModelPartInner) + , mpSurrogateModelPartOuter(&rSurrogateModelPartOuter) { mIsTrimmed = false; } @@ -611,8 +611,8 @@ class BrepSurface BrepCurveOnSurfaceArrayType mEmbeddedEdgesArray; - ModelPart* mpSurrogateModelPart_inner = nullptr; - ModelPart* mpSurrogateModelPart_outer = nullptr; + ModelPart* mpSurrogateModelPartInner = nullptr; + ModelPart* mpSurrogateModelPartOuter = nullptr; /** IsTrimmed is used to optimize processes as * e.g. creation of integration domain. diff --git a/kratos/geometries/nurbs_shape_function_utilities/nurbs_utilities.h b/kratos/geometries/nurbs_shape_function_utilities/nurbs_utilities.h index 66ce6449c17d..0d6ba0a1cf2b 100644 --- a/kratos/geometries/nurbs_shape_function_utilities/nurbs_utilities.h +++ b/kratos/geometries/nurbs_shape_function_utilities/nurbs_utilities.h @@ -66,8 +66,17 @@ namespace NurbsUtilities const Vector& rKnots, const double ParameterT) { + // Check if the ParameterT (gauss point) is coincident to a knot (laying on an edge of a knot span) + double parameter_t_corrected = ParameterT; + for (unsigned i = 0; i> matrix; matrix.reserve(n_knot_spans_uv[1]); - for (int j = 0; j <= n_knot_spans_uv[1]; ++j) { + for (int j = 0; j <= n_knot_spans_uv[1]-1; ++j) { std::vector row(n_knot_spans_uv[0]); matrix.push_back(row); } @@ -219,12 +219,12 @@ namespace Kratos n_knot_spans_uv, knot_step_uv, starting_pos_uv); if (is_inner) { - CreateSurrogateBuondaryFromSnake_inner (knot_spans_available, idInnerLoop, surrogate_sub_model_part, + CreateSurrogateBuondaryFromSnakeInner (knot_spans_available, idInnerLoop, surrogate_sub_model_part, n_knot_spans_uv, knot_vector_u, knot_vector_v, starting_pos_uv ); KRATOS_INFO_IF("::[SnakeSBMUtilities]::", rEchoLevel > 0) << "Inner :: Snake process has finished" << std::endl; } else { - CreateSurrogateBuondaryFromSnake_outer (testBins, skin_sub_model_part, knot_spans_available, idInnerLoop, surrogate_sub_model_part, + CreateSurrogateBuondaryFromSnakeOuter (testBins, skin_sub_model_part, knot_spans_available, idInnerLoop, surrogate_sub_model_part, n_knot_spans_uv, knot_vector_u, knot_vector_v, starting_pos_uv); KRATOS_INFO_IF("::[SnakeSBMUtilities]::", rEchoLevel > 0) << "Outer :: Snake process has finished" << std::endl; } @@ -257,6 +257,9 @@ namespace Kratos int knot_span_u_point_split = (x_true_boundary_split-starting_pos_uv[0]) / knot_step_uv[0] ; int knot_span_v_point_split = (y_true_boundary_split-starting_pos_uv[1]) / knot_step_uv[1] ; + if (knot_span_u_point_split == knot_spans_available[idMatrix][0].size()) knot_span_u_point_split--; + if (knot_span_v_point_split == knot_spans_available[idMatrix].size()) knot_span_v_point_split--; + // update xy_coord for the first split segment std::vector> xy_coord_i_cond_split(2); xy_coord_i_cond_split[0].resize(2); xy_coord_i_cond_split[1].resize(2); @@ -296,7 +299,6 @@ namespace Kratos // Find the "knot_spans_available" using the intersection knot_spans_available[idMatrix][knot_spans_uv[1][0]][knot_spans_uv[0][0]] = 2; knot_spans_available[idMatrix][knot_spans_uv[1][0]][knot_spans_uv[0][1]] = 2; - } else if (knot_spans_uv[1][0] != knot_spans_uv[1][1]) { // v knot value is crossed // Find the "knot_spans_available" using the intersection (Snake_coordinate classic -> External Boundary) @@ -320,7 +322,7 @@ namespace Kratos } - bool SnakeSBMUtilities::isPointInsideSkinBoundary(Point& point1, DynamicBins& testBins, ModelPart& skin_model_part) + bool SnakeSBMUtilities::IsPointInsideSkinBoundary(Point& point1, DynamicBins& testBins, ModelPart& skin_model_part) { // Get the nearest point of the true boundary PointerType pointToSearch = PointerType(new PointType(1000000, point1.X(), point1.Y(), 0.0)); @@ -373,50 +375,50 @@ namespace Kratos if (i != n_knot_spans_uv[1]-1) if (knot_spans_available[idMatrix][i+1][j] == 0) { Point gaussPoint = Point((j+0.5) * knot_step_uv[0] + starting_pos_uv[0], (i+1+0.5) * knot_step_uv[1] +starting_pos_uv[1], 0); - if (isPointInsideSkinBoundary(gaussPoint, testBin, skin_model_part)) {knot_spans_available[idMatrix][i+1][j] = 1;} + if (IsPointInsideSkinBoundary(gaussPoint, testBin, skin_model_part)) {knot_spans_available[idMatrix][i+1][j] = 1;} } // left node if (i != 0) if (knot_spans_available[idMatrix][i-1][j] == 0) { Point gaussPoint = Point((j+0.5) * knot_step_uv[0]+starting_pos_uv[0], (i-1+0.5) * knot_step_uv[1] + starting_pos_uv[1], 0); - if (isPointInsideSkinBoundary(gaussPoint, testBin, skin_model_part)) {knot_spans_available[idMatrix][i-1][j] = 1;} + if (IsPointInsideSkinBoundary(gaussPoint, testBin, skin_model_part)) {knot_spans_available[idMatrix][i-1][j] = 1;} } // up node if (j != n_knot_spans_uv[0]-1) if (knot_spans_available[idMatrix][i][j+1] == 0) { Point gaussPoint = Point((j+1+0.5) * knot_step_uv[0]+starting_pos_uv[0], (i+0.5) * knot_step_uv[1]+starting_pos_uv[1], 0); - if (isPointInsideSkinBoundary(gaussPoint, testBin, skin_model_part)) {knot_spans_available[idMatrix][i][j+1] = 1;} + if (IsPointInsideSkinBoundary(gaussPoint, testBin, skin_model_part)) {knot_spans_available[idMatrix][i][j+1] = 1;} } //down node if (j != 0) if (knot_spans_available[idMatrix][i][j-1] == 0) { Point gaussPoint = Point((j-1+0.5) * knot_step_uv[0]+starting_pos_uv[0], (i+0.5) * knot_step_uv[1]+starting_pos_uv[1], 0); - if (isPointInsideSkinBoundary(gaussPoint, testBin, skin_model_part)) {knot_spans_available[idMatrix][i][j-1] = 1;} + if (IsPointInsideSkinBoundary(gaussPoint, testBin, skin_model_part)) {knot_spans_available[idMatrix][i][j-1] = 1;} } // corner right-down node if (j != 0 && i != n_knot_spans_uv[1]-1) if (knot_spans_available[idMatrix][i+1][j-1] == 0) { Point gaussPoint = Point((j-1+0.5) * knot_step_uv[0]+starting_pos_uv[0], (i+1+0.5) * knot_step_uv[1]+starting_pos_uv[1], 0); - if (isPointInsideSkinBoundary(gaussPoint, testBin, skin_model_part)) {knot_spans_available[idMatrix][i+1][j-1] = 1;} + if (IsPointInsideSkinBoundary(gaussPoint, testBin, skin_model_part)) {knot_spans_available[idMatrix][i+1][j-1] = 1;} } // corner left-down node if (j != 0 && i != 0) if (knot_spans_available[idMatrix][i-1][j-1] == 0) { Point gaussPoint = Point((j-1+0.5) * knot_step_uv[0]+starting_pos_uv[0], (i-1+0.5) * knot_step_uv[1]+starting_pos_uv[1], 0); - if (isPointInsideSkinBoundary(gaussPoint, testBin, skin_model_part)) {knot_spans_available[idMatrix][i-1][j-1] = 1;} + if (IsPointInsideSkinBoundary(gaussPoint, testBin, skin_model_part)) {knot_spans_available[idMatrix][i-1][j-1] = 1;} } // corner right-up node if (j != n_knot_spans_uv[0]-1 && i != n_knot_spans_uv[1]-1) if (knot_spans_available[idMatrix][i+1][j+1] == 0) { Point gaussPoint = Point((j+1+0.5) * knot_step_uv[0]+starting_pos_uv[0], (i+1+0.5) * knot_step_uv[1]+starting_pos_uv[1], 0); - if (isPointInsideSkinBoundary(gaussPoint, testBin, skin_model_part)) {knot_spans_available[idMatrix][i+1][j+1] = 1;} + if (IsPointInsideSkinBoundary(gaussPoint, testBin, skin_model_part)) {knot_spans_available[idMatrix][i+1][j+1] = 1;} } // corner left-up node if (j != n_knot_spans_uv[0]-1 && i != 0) if (knot_spans_available[idMatrix][i-1][j+1] == 0) { Point gaussPoint = Point((j+1+0.5) * knot_step_uv[0]+starting_pos_uv[0], (i-1+0.5) * knot_step_uv[1]+starting_pos_uv[1], 0); - if (isPointInsideSkinBoundary(gaussPoint, testBin, skin_model_part)) {knot_spans_available[idMatrix][i-1][j+1] = 1;} + if (IsPointInsideSkinBoundary(gaussPoint, testBin, skin_model_part)) {knot_spans_available[idMatrix][i-1][j+1] = 1;} } // Create 25 "fake" GaussPoints to check if the majority are inside or outside @@ -430,7 +432,7 @@ namespace Kratos { double y_coord = i*knot_step_uv[1] + knot_step_uv[1]/(numFakeGaussPoints+1)*(i_GPy+1) + starting_pos_uv[1]; Point gaussPoint = Point(x_coord, y_coord, 0); // GAUSSIAN POINT - if (isPointInsideSkinBoundary(gaussPoint, testBin, skin_model_part)) { + if (IsPointInsideSkinBoundary(gaussPoint, testBin, skin_model_part)) { // Sum over the number of numFakeGaussPoints per knot span numberOfInsideGaussianPoints++; } @@ -451,7 +453,7 @@ namespace Kratos } - void SnakeSBMUtilities::CreateSurrogateBuondaryFromSnake_inner (std::vector>> & knot_spans_available, int idMatrix, + void SnakeSBMUtilities::CreateSurrogateBuondaryFromSnakeInner (std::vector>> & knot_spans_available, int idMatrix, ModelPart& surrogate_model_part_inner, std::vector& n_knot_spans_uv, Vector& knot_vector_u, Vector& knot_vector_v, const Vector& starting_pos_uv){ @@ -485,7 +487,7 @@ namespace Kratos // Follow the clockwise loop bool end = false; - // We are going orizontally + // We are going horizontally int direction = 0; // 0 = up_vertical, 1 = right_orizontal, 2 = down_vertical, 3 = left_orizontal int i = start_i; int j = start_j; @@ -506,7 +508,7 @@ namespace Kratos if (direction == -1) {direction = 3;} } else { - // Need to try to go straight or right w.r.t. the current direction; risetto e muovo + // Need to try to go straight or right w.r.t. the current direction; reset and move if (direction == 0) {I++ ; J++ ;} if (direction == 1) {J-- ; I++ ;} if (direction == 2) {I-- ; J-- ;} @@ -524,7 +526,7 @@ namespace Kratos idSnakeNode++; } else { - // Need to search to hte right; rese e muovo + // Need to search to hte right; reset and move if (direction == 0) {J-- ; I++ ;} if (direction == 1) {I-- ; J-- ;} if (direction == 2) {J++ ; I-- ;} @@ -620,7 +622,7 @@ namespace Kratos } - void SnakeSBMUtilities::CreateSurrogateBuondaryFromSnake_outer (DynamicBins& testBin_out, ModelPart& initial_skin_model_part_out, + void SnakeSBMUtilities::CreateSurrogateBuondaryFromSnakeOuter (DynamicBins& testBin_out, ModelPart& initial_skin_model_part_out, std::vector>> & knot_spans_available, int idMatrix, ModelPart& surrogate_model_part_outer, std::vector& n_knot_spans_uv, Vector& knot_vector_u, Vector& knot_vector_v, @@ -633,36 +635,36 @@ namespace Kratos double knot_step_v = knot_vector_v[1]-knot_vector_v[0]; for (int i = 0; i<2; i++) { - for (int j = 0; j < (n_knot_spans_uv[1]); j++ ) { + for (int j = 0; j < (n_knot_spans_uv[0]); j++ ) { Point centroidKnotSpan = Point((j+0.5)*knot_step_u+starting_pos_uv[0], (i+0.5)*knot_step_v+starting_pos_uv[1], 0); - if (isPointInsideSkinBoundary(centroidKnotSpan, testBin_out, initial_skin_model_part_out) && knot_spans_available[idMatrix][i][j] != -1) { + if (IsPointInsideSkinBoundary(centroidKnotSpan, testBin_out, initial_skin_model_part_out) && knot_spans_available[idMatrix][i][j] != -1) { knot_spans_available[idMatrix][i][j] = 1; } } } // TOP BOUNDARY - for (int j = int(knot_spans_available[idMatrix][0].size())-1; j > int(knot_spans_available[idMatrix][0].size())-3; j--) { - for (int i = 0; i < (n_knot_spans_uv[0]); i++) { + for (int j = knot_spans_available[idMatrix][0].size()-1; j > knot_spans_available[idMatrix][0].size()-3; j--) { + for (int i = 0; i < (n_knot_spans_uv[1]); i++) { Point centroidKnotSpan = Point((j+0.5)*knot_step_u+starting_pos_uv[0], (i+0.5)*knot_step_v+starting_pos_uv[1], 0); - if (isPointInsideSkinBoundary(centroidKnotSpan, testBin_out, initial_skin_model_part_out) && knot_spans_available[idMatrix][i][j] != -1) { + if (IsPointInsideSkinBoundary(centroidKnotSpan, testBin_out, initial_skin_model_part_out) && knot_spans_available[idMatrix][i][j] != -1) { knot_spans_available[idMatrix][i][j] = 1; } } } // RIGHT BOUNDARY - for (int i = int(knot_spans_available[idMatrix][0].size())-1; i > int(knot_spans_available[idMatrix][0].size())-3; i--) { - for (int j = n_knot_spans_uv[1]-1; j > -1; j-- ) { + for (int i = knot_spans_available[idMatrix].size()-1; i > knot_spans_available[idMatrix].size()-3; i--) { + for (int j = n_knot_spans_uv[0]-1; j > -1; j-- ) { Point centroidKnotSpan = Point((j+0.5)*knot_step_u+starting_pos_uv[0], (i+0.5)*knot_step_v+starting_pos_uv[1], 0); - if (isPointInsideSkinBoundary(centroidKnotSpan, testBin_out, initial_skin_model_part_out) && knot_spans_available[idMatrix][i][j] != -1) { + if (IsPointInsideSkinBoundary(centroidKnotSpan, testBin_out, initial_skin_model_part_out) && knot_spans_available[idMatrix][i][j] != -1) { knot_spans_available[idMatrix][i][j] = 1; } } } // BOTTOM BOUNDARY for (int j = 0; j<2; j++) { - for (int i = n_knot_spans_uv[0]-1; i > -1 ; i--) { + for (int i = n_knot_spans_uv[1]-1; i > -1 ; i--) { Point centroidKnotSpan = Point((j+0.5)*knot_step_u+starting_pos_uv[0], (i+0.5)*knot_step_v+starting_pos_uv[1], 0); - if (isPointInsideSkinBoundary(centroidKnotSpan, testBin_out, initial_skin_model_part_out) && knot_spans_available[idMatrix][i][j] != -1) { + if (IsPointInsideSkinBoundary(centroidKnotSpan, testBin_out, initial_skin_model_part_out) && knot_spans_available[idMatrix][i][j] != -1) { knot_spans_available[idMatrix][i][j] = 1; } } diff --git a/kratos/utilities/nurbs_utilities/snake_sbm_utilities.h b/kratos/utilities/nurbs_utilities/snake_sbm_utilities.h index 008c643a76c6..1763643982e7 100644 --- a/kratos/utilities/nurbs_utilities/snake_sbm_utilities.h +++ b/kratos/utilities/nurbs_utilities/snake_sbm_utilities.h @@ -108,7 +108,7 @@ namespace Kratos * @param skin_model_part_in Input skin model part. * @return True if the point is inside the boundary, false otherwise. */ - static bool isPointInsideSkinBoundary(Point& pointToSearch, + static bool IsPointInsideSkinBoundary(Point& pointToSearch, DynamicBins& testBins, ModelPart& skin_model_part_in); @@ -141,7 +141,7 @@ namespace Kratos * @param knot_vector_v Knot vector in the V-direction. * @param starting_pos_uv Starting position in UV space. */ - static void CreateSurrogateBuondaryFromSnake_inner (std::vector>> & knot_spans_available, + static void CreateSurrogateBuondaryFromSnakeInner (std::vector>> & knot_spans_available, int idMatrix, ModelPart& surrogate_model_part_inner, std::vector& n_knot_spans_uv, @@ -161,7 +161,7 @@ namespace Kratos * @param knot_vector_v Knot vector in the V-direction. * @param starting_pos_uv Starting position in UV space. */ - static void CreateSurrogateBuondaryFromSnake_outer (DynamicBins& testBin_out, + static void CreateSurrogateBuondaryFromSnakeOuter (DynamicBins& testBin_out, ModelPart& initial_skin_model_part_out, std::vector>> & knot_spans_available, int idMatrix, From cfb7ba453cb8daeffde47b048af21eb3a07d3428 Mon Sep 17 00:00:00 2001 From: andrewgorgi Date: Wed, 29 Jan 2025 18:05:23 +0100 Subject: [PATCH 7/8] IndexType error solved --- kratos/utilities/nurbs_utilities/snake_sbm_utilities.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/kratos/utilities/nurbs_utilities/snake_sbm_utilities.cpp b/kratos/utilities/nurbs_utilities/snake_sbm_utilities.cpp index d040e3166d13..5a26c4460dcd 100644 --- a/kratos/utilities/nurbs_utilities/snake_sbm_utilities.cpp +++ b/kratos/utilities/nurbs_utilities/snake_sbm_utilities.cpp @@ -643,7 +643,7 @@ namespace Kratos } } // TOP BOUNDARY - for (int j = knot_spans_available[idMatrix][0].size()-1; j > knot_spans_available[idMatrix][0].size()-3; j--) { + for (int j = int (knot_spans_available[idMatrix][0].size()-1); j > int (knot_spans_available[idMatrix][0].size()-3); j--) { for (int i = 0; i < (n_knot_spans_uv[1]); i++) { Point centroidKnotSpan = Point((j+0.5)*knot_step_u+starting_pos_uv[0], (i+0.5)*knot_step_v+starting_pos_uv[1], 0); if (IsPointInsideSkinBoundary(centroidKnotSpan, testBin_out, initial_skin_model_part_out) && knot_spans_available[idMatrix][i][j] != -1) { @@ -652,7 +652,7 @@ namespace Kratos } } // RIGHT BOUNDARY - for (int i = knot_spans_available[idMatrix].size()-1; i > knot_spans_available[idMatrix].size()-3; i--) { + for (int i = int (knot_spans_available[idMatrix].size()-1); i > int (knot_spans_available[idMatrix].size()-3); i--) { for (int j = n_knot_spans_uv[0]-1; j > -1; j-- ) { Point centroidKnotSpan = Point((j+0.5)*knot_step_u+starting_pos_uv[0], (i+0.5)*knot_step_v+starting_pos_uv[1], 0); if (IsPointInsideSkinBoundary(centroidKnotSpan, testBin_out, initial_skin_model_part_out) && knot_spans_available[idMatrix][i][j] != -1) { From 2abf0c3ad91cd1397b0c79cc3b421bc12088b402 Mon Sep 17 00:00:00 2001 From: andrewgorgi Date: Wed, 29 Jan 2025 22:27:05 +0100 Subject: [PATCH 8/8] Another int/IndexType correction --- kratos/utilities/nurbs_utilities/snake_sbm_utilities.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/kratos/utilities/nurbs_utilities/snake_sbm_utilities.cpp b/kratos/utilities/nurbs_utilities/snake_sbm_utilities.cpp index 5a26c4460dcd..6a3002c620c6 100644 --- a/kratos/utilities/nurbs_utilities/snake_sbm_utilities.cpp +++ b/kratos/utilities/nurbs_utilities/snake_sbm_utilities.cpp @@ -257,8 +257,8 @@ namespace Kratos int knot_span_u_point_split = (x_true_boundary_split-starting_pos_uv[0]) / knot_step_uv[0] ; int knot_span_v_point_split = (y_true_boundary_split-starting_pos_uv[1]) / knot_step_uv[1] ; - if (knot_span_u_point_split == knot_spans_available[idMatrix][0].size()) knot_span_u_point_split--; - if (knot_span_v_point_split == knot_spans_available[idMatrix].size()) knot_span_v_point_split--; + if (knot_span_u_point_split == int (knot_spans_available[idMatrix][0].size())) knot_span_u_point_split--; + if (knot_span_v_point_split == int (knot_spans_available[idMatrix].size())) knot_span_v_point_split--; // update xy_coord for the first split segment std::vector> xy_coord_i_cond_split(2);