Merge branch 'bbar_impr_more_tests' into 'master'

a Jupyter Notebook benchmark for evaluating the B-bar method with simple tests and some others

See merge request ogs/ogs!5105

ProcessLib/SmallDeformation/Tests.cmake

@@ -351,8 +351,9 @@ AddTest(
 )
 
 if(NOT OGS_USE_PETSC)
-    NotebookTest(NOTEBOOKFILE Mechanics/CooksMembrane/CooksMembraneBbar.ipynb RUNTIME 1)
+    NotebookTest(NOTEBOOKFILE Mechanics/CooksMembrane/CooksMembraneBbar.py RUNTIME 1)
     NotebookTest(NOTEBOOKFILE Mechanics/Linear/SimpleMechanics.ipynb RUNTIME 5)
+    NotebookTest(NOTEBOOKFILE Mechanics/EvaluatingBbarWithSimpleExamples/evaluating_bbbar_with_simple_examples.py RUNTIME 5)
     NotebookTest(NOTEBOOKFILE Mechanics/Linear/DiscWithHole/Linear_Disc_with_hole.md RUNTIME 15)
     if(NOT WIN32)
         NotebookTest(NOTEBOOKFILE Mechanics/Linear/DiscWithHole/Linear_Disc_with_hole_convergence_analysis.ipynb RUNTIME 40)

Tests/Data/Mechanics/CooksMembrane/CooksMembraneBbar.py

