[IgaApplication] Add Nurbs modeler for SBM #13009

NickNick9 · 2025-01-14T10:43:45Z

📝 Description
We have added the NURBS Modeler SBM:

It is a derived class from the existing NURBS modeler, which creates a rectangular domain.
The first addition is the identification of surrogate boundaries for the inner and outer loops.
The second addition is the creation of the surface and BREP curves from the surrogate boundaries.

In particular, we have developed two utilities: one to identify the surrogate boundaries and another to create the BREP entities. These utilities were created to enhance clarity and because they might be used separately in other applications.

Please mark the PR with appropriate tags:

IgaApplication

…rSbm Co-authored-by: Nicolò Antonelli <[email protected]>

AlejandroCornejo · 2025-01-15T11:15:04Z

applications/IgaApplication/custom_modelers/nurbs_geometry_modeler_sbm.cpp

+        if (mParameters.Has("skin_model_part_inner_initial_name")) {
+            skin_model_part_inner_initial_name = mParameters["skin_model_part_inner_initial_name"].GetString();
+
+            if (!mpModel->HasModelPart(skin_model_part_inner_initial_name)) 


you can use KRATOS_ERROR_IF()

AlejandroCornejo · 2025-01-15T11:15:25Z

applications/IgaApplication/custom_modelers/nurbs_geometry_modeler_sbm.cpp

+                                    : 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";


Wouls it be useful to set this name from the json file?

This is actually the "default" name which is overwritten few lines below if the "skin_model_part_inner_inner_initial_name" is provided in the json file.

Json defaults should automatically handle this.

Can you suggest us an example to start form for the json defaults? @rubenzorrilla

AlejandroCornejo · 2025-01-15T11:16:41Z

applications/IgaApplication/custom_modelers/nurbs_geometry_modeler_sbm.h

+//
+
+#if !defined(KRATOS_NURBS_GEOMETRY_MODELER_SBM_H_INCLUDED )
+#define  KRATOS_NURBS_GEOMETRY_MODELER_SBM_H_INCLUDED


pragma once

AlejandroCornejo · 2025-01-15T11:17:10Z

applications/IgaApplication/custom_modelers/nurbs_geometry_modeler_sbm.h

+#include "includes/model_part.h"
+#include "nurbs_geometry_modeler.h"
+#include "geometries/nurbs_volume_geometry.h"
+#include "geometries/nurbs_surface_geometry.h"


can those includes go to the .cpp?

AlejandroCornejo · 2025-01-15T11:17:32Z

applications/IgaApplication/custom_modelers/nurbs_geometry_modeler_sbm.h

+    typedef Node NodeType;
+
+    typedef Geometry<NodeType> GeometryType;
+    typedef typename GeometryType::Pointer GeometryPointerType;


now we use using GeometryType = Geometry<NodeType>;

AlejandroCornejo · 2025-01-15T11:23:13Z

kratos/utilities/nurbs_utilities/snake_sbm_utilities.cpp

+
+                xy_coord_i_cond[0][0] = i_cond.GetGeometry()[0].X(); // x_true_boundary1
+                xy_coord_i_cond[1][0] = i_cond.GetGeometry()[0].Y(); // y_true_boundary1
+                xy_coord_i_cond[0][1] = i_cond.GetGeometry()[1].X(); // x_true_boundary2


AlejandroCornejo · 2025-01-15T11:24:03Z

kratos/utilities/nurbs_utilities/snake_sbm_utilities.cpp

+                        knot_spans_available[idMatrix][i][j] = 1; // The knot span is considered DEACTIVE
+                    }
+
+                    // exit(0);


AlejandroCornejo · 2025-01-15T11:24:16Z

kratos/utilities/nurbs_utilities/snake_sbm_utilities.cpp

+        IndexType idSnakeNode = id_first_node+1;
+
+        // Follow the clockwise loop
+        int end = 0;


AlejandroCornejo · 2025-01-15T11:24:44Z

kratos/utilities/nurbs_utilities/snake_sbm_utilities.cpp

