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..ef2738013065 --- /dev/null +++ b/applications/IgaApplication/custom_modelers/nurbs_geometry_modeler_sbm.cpp @@ -0,0 +1,126 @@ +// | / | +// ' / __| _` | __| _ \ __| +// . \ | ( | | ( |\__ ` +// _|\_\_| \__,_|\__|\___/ ____/ +// Multi-Physics +// +// License: BSD License +// Kratos default license: kratos/license.txt +// +// Main authors: Nicolo' Antonelli +// Andrea Gorgi +// + +// Project includes +#include "nurbs_geometry_modeler_sbm.h" +#include "custom_utilities/create_breps_sbm_utilities.h" +#include "utilities/nurbs_utilities/snake_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(); + + 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(); + + 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 + if (!(mParameters.Has("skin_model_part_inner_initial_name") || mParameters.Has("skin_model_part_outer_initial_name"))){ + + // 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")) + 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); + skin_model_part.CreateSubModelPart("inner"); + 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..3d3e2a1a643d --- /dev/null +++ b/applications/IgaApplication/custom_modelers/nurbs_geometry_modeler_sbm.h @@ -0,0 +1,121 @@ +// | / | +// ' / __| _` | __| _ \ __| +// . \ | ( | | ( |\__ ` +// _|\_\_| \__,_|\__|\___/ ____/ +// Multi-Physics +// +// License: BSD License +// Kratos default license: kratos/license.txt +// +// Main authors: Nicolo' Antonelli +// Andrea Gorgi +// + +# pragma once + +// System includes + +// External includes + +// Project includes +#include "nurbs_geometry_modeler.h" + +namespace Kratos { + +class KRATOS_API(IGA_APPLICATION) NurbsGeometryModelerSbm + : public NurbsGeometryModeler +{ +public: + ///@name Type Definitions + ///@{ + KRATOS_CLASS_POINTER_DEFINITION( NurbsGeometryModelerSbm ); + + using IndexType = std::size_t; + using SizeType = std::size_t; + using NodeType = Node; + + using GeometryType = Geometry; + using GeometryPointerType = GeometryType::Pointer; + + using NurbsSurfaceGeometryType = NurbsSurfaceGeometry<3, PointerVector>; + using NurbsSurfaceGeometryPointerType = NurbsSurfaceGeometryType::Pointer; + + using NurbsVolumeGeometryType = NurbsVolumeGeometry>; + using NurbsVolumeGeometryPointerType = NurbsVolumeGeometryType::Pointer; + + using ContainerNodeType = PointerVector; + using ContainerEmbeddedNodeType = PointerVector; + + ///@} + ///@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 \ 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..5c41ab4c6128 --- /dev/null +++ b/applications/IgaApplication/custom_utilities/create_breps_sbm_utilities.h @@ -0,0 +1,443 @@ +// | / | +// ' / __| _` | __| _ \ __| +// . \ | ( | | ( |\__ ` +// _|\_\_| \__,_|\__|\___/ ____/ +// Multi-Physics +// +// License: BSD License +// Kratos default license: kratos/license.txt +// +// Main authors: Nicolo' Antonelli +// Andrea Gorgi +// + +# pragma once + +// System includes + +// External includes + +// Project includes +#include "includes/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); + + using SizeType = std::size_t; + using IndexType = std::size_t; + + using GeometryType = Geometry; + using GeometryPointerType = typename GeometryType::Pointer; + + using ContainerNodeType = PointerVector; + using ContainerEmbeddedNodeType = PointerVector; + + using CouplingGeometryType = CouplingGeometry; + + using NurbsSurfaceType = NurbsSurfaceGeometry<3, ContainerNodeType>; + using NurbsTrimmingCurveType = NurbsCurveGeometry<2, ContainerEmbeddedNodeType>; + + using NurbsSurfacePointerType = typename NurbsSurfaceType::Pointer; + using NurbsTrimmingCurvePointerType = typename NurbsTrimmingCurveType::Pointer; + + using NurbsSurfaceGeometryType = NurbsSurfaceGeometry<3, PointerVector>; + using NurbsSurfaceGeometryPointerType = typename NurbsSurfaceGeometryType::Pointer; + + using BrepSurfaceType = BrepSurface; + using BrepCurveOnSurfaceType = BrepCurveOnSurface; + + using BrepCurveType = BrepCurve; + using PointOnGeometryOnSurfaceType = PointOnGeometry; + using PointOnGeometryOnCurveType = PointOnGeometry; + + using BrepCurveOnSurfaceArrayType = DenseVector; + using BrepCurveOnSurfaceLoopType = DenseVector; + using BrepCurveOnSurfaceLoopArrayType = DenseVector>; + + ///@} + ///@name Life Cycle + ///@{ + + /// Constructor with path to input file. + CreateBrepsSBMUtilities( + SizeType EchoLevel = 0) + : mEchoLevel(EchoLevel) + { + } + + /// Destructor. + ~CreateBrepsSBMUtilities() = default; + + ///@} + ///@name Python exposed Functions + ///@{ + + /// 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: + + /** + * @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, + ModelPart& rSurrogateModelPart_inner, + ModelPart& rSurrogateModelPart_outer, + 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, + rSurrogateModelPart_inner, + rSurrogateModelPart_outer); + + // 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); + } + + + 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 (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; + + 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. 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..18e8969600c7 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 TestModelersSbm 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..a72a1d429a89 --- /dev/null +++ b/applications/IgaApplication/tests/test_modelers_sbm.py @@ -0,0 +1,342 @@ +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 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() + 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..d77b2368a967 100644 --- a/kratos/utilities/geometry_utilities/brep_trimming_utilities.cpp +++ b/kratos/utilities/geometry_utilities/brep_trimming_utilities.cpp @@ -14,8 +14,8 @@ namespace Kratos { ///@name Kratos Classes ///@{ - - void BrepTrimmingUtilities::CreateBrepSurfaceTrimmingIntegrationPoints( + template + void BrepTrimmingUtilities::CreateBrepSurfaceTrimmingIntegrationPoints( IntegrationPointsArrayType& rIntegrationPoints, const DenseVector>& rOuterLoops, const DenseVector>& rInnerLoops, @@ -152,4 +152,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..dcf4fc6086a4 100644 --- a/kratos/utilities/geometry_utilities/brep_trimming_utilities.h +++ b/kratos/utilities/geometry_utilities/brep_trimming_utilities.h @@ -7,8 +7,7 @@ // License: BSD License // Kratos default license: kratos/license.txt -#if !defined(KRATOS_BREP_TRIMMING_UTILITIES_H_INCLUDED) -#define KRATOS_BREP_TRIMMING_UTILITIES_H_INCLUDED +#pragma once // Std includes #include @@ -32,7 +31,7 @@ namespace Kratos ///@{ //using namespace ClipperLib; - + template class KRATOS_API(KRATOS_CORE) BrepTrimmingUtilities { public: @@ -52,7 +51,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( @@ -326,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 new file mode 100644 index 000000000000..c774bcf3c1b1 --- /dev/null +++ b/kratos/utilities/nurbs_utilities/snake_sbm_utilities.cpp @@ -0,0 +1,834 @@ +// | / | +// ' / __| _` | __| _ \ __| +// . \ | ( | | ( |\__ ` +// _|\_\_| \__,_|\__|\___/ ____/ +// 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); + + 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)] ) ; + + 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); + + array_1d starting_pos_uv; + 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 (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) { + 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; + IndexType 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 + 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] = 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] = (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 && + (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 (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); + + 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 + 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 + IndexType 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, + 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) { + // 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 (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 (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 + 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 + } + } + } + } + } + + + 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 + bool end = false; + // 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 && 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 = true; + 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 = 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) { + 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-- ) { + 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 (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]; + } + } + + 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 + bool end = false; + // 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 && 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 = true; + } + } + // Create "fictituos element" to memorize starting and ending node id for each surrogate boundary loop + std::vector elem_nodes{1, idSnakeNode-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. 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..008c643a76c6 --- /dev/null +++ b/kratos/utilities/nurbs_utilities/snake_sbm_utilities.h @@ -0,0 +1,176 @@ +// | / | +// ' / __| _` | __| _ \ __| +// . \ | ( | | ( |\__ ` +// _|\_\_| \__,_|\__|\___/ ____/ +// Multi-Physics +// +// License: BSD License +// Kratos default license: kratos/license.txt +// +// Main authors: Nicolo' Antonelli +// Andrea Gorgi +// + +#pragma once + +// 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, + array_1d& 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.