+# ---
+# jupyter:
+#   jupytext:
+#     text_representation:
+#       extension: .py
+#       format_name: percent
+#       format_version: '1.3'
+#       jupytext_version: 1.16.2
+#   kernelspec:
+#     display_name: Python 3 (ipykernel)
+#     language: python
+#     name: python3
+# ---
+
+# %% [raw]
+# +++
+# title = "Cook's membrane example"
+# date = "2024-06-11"
+# author = "Wenqing Wang"
+# image = "figures/cooks_membrane.png"
+# web_subsection = "small-deformations"
+# weight = 3
+# +++
+
+# %% [markdown]
+# $$
+# \newcommand{\B}{\text{B}}
+# \newcommand{\F}{\text{F}}
+# \newcommand{\I}{\mathbf I}
+# \newcommand{\intD}[1]{\int_{\Omega_e}#1\mathrm{d}\Omega}
+# $$
+#
+# # Cook's membrane example for nearly icompressible solid
+#
+# ## B bar method
+#  Considering a strain decomposition: $\mathbf\epsilon = \underbrace{\mathbf\epsilon- \frac{1}{3}(\epsilon:\mathbf I)}_{\text{deviatoric}}\I +  \underbrace{\frac{1}{3}(\epsilon:\mathbf I)}_{\text{dilatational}} \I$.
+# 	   The idea of the B bar method is to use another quadrature rule to interpolate the dilatational part, which leads to a modified B matrix [1]:
+# $$
+# 	     \bar\B = \underbrace{\B - \B^{\text{dil}}}_{\text{original B elements}}+ \underbrace{{\bar\B}^{\text{dil}}}_{\text{by another quadrature rule} }
+# $$
+# There are several methods to form ${\bar\B}^{\text{dil}}$  such as selective integration, generalization of the mean-dilatation formulation. In the current OGS, we use the latter, which reads
+# $$
+# 		  {\bar\B}^{\text{dil}} = \frac{\intD{\B^{\text{dil}}(\xi)}}{\intD{}}
+# $$
+#
+# ## Example
+# To verify the implementation of the B bar method, the so called Cook's membrane is used as a benchmark.
+# Illustrated in the following figure, this example simulates a tapered
+# and swept panel of unit thickness. The left edge is clamped and the right edge is applied with  a distributed shearing load $F$ = 100 N/mm. The plane strain condition is considered. This numerical model is exactly the same as that is presented in the paper by T. Elguedj et al [1,2].
+#
+# <img src="figures/cooks_membrane.png" alt="Cook's membrane" width="320" height="320" />
+#
+# ## Reference
+#
+# [1] T.J.R. Hughes (1980). Generalization of selective integration procedures to anisotropic and nonlinear media. International Journal for Numerical Methods in Engineering, 15(9), 1413-1418.
+#
+# [2] T. Elguedj, Y. Bazilevs, V.M. Calo, T.J.R. Hughes (2008),
+#  $\bar\B$ and $\bar\F$ projection methods for nearly incompressible linear and non-linear elasticity and plasticity using higher-order NURBS elements, Computer Methods in Applied Mechanics and Engineering, 197(33--40), 2732-2762.
+#
+
+# %%
+import os
+import xml.etree.ElementTree as ET
+from pathlib import Path
+
+import matplotlib.pyplot as plt
+import matplotlib.tri as tri
+import numpy as np
+import pyvista as pv
+import vtuIO
+from ogs6py.ogs import OGS
+
+out_dir = Path(os.environ.get("OGS_TESTRUNNER_OUT_DIR", "_out"))
+if not out_dir.exists():
+    out_dir.mkdir(parents=True)
+
+
+# %%
+def get_last_vtu_file_name(pvd_file_name):
+    tree = ET.parse(Path(out_dir) / pvd_file_name)
+    root = tree.getroot()
+    # Get the last DataSet tag
+    last_dataset = root.findall(".//DataSet")[-1]
+
+    # Get the 'file' attribute of the last DataSet tag
+    file_attribute = last_dataset.attrib["file"]
+    return f"{out_dir}/" + file_attribute
+
+
+def get_top_uy(pvd_file_name):
+    top_point = (48.0e-3, 60.0e-3, 0)
+    file_name = get_last_vtu_file_name(pvd_file_name)
+    mesh = pv.read(file_name)
+    p_id = mesh.find_closest_point(top_point)
+    u = mesh.point_data["displacement"][p_id]
+    return u[1]
+
+
+# %%
+def run_single_test(mesh_name, output_prefix, use_bbar="false"):
+    model = OGS(INPUT_FILE="CooksMembrane.prj", PROJECT_FILE=f"{out_dir}/modified.prj")
+    model.replace_text(mesh_name, xpath="./mesh")
+    model.replace_text(use_bbar, xpath="./processes/process/use_b_bar")
+    model.replace_text(output_prefix, xpath="./time_loop/output/prefix")
+    model.replace_text(
+        "BiCGSTAB", xpath="./linear_solvers/linear_solver/eigen/solver_type"
+    )
+    model.replace_text("ILUT", xpath="./linear_solvers/linear_solver/eigen/precon_type")
+    vtu_file_name = output_prefix + "_ts_1_t_1.000000.vtu"
+    model.replace_text(vtu_file_name, xpath="./test_definition/vtkdiff[1]/file")
+    model.replace_text(vtu_file_name, xpath="./test_definition/vtkdiff[2]/file")
+    model.replace_text(vtu_file_name, xpath="./test_definition/vtkdiff[3]/file")
+
+    model.write_input()
+
+    # Run OGS
+    model.run_model(logfile=f"{out_dir}/out.txt", args=f"-o {out_dir} -m .")
+
+    # Get uy at the top
+    return get_top_uy(output_prefix + ".pvd")
+
+
+# %%
+mesh_names = [
+    "mesh.vtu",
+    "mesh_n10.vtu",
+    "mesh_n15.vtu",
+    "mesh_n20.vtu",
+    "mesh_n25.vtu",
+    "mesh_n30.vtu",
+]
+output_prefices_non_bbar = [
+    "cooks_membrane_sd_edge_div_4_non_bbar",
+    "cooks_membrane_sd_refined_mesh_10_non_bbar",
+    "cooks_membrane_sd_refined_mesh_15_non_bbar",
+    "cooks_membrane_sd_refined_mesh_20_non_bbar",
+    "cooks_membrane_sd_refined_mesh_25_non_bbar",
+    "cooks_membrane_sd_refined_mesh_30_non_bbar",
+]
+
+uys_at_top_non_bbar = []
+for mesh_name, output_prefix in zip(mesh_names, output_prefices_non_bbar):
+    uy_at_top = run_single_test(mesh_name, output_prefix)
+    uys_at_top_non_bbar.append(uy_at_top)
+
+expected_uys_at_top_non_bbar = np.array(
+    [
+        0.002164586784123102,
+        0.0022603329644579383,
+        0.002375295856067169,
+        0.002519725590136146,
+        0.0026515294133790837,
+        0.002868289617025223,
+    ]
+)
+np.testing.assert_allclose(
+    actual=uys_at_top_non_bbar, desired=expected_uys_at_top_non_bbar, atol=1e-10
+)
+
+# %%
+output_prefices = [
+    "cooks_membrane_sd_edge_div_4",
+    "cooks_membrane_sd_refined_mesh_10",
+    "cooks_membrane_sd_refined_mesh_15",
+    "cooks_membrane_sd_refined_mesh_20",
+    "cooks_membrane_sd_refined_mesh_25",
+    "cooks_membrane_sd_refined_mesh_30",
+]
+
+uys_at_top_bbar = []
+for mesh_name, output_prefix in zip(mesh_names, output_prefices):
+    uy_at_top = run_single_test(mesh_name, output_prefix, "true")
+    uys_at_top_bbar.append(uy_at_top)
+
+expected_uys_at_top_bbar = np.array(
+    [
+        0.0069574713856979,
+        0.007772616910217863,
+        0.007897597955618913,
+        0.007951479575082158,
+        0.007976349858390623,
+        0.007999718483861992,
+    ]
+)
+
+np.testing.assert_allclose(
+    actual=uys_at_top_bbar, desired=expected_uys_at_top_bbar, atol=2e-4
+)
+
+# %%
+ne = [4, 10, 15, 20, 25, 30]
+
+
+def plot_data(ne, u_y_bbar, uy_non_bbar, file_name=""):
+    # Plotting
+    plt.rcParams["figure.figsize"] = [5, 5]
+
+    if len(u_y_bbar) != 0:
+        plt.plot(
+            ne, np.array(u_y_bbar) * 1e3, marker="o", linestyle="dashed", label="B bar"
+        )
+    if len(uy_non_bbar) != 0:
+        plt.plot(
+            ne,
+            np.array(uy_non_bbar) * 1e3,
+            marker="x",
+            linestyle="dashed",
+            label="non B bar",
+        )
+
+    plt.xlabel("Number of elements per side")
+    plt.ylabel("Top right corner displacement /mm")
+    plt.legend()
+
+    plt.tight_layout()
+    if file_name != "":
+        plt.savefig(file_name)
+    plt.show()
+
+
+# %% [markdown]
+# ## Result
+#
+# ### Vertical diplacement at the top point
+#
+# The following figure shows that the convergence of the solutions obtained by using the B bar method follows the one presented in the paper by T. Elguedj et al [1]. However, the results obtained without the B bar method are quit far from the converged solution with the finest mesh.
+
+# %%
+plot_data(ne, uys_at_top_bbar, uys_at_top_non_bbar, "b_bar_linear.png")
+
+# %% [markdown]
+# ### Contour plot
+
+# %%
+nedges = ["4", "10", "15", "20", "25", "30"]
+
+
+def contour_plot(pvd_file_name, title):
+    file_name = get_last_vtu_file_name(pvd_file_name)
+    m_plot = vtuIO.VTUIO(file_name, dim=2)
+    triang = tri.Triangulation(m_plot.points[:, 0], m_plot.points[:, 1])
+    triang = tri.Triangulation(m_plot.points[:, 0], m_plot.points[:, 1])
+    s_plot = m_plot.get_point_field("sigma")
+    s_trace = s_plot[:, 0] + s_plot[:, 1] + s_plot[:, 2]
+    u_plot = m_plot.get_point_field("displacement")
+
+    fig, ax = plt.subplots(ncols=2, figsize=(8, 3))
+    ax[0].set_title(title, loc="left", y=1.12)
+    plt.subplots_adjust(wspace=0.5)
+
+    contour_stress = ax[0].tricontourf(triang, s_trace, cmap="viridis")
+    contour_displacement = ax[1].tricontourf(triang, u_plot[:, 1], cmap="gist_rainbow")
+    fig.colorbar(contour_stress, ax=ax[0], label="Stress trace / MPa")
+    fig.colorbar(contour_displacement, ax=ax[1], label="Dispplacement / m")
+    fig.tight_layout()
+    plt.savefig(pvd_file_name + ".png")
+    plt.show()
+
+
+# %% [markdown]
+# #### Results obtained without the B bar method:
+
+# %%
+for nedge, output_prefix in zip(nedges, output_prefices_non_bbar):
+    contour_plot(output_prefix + ".pvd", "Number of elements per side: " + nedge)
+
+# %% [markdown]
+# #### Results obtained with the B bar method:
+
+# %%
+for nedge, output_prefix in zip(nedges, output_prefices):
+    contour_plot(output_prefix + ".pvd", "Number of elements per side: " + nedge)
+
+# %% [markdown]
+# The contour plots show that even with the coarsest mesh, the B bar method still gives reasonable stress result.