+        // 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; 


AlejandroCornejo · 2025-01-15T11:24:55Z

kratos/utilities/nurbs_utilities/snake_sbm_utilities.h

+//
+
+#if !defined(KRATOS_SNAKE_SBM_UTILITIES_H_INCLUDED )
+#define KRATOS_SNAKE_SBM_UTILITIES_H_INCLUDED


pragma once

Co-authored-by: Nicolò Antonelli <[email protected]>

rickyaristio

Sorry for my late review. I have minor comments related to the writing convention style.

rickyaristio · 2025-01-28T14:34:14Z

kratos/utilities/nurbs_utilities/snake_sbm_utilities.cpp

+
+        // Follow the clockwise loop
+        bool end = false;
+        // We are going orizontally


rickyaristio · 2025-01-28T14:35:38Z

kratos/utilities/nurbs_utilities/snake_sbm_utilities.cpp

+                if (direction == -1) {direction = 3;}
+            }
+            else {
+                // Need to try to go straight or right w.r.t. the current direction; risetto e muovo


check comment :)

rickyaristio · 2025-01-28T14:41:52Z

kratos/utilities/nurbs_utilities/snake_sbm_utilities.h

+         * @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, 


IsPointInsideSkinBoundary
(naming convention)

rickyaristio · 2025-01-28T14:42:27Z

kratos/utilities/nurbs_utilities/snake_sbm_utilities.h

+         * @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<std::vector<std::vector<int>>> & knot_spans_available, 


CreateSurrogateBuondaryFromSnakeInner

🐍❓🤣

rickyaristio · 2025-01-28T14:42:38Z

kratos/utilities/nurbs_utilities/snake_sbm_utilities.h

+         * @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, 


CreateSurrogateBuondaryFromSnakeOuter

rickyaristio · 2025-01-28T14:43:16Z

kratos/utilities/nurbs_utilities/snake_sbm_utilities.cpp

+                    idSnakeNode++;
+                }
+                else {
+                    // Need to search to hte right; rese e muovo


check comment :)

rickyaristio · 2025-01-28T14:44:11Z

kratos/geometries/brep_surface.h

+    ModelPart* mpSurrogateModelPart_inner = nullptr;
+    ModelPart* mpSurrogateModelPart_outer = nullptr;


mpSurrogateModelPartInner
mpSurrogateModelPartOuter

rickyaristio · 2025-01-28T14:46:13Z

applications/IgaApplication/iga_application.h

@@ -141,6 +142,7 @@ class KRATOS_API(IGA_APPLICATION) KratosIgaApplication : public KratosApplicatio
    const IgaModeler mIgaModeler;
    const RefinementModeler mRefinementModeler;
    const NurbsGeometryModeler mNurbsGeometryModeler;
+    const NurbsGeometryModelerSbm mNurbsGeometryModelerSbm;


mNurbsGeometryModelerSBM (make it consistent with mSBMLaplacianCondition)

rickyaristio · 2025-01-28T14:46:33Z

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);


rickyaristio · 2025-01-28T14:59:36Z

Also, I have a few remarks:
-In this workflow, where do you assign the physics (elements or conditions, such as the sbm_laplacian_condition)?
-It would be great if you can provide one full test example (like this)

NickNick9 · 2025-01-28T16:37:19Z

Thanks @rickyaristio for your comments!
We still use the iga_modeler with many modifications, the plan is to create an iga_sbm_modeler. Next "big" PR will be on that and then we can provide a full test of the sbm workflow.

rickyaristio · 2025-01-29T07:10:15Z

Thanks @rickyaristio for your comments! We still use the iga_modeler with many modifications, the plan is to create an iga_sbm_modeler. Next "big" PR will be on that and then we can provide a full test of the sbm workflow.

I see. Looking forward to the next PR about iga_sbm_modeler.

rubenzorrilla · 2025-01-29T10:19:30Z

I'll review this after @rickyaristio's and @AlejandroCornejo's suggestions are applied.

matekelemen · 2025-02-03T20:31:43Z

