Skip to content

Commit f6be8d8

Browse files
NickNick9NickNick9
authored
Merge pull request #147 from KratosMultiphysics/PR_example1
[IgaApplication] example with SBM in 2d with convergence study importing a nurbs
2 parents 12e2e3c + 4e97e65 commit f6be8d8

File tree

9 files changed

+827
-0
lines changed

9 files changed

+827
-0
lines changed

‎iga/README.md‎

Lines changed: 8 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,8 @@
1+
# Iga Examples
2+
3+
## Use Cases
4+
- [External Boundary Circle (NURBS) ](https://github.com/KratosMultiphysics/Examples/blob/master/iga/use_cases/README.md)
5+
6+
## Validation Cases
7+
8+
# Use Cases

‎iga/use_cases/README.md‎

Lines changed: 66 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,66 @@
1+
# Use Case IGA – External Boundary Circle (NURBS)
2+
3+
Short description
4+
This example runs a 2D IGA structural model whose outer boundary is a NURBS circle. A manufactured solution is enforced to validate the setup.
5+
6+
Files
7+
- ProjectParameters.json
8+
- materials.json
9+
- nurbs_files/circle_nurbs.json
10+
- run_and_post.py (optional)
11+
12+
Geometry
13+
- ImportNurbsSbmModeler loads the circle NURBS into “initial_skin_model_part_out”.
14+
- NurbsGeometryModelerSbm builds the analysis surface (order [2,2], knot spans [20,20]) and creates “skin_model_part”.
15+
- IgaModelerSbm creates:
16+
• SolidElement on IgaModelPart.StructuralAnalysisDomain
17+
• SbmSolidCondition on SurfaceEdge (brep_ids: [2]) → IgaModelPart.SBM_Support_outer
18+
19+
Manufactured BCs and loads
20+
- Displacement on outer boundary:
21+
u = [ -cos(x)*sinh(y), sin(x)*cosh(y), 0.0 ]
22+
- Body force in the domain (consistent with the above):
23+
BODY_FORCE = [
24+
-1000*(1+0.3)/(1-0.3*0.3)*cos(x)*sinh(y),
25+
-1000*(1+0.3)/(1-0.3*0.3)*sin(x)*cosh(y),
26+
0.0
27+
]
28+
29+
Solver (key settings)
30+
- Static, nonlinear; OpenMP
31+
- Single time step: dt = 0.1, end_time = 0.1
32+
- Linear solver: LinearSolversApplication.sparse_lu
33+
34+
Minimal materials.json (example)
35+
Use plane strain (or plane stress) linear elastic; keep it consistent with ν=0.3:
36+
{
37+
"properties": [{
38+
"model_part_name": "IgaModelPart.StructuralAnalysisDomain",
39+
"properties_id": 1,
40+
"Material": {
41+
"constitutive_law": { "name": "LinearElasticPlaneStrain2DLaw" },
42+
"Variables": {
43+
"YOUNG_MODULUS": 1000.0,
44+
"POISSON_RATIO": 0.3,
45+
"DENSITY": 1.0
46+
}
47+
}
48+
}]
49+
}
50+
51+
Run:
52+
python3 run_and_post.py # to reproduce the error
53+
python3 convergence.py # for the convergence analysis
54+
python3 plot_surrogate_boundaries.py # to visualize the surrogate and true boundaries
55+
56+
<img width="999" height="786" alt="image" src="https://github.com/user-attachments/assets/ea1c01a6-806c-486a-b11e-d3846bac3b85" />
57+
58+
<img width="588" height="430" alt="image" src="https://github.com/user-attachments/assets/9401a33e-4c81-47de-bc9b-a0323fdb89a9" />
59+
60+
results of the convergence:
61+
62+
h = [1.0, 0.5, 0.25, 0.125, 0.0625]
63+
64+
L2_error = [0.04245443538478202, 0.00612519522524315, 0.0006828403386498715, 7.033648516213708e-05, 1.0075704685614167e-05]
65+
66+

‎iga/use_cases/external_boundary_circle_with_nurbs/ProjectParameters.json‎

Lines changed: 152 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,152 @@
1+
{
2+
"problem_data": {
3+
"problem_name": "example_1_kratos",
4+
"echo_level": 0,
5+
"parallel_type": "OpenMP",
6+
"start_time": 0,
7+
"end_time": 0.1
8+
},
9+
"solver_settings": {
10+
"model_part_name": "IgaModelPart",
11+
"domain_size": 2,
12+
"echo_level": 1,
13+
"buffer_size": 2,
14+
"analysis_type": "non_linear",
15+
"model_import_settings": { "input_type": "use_input_model_part" },
16+
"material_import_settings": { "materials_filename": "materials.json" },
17+
"time_stepping": { "time_step": 0.1 },
18+
"rotation_dofs": false,
19+
"reform_dofs_at_each_step": false,
20+
"line_search": false,
21+
"compute_reactions": true,
22+
"block_builder": true,
23+
"clear_storage": false,
24+
"move_mesh_flag": false,
25+
"convergence_criterion": "residual_criterion",
26+
"displacement_relative_tolerance": 0.0000000001,
27+
"displacement_absolute_tolerance": 0.0000000001,
28+
"residual_relative_tolerance": 0.00000000001,
29+
"residual_absolute_tolerance": 0.00000000001,
30+
"max_iteration": 5,
31+
"solver_type": "static",
32+
"linear_solver_settings": {
33+
"solver_type": "LinearSolversApplication.sparse_lu",
34+
"write_system_matrix_to_file": true,
35+
"max_iteration": 500,
36+
"tolerance": 1E-09,
37+
"scaling": false,
38+
"verbosity": 0
39+
},
40+
"auxiliary_variables_list": [],
41+
"auxiliary_dofs_list": [],
42+
"auxiliary_reaction_list": []
43+
},
44+
"modelers": [
45+
{
46+
"modeler_name" : "ImportNurbsSbmModeler",
47+
"Parameters" : {
48+
"input_filename" : "nurbs_files/circle_nurbs.json",
49+
"model_part_name" : "initial_skin_model_part_out",
50+
"link_layer_to_condition_name": [
51+
{
52+
"layer_name" : "Layer0",
53+
"condition_name" : "SbmSolidCondition"
54+
}
55+
]
56+
}
57+
},
58+
{
59+
"modeler_name": "NurbsGeometryModelerSbm",
60+
"Parameters": {
61+
"echo_level": 1,
62+
"model_part_name" : "IgaModelPart",
63+
"lower_point_xyz": [-2.0,-2.0,0.0],
64+
"upper_point_xyz": [2.0,2.0,0.0],
65+
"lower_point_uvw": [-2.0,-2.0,0.0],
66+
"upper_point_uvw": [2.0,2.0,0.0],
67+
"polynomial_order" : [2, 2],
68+
"number_of_knot_spans" : [20,20],
69+
"lambda_outer": 0.5,
70+
"number_of_inner_loops": 0,
71+
"skin_model_part_outer_initial_name": "initial_skin_model_part_out",
72+
"skin_model_part_name": "skin_model_part"
73+
}
74+
},
75+
{
76+
"modeler_name": "IgaModelerSbm",
77+
"Parameters": {
78+
"echo_level": 0,
79+
"skin_model_part_name": "skin_model_part",
80+
"analysis_model_part_name": "IgaModelPart",
81+
"element_condition_list": [
82+
{
83+
"geometry_type": "GeometrySurface",
84+
"iga_model_part": "StructuralAnalysisDomain",
85+
"type": "element",
86+
"name": "SolidElement",
87+
"shape_function_derivatives_order": 2
88+
},
89+
{
90+
"geometry_type": "SurfaceEdge",
91+
"brep_ids" : [2],
92+
"iga_model_part": "SBM_Support_outer",
93+
"type": "condition",
94+
"name": "SbmSolidCondition",
95+
"shape_function_derivatives_order": 6,
96+
"sbm_parameters": {
97+
"is_inner" : false
98+
}
99+
}
100+
]
101+
}
102+
}
103+
],
104+
"processes": {
105+
"additional_processes": [
106+
{
107+
"python_module": "assign_iga_external_conditions_process",
108+
"kratos_module" : "KratosMultiphysics.IgaApplication",
109+
"process_name" : "AssignIgaExternalConditionsProcess",
110+
"Parameters": {
111+
"echo_level": 4,
112+
"element_condition_list": [
113+
{
114+
"iga_model_part": "IgaModelPart.StructuralAnalysisDomain",
115+
"variables": [
116+
{
117+
"variable_name": "BODY_FORCE",
118+
"value": ["-1000*(1+0.3)/(1-0.3*0.3)*cos(x)*sinh(y)",
119+
"-1000*(1+0.3)/(1-0.3*0.3)*sin(x)*cosh(y)", "0.0"]
120+
}]
121+
},
122+
{
123+
"iga_model_part": "IgaModelPart.SBM_Support_outer",
124+
"variables": [
125+
{
126+
"variable_name": "DISPLACEMENT",
127+
"value" : ["-cos(x)*sinh(y)","sin(x)*cosh(y)", "0.0"]
128+
}]
129+
}
130+
]
131+
}
132+
}
133+
],
134+
"dirichlet_process_list": [
135+
{
136+
"kratos_module": "KratosMultiphysics",
137+
"python_module": "assign_vector_variable_to_nodes_process",
138+
"Parameters": {
139+
"model_part_name": "skin_model_part.outer",
140+
"variable_name": "DISPLACEMENT",
141+
"value" : ["-cos(x)*sinh(y)","sin(x)*cosh(y)", "0.0"]
142+
}
143+
}
144+
],
145+
"constraints_process_list" : []
146+
},
147+
"output_processes": {
148+
"output_process_list": []
149+
}
150+
}
151+
152+

‎iga/use_cases/external_boundary_circle_with_nurbs/convergence.py‎

Lines changed: 126 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,126 @@
1+
2+
import KratosMultiphysics
3+
import KratosMultiphysics.IgaApplication
4+
from KratosMultiphysics.StructuralMechanicsApplication.structural_mechanics_analysis import StructuralMechanicsAnalysis
5+
6+
import matplotlib.pyplot as plt
7+
import sympy as sp
8+
import numpy as np
9+
import math
10+
import os
11+
12+
13+
if __name__ == "__main__":
14+
15+
# Read original parameters once
16+
with open('ProjectParameters.json', 'r') as f:
17+
parameters = KratosMultiphysics.Parameters(f.read())
18+
19+
L2error_vector = []
20+
computational_area_vec = []
21+
h = []
22+
23+
insertion = [4, 8, 16, 32, 64]
24+
degree = 2
25+
26+
tot = len(insertion)
27+
for i in range(0,tot) :
28+
insertions = insertion[i]
29+
print("insertions: ", insertions)
30+
31+
for i in range(parameters["modelers"].size()):
32+
if parameters["modelers"][i]["modeler_name"].GetString() == "NurbsGeometryModelerSbm":
33+
modeler_number = i
34+
break
35+
36+
parameters["modelers"][modeler_number]["Parameters"]["number_of_knot_spans"] = KratosMultiphysics.Parameters(f"[{insertions}, {insertions}]")
37+
parameters["modelers"][modeler_number]["Parameters"]["polynomial_order"] = KratosMultiphysics.Parameters(f"[{degree}, {degree}]")
38+
39+
model = KratosMultiphysics.Model()
40+
simulation = StructuralMechanicsAnalysis(model,parameters)
41+
simulation.Run()
42+
43+
# Exact solution as function handle:
44+
x = sp.symbols('x')
45+
y = sp.symbols('y')
46+
u_exact = -sp.cos(x)*sp.sinh(y)
47+
u_exact_handle = sp.lambdify((x, y), u_exact)
48+
49+
# Computation of the error:
50+
output = []
51+
x_coord = []
52+
y_coord = []
53+
weight = []
54+
55+
mp = model["IgaModelPart.StructuralAnalysisDomain"]
56+
L2_err_temp = 0.0
57+
L2_norm_solution = 0.0
58+
59+
computational_area = 0.0
60+
for elem in mp.Elements:
61+
geom = elem.GetGeometry()
62+
63+
N = geom.ShapeFunctionsValues()
64+
weight = elem.GetValue(KratosMultiphysics.INTEGRATION_WEIGHT)
65+
66+
x = 0.0
67+
y = 0.0
68+
i = 0
69+
curr_disp_x = 0
70+
curr_disp_y = 0
71+
72+
for n in geom:
73+
curr_disp_x += N[0, i] * n.GetSolutionStepValue(KratosMultiphysics.DISPLACEMENT_X)
74+
curr_disp_y += N[0, i] * n.GetSolutionStepValue(KratosMultiphysics.DISPLACEMENT_Y)
75+
x += N[0, i] * n.X
76+
y += N[0, i] * n.Y
77+
i += 1
78+
79+
x_coord.append(x)
80+
y_coord.append(y)
81+
82+
L2_err_temp += (curr_disp_x - u_exact_handle(x,y))**2 * weight
83+
L2_norm_solution += (u_exact_handle(x,y))**2 * weight
84+
85+
computational_area += weight
86+
87+
L2_err_temp = np.sqrt(L2_err_temp/L2_norm_solution)
88+
89+
L2error_vector.append(L2_err_temp)
90+
computational_area_vec.append(computational_area)
91+
92+
lower_corner_coords = parameters["modelers"][1]["Parameters"]["lower_point_xyz"].GetVector()
93+
upper_corner_coords = parameters["modelers"][1]["Parameters"]["upper_point_xyz"].GetVector()
94+
95+
h_canditate = max(upper_corner_coords[0] - lower_corner_coords[0], upper_corner_coords[1] - lower_corner_coords[1]) / insertions
96+
h.append(h_canditate)
97+
98+
simulation.Clear()
99+
100+
print("\n h = ", h )
101+
print('\n L2 error', L2error_vector)
102+
103+
plt.xscale('log')
104+
plt.yscale('log')
105+
plt.grid(True, which="both", linestyle='--')
106+
107+
stored = np.zeros(6)
108+
stored[0] = (-(1)*(np.log(h[0])) + (np.log(L2error_vector[0])))
109+
stored[1] = (-(2)*(np.log(h[0])) + (np.log(L2error_vector[0])))
110+
stored[2] = (-(3)*(np.log(h[0])) + (np.log(L2error_vector[0])))
111+
stored[3] = (-(4)*(np.log(h[0])) + (np.log(L2error_vector[0])))
112+
stored[4] = (-(5)*(np.log(h[0])) + (np.log(L2error_vector[0])))
113+
stored[5] = (-(6)*(np.log(h[0])) + (np.log(L2error_vector[0])))
114+
degrees = np.arange(1, 7)
115+
yDegrees = np.zeros((degrees.size, len(h)))
116+
for i in range(0, degrees.size):
117+
for jtest in range(0, len(h)):
118+
yDegrees[i][jtest] = np.exp(degrees[i] * np.log(h[jtest]) + stored[i])
119+
plt.plot(h, yDegrees[i], "--", label='vel %d' % tuple([degrees[i]]))
120+
121+
plt.plot(h, L2error_vector, 's-', markersize=1.5, linewidth=2, label='L^2 error')
122+
plt.ylabel('L2 err',fontsize=14, color='blue')
123+
plt.xlabel('h',fontsize=14, color='blue')
124+
plt.legend(loc='lower right')
125+
plt.show()
126+

‎iga/use_cases/external_boundary_circle_with_nurbs/materials.json‎

Lines changed: 27 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,27 @@
1+
{
2+
"properties": [
3+
{
4+
"model_part_name": "IgaModelPart",
5+
"properties_id": 2,
6+
"Material": {
7+
"Variables": { "PENALTY_FACTOR": -1},
8+
"Tables": {}
9+
}
10+
},
11+
{
12+
"model_part_name": "IgaModelPart.StructuralAnalysisDomain",
13+
"properties_id": 2,
14+
"Material": {
15+
"name": "lin_el",
16+
"constitutive_law": { "name": "LinearElasticPlaneStress2DLaw" },
17+
"Variables": {
18+
"THICKNESS": 1.0,
19+
"YOUNG_MODULUS": 1000,
20+
"POISSON_RATIO": 0.3,
21+
"DENSITY": 1
22+
},
23+
"Tables": {}
24+
}
25+
}
26+
]
27+
}

0 commit comments

Comments
 (0)