Tests/Data/Mechanics/EvaluatingBbarWithSimpleExamples/quad.gml

+<?xml version="1.0" encoding="utf-8"?>
+<OpenGeoSysGLI>
+    <name>geometry</name>
+    <points>
+         <point id="0"  x="0." y="0."  z="0." name = "origin"/>
+         <point id="1"  x="1.0" y="0.0"  z="0." />
+         <point id="2"  x="1.0" y="1.0"  z="0." />
+         <point id="3"  x="0." y="1.0"  z="0." />
+    </points>
+
+    <polylines>
+        <polyline id="1" name="bottom">
+            <pnt>0</pnt>
+            <pnt>1</pnt>
+        </polyline>
+        <polyline id="2" name="top">
+            <pnt>2</pnt>
+            <pnt>3</pnt>
+        </polyline>
+    </polylines>
+</OpenGeoSysGLI>

Tests/Data/Mechanics/EvaluatingBbarWithSimpleExamples/quad_edge_div_2.vtu

+<?xml version="1.0"?>
+<VTKFile type="UnstructuredGrid" version="1.0" byte_order="LittleEndian" header_type="UInt64">
+  <UnstructuredGrid>
+    <Piece NumberOfPoints="21"                   NumberOfCells="4"                   >
+      <PointData>
+      </PointData>
+      <CellData>
+        <DataArray type="Int32" Name="MaterialIDs" format="appended" RangeMin="0"                    RangeMax="0"                    offset="0"                   />
+      </CellData>
+      <Points>
+        <DataArray type="Float64" Name="Points" NumberOfComponents="3" format="appended" RangeMin="0"                    RangeMax="1.4142135624"         offset="32"                  />
+      </Points>
+      <Cells>
+        <DataArray type="Int64" Name="connectivity" format="appended" RangeMin=""                     RangeMax=""                     offset="716"                 />
+        <DataArray type="Int64" Name="offsets" format="appended" RangeMin=""                     RangeMax=""                     offset="1068"                />
+        <DataArray type="UInt8" Name="types" format="appended" RangeMin=""                     RangeMax=""                     offset="1124"                />
+      </Cells>
+    </Piece>
+  </UnstructuredGrid>
+  <AppendedData encoding="base64">
+   _EAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA+AEAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAADgPwAAAAAAAAAAAAAAAAAAAAAAAAAAAADwPwAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAOA/AAAAAAAAAAAAAAAAAADgPwAAAAAAAOA/AAAAAAAAAAAAAAAAAADwPwAAAAAAAOA/AAAAAAAAAAAAAAAAAAAAAAAAAAAAAPA/AAAAAAAAAAAAAAAAAADgPwAAAAAAAPA/AAAAAAAAAAAAAAAAAADwPwAAAAAAAPA/AAAAAAAAAAAAAAAAAAAAAAAAAAAAANA/AAAAAAAAAAAAAAAAAAAAAAAAAAAAAOg/AAAAAAAAAAAAAAAAAADQPwAAAAAAAAAAAAAAAAAAAAAAAAAAAADQPwAAAAAAAOA/AAAAAAAAAAAAAAAAAADQPwAAAAAAAPA/AAAAAAAAAAAAAAAAAADgPwAAAAAAANA/AAAAAAAAAAAAAAAAAADgPwAAAAAAAOg/AAAAAAAAAAAAAAAAAADoPwAAAAAAAAAAAAAAAAAAAAAAAAAAAADoPwAAAAAAAOA/AAAAAAAAAAAAAAAAAADoPwAAAAAAAPA/AAAAAAAAAAAAAAAAAADwPwAAAAAAANA/AAAAAAAAAAAAAAAAAADwPwAAAAAAAOg/AAAAAAAAAAA=AAEAAAAAAAAAAAAAAAAAAAEAAAAAAAAABAAAAAAAAAADAAAAAAAAAAsAAAAAAAAADgAAAAAAAAAMAAAAAAAAAAkAAAAAAAAAAQAAAAAAAAACAAAAAAAAAAUAAAAAAAAABAAAAAAAAAAQAAAAAAAAABMAAAAAAAAAEQAAAAAAAAAOAAAAAAAAAAMAAAAAAAAABAAAAAAAAAAHAAAAAAAAAAYAAAAAAAAADAAAAAAAAAAPAAAAAAAAAA0AAAAAAAAACgAAAAAAAAAEAAAAAAAAAAUAAAAAAAAACAAAAAAAAAAHAAAAAAAAABEAAAAAAAAAFAAAAAAAAAASAAAAAAAAAA8AAAAAAAAAIAAAAAAAAAAIAAAAAAAAABAAAAAAAAAAGAAAAAAAAAAgAAAAAAAAAA==BAAAAAAAAAAXFxcX
+  </AppendedData>
+</VTKFile>

