diff --git a/iga/README.md b/iga/README.md
new file mode 100644
index 00000000..7e14057f
--- /dev/null
+++ b/iga/README.md
@@ -0,0 +1,8 @@
+# Iga Examples
+
+## Use Cases
+- [External Boundary Circle (NURBS) ](https://github.com/KratosMultiphysics/Examples/blob/master/iga/use_cases/README.md)
+
+## Validation Cases
+
+# Use Cases
diff --git a/iga/use_cases/README.md b/iga/use_cases/README.md
new file mode 100644
index 00000000..5216f315
--- /dev/null
+++ b/iga/use_cases/README.md
@@ -0,0 +1,66 @@
+# Use Case IGA – External Boundary Circle (NURBS)
+
+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
+- 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
+
+
+
+
+
+results of the convergence:
+
+ h = [1.0, 0.5, 0.25, 0.125, 0.0625]
+
+ L2_error = [0.04245443538478202, 0.00612519522524315, 0.0006828403386498715, 7.033648516213708e-05, 1.0075704685614167e-05]
+
+
diff --git a/iga/use_cases/external_boundary_circle_with_nurbs/ProjectParameters.json b/iga/use_cases/external_boundary_circle_with_nurbs/ProjectParameters.json
new file mode 100644
index 00000000..792b2f32
--- /dev/null
+++ b/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": []
+ }
+}
+
+
diff --git a/iga/use_cases/external_boundary_circle_with_nurbs/convergence.py b/iga/use_cases/external_boundary_circle_with_nurbs/convergence.py
new file mode 100644
index 00000000..29cfcc54
--- /dev/null
+++ b/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')
+ plt.ylabel('L2 err',fontsize=14, color='blue')
+ plt.xlabel('h',fontsize=14, color='blue')
+ plt.legend(loc='lower right')
+ plt.show()
+
diff --git a/iga/use_cases/external_boundary_circle_with_nurbs/materials.json b/iga/use_cases/external_boundary_circle_with_nurbs/materials.json
new file mode 100644
index 00000000..a64c8840
--- /dev/null
+++ b/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": {}
+ }
+ }
+ ]
+}
\ No newline at end of file
diff --git a/iga/use_cases/external_boundary_circle_with_nurbs/nurbs_files/circle_nurbs.json b/iga/use_cases/external_boundary_circle_with_nurbs/nurbs_files/circle_nurbs.json
new file mode 100644
index 00000000..0a604303
--- /dev/null
+++ b/iga/use_cases/external_boundary_circle_with_nurbs/nurbs_files/circle_nurbs.json
@@ -0,0 +1,116 @@
+{
+ "Surfaces": [
+ {
+ "Weights": [
+ 0,
+ 0,
+ 0,
+ 0
+ ],
+ "CPCoordinates": [
+ [
+ -1.0991860894187113,
+ -1.1000035979302478,
+ 0.0
+ ],
+ [
+ 1.0991860894187109,
+ -1.1000035979302478,
+ 0.0
+ ],
+ [
+ -1.0991860894187113,
+ 1.1000755565352007,
+ 0.0
+ ],
+ [
+ 1.0991860894187109,
+ 1.1000755565352007,
+ 0.0
+ ]
+ ],
+ "knotVectors": {
+ "Xi": [
+ 0.0,
+ 0.0,
+ 1.0,
+ 1.0
+ ],
+ "Eta": [
+ 0.0,
+ 0.0,
+ 1.0,
+ 1.0
+ ]
+ },
+ "id": 1,
+ "pDegreeU": 1,
+ "pDegreeV": 1
+ }
+ ],
+ "Lines": [
+ {
+ "Layer": "Layer0",
+ "Weights": [
+ 1.0,
+ 0.5000000000000001,
+ 1.0,
+ 0.5000000000000001,
+ 1.0,
+ 0.5000000000000001,
+ 1.0
+ ],
+ "CPCoordinates": [
+ [
+ 0.0,
+ 0.9999999999999984,
+ 0.0
+ ],
+ [
+ -1.7320508075688767,
+ 0.9999999999999984,
+ 0.0
+ ],
+ [
+ -0.8660254037844387,
+ -0.4999999999999998,
+ 0.0
+ ],
+ [
+ -2.449293598294706e-16,
+ -1.9999999999999996,
+ 0.0
+ ],
+ [
+ 0.8660254037844384,
+ -0.5000000000000004,
+ 0.0
+ ],
+ [
+ 1.7320508075688776,
+ 0.9999999999999984,
+ 0.0
+ ],
+ [
+ 2.4492935982947064e-16,
+ 0.9999999999999984,
+ 0.0
+ ]
+ ],
+ "knotVector": [
+ 0.0,
+ 0.0,
+ 0.0,
+ 0.3333333333333333,
+ 0.3333333333333333,
+ 0.6666666666666666,
+ 0.6666666666666666,
+ 1.0,
+ 1.0,
+ 1.0
+ ],
+ "id": 1,
+ "pDegree": 2
+ }
+ ]
+}
\ No newline at end of file
diff --git a/iga/use_cases/external_boundary_circle_with_nurbs/plot_surrogate_boundaries.py b/iga/use_cases/external_boundary_circle_with_nurbs/plot_surrogate_boundaries.py
new file mode 100644
index 00000000..f31dccc0
--- /dev/null
+++ b/iga/use_cases/external_boundary_circle_with_nurbs/plot_surrogate_boundaries.py
@@ -0,0 +1,165 @@
+import KratosMultiphysics
+import KratosMultiphysics.IgaApplication
+from KratosMultiphysics.StructuralMechanicsApplication.structural_mechanics_analysis import StructuralMechanicsAnalysis
+import numpy as np
+
+import matplotlib.pyplot as plt
+import matplotlib.patches as patches
+
+from matplotlib.lines import Line2D
+from matplotlib.patches import Patch
+
+if __name__ == "__main__":
+
+ with open("ProjectParameters.json",'r') as parameter_file:
+ parameters = KratosMultiphysics.Parameters(parameter_file.read())
+
+ model = KratosMultiphysics.Model()
+
+ simulation = StructuralMechanicsAnalysis(model, parameters)
+ simulation.Initialize()
+ simulation.InitializeSolutionStep()
+
+ index_param_space= 1
+
+ n_knot_span_u = parameters["modelers"][index_param_space]["Parameters"]["number_of_knot_spans"][0].GetInt()
+ n_knot_span_v = parameters["modelers"][index_param_space]["Parameters"]["number_of_knot_spans"][1].GetInt()
+ initial_u = parameters["modelers"][index_param_space]["Parameters"]["lower_point_uvw"][0].GetDouble()
+ initial_v = parameters["modelers"][index_param_space]["Parameters"]["lower_point_uvw"][1].GetDouble()
+ final_u = parameters["modelers"][index_param_space]["Parameters"]["upper_point_uvw"][0].GetDouble()
+ final_v = parameters["modelers"][index_param_space]["Parameters"]["upper_point_uvw"][1].GetDouble()
+
+ knots_u_1 = [initial_u]
+ knots_v_1 = [initial_v]
+
+ for j in range(1, n_knot_span_u):
+ knots_u_1.append(initial_u + (final_u-initial_u) / (n_knot_span_u) * j)
+ for j in range(1, n_knot_span_v):
+ knots_v_1.append(initial_v + (final_v-initial_v) / (n_knot_span_v) * j)
+ knots_u_1.append(final_u)
+ knots_v_1.append(final_v)
+
+ # Create figure and axes
+ fig, ax = plt.subplots()
+ ax.set_aspect('equal', adjustable='box')
+ ax.grid(False)
+ label_used = set()
+
+ # === PARAMETER SPACE GRID: BODY 1 ===
+ for u in knots_u_1:
+ ax.plot([u, u], [knots_v_1[0], knots_v_1[-1]], linestyle=':', color='gray', linewidth=1.5, label='Parameter space' if u == knots_u_1[0] else "")
+ for v in knots_v_1:
+ ax.plot([knots_u_1[0], knots_u_1[-1]], [v, v], linestyle=':', color='gray', linewidth=1.5)
+
+ label_used.add("Parameter space")
+
+ for side in ["outer", "inner"]:
+ surrogate_part_name = f"IgaModelPart.surrogate_{side}"
+ if model.HasModelPart(surrogate_part_name):
+ mp = model.GetModelPart(surrogate_part_name)
+ if len(mp.Conditions) > 0:
+ for cond in mp.Conditions:
+ x_surr, y_surr = [], []
+ geom = cond.GetGeometry()
+ x_surr.append(geom[0].X)
+ y_surr.append(geom[0].Y)
+ x_surr.append(geom[1].X)
+ y_surr.append(geom[1].Y)
+ ax.plot(x_surr, y_surr, color='darkred', linewidth=5, label="Surrogate boundary")
+ label_used.add("Surrogate boundary")
+
+ skin_part_name = f"skin_model_part.{side}"
+ if model.HasModelPart(skin_part_name):
+ mp = model.GetModelPart(skin_part_name)
+ if len(mp.Conditions) > 0:
+ x_skin, y_skin = [], []
+ for i, cond in enumerate(mp.Conditions):
+ geom = cond.GetGeometry()
+ if i == 0:
+ x_skin.append(geom[0].X)
+ y_skin.append(geom[0].Y)
+ x_skin.append(geom[1].X)
+ y_skin.append(geom[1].Y)
+ x_skin.append(x_skin[0])
+ y_skin.append(y_skin[0])
+ ax.plot(x_skin, y_skin, color='deepskyblue', linewidth=2, label="True boundary")
+ label_used.add("True boundary")
+
+ gp_part_name = f"IgaModelPart.SBM_Support_{side}"
+ if model.HasModelPart(gp_part_name):
+ mp = model.GetModelPart(gp_part_name)
+ if len(mp.Conditions) > 0:
+ for cond in mp.Conditions:
+ geom = cond.GetGeometry()
+ N = geom.ShapeFunctionsValues()
+ gp = np.zeros(2)
+ for i, node in enumerate(geom):
+ gp[0] += N[0, i] * node.X
+ gp[1] += N[0, i] * node.Y
+ # plot gp, arrows, etc...
+
+
+ # === ACTIVE GAUSS POINTS: FROM ELEMENT CENTERS ===
+ mp = model["IgaModelPart.StructuralAnalysisDomain"]
+ active_gauss_points = []
+ for elem in mp.Elements:
+ center = elem.GetGeometry().Center()
+ x = center.X
+ y = center.Y
+ active_gauss_points.append((x, y))
+ ax.plot(x, y, "bo", markersize=3, label="Elemental Gauss Points" if "Elemental Gauss Points" not in label_used else "")
+ label_used.add("Elemental Gauss Points")
+
+ # === FIND ACTIVE CELLS AND COLOR THEM LIGHT YELLOW ===
+ active_cells = set()
+ for (x_gp, y_gp) in active_gauss_points:
+ for i in range(len(knots_u_1) - 1):
+ for j in range(len(knots_v_1) - 1):
+ u_min, u_max = knots_u_1[i], knots_u_1[i+1]
+ v_min, v_max = knots_v_1[j], knots_v_1[j+1]
+ if u_min <= x_gp <= u_max and v_min <= y_gp <= v_max:
+ active_cells.add((i, j))
+
+ for (i, j) in active_cells:
+ u_min, u_max = knots_u_1[i], knots_u_1[i+1]
+ v_min, v_max = knots_v_1[j], knots_v_1[j+1]
+ rect = patches.Rectangle((u_min, v_min), u_max - u_min, v_max - v_min,
+ linewidth=0, edgecolor=None, facecolor="#FFFACD", zorder=0,
+ label="Active elements" if "Active elements" not in label_used else "")
+ ax.add_patch(rect)
+ label_used.add("Active elements")
+
+ # Custom handles
+ custom_handles = []
+
+ if "Parameter space" in label_used:
+ custom_handles.append(Line2D([0], [0], color='gray', linestyle=':', label='Parameter space'))
+
+ if "Surrogate boundary" in label_used:
+ custom_handles.append(Line2D([0], [0], color='darkred', linewidth=2, label='Surrogate boundary'))
+
+ if "True boundary" in label_used:
+ custom_handles.append(Line2D([0], [0], color='deepskyblue', linewidth=2, label='True boundary'))
+
+ if "Elemental Gauss Points" in label_used:
+ custom_handles.append(Line2D([0], [0], marker='o', color='blue', linestyle='None', markersize=5, label='Elemental Gauss Points'))
+
+ if "Boundary Gauss Points" in label_used:
+ custom_handles.append(Line2D([0], [0], marker='o', color='red', linestyle='None', markersize=5, label='Boundary Gauss Points'))
+
+ if "Active elements" in label_used:
+ custom_handles.append(Patch(facecolor='#FFFACD', edgecolor='goldenrod', linestyle=':', linewidth=1.5, label='Active elements'))
+
+ # Clean layout for publication
+ ax.set_aspect('equal')
+ ax.set_xlabel("") # or add label if meaningful
+ ax.set_ylabel("")
+ ax.tick_params(left=False, bottom=False, labelleft=False, labelbottom=False) # hide ticks
+ plt.tight_layout()
+ plt.subplots_adjust(right=0.8) # extra space for legend
+
+ ax.legend(handles=custom_handles, loc='center left', bbox_to_anchor=(1, 0.5))
+ plt.show()
+
+
+
diff --git a/iga/use_cases/external_boundary_circle_with_nurbs/run_and_post.py b/iga/use_cases/external_boundary_circle_with_nurbs/run_and_post.py
new file mode 100644
index 00000000..ff99e7d1
--- /dev/null
+++ b/iga/use_cases/external_boundary_circle_with_nurbs/run_and_post.py
@@ -0,0 +1,163 @@
+import KratosMultiphysics
+import KratosMultiphysics.IgaApplication
+from KratosMultiphysics.StructuralMechanicsApplication.structural_mechanics_analysis import StructuralMechanicsAnalysis
+
+import numpy as np
+import matplotlib.pyplot as plt
+import matplotlib.tri as mtri
+
+class CustomAnalysisStage(StructuralMechanicsAnalysis):
+ time_step = []
+ nSTEP = 0
+
+ def InitializeSolutionStep(self):
+ super().InitializeSolutionStep()
+ current_time = self._GetSolver().GetComputingModelPart().ProcessInfo[KratosMultiphysics.TIME]
+ self.time_step.append(current_time)
+ self.nSTEP += 1
+
+# Define the analytical solution
+def analytical_temperature(x, y, t):
+ return -np.cos(x)*np.sinh(y)
+
+# Main simulation function
+if __name__ == "__main__":
+ with open("ProjectParameters.json", 'r') as parameter_file:
+ parameters = KratosMultiphysics.Parameters(parameter_file.read())
+
+ model = KratosMultiphysics.Model()
+ simulation = CustomAnalysisStage(model, parameters)
+ simulation.Run()
+
+ # Extract the computed solution at a specific time step
+ mp = model["IgaModelPart"]
+
+ x_coord = []
+ y_coord = []
+ computed_temperature = []
+ weights = []
+
+ # Loop over elements to gather computed solution
+ for elem in mp.Elements:
+ if (elem.Id == 1):
+ continue
+ geom = elem.GetGeometry()
+ N = geom.ShapeFunctionsValues()
+ center = geom.Center()
+ weight = elem.GetValue(KratosMultiphysics.INTEGRATION_WEIGHT)
+ weights.append(weight)
+ # Extract Gauss point (center) coordinates
+ x_coord.append(center.X)
+ y_coord.append(center.Y)
+
+ # Initialize solution values at the center
+ curr_temperature = 0
+ # Compute nodal contributions using shape functions
+ for i, node in enumerate(geom):
+ curr_temperature += N[0, i] * node.GetSolutionStepValue(KratosMultiphysics.DISPLACEMENT_X, 0)
+ computed_temperature.append(curr_temperature)
+
+ # Get the current time after simulation run
+ current_time = simulation._GetSolver().GetComputingModelPart().ProcessInfo[KratosMultiphysics.TIME]
+ print("Current time after simulation:", current_time)
+
+ # Calculate errors and analytical solution
+ temperature_error = []
+ temperature_analytical_values = []
+
+ for i in range(len(x_coord)):
+ x = x_coord[i]
+ y = y_coord[i]
+
+ # Analytical solution at the point
+ temperature_analytical = analytical_temperature(x, y, current_time)
+
+ # Append analytical values for plotting
+ temperature_analytical_values.append(temperature_analytical)
+
+ # Error calculations
+ temperature_error.append(abs(computed_temperature[i] - temperature_analytical))
+
+ # Compute the L2 norm of the error for temperature
+ temperature_l2_norm = 0
+ for i in range(len(weights)):
+ temperature_l2_norm += weights[i] * temperature_error[i]**2
+ # Take the square root to finalize the L2 norm
+ temperature_l2_norm = (temperature_l2_norm)**0.5
+
+ # Compute the L2 norm of the analytical solution
+ analytical_l2_norm = 0
+ for i in range(len(weights)):
+ analytical_l2_norm += weights[i] * temperature_analytical_values[i]**2
+ analytical_l2_norm = (analytical_l2_norm)**0.5
+
+ # Compute the normalized L2 error
+ normalized_l2_error = temperature_l2_norm / analytical_l2_norm
+
+ # Print the results
+ print("L2 norm of temperature error:", temperature_l2_norm)
+ # print("L2 norm of analytical solution:", analytical_l2_norm)
+ print("Normalized L2 error:", normalized_l2_error)
+
+ # Plot computed solution and analytical solution
+ fig, ax = plt.subplots(1, 2, figsize=(12, 6))
+
+ # triangulation = mtri.Triangulation(x_coord, y_coord)
+ # Create a triangulation of the data points
+ triangulation = mtri.Triangulation(x_coord, y_coord)
+ areas = []
+ for triangle in triangulation.triangles:
+ x1, y1 = x_coord[triangle[0]], y_coord[triangle[0]]
+ x2, y2 = x_coord[triangle[1]], y_coord[triangle[1]]
+ x3, y3 = x_coord[triangle[2]], y_coord[triangle[2]]
+ # Calculate the area of the triangle using the shoelace formula
+ area = 0.5 * abs(x1 * (y2 - y3) + x2 * (y3 - y1) + x3 * (y1 - y2))
+ areas.append(area)
+ average_area = sum(areas) / len(areas)
+ new_triangles = []
+ for i in range(len(areas)):
+ # if areas[i] <= average_area*2.09:
+ # if areas[i] <= average_area*1.7:
+ if areas[i] <= average_area*1.5:
+ new_triangles.append(triangulation.triangles[i])
+ new_triangulation = mtri.Triangulation(x_coord, y_coord, new_triangles)
+
+ # Compute shared color limits (optional but keeps colormap consistent)
+ vmin = min(np.min(computed_temperature), np.min(temperature_analytical_values))
+ vmax = max(np.max(computed_temperature), np.max(temperature_analytical_values))
+
+ # Computed solution
+ tcf0 = ax[0].tricontourf(new_triangulation, computed_temperature, cmap='jet', vmin=vmin, vmax=vmax)
+ ax[0].set_title("Computed Temperature")
+ ax[0].set_xlabel("x")
+ ax[0].set_ylabel("y")
+ ax[0].set_aspect('equal')
+
+ # Analytical solution
+ tcf1 = ax[1].tricontourf(new_triangulation, temperature_analytical_values, cmap='jet', vmin=vmin, vmax=vmax)
+ ax[1].set_title("Analytical Temperature")
+ ax[1].set_xlabel("x")
+ ax[1].set_aspect('equal')
+
+ # Correct way to add colorbars
+ plt.colorbar(tcf0, ax=ax[0], orientation='vertical')
+ plt.colorbar(tcf1, ax=ax[1], orientation='vertical')
+
+ plt.show()
+
+
+ # Plot error distribution
+ fig, ax = plt.subplots(1, 1, figsize=(6, 6))
+
+ # Plot temperature error
+ tcf = ax.tricontourf(new_triangulation, temperature_error, cmap='jet')
+ ax.set_title("Error in Temperature")
+ ax.set_xlabel("x")
+ ax.set_ylabel("y")
+ ax.set_aspect('equal')
+ # ax.plot(x_coord, y_coord, 'y*', markersize=4, label='Gauss Points')
+
+ # Adding colorbar using the tricontourf return
+ plt.colorbar(tcf, ax=ax, orientation='vertical')
+ plt.show()
+
diff --git a/iga/validation/README.md b/iga/validation/README.md
new file mode 100644
index 00000000..78c1f8ac
--- /dev/null
+++ b/iga/validation/README.md
@@ -0,0 +1,4 @@
+# Validation Cases
+
+This folder contains the validation cases:
+