[IgaApplication] example with SBM in 2d with convergence study importing a nurbs

iga/README.md

@@ -0,0 +1,7 @@
+# Iga Examples
+
+## Use Cases
+- External Boundary Circle (NURBS) 
+
+## Validation Cases
+- External Boundary Circle (NURBS) 

iga/use_cases/README.md

@@ -0,0 +1,56 @@
+# Use Cases IGA – External Boundary Circle (NURBS) — README.txt
+
+Short description
+This example runs a 2D IGA structural model whose outer boundary is a NURBS circle. A manufactured solution is enforced to validate the setup.
+
+Files
+- ProjectParameters.json
+- materials.json
+- nurbs_files/circle_nurbs.json
+- run_and_post.py (optional)
+
+Geometry / modeling (very brief)
+- ImportNurbsSbmModeler loads the circle NURBS into “initial_skin_model_part_out”.
+- NurbsGeometryModelerSbm builds the analysis surface (order [2,2], knot spans [20,20]) and creates “skin_model_part”.
+- IgaModelerSbm creates:
+  • SolidElement on IgaModelPart.StructuralAnalysisDomain
+  • SbmSolidCondition on SurfaceEdge (brep_ids: [2]) → IgaModelPart.SBM_Support_outer
+
+Manufactured BCs and loads
+- Displacement on outer boundary:
+  u = [ -cos(x)*sinh(y),  sin(x)*cosh(y), 0.0 ]
+- Body force in the domain (consistent with the above):
+  BODY_FORCE = [
+    -1000*(1+0.3)/(1-0.3*0.3)*cos(x)*sinh(y),
+    -1000*(1+0.3)/(1-0.3*0.3)*sin(x)*cosh(y),
+    0.0
+  ]
+
+Solver (key settings)
+- Static, nonlinear; OpenMP
+- Single time step: dt = 0.1, end_time = 0.1
+- Linear solver: LinearSolversApplication.sparse_lu
+
+Minimal materials.json (example)
+Use plane strain (or plane stress) linear elastic; keep it consistent with ν=0.3:
+{
+  "properties": [{
+    "model_part_name": "IgaModelPart.StructuralAnalysisDomain",
+    "properties_id": 1,
+    "Material": {
+      "constitutive_law": { "name": "LinearElasticPlaneStrain2DLaw" },
+      "Variables": {
+        "YOUNG_MODULUS": 1000.0,
+        "POISSON_RATIO": 0.3,
+        "DENSITY": 1.0
+      }
+    }
+  }]
+}
+
+Run:
+  python3 run_and_post.py        # to reproduce the error
+  python3 convergence.py         # for the convergence analysis
+  python3 plot_surrogate_boundaries.py  # to visualize the surrogate and true boundaries
+
+

iga/use_cases/external_boundary_circle_with_nurbs/ProjectParameters.json