Tests/Data/Mechanics/EvaluatingBbarWithSimpleExamples/simple_b_bar_test.prj

+<?xml version="1.0" encoding="ISO-8859-1"?>
+<OpenGeoSysProject>
+    <mesh>quad_edge_div_2.vtu</mesh>
+    <geometry>quad.gml</geometry>
+    <search_length_algorithm>
+        <type>fixed</type>
+        <value>1e-5</value>
+    </search_length_algorithm>    
+    <processes>
+        <process>
+            <name>SD</name>
+            <type>SMALL_DEFORMATION</type>
+            <integration_order>3</integration_order>
+            <use_b_bar>true</use_b_bar>
+            <constitutive_relation>
+                <type>LinearElasticIsotropic</type>
+                <youngs_modulus>E</youngs_modulus>
+                <poissons_ratio>nu</poissons_ratio>
+            </constitutive_relation>
+            <specific_body_force>0 0</specific_body_force>
+            <process_variables>
+                <process_variable>displacement</process_variable>
+            </process_variables>
+            <secondary_variables>
+                <secondary_variable internal_name="epsilon" output_name="epsilon"/>
+                <secondary_variable internal_name="sigma" output_name="sigma"/>
+            </secondary_variables>
+        </process>
+    </processes>
+    <time_loop>
+        <processes>
+            <process ref="SD">
+                <nonlinear_solver>basic_newton</nonlinear_solver>
+                <convergence_criterion>
+                    <type>PerComponentDeltaX</type>
+                    <norm_type>NORM2</norm_type>
+                    <reltols> 1.0e-9 1.0e-9 </reltols>
+                </convergence_criterion>
+                <time_discretization>
+                    <type>BackwardEuler</type>
+                </time_discretization>
+                <time_stepping>
+                    <type>FixedTimeStepping</type>
+                    <t_initial> 0.0 </t_initial>
+                    <t_end>1</t_end>
+                    <timesteps>
+                        <pair>
+                            <repeat>1</repeat>
+                            <delta_t>1</delta_t>
+                        </pair>
+                    </timesteps>
+                </time_stepping>
+            </process>
+        </processes>
+        <output>
+            <type>VTK</type>
+            <prefix>simple_b_bar_test</prefix>
+            <timesteps>
+                <pair>
+                    <repeat> 1 </repeat>
+                    <each_steps> 1 </each_steps>
+                </pair>
+            </timesteps>
+            <variables>
+                <variable>displacement</variable>
+                <variable>epsilon</variable>
+                <variable>sigma</variable>
+            </variables>
+            <suffix>_ts_{:timestep}_t_{:time}</suffix>
+        </output>
+    </time_loop>
+    <media>
+        <medium>
+            <phases>
+                <phase>
+                    <type>Solid</type>
+                    <properties>
+                        <property>
+                            <name>density</name>
+                            <type>Parameter</type>
+                            <parameter_name>rho</parameter_name>
+                        </property>
+                    </properties>
+                </phase>
+            </phases>
+        </medium>
+    </media>
+    <parameters>
+        <parameter>
+            <name>E</name>
+            <type>Constant</type>
+            <value>1.0E+10</value>
+        </parameter>
+        <parameter>
+            <name>nu</name>
+            <type>Constant</type>
+            <value>0.2</value>
+        </parameter>
+        <parameter>
+            <name>rho</name>
+            <type>Constant</type>
+            <value>0.0</value>
+        </parameter>
+        <parameter>
+            <name>zero</name>
+            <type>Constant</type>
+            <value>0.0</value>
+        </parameter>
+        <parameter>
+            <name>u0</name>
+            <type>Constant</type>
+            <values>0 0</values>
+        </parameter>
+        <parameter>
+            <name>load</name>
+            <type>Constant</type>
+            <values>-10e+6</values>
+        </parameter>
+    </parameters>
+    <process_variables>
+        <process_variable>
+            <name>displacement</name>
+            <components>2</components>
+            <order>2</order>
+            <initial_condition>u0</initial_condition>
+            <boundary_conditions>
+                <boundary_condition>
+                    <geometrical_set>geometry</geometrical_set>
+                    <geometry>bottom</geometry>
+                    <type>Dirichlet</type>
+                    <component>0</component>
+                    <parameter>zero</parameter>
+                </boundary_condition>
+                <boundary_condition>
+                    <geometrical_set>geometry</geometrical_set>
+                    <geometry>bottom</geometry>
+                    <type>Dirichlet</type>
+                    <component>1</component>
+                    <parameter>zero</parameter>
+                </boundary_condition>
+                <boundary_condition>
+                    <geometrical_set>geometry</geometrical_set>
+                    <geometry>top</geometry>
+                    <type>Neumann</type>
+                    <component>1</component>
+                    <parameter>load</parameter>
+                </boundary_condition>
+            </boundary_conditions>
+        </process_variable>
+    </process_variables>
+    <nonlinear_solvers>
+        <nonlinear_solver>
+            <name>basic_newton</name>
+            <type>Newton</type>
+            <max_iter>10</max_iter>
+            <linear_solver>general_linear_solver</linear_solver>
+        </nonlinear_solver>
+    </nonlinear_solvers>
+    <linear_solvers>
+        <linear_solver>
+            <name>general_linear_solver</name>
+            <eigen>
+                <solver_type>SparseLU</solver_type>
+                <scaling>true</scaling>
+            </eigen>
+        </linear_solver>
+    </linear_solvers>
+</OpenGeoSysProject>