new solution metrics added

Author: div-tyg

docs/benchmarks/linear elasticity/plate-with-hole.md

@@ -96,31 +96,40 @@ $$
 
 In order to solve the weak formulation, the finite-element method (FEM) can be used. This method discretizes the domain $\Omega$ into so called finite elements that can for example be triangles or quadrilaterals in 2D. On these elements, ansatz functions are defined such that they are continous on the boundaries between elements. These functions form a basis for the solution space for an approximate solution $\boldsymbol{u}_h$ of the problem.
 
-## Comparison of approximate solution with analytical solution
+## Solution metrics
 
-The approxiamte solution and the analytical solution can be compared with the $L_2$ norm which is defined as
+The approximate solution is compared with the analytical solution using the $L_2$ norm which is defined as
 
 $$
-\Vert \boldsymbol{u}\Vert_{L_2} = \sqrt{\int_\Omega \Vert\boldsymbol{u}(\boldsymbol{x})\Vert_2^2 \mathrm d \boldsymbol x}
+\Vert \boldsymbol{f}\Vert_{L_2} = \sqrt{\int_\Omega \left|\boldsymbol{f}(\boldsymbol{x})\right|^2 \mathrm d \boldsymbol x}
 $$
 
-and the error in the $L_2$ norm
+This norm is computed for the error between the finite element solution and the analytical solution of displacements i.e., $\Vert \boldsymbol{u}-\boldsymbol{u}_h\Vert_{L_2}$.
+
+The maximum displacement error is computed at the nodes of the finite element mesh with respect to the analytical solution.
 
 $$
-e_{L_2} = \Vert \boldsymbol{u}-\boldsymbol{u}_h\Vert_{L_2}.
+e_{\max}
+=
+\max_{i \in \mathcal{N}}
+\left\|
+\mathbf{u}(\mathbf{x}_i)
+-
+\mathbf{u}_h(\mathbf{x}_i)
+\right\|
 $$
 
-Alternatively, the sup norm can be used which is defined as
+where $\mathcal{N}$ denotes the set of nodes of the finite element mesh.
 
-$$
-\Vert \boldsymbol{u}\Vert_{\inf} = \sup_{\boldsymbol x} \Vert \boldsymbol{u}(\boldsymbol{x}) \Vert
-$$
+Max Von-Mises stress
+
+Displacement at the top-right corner of the plate.
+
+The reaction force at the left boundary.
+
+Number of Degrees of Freedom which supports in comparing the outputs between different meshes.
 
-and the error in the sup norm
 
-$$
-e_{\inf} = \Vert \boldsymbol{u}-\boldsymbol{u}_h\Vert_{\inf}.
-$$
 
 
 With these metrices, we can perform a convergence analysis for different approximations $\boldsymbol{u}_h$ which differ in the element size $h$. Plotting the error over the used element-size in a log-log plot lets us determine the convergence order of the approximation.

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

@@ -170,6 +170,7 @@ def traction_top_expr(x: np.ndarray) -> np.ndarray:
     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)
+    num_dofs = V.dofmap.index_map.size_global * V.dofmap.index_map_bs
 
 
     # Compute L2 error norm between FE displacement and analytical displacement.
@@ -187,6 +188,14 @@ def analytical_displacement_expr(x: np.ndarray) -> np.ndarray:
     l2_error_sq_global = MPI.COMM_WORLD.allreduce(l2_error_sq_local, op=MPI.SUM)
     l2_error_displacement = np.sqrt(l2_error_sq_global)
 
+    # Compute max nodal displacement error magnitude (global across MPI)
+    block_size = V.dofmap.index_map_bs
+    nodal_error = (u.x.array - u_analytical.x.array).reshape(-1, block_size)
+    max_displacement_error_nodes_local = np.max(np.linalg.norm(nodal_error, axis=1))
+    max_displacement_error_nodes = MPI.COMM_WORLD.allreduce(
+        max_displacement_error_nodes_local, op=MPI.MAX
+    )
+
     # Evaluate displacement at the specified evaluation point
     displacement_eval_point = np.array(
         [[displacement_evaluation_x, displacement_evaluation_y, 0.0]],
@@ -281,8 +290,6 @@ def mises_stress(u):
     ) as vtk:
         vtk.write_function(mises_stress_nodes, 0.0)
 
-    # extract maximum von Mises stress
-    max_mises_stress_nodes = np.max(mises_stress_nodes.x.array)
 
     # Compute von Mises stress at quadrature (Gauss) points and extract maximum (global across MPI)
     quad_element = basix.ufl.quadrature_element(
@@ -319,9 +326,10 @@ def mises_stress(u):
 
     # Save metrics
     metrics = {
-        "max_von_mises_stress_nodes[Pa]": max_mises_stress_nodes,
-        "max_von_mises_stress_gauss_points[Pa]": max_mises_stress_gauss_points,
+        "number_of_dofs[-]": num_dofs,
+        "max_von_mises_stress[Pa]": max_mises_stress_gauss_points,
         "l2_error_displacement[m]": l2_error_displacement,
+        "max_displacement_error[m]": max_displacement_error_nodes,
         "reaction_force_left_boundary_x[N]": reaction_left_x,
         "reaction_force_left_boundary_y[N]": reaction_left_y,
         f"displacement_x_at_evaluation_point (x={displacement_evaluation_x}, y={displacement_evaluation_y})[m]": displacement_x_at_evaluation_point,

examples/linear-elastic-plate-with-hole/kratos/postprocess_results.py

@@ -77,10 +77,8 @@ def postprocess_results(input_parameter_file, input_result_vtk, output_metrics_f
         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
+    # Compute maximum von Mises stress at Gauss points.
+    max_von_mises_stress_gauss_points = 0
     for key, values in mesh.cell_data.items():
         if "VON_MISES_STRESS" in key:
             max_von_mises_stress_gauss_points = float(np.max(values))
@@ -123,10 +121,18 @@ def postprocess_results(input_parameter_file, input_result_vtk, output_metrics_f
     else:
         displacement_x_at_evaluation_point = float(displacement_sampled[0, 0])
 
+    # Compute nodal displacement error (Euclidean norm of error vector at each node)
+    nodal_displacement_error = np.linalg.norm(displacement - u_ref, axis=1)
+    max_displacement_error_nodes = float(np.max(nodal_displacement_error))
+    
+    # Compute the number of dofs
+    num_dofs = int(mesh.n_points * 2)
+
     metrics = {
-        "max_von_mises_stress_nodes[Pa]": max_von_mises_stress_nodes,
-        "max_von_mises_stress_gauss_points[Pa]": max_von_mises_stress_gauss_points,
+        "number_of_dofs[-]": num_dofs,
+        "max_von_mises_stress[Pa]": max_von_mises_stress_gauss_points,
         "l2_error_displacement[m]": l2_error_displacement,
+        "max_displacement_error[m]": max_displacement_error_nodes,
         "reaction_force_left_boundary_x[N]": reaction_force_left_boundary_x,
         "reaction_force_left_boundary_y[N]": reaction_force_left_boundary_y,
         f"displacement_x_at_evaluation_point (x={displacement_evaluation_x}, y={displacement_evaluation_y})[m]": displacement_x_at_evaluation_point,