@@ -0,0 +1,152 @@
+{
+  "problem_data": {
+    "problem_name": "example_1_kratos",
+    "echo_level": 0,
+    "parallel_type": "OpenMP",
+    "start_time": 0,
+    "end_time": 0.1
+  },
+  "solver_settings": {
+    "model_part_name": "IgaModelPart",
+    "domain_size": 2,
+    "echo_level": 1,
+    "buffer_size": 2,
+    "analysis_type": "non_linear",
+    "model_import_settings": { "input_type": "use_input_model_part" },
+    "material_import_settings": { "materials_filename": "materials.json" },
+    "time_stepping": { "time_step": 0.1 },
+    "rotation_dofs": false,
+    "reform_dofs_at_each_step": false,
+    "line_search": false,
+    "compute_reactions": true,
+    "block_builder": true,
+    "clear_storage": false,
+    "move_mesh_flag": false,
+    "convergence_criterion": "residual_criterion",
+    "displacement_relative_tolerance": 0.0000000001,
+    "displacement_absolute_tolerance": 0.0000000001,
+    "residual_relative_tolerance": 0.00000000001,
+    "residual_absolute_tolerance": 0.00000000001,
+    "max_iteration": 5,
+    "solver_type": "static",
+    "linear_solver_settings": {
+      "solver_type": "LinearSolversApplication.sparse_lu",
+      "write_system_matrix_to_file": true,
+      "max_iteration": 500,
+      "tolerance": 1E-09,
+      "scaling": false,
+      "verbosity": 0
+    },
+    "auxiliary_variables_list": [],
+    "auxiliary_dofs_list": [],
+    "auxiliary_reaction_list": []
+  },
+  "modelers": [
+    {
+      "modeler_name" : "ImportNurbsSbmModeler",
+      "Parameters" : {
+        "input_filename" : "nurbs_files/circle_nurbs.json",
+        "model_part_name" : "initial_skin_model_part_out",
+        "link_layer_to_condition_name": [
+          {
+            "layer_name" : "Layer0",
+            "condition_name" : "SbmSolidCondition"
+          }
+        ]
+        }
+    },
+    {
+      "modeler_name": "NurbsGeometryModelerSbm",
+      "Parameters": {
+          "echo_level": 1,
+          "model_part_name" : "IgaModelPart",
+          "lower_point_xyz": [-2.0,-2.0,0.0],
+          "upper_point_xyz": [2.0,2.0,0.0],
+          "lower_point_uvw": [-2.0,-2.0,0.0],
+          "upper_point_uvw": [2.0,2.0,0.0],
+          "polynomial_order" : [2, 2],
+          "number_of_knot_spans" : [20,20],
+          "lambda_outer": 0.5,
+          "number_of_inner_loops": 0,
+          "skin_model_part_outer_initial_name": "initial_skin_model_part_out",    
+          "skin_model_part_name": "skin_model_part"
+      }
+    },
+    {
+      "modeler_name": "IgaModelerSbm",
+      "Parameters": {
+      "echo_level": 0,
+      "skin_model_part_name": "skin_model_part",
+      "analysis_model_part_name": "IgaModelPart",
+      "element_condition_list": [
+          {
+              "geometry_type": "GeometrySurface",
+              "iga_model_part": "StructuralAnalysisDomain",
+              "type": "element",
+              "name": "SolidElement",
+              "shape_function_derivatives_order": 2
+          },
+          {
+            "geometry_type": "SurfaceEdge",
+            "brep_ids" : [2],
+            "iga_model_part": "SBM_Support_outer",
+            "type": "condition",
+            "name": "SbmSolidCondition",
+            "shape_function_derivatives_order": 6,
+            "sbm_parameters": {
+                "is_inner" : false
+            }
+          }
+      ]
+      }
+  }
+  ],
+  "processes": {
+    "additional_processes": [
+      {
+        "python_module": "assign_iga_external_conditions_process",
+        "kratos_module" : "KratosMultiphysics.IgaApplication",
+        "process_name"  : "AssignIgaExternalConditionsProcess",
+        "Parameters": {
+          "echo_level": 4,
+          "element_condition_list": [
+            {
+                "iga_model_part": "IgaModelPart.StructuralAnalysisDomain",
+                "variables": [
+                {
+                    "variable_name": "BODY_FORCE",
+                    "value": ["-1000*(1+0.3)/(1-0.3*0.3)*cos(x)*sinh(y)", 
+                             "-1000*(1+0.3)/(1-0.3*0.3)*sin(x)*cosh(y)", "0.0"]
+                }]
+            },
+            {
+                "iga_model_part": "IgaModelPart.SBM_Support_outer",
+                "variables": [
+                {
+                    "variable_name": "DISPLACEMENT",
+                    "value"           : ["-cos(x)*sinh(y)","sin(x)*cosh(y)", "0.0"]
+                }]
+            }
+          ]
+        }
+      }
+    ],
+    "dirichlet_process_list": [
+      {
+        "kratos_module": "KratosMultiphysics",
+        "python_module": "assign_vector_variable_to_nodes_process",
+        "Parameters": {
+          "model_part_name": "skin_model_part.outer",
+          "variable_name": "DISPLACEMENT",
+          "value"           : ["-cos(x)*sinh(y)","sin(x)*cosh(y)", "0.0"]
+        }
+      }
+    ],
+    "constraints_process_list"        : []
+},
+  "output_processes": {
+    "output_process_list": []
+  }
+}
+
+

iga/use_cases/external_boundary_circle_with_nurbs/convergence.py