applications/IgaApplication/custom_modelers/nurbs_geometry_modeler.cpp

+                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();
+                    }
                }


Suggested change

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();

}

}

const SizeType last_geometry_id = r_model_part.NumberOfGeometries()

? r_model_part.Geometries().Geometries().back().Id()

: 0;

I know, it looks stupid. In fact, it is stupid, but only because the geometry container is still weird.

matekelemen · 2025-02-03T20:38:32Z

applications/IgaApplication/custom_modelers/nurbs_geometry_modeler.h

+    ///@name Private Member Variables
+    ///@{


Just an FYI on doxygen annotations:

no doxygen docs are generated for private members, so there's no point marking these sections. Of course commenting the purpose of members is still good practice but it will only be visible in the code, and not the docs.

doxygen automatically detects and highlights the visibility (public/protected/private) of members so again: there's no point in doing section annontations.

matekelemen · 2025-02-03T20:44:14Z

applications/IgaApplication/custom_modelers/nurbs_geometry_modeler_sbm.cpp

+    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)


I assume this is an override of a virtual member. If so, mark it as such. Also, you can use the @copydoc doxygen command to pull in the annotations from the base class.

Suggested change

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)

/// @copydoc NurbsGeometryModeler::CreateAndAddRegularGrid2D

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) override

+ same comment about variable naming in general.

matekelemen · 2025-02-03T20:53:54Z

applications/IgaApplication/custom_modelers/nurbs_geometry_modeler.h

+    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);
+


That's a lot of arguments, and the docstring seems to be wildly out of date. Please update it and write a @param entry for each argument with some explanations.

Oh and please pick a naming format and stick with it. I suggest the one in the kratos style guide if you don't want the technical committee get mad at you. It would also help you avoid @loumalouomega reformatting your code.

matekelemen · 2025-02-03T21:02:04Z

applications/IgaApplication/custom_modelers/nurbs_geometry_modeler_sbm.h

+    ///@name Private Member Variables
+    ///@{
+
+    Model* mpModel;


Wait, doesn't the base class already have a pointer to the Model?

You could just expose that through a protected member of the base class to avoid duplication/ambiguity.

matekelemen · 2025-02-03T21:27:27Z

kratos/geometries/nurbs_shape_function_utilities/nurbs_utilities.h

        const auto span = std::lower_bound(std::begin(rKnots) + PolynomialDegree,
-            std::end(rKnots) - PolynomialDegree, ParameterT) - std::begin(rKnots) - 1;
+            std::end(rKnots) - PolynomialDegree, parameter_t_corrected) - std::begin(rKnots) - 1;


Can parameter_t_corrected <= rKnots.front() ever happen? If so, that'd be scary.

As a general comment for this function: you don't really need that first loop. As I understand, it's just handling cases that are close to knots. You could just use the output of std::lower_bound and check its neighbors.

e.g.:

auto it = std::lower_bound(...); if (it + 1 != rKnots.end()) { it = std::abs(ParameterT - it[1]) < tolerance ? it + 1 : it; } return std::distance(rKnots.begin(), it - 1);

Also keep in mind that you might have to clamp your output to 0 if you don't handle potentially negative spans later on.

matekelemen · 2025-02-03T21:29:49Z

kratos/utilities/nurbs_utilities/snake_sbm_utilities.h

+                                                DynamicBins& testBins,
+                                                ModelPart& skin_model_part_in);


my fav ❤️:

DynamicBins& testBins, ModelPart& skinModelPartIn);

ancient C dev (or annoying 🦀 crusader):

DynamicBins& test_bins, ModelPart& skin_model_part_in);

kratos style:

DynamicBins& TestBins, ModelPart& SkinModelPartIn);

matekelemen · 2025-02-03T21:31:39Z

kratos/utilities/nurbs_utilities/snake_sbm_utilities.h

+         * @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.


I appreciate the docs 👍.

You need to explain what you expect as parameters here though.

matekelemen · 2025-02-03T21:40:04Z

