From f9bc5aa619a7826aad5e6c1915186b067facba87 Mon Sep 17 00:00:00 2001 From: NickNick9 Date: Wed, 18 Jun 2025 14:45:07 +0200 Subject: [PATCH 1/7] first commit --- iga/README.md | 5 +++++ iga/use_cases/README.md | 4 ++++ iga/validation/README.md | 4 ++++ 3 files changed, 13 insertions(+) create mode 100644 iga/README.md create mode 100644 iga/use_cases/README.md create mode 100644 iga/validation/README.md diff --git a/iga/README.md b/iga/README.md new file mode 100644 index 00000000..6ea459d9 --- /dev/null +++ b/iga/README.md @@ -0,0 +1,5 @@ +# Fluid Dynamics Examples + +## Use Cases + +## Validation Cases diff --git a/iga/use_cases/README.md b/iga/use_cases/README.md new file mode 100644 index 00000000..095f34bc --- /dev/null +++ b/iga/use_cases/README.md @@ -0,0 +1,4 @@ +# Use Cases + +This folder contains the validation cases: + 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: + From 8cad31be509eecf834fd51a1d98baec7a7054bb9 Mon Sep 17 00:00:00 2001 From: NickNick9 Date: Wed, 18 Jun 2025 14:49:43 +0200 Subject: [PATCH 2/7] example 1 data --- .../ProjectParameters.json | 152 ++++++++++++++ .../convergence.py | 126 ++++++++++++ .../materials.json | 27 +++ .../nurbs_files/circle_nurbs.json | 116 +++++++++++ .../plot_surrogate_boundaries.py | 167 ++++++++++++++++ .../run_and_post.py | 186 ++++++++++++++++++ 6 files changed, 774 insertions(+) create mode 100644 iga/use_cases/external_boundary_circle_with_nurbs/ProjectParameters.json create mode 100644 iga/use_cases/external_boundary_circle_with_nurbs/convergence.py create mode 100644 iga/use_cases/external_boundary_circle_with_nurbs/materials.json create mode 100644 iga/use_cases/external_boundary_circle_with_nurbs/nurbs_files/circle_nurbs.json create mode 100644 iga/use_cases/external_boundary_circle_with_nurbs/plot_surrogate_boundaries.py create mode 100644 iga/use_cases/external_boundary_circle_with_nurbs/run_and_post.py 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..3e9fb78b --- /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 with Weight') + 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..0875fac9 --- /dev/null +++ b/iga/use_cases/external_boundary_circle_with_nurbs/plot_surrogate_boundaries.py @@ -0,0 +1,167 @@ +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..b17f4cc3 --- /dev/null +++ b/iga/use_cases/external_boundary_circle_with_nurbs/run_and_post.py @@ -0,0 +1,186 @@ +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 x**2*y + return -np.cos(x)*np.sinh(y) + # return x**3+y**3 + # return x**2+y**2 + # return x+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) + + + # # Loop over elements to gather computed solution + # for cond in mp.Conditions: + # geom = cond.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.TEMPERATURE, 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() + From cc6631a34a7d018235932e91d6ce9cd43e1fc948 Mon Sep 17 00:00:00 2001 From: NickNick9 Date: Fri, 5 Sep 2025 16:23:17 +0200 Subject: [PATCH 3/7] minor stuff --- .../plot_surrogate_boundaries.py | 2 -- .../run_and_post.py | 23 ------------------- 2 files changed, 25 deletions(-) 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 index 0875fac9..f31dccc0 100644 --- 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 @@ -99,8 +99,6 @@ # plot gp, arrows, etc... - - # === ACTIVE GAUSS POINTS: FROM ELEMENT CENTERS === mp = model["IgaModelPart.StructuralAnalysisDomain"] active_gauss_points = [] 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 index b17f4cc3..ff99e7d1 100644 --- 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 @@ -18,11 +18,7 @@ def InitializeSolutionStep(self): # Define the analytical solution def analytical_temperature(x, y, t): - # return x**2*y return -np.cos(x)*np.sinh(y) - # return x**3+y**3 - # return x**2+y**2 - # return x+y # Main simulation function if __name__ == "__main__": @@ -61,25 +57,6 @@ def analytical_temperature(x, y, t): curr_temperature += N[0, i] * node.GetSolutionStepValue(KratosMultiphysics.DISPLACEMENT_X, 0) computed_temperature.append(curr_temperature) - - # # Loop over elements to gather computed solution - # for cond in mp.Conditions: - # geom = cond.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.TEMPERATURE, 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) From 74dc3731250189e283e204ad2ec22998254dc9e1 Mon Sep 17 00:00:00 2001 From: NickNick9 Date: Fri, 5 Sep 2025 16:34:10 +0200 Subject: [PATCH 4/7] readme --- iga/README.md | 4 ++- iga/use_cases/README.md | 56 +++++++++++++++++++++++++++++++++++++++-- 2 files changed, 57 insertions(+), 3 deletions(-) diff --git a/iga/README.md b/iga/README.md index 6ea459d9..43c60af9 100644 --- a/iga/README.md +++ b/iga/README.md @@ -1,5 +1,7 @@ -# Fluid Dynamics Examples +# Iga Examples ## Use Cases +- External Boundary Circle (NURBS) ## Validation Cases +- External Boundary Circle (NURBS) \ No newline at end of file diff --git a/iga/use_cases/README.md b/iga/use_cases/README.md index 095f34bc..6ecb3fd8 100644 --- a/iga/use_cases/README.md +++ b/iga/use_cases/README.md @@ -1,4 +1,56 @@ -# Use Cases +# 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 -This folder contains the validation cases: From a3e9e5a78e75898ee634ba4acf7c6762a93cc43a Mon Sep 17 00:00:00 2001 From: NickNick9 Date: Fri, 5 Sep 2025 17:22:38 +0200 Subject: [PATCH 5/7] correct readme --- iga/README.md | 5 +++-- iga/use_cases/README.md | 4 ++-- 2 files changed, 5 insertions(+), 4 deletions(-) diff --git a/iga/README.md b/iga/README.md index 43c60af9..7e14057f 100644 --- a/iga/README.md +++ b/iga/README.md @@ -1,7 +1,8 @@ # Iga Examples ## Use Cases -- External Boundary Circle (NURBS) +- [External Boundary Circle (NURBS) ](https://github.com/KratosMultiphysics/Examples/blob/master/iga/use_cases/README.md) ## Validation Cases -- External Boundary Circle (NURBS) \ No newline at end of file + +# Use Cases diff --git a/iga/use_cases/README.md b/iga/use_cases/README.md index 6ecb3fd8..fec37dfe 100644 --- a/iga/use_cases/README.md +++ b/iga/use_cases/README.md @@ -1,4 +1,4 @@ -# Use Cases IGA – External Boundary Circle (NURBS) — README.txt +# 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. @@ -9,7 +9,7 @@ Files - nurbs_files/circle_nurbs.json - run_and_post.py (optional) -Geometry / modeling (very brief) +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: From f2e887b8ee3142577cc857775b6283fd5dc5c1a1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Nicol=C3=B2=20Antonelli?= <114093318+NickNick9@users.noreply.github.com> Date: Fri, 5 Sep 2025 17:31:03 +0200 Subject: [PATCH 6/7] images --- iga/use_cases/README.md | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/iga/use_cases/README.md b/iga/use_cases/README.md index fec37dfe..5216f315 100644 --- a/iga/use_cases/README.md +++ b/iga/use_cases/README.md @@ -53,4 +53,14 @@ Run: python3 convergence.py # for the convergence analysis python3 plot_surrogate_boundaries.py # to visualize the surrogate and true boundaries +image + +image + +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] + From 4e97e65086abda6f83b4ba4d6ec3b6435b393a52 Mon Sep 17 00:00:00 2001 From: NickNick9 Date: Fri, 5 Sep 2025 17:32:19 +0200 Subject: [PATCH 7/7] last commiy --- .../external_boundary_circle_with_nurbs/convergence.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/iga/use_cases/external_boundary_circle_with_nurbs/convergence.py b/iga/use_cases/external_boundary_circle_with_nurbs/convergence.py index 3e9fb78b..29cfcc54 100644 --- a/iga/use_cases/external_boundary_circle_with_nurbs/convergence.py +++ b/iga/use_cases/external_boundary_circle_with_nurbs/convergence.py @@ -118,7 +118,7 @@ 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.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')