@@ -0,0 +1,126 @@
+
+import KratosMultiphysics
+import KratosMultiphysics.IgaApplication
+from KratosMultiphysics.StructuralMechanicsApplication.structural_mechanics_analysis import StructuralMechanicsAnalysis
+
+import matplotlib.pyplot as plt
+import sympy as sp
+import numpy as np
+import math
+import os
+
+
+if __name__ == "__main__":
+
+    # Read original parameters once
+    with open('ProjectParameters.json', 'r') as f:
+        parameters = KratosMultiphysics.Parameters(f.read())
+
+    L2error_vector = []
+    computational_area_vec = []
+    h = []
+
+    insertion = [4, 8, 16, 32, 64]
+    degree = 2
+
+    tot = len(insertion)
+    for i in range(0,tot) :
+        insertions = insertion[i]
+        print("insertions: ", insertions) 
+
+        for i in range(parameters["modelers"].size()):
+            if parameters["modelers"][i]["modeler_name"].GetString() == "NurbsGeometryModelerSbm":
+                modeler_number = i
+                break
+
+        parameters["modelers"][modeler_number]["Parameters"]["number_of_knot_spans"] = KratosMultiphysics.Parameters(f"[{insertions}, {insertions}]")
+        parameters["modelers"][modeler_number]["Parameters"]["polynomial_order"] = KratosMultiphysics.Parameters(f"[{degree}, {degree}]")
+
+        model = KratosMultiphysics.Model()
+        simulation = StructuralMechanicsAnalysis(model,parameters)
+        simulation.Run()
+
+        # Exact solution as function handle:
+        x = sp.symbols('x')
+        y = sp.symbols('y')
+        u_exact = -sp.cos(x)*sp.sinh(y)
+        u_exact_handle = sp.lambdify((x, y), u_exact)
+
+        # Computation of the error:
+        output = []
+        x_coord = []
+        y_coord = []
+        weight = []
+
+        mp = model["IgaModelPart.StructuralAnalysisDomain"]
+        L2_err_temp = 0.0
+        L2_norm_solution = 0.0
+
+        computational_area = 0.0
+        for elem in mp.Elements:
+            geom = elem.GetGeometry()
+
+            N = geom.ShapeFunctionsValues()
+            weight = elem.GetValue(KratosMultiphysics.INTEGRATION_WEIGHT)
+
+            x = 0.0
+            y = 0.0
+            i = 0
+            curr_disp_x = 0
+            curr_disp_y = 0
+
+            for n in geom:
+                curr_disp_x += N[0, i] * n.GetSolutionStepValue(KratosMultiphysics.DISPLACEMENT_X)
+                curr_disp_y += N[0, i] * n.GetSolutionStepValue(KratosMultiphysics.DISPLACEMENT_Y)
+                x += N[0, i] * n.X
+                y += N[0, i] * n.Y
+                i += 1
+
+            x_coord.append(x)
+            y_coord.append(y)
+
+            L2_err_temp += (curr_disp_x - u_exact_handle(x,y))**2 * weight
+            L2_norm_solution += (u_exact_handle(x,y))**2 * weight
+
+            computational_area += weight
+
+        L2_err_temp = np.sqrt(L2_err_temp/L2_norm_solution)
+
+        L2error_vector.append(L2_err_temp)
+        computational_area_vec.append(computational_area)
+
+        lower_corner_coords = parameters["modelers"][1]["Parameters"]["lower_point_xyz"].GetVector()
+        upper_corner_coords = parameters["modelers"][1]["Parameters"]["upper_point_xyz"].GetVector()
+
+        h_canditate = max(upper_corner_coords[0] - lower_corner_coords[0], upper_corner_coords[1] - lower_corner_coords[1]) / insertions
+        h.append(h_canditate)
+
+        simulation.Clear()
+
+    print("\n h = ", h )
+    print('\n L2 error', L2error_vector)
+
+    plt.xscale('log')
+    plt.yscale('log')
+    plt.grid(True, which="both", linestyle='--')
+
+    stored = np.zeros(6)
+    stored[0] = (-(1)*(np.log(h[0])) + (np.log(L2error_vector[0])))
+    stored[1] = (-(2)*(np.log(h[0])) + (np.log(L2error_vector[0])))
+    stored[2] = (-(3)*(np.log(h[0])) + (np.log(L2error_vector[0])))
+    stored[3] = (-(4)*(np.log(h[0])) + (np.log(L2error_vector[0])))
+    stored[4] = (-(5)*(np.log(h[0])) + (np.log(L2error_vector[0])))
+    stored[5] = (-(6)*(np.log(h[0])) + (np.log(L2error_vector[0])))
+    degrees = np.arange(1, 7)
+    yDegrees = np.zeros((degrees.size, len(h)))
+    for i in range(0, degrees.size):
+        for jtest in range(0, len(h)):
+            yDegrees[i][jtest] = np.exp(degrees[i] * np.log(h[jtest]) + stored[i])
+        plt.plot(h, yDegrees[i], "--", label='vel %d' % tuple([degrees[i]]))
+
+    plt.plot(h, L2error_vector, 's-', markersize=1.5, linewidth=2, label='L^2 error with Weight')
+    plt.ylabel('L2 err',fontsize=14, color='blue')
+    plt.xlabel('h',fontsize=14, color='blue')
+    plt.legend(loc='lower right')
+    plt.show()
+

iga/use_cases/external_boundary_circle_with_nurbs/materials.json

@@ -0,0 +1,27 @@
+{
+  "properties": [
+    {
+      "model_part_name": "IgaModelPart",
+      "properties_id": 2,
+      "Material": {
+        "Variables": { "PENALTY_FACTOR": -1},
+        "Tables": {}
+      }
+    },
+    {
+      "model_part_name": "IgaModelPart.StructuralAnalysisDomain",
+      "properties_id": 2,
+      "Material": {
+        "name": "lin_el",
+        "constitutive_law": { "name": "LinearElasticPlaneStress2DLaw" },
+        "Variables": {
+          "THICKNESS": 1.0,
+          "YOUNG_MODULUS": 1000,
+          "POISSON_RATIO": 0.3,
+          "DENSITY": 1
+        },
+        "Tables": {}
+      }
+    }
+  ]
+}