kratos/utilities/nurbs_utilities/snake_sbm_utilities.cpp

+        for (IndexType i = 0; i < numberOfLoops; ++i) {
+            std::vector<std::vector<int>> matrix; 
+            matrix.reserve(n_knot_spans_uv[1]);
+            for (int j = 0; j <= n_knot_spans_uv[1]-1; ++j) {
+                std::vector<int> row(n_knot_spans_uv[0]); 
+                matrix.push_back(row); 
+            }
+            knot_spans_available.push_back(matrix);
+        }


std::vector<std::vector<int>> matrix; ouch

How big can numberOfLoops get? If it's on the order of magnitude as elements/nodes or larger, this is going to be very pricey.

I see vector nesting is a running theme in this utility but I'd strongly advise against it. Consider using a dense matrix instead.

I know transitioning to a more efficient data structure would be a big change so I'm not suggesting that you do it in this PR. Just keep it in mind and do it when you have the time.

matekelemen · 2025-02-03T21:45:01Z

kratos/utilities/nurbs_utilities/snake_sbm_utilities.cpp

+        }
+        // Create "fictituos element" to memorize starting and ending node id for each surrogate boundary loop
+        std::vector<ModelPart::IndexType> elem_nodes{1, idSnakeNode-1};
+        IndexType elem_id = surrogate_model_part_outer.GetRootModelPart().Elements().size()+1;


Watch out: index != identifier.

Suggested change

IndexType elem_id = surrogate_model_part_outer.GetRootModelPart().Elements().size()+1;

const ModelPart& r_root = surrogate_model_part_outer.GetRootModelPart();

IndexType elem_id = r_root.Elements().empty()

? 1

: r_root.Elements().back().Id() + 1;

Ready for PR

92fd157

NickNick9 assigned andrewgorgi Jan 14, 2025

NickNick9 requested review from a team as code owners January 14, 2025 10:43

NickNick9 changed the title ~~Add Nurbs modeler for SBM for IgaApplication~~ [IgaApplication] Add Nurbs modeler for SBM Jan 14, 2025

andrewgorgi and others added 2 commits January 14, 2025 16:00

Modified int->IndexType + Body fitted case in the NurbsGeometryModele…

859abbd

…rSbm Co-authored-by: Nicolò Antonelli <[email protected]>

Correct the typo in test_IgaApplication.py

7a518dd

rickyaristio self-assigned this Jan 15, 2025

AlejandroCornejo reviewed Jan 15, 2025

View reviewed changes

andrewgorgi and others added 2 commits January 16, 2025 18:05

Suggestions by Alejandro

9f49410

Co-authored-by: Nicolò Antonelli <[email protected]>

direction variable converted to integer (not SizeType)

625d3af

rickyaristio reviewed Jan 28, 2025

View reviewed changes

andrewgorgi added 3 commits January 29, 2025 16:05

Correction to SnakeSbmUtilities + PR suggestions by Ricky

05c0290

IndexType error solved

cfb7ba4

Another int/IndexType correction

2abf0c3

matekelemen reviewed Feb 3, 2025

View reviewed changes

		ModelPart* mpSurrogateModelPart_inner = nullptr;
		ModelPart* mpSurrogateModelPart_outer = nullptr;

		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)

		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);

-        IndexType elem_id = surrogate_model_part_outer.GetRootModelPart().Elements().size()+1;
+        const ModelPart& r_root = surrogate_model_part_outer.GetRootModelPart();
+        IndexType elem_id = r_root.Elements().empty()
+                            ? 1
+                            : r_root.Elements().back().Id() + 1;

[IgaApplication] Add Nurbs modeler for SBM #13009

Are you sure you want to change the base?

[IgaApplication] Add Nurbs modeler for SBM #13009

Conversation

NickNick9 commented Jan 14, 2025

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

rickyaristio left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

rickyaristio commented Jan 28, 2025

NickNick9 commented Jan 28, 2025

rickyaristio commented Jan 29, 2025

rubenzorrilla commented Jan 29, 2025

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment