Neumann BC for linear elastic plate example in fenics and kratos

.github/workflows/run-benchmark-examples.yml

@@ -14,8 +14,6 @@ on:
 
 env:
   CACHE_NUMBER: 1  # increase to reset cache manually
-  SNAKEMAKE_RESULT_FILE: metadata4ing_provenance
-  PROVENANACE_FILE_NAME: element_size_vs_max_mises_stress.pdf
 
 jobs:
   run-examples:
@@ -43,12 +41,25 @@ jobs:
           cd $GITHUB_WORKSPACE/examples/linear-elastic-plate-with-hole/fenics/
           python run_benchmark.py
 
+      - name: Run linear-elastic-plate-with-hole using kratos
+        shell: bash -l {0}
+        run: |
+          cd $GITHUB_WORKSPACE/examples/linear-elastic-plate-with-hole/kratos/
+          python run_benchmark.py
+
       - name: Archive results of the fenics run of linear-elastic-plate-with-hole
         uses: actions/upload-artifact@v4
         with:
-            name: snakemake_results_linear-elastic-plate-with-hole
+            name: fenics_results_linear-elastic-plate-with-hole
             path: |
               examples/linear-elastic-plate-with-hole/fenics/results/
 
+      - name: Archive results of the kratos run of linear-elastic-plate-with-hole
+        uses: actions/upload-artifact@v4
+        with:
+            name: kratos_results_linear-elastic-plate-with-hole
+            path: |
+              examples/linear-elastic-plate-with-hole/kratos/results/
+
 
 

environment_benchmarks.yml

@@ -11,4 +11,4 @@ dependencies:
   - conda
   - pip
   - pip:
-      - "git+https://github.com/izus-fokus/snakemake-report-plugin-metadata4ing@v1.2.5#egg=snakemake-report-plugin-metadata4ing"
+      - "git+https://github.com/izus-fokus/snakemake-report-plugin-metadata4ing@v1.2.6#egg=snakemake-report-plugin-metadata4ing"

examples/linear-elastic-plate-with-hole/fenics/run_benchmark.py

@@ -44,7 +44,7 @@
     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) 

examples/linear-elastic-plate-with-hole/fenics/run_fenics_simulation.py → examples/linear-elastic-plate-with-hole/fenics/run_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_

examples/linear-elastic-plate-with-hole/kratos/Snakefile_kratos.smk

@@ -0,0 +1,77 @@
+import json
+import os
+
+kratos_input_template = "input_template.json"
+kratos_material_template = "StructuralMaterials_template.json"
+
+rule all:
+    input:
+        "solution_metrics.json",
+        "solution_field_data.zip"
+
+rule mesh_to_mdpa:
+    input:
+        parameters = "parameters.json",
+        mesh = "mesh.msh",
+        script = "msh_to_mdpa.py",
+    output:
+        mdpa = "mesh.mdpa",
+    conda:
+        "environment_simulation.yml",
+    shell:
+        """
+        python3 {input.script} \
+            --input_parameter_file {input.parameters} \
+            --input_mesh_file {input.mesh} \
+            --output_mdpa_file {output.mdpa}
+        """
+
+rule create_kratos_input_and_run_simulation:
+    input:
+        parameters = "parameters.json",
+        mdpa = "mesh.mdpa",
+        kratos_input_template = kratos_input_template,
+        kratos_material_template = kratos_material_template,
+        script_create_kratos_input = "create_kratos_input.py",
+        script_run_kratos_simulation = "run_simulation.py",
+    output:
+        kratos_inputfile = "ProjectParameters.json",
+        kratos_materialfile = "MaterialParameters.json",
+        result_vtk = "vtk/Structure_0_1.vtk",
+    conda:
+        "environment_simulation.yml",
+    shell:
+        """
+        python3 {input.script_create_kratos_input} \
+            --input_parameter_file {input.parameters} \
+            --input_mdpa_file {input.mdpa} \
+            --input_kratos_input_template {input.kratos_input_template} \
+            --input_material_template {input.kratos_material_template} \
+            --output_kratos_inputfile {output.kratos_inputfile} \
+            --output_kratos_materialfile {output.kratos_materialfile}
+
+        python3 {input.script_run_kratos_simulation} \
+            --input_parameter_file {input.parameters} \
+            --input_kratos_inputfile {output.kratos_inputfile} \
+            --input_kratos_materialfile {output.kratos_materialfile} \
+        """
+
+rule postprocess_kratos_results:
+    input:
+        parameters = "parameters.json",
+        result_vtk = "vtk/Structure_0_1.vtk",
+        script = "postprocess_results.py",
+    output:
+        zip = "solution_field_data.zip",
+        metrics = "solution_metrics.json",
+    conda:
+        "environment_simulation.yml",
+    shell:
+        """
+        python3 {input.script} \
+            --input_parameter_file {input.parameters} \
+            --input_result_vtk {input.result_vtk} \
+            --output_solution_file_zip {output.zip} \
+            --output_metrics_file {output.metrics}
+        """
+

examples/linear-elastic-plate-with-hole/kratos/StructuralMaterials_template.json

@@ -0,0 +1,18 @@
+{
+    "properties": [
+        {
+            "model_part_name": "Structure",
+            "properties_id": 1,
+            "Material": {
+                "constitutive_law": {
+                    "name": "LinearElasticPlaneStress2DLaw"
+                },
+                "Variables": {
+                    "YOUNG_MODULUS": "{{YOUNG_MODULUS}}",
+                    "POISSON_RATIO": "{{POISSON_RATIO}}"
+                },
+                "Tables": {}
+            }
+        }
+    ]
+}