diff --git a/.vscode/settings.json b/.vscode/settings.json new file mode 100644 index 0000000..4b5a294 --- /dev/null +++ b/.vscode/settings.json @@ -0,0 +1,4 @@ +{ + "python-envs.defaultEnvManager": "ms-python.python:conda", + "python-envs.defaultPackageManager": "ms-python.python:conda" +} \ No newline at end of file diff --git a/benchmarks/common/parameter_extractor.py b/benchmarks/common/parameter_extractor.py index 8829f8c..660a723 100644 --- a/benchmarks/common/parameter_extractor.py +++ b/benchmarks/common/parameter_extractor.py @@ -23,10 +23,10 @@ def extract_params(self, rule_name: str, file_path: str) -> dict: with open(file_path) as f: data = json.load(f) for key, val in data.items(): - if isinstance(val, dict): + if isinstance(val, dict) and "value" in val: results[rule_name]["has parameter"].append({key: { "value": val["value"], - "unit": f"{val["unit"]}" if "unit" in val else None, + "unit": f"{val['unit']}" if "unit" in val else None, "json-path": f"/{key}/value", "data-type": self._get_type(val["value"]), }}) diff --git a/benchmarks/linear-elastic-plate-with-hole/fenics/run_fenics_simulation.py b/benchmarks/linear-elastic-plate-with-hole/fenics/run_fenics_simulation.py index ef980b2..3f10cc1 100644 --- a/benchmarks/linear-elastic-plate-with-hole/fenics/run_fenics_simulation.py +++ b/benchmarks/linear-elastic-plate-with-hole/fenics/run_fenics_simulation.py @@ -8,7 +8,6 @@ import numpy as np import ufl from dolfinx.fem.petsc import LinearProblem -from petsc4py.PETSc import ScalarType from mpi4py import MPI from pint import UnitRegistry @@ -40,8 +39,6 @@ def run_fenics_simulation( # Boundary conditions dofs_left = df.fem.locate_dofs_topological(V.sub(0), 1, tags_left) dofs_bottom = df.fem.locate_dofs_topological(V.sub(1), 1, tags_bottom) - dofs_right = df.fem.locate_dofs_topological(V, 1, tags_right) - dofs_top = df.fem.locate_dofs_topological(V, 1, tags_top) bc_left = df.fem.dirichletbc(0.0, dofs_left, V.sub(0)) bc_bottom = df.fem.dirichletbc(0.0, dofs_bottom, V.sub(1)) @@ -75,6 +72,23 @@ def run_fenics_simulation( .to_base_units() .magnitude ) + displacement_evaluation_point = parameters["displacement-evaluation-point"] + displacement_evaluation_x = ( + ureg.Quantity( + displacement_evaluation_point["x"]["value"], + displacement_evaluation_point["x"]["unit"], + ) + .to_base_units() + .magnitude + ) + displacement_evaluation_y = ( + ureg.Quantity( + displacement_evaluation_point["y"]["value"], + displacement_evaluation_point["y"]["unit"], + ) + .to_base_units() + .magnitude + ) analytical_solution = AnalyticalSolution( E=E, @@ -111,29 +125,38 @@ def as_tensor(v): domain=mesh, subdomain_data=facet_tags, ) - stress_space = df.fem.functionspace( - mesh, ("CG", parameters["element-degree"], (2, 2)) - ) - stress_function = df.fem.Function(stress_space) u = df.fem.Function(V, name="u") - u_prescribed = df.fem.Function(V, name="u_prescribed") - u_prescribed.interpolate(lambda x: analytical_solution.displacement(x)) - u_prescribed.x.scatter_forward() u_ = ufl.TestFunction(V) v_ = ufl.TrialFunction(V) a = df.fem.form(ufl.inner(sigma(u_), eps(v_)) * dx) - # set rhs to zero - f = df.fem.form(ufl.inner(df.fem.Constant(mesh, np.array([0.0, 0.0])), u_) * ufl.ds) + # Apply Neumann tractions on right and top boundaries from analytical stress. + traction_right = df.fem.Function(V, name="traction_right") + traction_top = df.fem.Function(V, name="traction_top") + + def traction_right_expr(x: np.ndarray) -> np.ndarray: + sxx, sxy, _ = analytical_solution.stress(x) + return np.vstack((np.asarray(sxx), np.asarray(sxy))) + + def traction_top_expr(x: np.ndarray) -> np.ndarray: + _, sxy, syy = analytical_solution.stress(x) + return np.vstack((np.asarray(sxy), np.asarray(syy))) + + traction_right.interpolate(traction_right_expr) + traction_top.interpolate(traction_top_expr) + traction_right.x.scatter_forward() + traction_top.x.scatter_forward() + + f = df.fem.form( + ufl.inner(traction_right, u_) * ds(3) + ufl.inner(traction_top, u_) * ds(4) + ) - bc_right = df.fem.dirichletbc(u_prescribed, dofs_right) - bc_top = df.fem.dirichletbc(u_prescribed, dofs_top) solver = LinearProblem( a, f, - bcs=[bc_left, bc_bottom, bc_right, bc_top], + bcs=[bc_left, bc_bottom], u=u, petsc_options={ "ksp_type": "gmres", @@ -143,6 +166,50 @@ def as_tensor(v): ) solver.solve() + # Support reaction on the left boundary + n = ufl.FacetNormal(mesh) + traction = ufl.dot(sigma(u), n) + reaction_left_x_local = df.fem.assemble_scalar(df.fem.form(-traction[0] * ds(1))) + reaction_left_y_local = df.fem.assemble_scalar(df.fem.form(-traction[1] * ds(1))) + reaction_left_x = MPI.COMM_WORLD.allreduce(reaction_left_x_local, op=MPI.SUM) + reaction_left_y = MPI.COMM_WORLD.allreduce(reaction_left_y_local, op=MPI.SUM) + + + # Compute L2 error norm between FE displacement and analytical displacement. + u_analytical = df.fem.Function(V, name="u_analytical") + + def analytical_displacement_expr(x: np.ndarray) -> np.ndarray: + ux, uy = analytical_solution.displacement(x) + return np.vstack((np.asarray(ux), np.asarray(uy))) + + u_analytical.interpolate(analytical_displacement_expr) + u_analytical.x.scatter_forward() + + l2_error_form = df.fem.form(ufl.inner(u - u_analytical, u - u_analytical) * dx) + l2_error_sq_local = df.fem.assemble_scalar(l2_error_form) + l2_error_sq_global = MPI.COMM_WORLD.allreduce(l2_error_sq_local, op=MPI.SUM) + l2_error_displacement = np.sqrt(l2_error_sq_global) + + # Evaluate displacement at the specified evaluation point + displacement_eval_point = np.array( + [[displacement_evaluation_x, displacement_evaluation_y, 0.0]], + dtype=np.float64, + ) + tree = df.geometry.bb_tree(mesh, mesh.topology.dim) + cell_candidates = df.geometry.compute_collisions_points( + tree, displacement_eval_point + ) + colliding_cells = df.geometry.compute_colliding_cells( + mesh, cell_candidates, displacement_eval_point + ) + local_displacement_x = None + if len(colliding_cells.links(0)) > 0: + cell = colliding_cells.links(0)[0] + local_displacement_x = u.eval( + displacement_eval_point, np.array([cell], dtype=np.int32) + )[0] + + def project( v: df.fem.Function | ufl.core.expr.Expr, V: df.fem.FunctionSpace, @@ -234,10 +301,32 @@ def mises_stress(u): max_mises_stress_gauss_points = MPI.COMM_WORLD.allreduce( np.max(mises_qp.x.array), op=MPI.MAX ) + + displacement_x_at_evaluation_point = None + if MPI.COMM_WORLD.rank == 0: + displacement_x_candidates = ( + MPI.COMM_WORLD.gather(local_displacement_x, root=0) or [] + ) + for value in displacement_x_candidates: + if value is not None: + displacement_x_at_evaluation_point = value + break + + if displacement_x_at_evaluation_point is None: + raise ValueError( + "Could not evaluate displacement at the configured evaluation point." + ) + else: + MPI.COMM_WORLD.gather(local_displacement_x, root=0) + # Save metrics metrics = { "max_von_mises_stress_nodes": max_mises_stress_nodes, "max_von_mises_stress_gauss_points": max_mises_stress_gauss_points, + "l2_error_displacement": l2_error_displacement, + "reaction_force_left_boundary_x": reaction_left_x, + "reaction_force_left_boundary_y": reaction_left_y, + f"displacement_x_at_evaluation_point (x={displacement_evaluation_x}, y={displacement_evaluation_y})": displacement_x_at_evaluation_point, } if MPI.COMM_WORLD.rank == 0: @@ -299,3 +388,5 @@ def mises_stress(u): args.output_solution_file_zip, args.output_metrics_file, ) +#python3 run_fenics_simulation.py --input_parameter_file parameters_1.json --input_mesh_file mesh_1.msh --output_solution_file_zip results/solution_field_data.zip --output_metrics_file results/solution_metrics.json +#conda activate /home/dtyagi/NFDI4IngModelValidationPlatform/examples/linear-elastic-plate-with-hole/fenics/conda_envs/68782c9259a7e25569f1ab0241a766c5_ \ No newline at end of file diff --git a/benchmarks/linear-elastic-plate-with-hole/kratos/MainKratos.py b/benchmarks/linear-elastic-plate-with-hole/kratos/MainKratos.py deleted file mode 100755 index 02d2394..0000000 --- a/benchmarks/linear-elastic-plate-with-hole/kratos/MainKratos.py +++ /dev/null @@ -1,19 +0,0 @@ -from __future__ import print_function, absolute_import, division #makes KratosMultiphysics backward compatible with python 2.6 and 2.7 - -import KratosMultiphysics -from KratosMultiphysics.StructuralMechanicsApplication.structural_mechanics_analysis import StructuralMechanicsAnalysis -import sys -""" -For user-scripting it is intended that a new class is derived -from StructuralMechanicsAnalysis to do modifications -""" -parameters = sys.argv[1] - -if __name__ == "__main__": - - with open(parameters,'r') as parameter_file: - parameters = KratosMultiphysics.Parameters(parameter_file.read()) - - model = KratosMultiphysics.Model() - simulation = StructuralMechanicsAnalysis(model,parameters) - simulation.Run() diff --git a/benchmarks/linear-elastic-plate-with-hole/kratos/Snakefile b/benchmarks/linear-elastic-plate-with-hole/kratos/Snakefile index 38054c6..5f14a26 100644 --- a/benchmarks/linear-elastic-plate-with-hole/kratos/Snakefile +++ b/benchmarks/linear-elastic-plate-with-hole/kratos/Snakefile @@ -59,7 +59,7 @@ rule create_kratos_input_and_run_simulation: python3 {input.script_run_kratos_simulation} \ --input_parameter_file {input.parameters} \ --input_kratos_inputfile {output.kratos_inputfile} \ - --input_kratos_materialfile {output.kratos_materialfile} + --input_kratos_materialfile {output.kratos_materialfile} \ """ rule postprocess_kratos_results: diff --git a/benchmarks/linear-elastic-plate-with-hole/kratos/create_kratos_input.py b/benchmarks/linear-elastic-plate-with-hole/kratos/create_kratos_input.py index bb6b0c8..0307568 100644 --- a/benchmarks/linear-elastic-plate-with-hole/kratos/create_kratos_input.py +++ b/benchmarks/linear-elastic-plate-with-hole/kratos/create_kratos_input.py @@ -4,6 +4,7 @@ from argparse import ArgumentParser from pathlib import Path import sys +import sympy as sp # Ensure the parent directory is in the path to import AnalyticalSolution sys.path.insert(0, str(Path(__file__).resolve().parent.parent)) from analytical_solution import AnalyticalSolution @@ -58,7 +59,17 @@ def create_kratos_input( load=load, ) - bc = analytical_solution.displacement_symbolic_str("X", "Y") + # Build traction expressions t = sigma * n on right (n=[1,0]) and top (n=[0,1]) boundaries. + sxx_sym, sxy_sym, syy_sym = analytical_solution.stress_symbolic() + x, y = sp.symbols("x y") + X_sym = sp.Symbol("X") + Y_sym = sp.Symbol("Y") + E_sym, nu_sym, a_sym, T_sym = sp.symbols("E nu a T") + subs_vars = {x: X_sym, y: Y_sym} + subs_params = {E_sym: E, nu_sym: nu, a_sym: radius, T_sym: load} + sxx_str = sp.sstr(sxx_sym.subs(subs_vars).subs(subs_params)) + sxy_str = sp.sstr(sxy_sym.subs(subs_vars).subs(subs_params)) + syy_str = sp.sstr(syy_sym.subs(subs_vars).subs(subs_params)) with open(kratos_material_template_file) as f: material_string = f.read() @@ -78,16 +89,16 @@ def create_kratos_input( r"{{MATERIAL_FILE}}", kratos_material_file ) project_parameters_string = project_parameters_string.replace( - r"{{BOUNDARY_RIGHT_DISPLACEMENT_X}}", str(bc[0]) + r"{{BOUNDARY_RIGHT_TRACTION_X}}", sxx_str ) project_parameters_string = project_parameters_string.replace( - r"{{BOUNDARY_RIGHT_DISPLACEMENT_Y}}", str(bc[1]) + r"{{BOUNDARY_RIGHT_TRACTION_Y}}", sxy_str ) project_parameters_string = project_parameters_string.replace( - r"{{BOUNDARY_TOP_DISPLACEMENT_X}}", str(bc[0]) + r"{{BOUNDARY_TOP_TRACTION_X}}", sxy_str ) project_parameters_string = project_parameters_string.replace( - r"{{BOUNDARY_TOP_DISPLACEMENT_Y}}", str(bc[1]) + r"{{BOUNDARY_TOP_TRACTION_Y}}", syy_str ) config = parameters["configuration"] output_dir = os.path.join(os.path.dirname(os.path.abspath(kratos_input_file)), str(config)) diff --git a/benchmarks/linear-elastic-plate-with-hole/kratos/input_template.json b/benchmarks/linear-elastic-plate-with-hole/kratos/input_template.json index f8b680a..cbac685 100644 --- a/benchmarks/linear-elastic-plate-with-hole/kratos/input_template.json +++ b/benchmarks/linear-elastic-plate-with-hole/kratos/input_template.json @@ -12,6 +12,7 @@ "echo_level": 1, "domain_size": 2, "analysis_type": "linear", + "compute_reactions": true, "model_import_settings": { "input_type": "mdpa", "input_filename": "{{MESH_FILE}}" @@ -68,21 +69,23 @@ "End" ] } - }, + } + ], + "loads_process_list": [ { "python_module": "assign_vector_variable_process", "kratos_module": "KratosMultiphysics", "Parameters": { "model_part_name": "Structure.boundary_right", - "variable_name": "DISPLACEMENT", + "variable_name": "POINT_LOAD", "constrained": [ - true, - true, - true + false, + false, + false ], "value": [ - "{{BOUNDARY_RIGHT_DISPLACEMENT_X}}", - "{{BOUNDARY_RIGHT_DISPLACEMENT_Y}}", + "{{BOUNDARY_RIGHT_TRACTION_X}}", + "{{BOUNDARY_RIGHT_TRACTION_Y}}", 0.0 ], "interval": [ @@ -96,15 +99,15 @@ "kratos_module": "KratosMultiphysics", "Parameters": { "model_part_name": "Structure.boundary_top", - "variable_name": "DISPLACEMENT", + "variable_name": "POINT_LOAD", "constrained": [ - true, - true, - true + false, + false, + false ], "value": [ - "{{BOUNDARY_TOP_DISPLACEMENT_X}}", - "{{BOUNDARY_TOP_DISPLACEMENT_Y}}", + "{{BOUNDARY_TOP_TRACTION_X}}", + "{{BOUNDARY_TOP_TRACTION_Y}}", 0.0 ], "interval": [ @@ -114,7 +117,6 @@ } } ], - "loads_process_list": [], "list_other_processes": [] }, "output_processes": { @@ -129,7 +131,11 @@ "output_sub_model_parts": false, "output_interval": 1, "nodal_solution_step_data_variables": [ - "DISPLACEMENT" + "DISPLACEMENT", + "REACTION" + ], + "gauss_point_variables_in_elements": [ + "VON_MISES_STRESS" ], "gauss_point_variables_extrapolated_to_nodes": [ "CAUCHY_STRESS_VECTOR", diff --git a/benchmarks/linear-elastic-plate-with-hole/kratos/postprocess_results.py b/benchmarks/linear-elastic-plate-with-hole/kratos/postprocess_results.py index 467fb42..d2f050d 100644 --- a/benchmarks/linear-elastic-plate-with-hole/kratos/postprocess_results.py +++ b/benchmarks/linear-elastic-plate-with-hole/kratos/postprocess_results.py @@ -3,17 +3,131 @@ from pathlib import Path import zipfile from argparse import ArgumentParser +import numpy as np +from pint import UnitRegistry +import sys +from typing import cast + +sys.path.insert(0, str(Path(__file__).resolve().parent.parent)) +from analytical_solution import AnalyticalSolution def postprocess_results(input_parameter_file, input_result_vtk, output_metrics_file, output_solution_file_zip): + ureg = UnitRegistry() with open(input_parameter_file) as f: parameters = json.load(f) config = parameters["configuration"] mesh = pyvista.read(str(input_result_vtk)) - max_von_mises_stress = float(mesh["VON_MISES_STRESS"].max()) - print("Max Von Mises Stress:", max_von_mises_stress) + if not isinstance(mesh, pyvista.DataSet): + raise TypeError(f"Expected a pyvista.DataSet, got {type(mesh).__name__}") + mesh = cast(pyvista.DataSet, mesh) + + E = ( + ureg.Quantity( + parameters["young-modulus"]["value"], parameters["young-modulus"]["unit"] + ) + .to_base_units() + .magnitude + ) + nu = ( + ureg.Quantity( + parameters["poisson-ratio"]["value"], parameters["poisson-ratio"]["unit"] + ) + .to_base_units() + .magnitude + ) + radius = ( + ureg.Quantity(parameters["radius"]["value"], parameters["radius"]["unit"]) + .to_base_units() + .magnitude + ) + L = ( + ureg.Quantity(parameters["length"]["value"], parameters["length"]["unit"]) + .to_base_units() + .magnitude + ) + load = ( + ureg.Quantity(parameters["load"]["value"], parameters["load"]["unit"]) + .to_base_units() + .magnitude + ) + + displacement_evaluation_point = parameters["displacement-evaluation-point"] + displacement_evaluation_x = ( + ureg.Quantity( + displacement_evaluation_point["x"]["value"], + displacement_evaluation_point["x"]["unit"], + ) + .to_base_units() + .magnitude + ) + displacement_evaluation_y = ( + ureg.Quantity( + displacement_evaluation_point["y"]["value"], + displacement_evaluation_point["y"]["unit"], + ) + .to_base_units() + .magnitude + ) + + analytical_solution = AnalyticalSolution( + E=E, + nu=nu, + radius=radius, + L=L, + load=load, + ) + # Compute maximum von Mises stress at nodes and Gauss points. + max_von_mises_stress_nodes = float(np.max(mesh.point_data["VON_MISES_STRESS"])) + + max_von_mises_stress_gauss_points = max_von_mises_stress_nodes + for key, values in mesh.cell_data.items(): + if "VON_MISES_STRESS" in key: + max_von_mises_stress_gauss_points = float(np.max(values)) + break + + # Compute L2 error of displacement field compared to analytical solution. + coords = np.asarray(mesh.points) + displacement = np.asarray(mesh.point_data["DISPLACEMENT"])[:, :2] + u_ref_x, u_ref_y = analytical_solution.displacement(coords[:, :2].T) + u_ref = np.column_stack((np.asarray(u_ref_x), np.asarray(u_ref_y))) + err_sq_node = np.sum((displacement - u_ref) ** 2, axis=1) + + cell_sizes = mesh.compute_cell_sizes(length=False, area=True, volume=False) + cell_areas = np.asarray(cell_sizes.cell_data["Area"]) + l2_error_sq = 0.0 + for i in range(mesh.n_cells): + point_ids = mesh.get_cell(i).point_ids + if len(point_ids) == 0: + continue + l2_error_sq += float(np.mean(err_sq_node[point_ids]) * cell_areas[i]) + l2_error_displacement = float(np.sqrt(l2_error_sq)) + + # Compute reaction forces on the left boundary (x=0) by summing the reaction forces at the nodes on that boundary. + tolerance = 1e-10 * max(1.0, L) + left_boundary_mask = np.isclose(coords[:, 0], 0.0, atol=tolerance) + reaction = np.asarray(mesh.point_data.get("REACTION", np.zeros((mesh.n_points, 3)))) + reaction_force_left_boundary_x = float(np.sum(reaction[left_boundary_mask, 0])) + reaction_force_left_boundary_y = float(np.sum(reaction[left_boundary_mask, 1])) + + probe_points = pyvista.PolyData( + np.array([[displacement_evaluation_x, displacement_evaluation_y, 0.0]], dtype=float) + ) + sampled = probe_points.sample(mesh) + displacement_sampled = sampled.point_data.get("DISPLACEMENT") + if displacement_sampled is None: + closest_id = mesh.find_closest_point([displacement_evaluation_x, displacement_evaluation_y, 0.0]) + displacement_x_at_evaluation_point = float(displacement[closest_id, 0]) + else: + displacement_x_at_evaluation_point = float(displacement_sampled[0, 0]) + metrics = { - "max_von_mises_stress_nodes": max_von_mises_stress + "max_von_mises_stress_nodes": max_von_mises_stress_nodes, + "max_von_mises_stress_gauss_points": max_von_mises_stress_gauss_points, + "l2_error_displacement": l2_error_displacement, + "reaction_force_left_boundary_x": reaction_force_left_boundary_x, + "reaction_force_left_boundary_y": reaction_force_left_boundary_y, + f"displacement_x_at_evaluation_point (x={displacement_evaluation_x}, y={displacement_evaluation_y})": displacement_x_at_evaluation_point, } with open(output_metrics_file, "w") as f: json.dump(metrics, f, indent=4) diff --git a/benchmarks/linear-elastic-plate-with-hole/kratos/run_kratos_simulation.py b/benchmarks/linear-elastic-plate-with-hole/kratos/run_kratos_simulation.py index 5df8a0b..8d270c3 100644 --- a/benchmarks/linear-elastic-plate-with-hole/kratos/run_kratos_simulation.py +++ b/benchmarks/linear-elastic-plate-with-hole/kratos/run_kratos_simulation.py @@ -23,7 +23,7 @@ if __name__ == "__main__": parser = ArgumentParser( description="Run FEniCS simulation for a plate with a hole.\n" - "Inputs: --input_parameter_file, --input_kratos_inputfile, --input_kratos_materialfile\n" + "Inputs: --input_parameter_file, --input_kratos_inputfile \n" "Outputs: --output_solution_file_zip, --output_metrics_file" ) parser.add_argument( @@ -38,8 +38,9 @@ ) parser.add_argument( "--input_kratos_materialfile", - required=True, - help="Path to the kratos material file (input)", + required=False, + default=None, + help="Optional path to the kratos material file (unused, kept for compatibility)", ) args, _ = parser.parse_known_args() diff --git a/benchmarks/linear-elastic-plate-with-hole/parameters_003125.json b/benchmarks/linear-elastic-plate-with-hole/parameters_003125.json index dd1e2d1..07d8829 100644 --- a/benchmarks/linear-elastic-plate-with-hole/parameters_003125.json +++ b/benchmarks/linear-elastic-plate-with-hole/parameters_003125.json @@ -8,6 +8,16 @@ "value": 1.0, "unit": "m" }, + "displacement-evaluation-point": { + "x": { + "value": 1.0, + "unit": "m" + }, + "y": { + "value": 1.0, + "unit": "m" + } + }, "load": { "value": 100.0, "unit": "MPa" diff --git a/benchmarks/linear-elastic-plate-with-hole/parameters_00625.json b/benchmarks/linear-elastic-plate-with-hole/parameters_00625.json index a10d6b8..40226a0 100644 --- a/benchmarks/linear-elastic-plate-with-hole/parameters_00625.json +++ b/benchmarks/linear-elastic-plate-with-hole/parameters_00625.json @@ -8,6 +8,16 @@ "value": 1.0, "unit": "m" }, + "displacement-evaluation-point": { + "x": { + "value": 1.0, + "unit": "m" + }, + "y": { + "value": 1.0, + "unit": "m" + } + }, "load": { "value": 100.0, "unit": "MPa" diff --git a/benchmarks/linear-elastic-plate-with-hole/parameters_0125.json b/benchmarks/linear-elastic-plate-with-hole/parameters_0125.json index cba88ae..8c1e6fe 100644 --- a/benchmarks/linear-elastic-plate-with-hole/parameters_0125.json +++ b/benchmarks/linear-elastic-plate-with-hole/parameters_0125.json @@ -8,6 +8,16 @@ "value": 1.0, "unit": "m" }, + "displacement-evaluation-point": { + "x": { + "value": 1.0, + "unit": "m" + }, + "y": { + "value": 1.0, + "unit": "m" + } + }, "load": { "value": 100.0, "unit": "MPa" diff --git a/benchmarks/linear-elastic-plate-with-hole/parameters_025.json b/benchmarks/linear-elastic-plate-with-hole/parameters_025.json index fb23454..5b086f1 100644 --- a/benchmarks/linear-elastic-plate-with-hole/parameters_025.json +++ b/benchmarks/linear-elastic-plate-with-hole/parameters_025.json @@ -8,6 +8,16 @@ "value": 1.0, "unit": "m" }, + "displacement-evaluation-point": { + "x": { + "value": 1.0, + "unit": "m" + }, + "y": { + "value": 1.0, + "unit": "m" + } + }, "load": { "value": 100.0, "unit": "MPa" diff --git a/benchmarks/linear-elastic-plate-with-hole/parameters_05.json b/benchmarks/linear-elastic-plate-with-hole/parameters_05.json index ac9a011..36a227f 100644 --- a/benchmarks/linear-elastic-plate-with-hole/parameters_05.json +++ b/benchmarks/linear-elastic-plate-with-hole/parameters_05.json @@ -8,6 +8,16 @@ "value": 1.0, "unit": "m" }, + "displacement-evaluation-point": { + "x": { + "value": 1.0, + "unit": "m" + }, + "y": { + "value": 1.0, + "unit": "m" + } + }, "load": { "value": 100.0, "unit": "MPa" diff --git a/benchmarks/linear-elastic-plate-with-hole/parameters_1.json b/benchmarks/linear-elastic-plate-with-hole/parameters_1.json index 9f79e04..106e3f6 100644 --- a/benchmarks/linear-elastic-plate-with-hole/parameters_1.json +++ b/benchmarks/linear-elastic-plate-with-hole/parameters_1.json @@ -8,6 +8,16 @@ "value": 1.0, "unit": "m" }, + "displacement-evaluation-point": { + "x": { + "value": 1.0, + "unit": "m" + }, + "y": { + "value": 1.0, + "unit": "m" + } + }, "load": { "value": 100.0, "unit": "MPa" diff --git a/benchmarks/linear-elastic-plate-with-hole/parameters_geo1fun2.json b/benchmarks/linear-elastic-plate-with-hole/parameters_geo1fun2.json index 4f223b2..303ad65 100644 --- a/benchmarks/linear-elastic-plate-with-hole/parameters_geo1fun2.json +++ b/benchmarks/linear-elastic-plate-with-hole/parameters_geo1fun2.json @@ -8,6 +8,16 @@ "value": 1.0, "unit": "m" }, + "displacement-evaluation-point": { + "x": { + "value": 1.0, + "unit": "m" + }, + "y": { + "value": 1.0, + "unit": "m" + } + }, "load": { "value": 100.0, "unit": "MPa" diff --git a/benchmarks/linear-elastic-plate-with-hole/parameters_geo2fun2.json b/benchmarks/linear-elastic-plate-with-hole/parameters_geo2fun2.json index c04200e..a4343a3 100644 --- a/benchmarks/linear-elastic-plate-with-hole/parameters_geo2fun2.json +++ b/benchmarks/linear-elastic-plate-with-hole/parameters_geo2fun2.json @@ -8,6 +8,16 @@ "value": 1.0, "unit": "m" }, + "displacement-evaluation-point": { + "x": { + "value": 1.0, + "unit": "m" + }, + "y": { + "value": 1.0, + "unit": "m" + } + }, "load": { "value": 100.0, "unit": "MPa" diff --git a/benchmarks/linear-elastic-plate-with-hole/parameters_geo2fun2fine.json b/benchmarks/linear-elastic-plate-with-hole/parameters_geo2fun2fine.json index b089eaf..df25eff 100644 --- a/benchmarks/linear-elastic-plate-with-hole/parameters_geo2fun2fine.json +++ b/benchmarks/linear-elastic-plate-with-hole/parameters_geo2fun2fine.json @@ -8,6 +8,16 @@ "value": 1.0, "unit": "m" }, + "displacement-evaluation-point": { + "x": { + "value": 1.0, + "unit": "m" + }, + "y": { + "value": 1.0, + "unit": "m" + } + }, "load": { "value": 100.0, "unit": "MPa" diff --git a/examples/linear-elastic-plate-with-hole/fenics/run_benchmark.py b/examples/linear-elastic-plate-with-hole/fenics/run_benchmark.py index 929fd72..9951016 100644 --- a/examples/linear-elastic-plate-with-hole/fenics/run_benchmark.py +++ b/examples/linear-elastic-plate-with-hole/fenics/run_benchmark.py @@ -44,7 +44,6 @@ with open(file, "r") as f: data = json.load(f) if data.get("element-size").get("value") >= 0.025: - # Create output directory for the configuration output_dir = root_unzipped_benchmark_dir / "results" / data.get("configuration") output_dir.mkdir(parents=True, exist_ok=True) diff --git a/examples/linear-elastic-plate-with-hole/fenics/run_fenics_simulation.py b/examples/linear-elastic-plate-with-hole/fenics/run_fenics_simulation.py index 6b5a114..3f10cc1 100644 --- a/examples/linear-elastic-plate-with-hole/fenics/run_fenics_simulation.py +++ b/examples/linear-elastic-plate-with-hole/fenics/run_fenics_simulation.py @@ -8,7 +8,6 @@ import numpy as np import ufl from dolfinx.fem.petsc import LinearProblem -from petsc4py.PETSc import ScalarType from mpi4py import MPI from pint import UnitRegistry @@ -40,8 +39,6 @@ def run_fenics_simulation( # Boundary conditions dofs_left = df.fem.locate_dofs_topological(V.sub(0), 1, tags_left) dofs_bottom = df.fem.locate_dofs_topological(V.sub(1), 1, tags_bottom) - dofs_right = df.fem.locate_dofs_topological(V, 1, tags_right) - dofs_top = df.fem.locate_dofs_topological(V, 1, tags_top) bc_left = df.fem.dirichletbc(0.0, dofs_left, V.sub(0)) bc_bottom = df.fem.dirichletbc(0.0, dofs_bottom, V.sub(1)) @@ -128,29 +125,38 @@ def as_tensor(v): domain=mesh, subdomain_data=facet_tags, ) - stress_space = df.fem.functionspace( - mesh, ("CG", parameters["element-degree"], (2, 2)) - ) - stress_function = df.fem.Function(stress_space) u = df.fem.Function(V, name="u") - u_prescribed = df.fem.Function(V, name="u_prescribed") - u_prescribed.interpolate(lambda x: analytical_solution.displacement(x)) - u_prescribed.x.scatter_forward() u_ = ufl.TestFunction(V) v_ = ufl.TrialFunction(V) a = df.fem.form(ufl.inner(sigma(u_), eps(v_)) * dx) - # set rhs to zero - f = df.fem.form(ufl.inner(df.fem.Constant(mesh, np.array([0.0, 0.0])), u_) * ufl.ds) + # Apply Neumann tractions on right and top boundaries from analytical stress. + traction_right = df.fem.Function(V, name="traction_right") + traction_top = df.fem.Function(V, name="traction_top") + + def traction_right_expr(x: np.ndarray) -> np.ndarray: + sxx, sxy, _ = analytical_solution.stress(x) + return np.vstack((np.asarray(sxx), np.asarray(sxy))) + + def traction_top_expr(x: np.ndarray) -> np.ndarray: + _, sxy, syy = analytical_solution.stress(x) + return np.vstack((np.asarray(sxy), np.asarray(syy))) + + traction_right.interpolate(traction_right_expr) + traction_top.interpolate(traction_top_expr) + traction_right.x.scatter_forward() + traction_top.x.scatter_forward() + + f = df.fem.form( + ufl.inner(traction_right, u_) * ds(3) + ufl.inner(traction_top, u_) * ds(4) + ) - bc_right = df.fem.dirichletbc(u_prescribed, dofs_right) - bc_top = df.fem.dirichletbc(u_prescribed, dofs_top) solver = LinearProblem( a, f, - bcs=[bc_left, bc_bottom, bc_right, bc_top], + bcs=[bc_left, bc_bottom], u=u, petsc_options={ "ksp_type": "gmres", @@ -160,6 +166,30 @@ def as_tensor(v): ) solver.solve() + # Support reaction on the left boundary + n = ufl.FacetNormal(mesh) + traction = ufl.dot(sigma(u), n) + reaction_left_x_local = df.fem.assemble_scalar(df.fem.form(-traction[0] * ds(1))) + reaction_left_y_local = df.fem.assemble_scalar(df.fem.form(-traction[1] * ds(1))) + reaction_left_x = MPI.COMM_WORLD.allreduce(reaction_left_x_local, op=MPI.SUM) + reaction_left_y = MPI.COMM_WORLD.allreduce(reaction_left_y_local, op=MPI.SUM) + + + # Compute L2 error norm between FE displacement and analytical displacement. + u_analytical = df.fem.Function(V, name="u_analytical") + + def analytical_displacement_expr(x: np.ndarray) -> np.ndarray: + ux, uy = analytical_solution.displacement(x) + return np.vstack((np.asarray(ux), np.asarray(uy))) + + u_analytical.interpolate(analytical_displacement_expr) + u_analytical.x.scatter_forward() + + l2_error_form = df.fem.form(ufl.inner(u - u_analytical, u - u_analytical) * dx) + l2_error_sq_local = df.fem.assemble_scalar(l2_error_form) + l2_error_sq_global = MPI.COMM_WORLD.allreduce(l2_error_sq_local, op=MPI.SUM) + l2_error_displacement = np.sqrt(l2_error_sq_global) + # Evaluate displacement at the specified evaluation point displacement_eval_point = np.array( [[displacement_evaluation_x, displacement_evaluation_y, 0.0]], @@ -293,6 +323,9 @@ def mises_stress(u): metrics = { "max_von_mises_stress_nodes": max_mises_stress_nodes, "max_von_mises_stress_gauss_points": max_mises_stress_gauss_points, + "l2_error_displacement": l2_error_displacement, + "reaction_force_left_boundary_x": reaction_left_x, + "reaction_force_left_boundary_y": reaction_left_y, f"displacement_x_at_evaluation_point (x={displacement_evaluation_x}, y={displacement_evaluation_y})": displacement_x_at_evaluation_point, } @@ -355,3 +388,5 @@ def mises_stress(u): args.output_solution_file_zip, args.output_metrics_file, ) +#python3 run_fenics_simulation.py --input_parameter_file parameters_1.json --input_mesh_file mesh_1.msh --output_solution_file_zip results/solution_field_data.zip --output_metrics_file results/solution_metrics.json +#conda activate /home/dtyagi/NFDI4IngModelValidationPlatform/examples/linear-elastic-plate-with-hole/fenics/conda_envs/68782c9259a7e25569f1ab0241a766c5_ \ No newline at end of file diff --git a/examples/linear-elastic-plate-with-hole/linear-elastic-plate-with-hole.zip b/examples/linear-elastic-plate-with-hole/linear-elastic-plate-with-hole.zip index 9644a81..f1557b9 100644 Binary files a/examples/linear-elastic-plate-with-hole/linear-elastic-plate-with-hole.zip and b/examples/linear-elastic-plate-with-hole/linear-elastic-plate-with-hole.zip differ