Neumann BC for linear elastic plate example in fenics and kratos

.vscode/settings.json

@@ -0,0 +1,4 @@
+{
+    "python-envs.defaultEnvManager": "ms-python.python:conda",
+    "python-envs.defaultPackageManager": "ms-python.python:conda"
+}

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"]),
                     }})

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_

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:

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))

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",