diff --git a/ProcessLib/PhaseField/Tests.cmake b/ProcessLib/PhaseField/Tests.cmake
index 118afacc655e492450683eba827f877d57a0ca42..dbe103f1b4c65cac74ae1ab984b61bfe2c4fec76 100644
--- a/ProcessLib/PhaseField/Tests.cmake
+++ b/ProcessLib/PhaseField/Tests.cmake
@@ -137,4 +137,9 @@ AddTest(
if(OGS_USE_PETSC)
NotebookTest(NOTEBOOKFILE PhaseField/surfing_jupyter_notebook/surfing_pyvista.ipynb RUNTIME 25 RESOURCE_LOCK PYVISTA)
+ NotebookTest(NOTEBOOKFILE PhaseField/beam_jupyter_notebook/beam.ipynb RUNTIME 500 RESOURCE_LOCK PYVISTA)
+endif()
+
+if(TEST nb-PhaseField/beam_jupyter_notebook/beam-LARGE)
+ set_tests_properties(nb-PhaseField/beam_jupyter_notebook/beam-LARGE PROPERTIES PROCESSORS 3)
endif()
diff --git a/Tests/Data/PhaseField/beam_jupyter_notebook/beam.ipynb b/Tests/Data/PhaseField/beam_jupyter_notebook/beam.ipynb
new file mode 100644
index 0000000000000000000000000000000000000000..8a1e42980a51e196f318bed2356ac9550f81727b
--- /dev/null
+++ b/Tests/Data/PhaseField/beam_jupyter_notebook/beam.ipynb
@@ -0,0 +1,373 @@
+{
+ "cells": [
+ {
+ "cell_type": "raw",
+ "id": "90c995c2",
+ "metadata": {},
+ "source": [
+ "notebook = \"PhaseField/beam_jupyter_notebook/beam_pyvista.ipynb\"\n",
+ "author = \"Matthes Kantzenbach, Keita Yoshioka, Mostafa Mollaali\"\n",
+ "date = \"2022-11-28\"\n",
+ "title = \"Beam\"\n",
+ "web_subsection = \"phase-field\"\n",
+ "<!--eofm-->"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "41d0a5bc",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "from ogs6py import ogs\n",
+ "import os\n",
+ "import ogs6py\n",
+ "import numpy as np\n",
+ "import matplotlib.pyplot as plt\n",
+ "import pyvista as pv\n",
+ "import time\n",
+ "from xml.dom import minidom\n",
+ "from types import MethodType"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "aa987d95",
+ "metadata": {},
+ "source": [
+ "\n",
+ "## Problem description\n",
+ "\n",
+ "In this example, it is shown, how phasefield models can be applied to simulate crack initiation.\n",
+ "Consider a bar $\\Omega=[0,1]\\times [0,0.05]\\times [0,0.05]$ where the left end $\\partial \\Omega_{Left}$ is fixed and a displacement $U_x(t)$ is applied to the right end $\\partial \\Omega_{Right}$, resulting in a uniaxial state of stress $\\boldsymbol{\\sigma} = \\sigma_{33} \\, \\mathbf{n}_3 \\otimes \\mathbf{n}_3$. At a certain point, fracture will initiate resulting in failure of the beam. $\\newline$ This will be examined for positive $U_x(t)$ (tensile stress) and for negative $U_x(t)$ (compressive stress). As energy-split models the Isotropic and the Volumetric Deviatoric are used and as phasefield models $\\texttt{AT}_1$ and $\\texttt{AT}_2$. So in total 8 different cases are tested.\n",
+ "The aim is to estimate the ultimate load before fracture initiates for each case, as this can be veryfied against analytical results.\n",
+ "\n",
+ "![Schematic view of beam](./figures/bar_.png#one-half \"Schematic view of beam.\")\n",
+ "\n",
+ "Also considering homogeneous damage development (i.e.$\\nabla v = 0$), the analytical solution for the strengths under uniaxial tension and compression are as follows.\n",
+ "\n",
+ "| | Isotropic | Vol-Dev |\n",
+ "|:---:|:---:|:---:|\n",
+ "|Tensile $$\\texttt{AT}_1$$| $$\\sqrt{ \\dfrac{3 G_\\mathrm{c} E} {8 \\ell}}$$| $$\\sqrt{\\dfrac{3 G_\\mathrm{c} E}{8 \\ell}}$$ |\n",
+ "|Tensile $$\\texttt{AT}_2$$| $$ \\dfrac{3}{16} \\sqrt{ \\dfrac{3 G_\\mathrm{c} E } { \\ell}}$$ | $$\\dfrac{3}{16} \\sqrt{\\dfrac{3 G_\\mathrm{c} E}{\\ell}}$$ |\n",
+ "|Compressive $$\\texttt{AT}_1$$| $$-\\sqrt{ \\dfrac{3 G_\\mathrm{c} E} {8 \\ell}}$$ | $$- \\sqrt{ \\dfrac{9 G_\\mathrm{c} E} {16 \\ell(1+\\nu)}}$$ |\n",
+ "|Compressive $$\\texttt{AT}_2$$| $$ -\\dfrac{3}{16} \\sqrt{ \\dfrac{3 G_\\mathrm{c} E } { \\ell}}$$ | $$ -\\dfrac{9}{32} \\sqrt{\\dfrac{2 G_\\mathrm{c} E}{ \\ell (1+\\nu)}}$$ |"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "7ac26dee",
+ "metadata": {},
+ "source": [
+ "## Define some helper functions"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "b137ce19",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "data_dir = os.environ.get('OGS_DATA_DIR', '../../..')\n",
+ "out_dir = os.environ.get('OGS_TESTRUNNER_OUT_DIR', '_out')\n",
+ "if not os.path.exists(out_dir):\n",
+ " os.makedirs(out_dir) \n",
+ "\n",
+ "output_dir= out_dir\n",
+ "\n",
+ "# define method to be assigned to model, to replace a specific curve, given by name \n",
+ "# (analogue to replace_parameter method)\n",
+ "def replace_curve(self, name=None, value=None, coords=None, parametertype=None, valuetag=\"values\", coordstag=\"coords\"):\n",
+ " root = self._get_root()\n",
+ " parameterpath = \"./curves/curve\"\n",
+ " parameterpointer = self._get_parameter_pointer(root, name, parameterpath)\n",
+ " self._set_type_value(parameterpointer, value, parametertype, valuetag=valuetag)\n",
+ " self._set_type_value(parameterpointer, coords, parametertype, valuetag=coordstag)\n",
+ "# define method to change timstepping in project file \n",
+ "def set_timestepping(model,repeat_list, delta_t_list):\n",
+ " model.remove_element(xpath='./time_loop/processes/process/time_stepping/timesteps/pair')\n",
+ " for i in range(len(repeat_list)):\n",
+ " model.add_block(blocktag = 'pair',parent_xpath='./time_loop/processes/process/time_stepping/timesteps', taglist = ['repeat', 'delta_t'], textlist = [repeat_list[i], delta_t_list[i]])"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "d7d09975",
+ "metadata": {},
+ "source": [
+ "## Define function generating mesh, modifying project file and running ogs with given parameters"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "a34e61bd",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "def ogs_beam(phasefield_model, energy_split_model, mesh_size = 0.01, length_scale = 0.02, bc_displacement = 5, ts_coords='0 0.05 1', values ='0 0.25 1', repeat_list=None, delta_t_list=None, hypre = False):\n",
+ " ##phasefield_model: 'AT1' or 'AT2' \n",
+ " ##energy_split_model: 'VolumetricDeviatoric' or 'Isotropic'\n",
+ " \n",
+ " without_hypre='-ksp_type cg -pc_type bjacobi -ksp_atol 1e-14 -ksp_rtol 1e-14'\n",
+ " with_hypre='-ksp_type cg -pc_type hypre -pc_hypre_type boomeramg -pc_hypre_boomeramg_strong_threshold 0.7 -ksp_atol 1e-8 -ksp_rtol 1e-8' \n",
+ " #file's name\n",
+ " prj_name = \"beam.prj\"\n",
+ " print(f\"> Running beam model {phasefield_model} – {energy_split_model} ... <\")\n",
+ " logfile = f\"{out_dir}/log_{phasefield_model}_{energy_split_model}.txt\"\n",
+ " #beam dimensions\n",
+ " beam_height = 0.05\n",
+ " beam_depth = beam_height\n",
+ " beam_length = 1.\n",
+ " #mesh properties\n",
+ " h = mesh_size # distance between nodes\n",
+ " ls = length_scale\n",
+ " #generate prefix from properties\n",
+ " if energy_split_model == 'VolumetricDeviatoric':\n",
+ " prefix = phasefield_model + '_vd'\n",
+ " elif energy_split_model == 'Isotropic':\n",
+ " prefix = phasefield_model + '_iso'\n",
+ " else :\n",
+ " raise ValueError('\"'+energy_split_model+'\"'+' is no valid input for energy_split_model, choose between \"VolumetricDeviatoric\" and \"Isotropic\"') \n",
+ " if bc_displacement > 0:\n",
+ " prefix = prefix + '_tensile'\n",
+ " else :\n",
+ " prefix = prefix + '_compressive'\n",
+ " #generate mesh \n",
+ " ! generateStructuredMesh -o {out_dir}/bar_.vtu -e hex --lx {beam_length} --nx {round(beam_length/h)} --ly {beam_height} --ny {round(beam_height/h)} --lz {beam_depth} --nz {round(beam_depth/h)} > {logfile}\n",
+ " ! NodeReordering -i {out_dir}/bar_.vtu -o {out_dir}/bar.vtu >> {logfile}\n",
+ " ! ExtractSurface -i {out_dir}/bar.vtu -o {out_dir}/bar_left.vtu -x 1 -y 0 -z 0 >> {logfile}\n",
+ " ! ExtractSurface -i {out_dir}/bar.vtu -o {out_dir}/bar_right.vtu -x -1 -y 0 -z 0 >> {logfile}\n",
+ " ! partmesh -s -o {out_dir} -i {out_dir}/bar.vtu >> {logfile}\n",
+ " ! partmesh -m -n 3 -o {out_dir} -i {out_dir}/bar.vtu -- {out_dir}/bar_right.vtu {out_dir}/bar_left.vtu >> {logfile}\n",
+ " #change properties in prj file\n",
+ " model = ogs.OGS(INPUT_FILE=prj_name, PROJECT_FILE=f\"{out_dir}/{prj_name}\", MKL=True)\n",
+ " model.replace_parameter_value(name=\"ls\", value=ls)\n",
+ " model.replace_text(phasefield_model, xpath=\"./processes/process/phasefield_model\")\n",
+ " model.replace_text(energy_split_model, xpath=\"./processes/process/energy_split_model\")\n",
+ " model.replace_text(prefix, xpath=\"./time_loop/output/prefix\")\n",
+ " model.replace_parameter_value(name=\"dirichlet_right\", value=bc_displacement)\n",
+ " model.replace_curve = MethodType(replace_curve, model)\n",
+ " model.replace_curve(name=\"dirichlet_time\",value=values, coords=ts_coords)\n",
+ " if repeat_list != None and delta_t_list != None: \n",
+ " set_timestepping(model,repeat_list, delta_t_list)\n",
+ " else:\n",
+ " set_timestepping(model,['1'], ['1e-2']) \n",
+ " if hypre == True:\n",
+ " model.replace_text(with_hypre, xpath='./linear_solvers/linear_solver/petsc/parameters',occurrence=1)\n",
+ " else:\n",
+ " model.replace_text(without_hypre, xpath='./linear_solvers/linear_solver/petsc/parameters', occurrence=1)\n",
+ " model.write_input()\n",
+ " #run ogs\n",
+ " t0 = time.time()\n",
+ " print(\" > OGS started execution ...\")\n",
+ " ! mpirun -n 3 ogs {out_dir}/{prj_name} -o {output_dir} >> {logfile}\n",
+ " tf = time.time()\n",
+ " print(\" > OGS terminated execution. Elapsed time: \", round(tf - t0, 2), \" s.\")"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "324f380a",
+ "metadata": {},
+ "source": [
+ "## Input data\n",
+ "\n",
+ "The values of the material properties are choosen as follows.\n",
+ "\n",
+ "\n",
+ "| **Name** | **Value** | **Unit** | **Symbol** |\n",
+ "|--------------------------------|--------------------|--------------|------------|\n",
+ "| _Young's modulus_ | 1 | Pa | $E$ |\n",
+ "| _Critical energy release rate_ | 1 | Pa$\\cdot$m | $G_{c}$ |\n",
+ "| _Poisson's ratio_ | 0.15 | $-$ | $\\nu$ |\n",
+ "| _Regularization parameter_ | 3$h$ | m | $\\ell$ |\n",
+ "| _Length_ | $1$ | m | $L$ |\n",
+ "| _Height_ | $0.05$ | m | $H$ |\n",
+ "| _Depth_ | $0.05$ | m | $D$ |\n"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "02482e8b",
+ "metadata": {},
+ "source": [
+ "## Run Simulations"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "71c6b710",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "pf_ms = ['AT1', 'AT2']\n",
+ "es_ms = ['VolumetricDeviatoric','Isotropic']\n",
+ "displ = [4.,-4.]\n",
+ "'''\n",
+ "for a in pf_ms:\n",
+ " for c in displ:\n",
+ " ogs_beam(a,es_ms[1],bc_displacement = c,mesh_size = 0.01, length_scale = 0.03)\n",
+ "ogs_beam(pf_ms[1],es_ms[0],bc_displacement = 4,mesh_size = 0.01, length_scale = 0.03)\n",
+ " \n",
+ "# run AT1_vd_tensile with smaller timesteps in critical time range\n",
+ "ogs_beam(pf_ms[0],es_ms[0],bc_displacement = 5,mesh_size = 0.01, length_scale = 0.03,repeat_list=['62','2','20','1'], delta_t_list=['1e-2','1e-3','1e-4','1e-2'])\n",
+ "\n",
+ "# run VolumetricDeviatoric in Compression with Hypre and smaller timesteps in critical time range\n",
+ "ogs_beam(pf_ms[1],es_ms[0],bc_displacement = -4.5,mesh_size = 0.01, length_scale = 0.03, hypre = True, repeat_list=['70','4','30','1'], delta_t_list=['1e-2','1e-3','1e-4','1e-2'])\n",
+ "\n",
+ "# loosen relative error tolerance for displacement process in order to get convergence for the AT1 case\n",
+ "prj_path='./'\n",
+ "prj_name = \"beam.prj\"\n",
+ "model = ogs.OGS(INPUT_FILE=prj_path+prj_name, PROJECT_FILE=prj_path+prj_name, MKL=True)\n",
+ "model.replace_text('1e-6', xpath=\"./time_loop/processes/process/convergence_criterion/reltol\",occurrence=0)\n",
+ "model.write_input()\n",
+ "ogs_beam(pf_ms[0],es_ms[0],bc_displacement = -4.95, mesh_size = 0.01, length_scale = 0.03, hypre= True, repeat_list=['66', '8','3','3','20','1'], delta_t_list=['1e-2','1e-3','1e-4','1e-5','1e-6','1e-2'],ts_coords='0 0.1 1', values ='0 0.5 1')\n",
+ "model = ogs.OGS(INPUT_FILE=prj_path+prj_name, PROJECT_FILE=prj_path+prj_name, MKL=True)\n",
+ "model.replace_text('1e-14', xpath=\"./time_loop/processes/process/convergence_criterion/reltol\",occurrence=0)\n",
+ "model.write_input()\n",
+ "''' \n",
+ "## run only cases easy to handle with coarse timestepping:\n",
+ "for a in pf_ms:\n",
+ " for b in es_ms:\n",
+ " for c in displ:\n",
+ " if a == 'AT1' and b == 'VolumetricDeviatoric':\n",
+ " continue\n",
+ " if a == 'AT2' and b == 'VolumetricDeviatoric' and c < 0:\n",
+ " ogs_beam(a,b,bc_displacement = c,mesh_size = 0.01,length_scale = 0.03, hypre = True, repeat_list=['1'],delta_t_list=['1e-1'])\n",
+ " else:\n",
+ " ogs_beam(a,b,bc_displacement = c,mesh_size = 0.01,length_scale = 0.03, repeat_list=['1'],delta_t_list=['1e-1'])\n",
+ " "
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "44016656",
+ "metadata": {},
+ "source": [
+ "## Results\n",
+ "\n",
+ "The final timestep of the AT1 iso tensile case is shown exemplary for the phasefield after fracture.\n",
+ "\n",
+ "![Phasefield after fracture](./figures/beam_final.png \"phasefield after fracture\")"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "d2b7386b",
+ "metadata": {},
+ "source": [
+ "## Post-processing\n",
+ "The force applied to the beam is compared to the change in length of the beam."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "b1357ab9",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# define function to obtain displacement applied at the right end of the beam from rvtu file\n",
+ "def displ_right(filename):\n",
+ " data = pv.read(filename)\n",
+ " data.point_data[\"displacement\"]\n",
+ " max_x = max(data.points[:,0])\n",
+ " return np.mean(data.point_data[\"displacement\"][:,0], where= np.transpose(data.points[:,0]==max_x))\n",
+ "# define fuction to obtain force acting on the right end of the beam from vtu file\n",
+ "def force_right(filename):\n",
+ " data = pv.read(filename)\n",
+ " data.point_data[\"NodalForces\"]\n",
+ " max_x = max(data.points[:,0])\n",
+ " return np.sum(data.point_data[\"NodalForces\"][:,0], where= np.transpose(data.points[:,0]==max_x))\n",
+ "# define function applying obove functions on all vtu file listet in pvd file, returning force-displacement curve\n",
+ "def force_displ_from_pvd(pvd):\n",
+ " doc = minidom.parse(pvd)\n",
+ " DataSets = doc.getElementsByTagName(\"DataSet\")\n",
+ " vtu_files = [x.getAttribute(\"file\") for x in DataSets]\n",
+ " forces_right = [force_right(f\"{out_dir}/{x}\") for x in vtu_files]\n",
+ " displs_right = [displ_right(f\"{out_dir}/{x}\") for x in vtu_files]\n",
+ " return [forces_right,displs_right]\n",
+ "\n"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "5375302a",
+ "metadata": {},
+ "source": [
+ "## Plot force-strain curves"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "5ec1718d",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "prefixes = ['AT1_vd_tensile','AT1_iso_tensile','AT2_vd_tensile','AT2_iso_tensile', 'AT1_vd_compressive', 'AT1_iso_compressive', 'AT2_vd_compressive', 'AT2_iso_compressive']\n",
+ "labels = [r'$\\texttt{AT}_1$ vol-dev tensile',r'$\\texttt{AT}_1$ iso tensile',r'$\\texttt{AT}_2$ vol-dev tensile',r'$\\texttt{AT}_2$ iso tensile', r'$\\texttt{AT}_1$ vol-dev compressive', r'$\\texttt{AT}_1$ iso compressive', r'$\\texttt{AT}_2$ vol-dev compressive', r'$\\texttt{AT}_2$ iso compressive']\n",
+ "ls=['-','--']\n",
+ "colors = ['#ffdf4d','#006ddb','#8f4e00','#ff6db6','#920000','#b66dff','#db6d00','#490092']\n",
+ "\n",
+ "fig, ax = plt.subplots()\n",
+ "plt.rc('text', usetex=True)\n",
+ "fig.set_size_inches(18.5, 10.5)\n",
+ "for i,pre in enumerate(prefixes):\n",
+ " pvd = f\"{output_dir}/{pre}.pvd\"\n",
+ " if os.path.isfile(pvd) :\n",
+ " curve = force_displ_from_pvd(pvd)\n",
+ " ax.plot(curve[1],curve[0],ls[i%2], label = labels[i],linewidth=5, color = colors[i], alpha= 1)\n",
+ " \n",
+ "plt.rcParams['xtick.labelsize'] = 16 \n",
+ "plt.rcParams['ytick.labelsize'] = 16 \n",
+ "ax.grid(linestyle='dashed') \n",
+ "ax.set_xlabel('$\\Delta [m]$',fontsize =18)\n",
+ "ax.set_ylabel('$F_y [N]$',fontsize =18)\n",
+ "plt.legend(fontsize =18, ncol = 2)\n",
+ "ax.axhline(y = 0, color = 'black',linewidth=1) \n",
+ "ax.axvline(x = 0, color = 'black',linewidth=1) \n",
+ "ax.set_xlim(-4.5, 4.5);\n"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "9346910f",
+ "metadata": {},
+ "source": [
+ "Running all the cases with finer timesteps yields the following figure. The failure of the beam in the simulation is observed, when the loading reaches the analytical strength for the corresponding case.\n",
+ "![Results](./figures/beam_result_with_analytical.png \"Results\")"
+ ]
+ }
+ ],
+ "metadata": {
+ "kernelspec": {
+ "display_name": "Python 3.10.8 ('.venv': venv)",
+ "language": "python",
+ "name": "python3"
+ },
+ "language_info": {
+ "codemirror_mode": {
+ "name": "ipython",
+ "version": 3
+ },
+ "file_extension": ".py",
+ "mimetype": "text/x-python",
+ "name": "python",
+ "nbconvert_exporter": "python",
+ "pygments_lexer": "ipython3",
+ "version": "3.10.8"
+ },
+ "vscode": {
+ "interpreter": {
+ "hash": "99f98b0875d01078cd28087d74fcb98ba3e814a2f158be930370408623897ead"
+ }
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
diff --git a/Tests/Data/PhaseField/beam_jupyter_notebook/beam.prj b/Tests/Data/PhaseField/beam_jupyter_notebook/beam.prj
new file mode 100644
index 0000000000000000000000000000000000000000..484c0e7e06adf25d69634bd46c06ca3686091de1
--- /dev/null
+++ b/Tests/Data/PhaseField/beam_jupyter_notebook/beam.prj
@@ -0,0 +1,284 @@
+<?xml version='1.0' encoding='ISO-8859-1'?>
+<OpenGeoSysProject>
+ <meshes>
+ <mesh>bar.vtu</mesh>
+ <mesh>bar_left.vtu</mesh>
+ <mesh>bar_right.vtu</mesh>
+ </meshes>
+ <processes>
+ <process>
+ <name>PhaseField</name>
+ <type>PHASE_FIELD</type>
+ <coupling_scheme>staggered</coupling_scheme>
+ <integration_order>2</integration_order>
+ <phasefield_model>AT1</phasefield_model>
+ <energy_split_model>Isotropic</energy_split_model>
+ <irreversible_threshold>0.05</irreversible_threshold>
+ <constitutive_relation>
+ <type>LinearElasticIsotropic</type>
+ <youngs_modulus>E</youngs_modulus>
+ <poissons_ratio>nu</poissons_ratio>
+ </constitutive_relation>
+ <phasefield_parameters>
+ <residual_stiffness>k</residual_stiffness>
+ <crack_resistance>gc</crack_resistance>
+ <crack_length_scale>ls</crack_length_scale>
+ </phasefield_parameters>
+ <solid_density>rho_sr</solid_density>
+ <process_variables>
+ <phasefield>phasefield</phasefield>
+ <displacement>displacement</displacement>
+ </process_variables>
+ <secondary_variables>
+ <secondary_variable type="static" internal_name="sigma" output_name="sigma"/>
+ <secondary_variable type="static" internal_name="epsilon" output_name="epsilon"/>
+ </secondary_variables>
+ <specific_body_force>0 0 0</specific_body_force>
+ </process>
+ </processes>
+ <time_loop>
+ <global_process_coupling>
+ <max_iter> 1000 </max_iter>
+ <convergence_criteria>
+ <!-- convergence criterion for the first process -->
+ <convergence_criterion>
+ <type>DeltaX</type>
+ <norm_type>INFINITY_N</norm_type>
+ <abstol>1.e-1</abstol>
+ <reltol>1.e-1</reltol>
+ </convergence_criterion>
+ <!-- convergence criterion for the second process -->
+ <convergence_criterion>
+ <type>DeltaX</type>
+ <norm_type>INFINITY_N</norm_type>
+ <abstol>1.e-4</abstol>
+ <reltol>1.e-10</reltol>
+ </convergence_criterion>
+ </convergence_criteria>
+ </global_process_coupling>
+ <processes>
+ <process ref="PhaseField">
+ <nonlinear_solver>basic_newton_u</nonlinear_solver>
+ <convergence_criterion>
+ <type>DeltaX</type>
+ <norm_type>NORM2</norm_type>
+ <reltol>1e-14</reltol>
+ </convergence_criterion>
+ <time_discretization>
+ <type>BackwardEuler</type>
+ </time_discretization>
+ <time_stepping>
+ <type>FixedTimeStepping</type>
+ <t_initial>0</t_initial>
+ <t_end>1.0</t_end>
+ <timesteps>
+ <pair>
+ <repeat>1</repeat>
+ <delta_t>1e-1</delta_t>
+ </pair>
+ </timesteps>
+ </time_stepping>
+ </process>
+ <process ref="PhaseField">
+ <nonlinear_solver>petsc_snes</nonlinear_solver>
+ <convergence_criterion>
+ <type>DeltaX</type>
+ <norm_type>NORM2</norm_type>
+ <reltol>1.e-14</reltol>
+ </convergence_criterion>
+ <time_discretization>
+ <type>BackwardEuler</type>
+ </time_discretization>
+ <time_stepping>
+ <type>FixedTimeStepping</type>
+ <t_initial>0</t_initial>
+ <t_end>1.0</t_end>
+ <timesteps>
+ <pair>
+ <repeat>1</repeat>
+ <delta_t>1e-1</delta_t>
+ </pair>
+ </timesteps>
+ </time_stepping>
+ </process>
+ </processes>
+ <output>
+ <variables>
+ <variable>displacement</variable>
+ <variable>phasefield</variable>
+ <variable>sigma</variable>
+ <variable>epsilon</variable>
+ </variables>
+ <type>VTK</type>
+ <prefix>AT1_iso_compressive</prefix>
+ <timesteps>
+ <pair>
+ <repeat>1</repeat>
+ <each_steps>1</each_steps>
+ </pair>
+ </timesteps>
+ </output>
+ </time_loop>
+ <parameters>
+ <!-- Mechanics -->
+ <parameter>
+ <name>E</name>
+ <type>Constant</type>
+ <value>1.0</value>
+ </parameter>
+ <parameter>
+ <name>nu</name>
+ <type>Constant</type>
+ <value>0.15</value>
+ </parameter>
+ <parameter>
+ <name>k</name>
+ <type>Constant</type>
+ <value>1e-8</value>
+ </parameter>
+ <parameter>
+ <name>gc</name>
+ <type>Constant</type>
+ <value>1.0</value>
+ </parameter>
+ <parameter>
+ <name>ls</name>
+ <type>Constant</type>
+ <value>0.03</value>
+ </parameter>
+ <parameter>
+ <name>rho_sr</name>
+ <type>Constant</type>
+ <value>0.0</value>
+ </parameter>
+ <parameter>
+ <name>displacement0</name>
+ <type>Constant</type>
+ <values>0 0 0</values>
+ </parameter>
+ <parameter>
+ <name>phasefield_ic</name>
+ <type>Constant</type>
+ <value>1</value>
+ </parameter>
+ <parameter>
+ <name>phasefield_bc</name>
+ <type>Constant</type>
+ <value>1</value>
+ </parameter>
+ <parameter>
+ <name>dirichlet0</name>
+ <type>Constant</type>
+ <value>0</value>
+ </parameter>
+ <parameter>
+ <name>Dirichlet_spatial</name>
+ <type>Constant</type>
+ <value>1</value>
+ </parameter>
+ <parameter>
+ <name>dirichlet_right_time</name>
+ <type>CurveScaled</type>
+ <curve>dirichlet_time</curve>
+ <parameter>dirichlet_right</parameter>
+ </parameter>
+ <parameter>
+ <name>dirichlet_right</name>
+ <type>Constant</type>
+ <value>-4.0</value>
+ </parameter>
+ </parameters>
+ <curves>
+ <curve>
+ <name>dirichlet_time</name>
+ <coords>0 0.05 1</coords>
+ <values>0 0.25 1</values>
+ </curve>
+ </curves>
+ <process_variables>
+ <process_variable>
+ <name>phasefield</name>
+ <components>1</components>
+ <order>1</order>
+ <initial_condition>phasefield_ic</initial_condition>
+ <boundary_conditions>
+ <boundary_condition>
+ <geometrical_set>bar</geometrical_set>
+ <geometry>left</geometry>
+ <type>Dirichlet</type>
+ <component>0</component>
+ <parameter>phasefield_bc</parameter>
+ </boundary_condition>
+ <boundary_condition>
+ <geometrical_set>bar</geometrical_set>
+ <geometry>right</geometry>
+ <type>Dirichlet</type>
+ <component>0</component>
+ <parameter>phasefield_bc</parameter>
+ </boundary_condition>
+ <boundary_condition>
+ <geometrical_set>bar</geometrical_set>
+ <geometry>right</geometry>
+ <type>PhaseFieldIrreversibleDamageOracleBoundaryCondition</type>
+ <component>0</component>
+ </boundary_condition>
+ </boundary_conditions>
+ </process_variable>
+ <process_variable>
+ <name>displacement</name>
+ <components>3</components>
+ <order>1</order>
+ <initial_condition>displacement0</initial_condition>
+ <boundary_conditions>
+ <boundary_condition>
+ <geometrical_set>bar</geometrical_set>
+ <geometry>left</geometry>
+ <type>Dirichlet</type>
+ <component>0</component>
+ <parameter>dirichlet0</parameter>
+ </boundary_condition>
+ <boundary_condition>
+ <geometrical_set>bar</geometrical_set>
+ <geometry>right</geometry>
+ <type>Dirichlet</type>
+ <component>0</component>
+ <parameter>dirichlet_right_time</parameter>
+ </boundary_condition>
+ </boundary_conditions>
+ </process_variable>
+ </process_variables>
+ <nonlinear_solvers>
+ <nonlinear_solver>
+ <name>petsc_snes</name>
+ <type>PETScSNES</type>
+ <max_iter>200</max_iter>
+ <linear_solver>linear_solver_d</linear_solver>
+ </nonlinear_solver>
+ <nonlinear_solver>
+ <name>basic_newton_u</name>
+ <type>Newton</type>
+ <max_iter>10000</max_iter>
+ <linear_solver>linear_solver_u</linear_solver>
+ </nonlinear_solver>
+ </nonlinear_solvers>
+ <linear_solvers>
+ <linear_solver>
+ <name>linear_solver_d</name>
+ <eigen>
+ <solver_type>BiCGSTAB</solver_type>
+ <precon_type>ILUT</precon_type>
+ <max_iteration_step>10000</max_iteration_step>
+ <error_tolerance>1e-16</error_tolerance>
+ </eigen>
+ <petsc>
+ <parameters>-ksp_type cg -pc_type bjacobi -ksp_atol 1e-10 -ksp_rtol 1e-10 -snes_type vinewtonrsls -snes_linesearch_type l2 -snes_atol 1.e-8 -snes_rtol 1.e-8 -snes_max_it 1000 -snes_monitor</parameters>
+ </petsc>
+ </linear_solver>
+ <linear_solver>
+ <name>linear_solver_u</name>
+ <petsc>
+ <parameters>-ksp_type cg -pc_type bjacobi -ksp_atol 1e-14 -ksp_rtol 1e-14</parameters>
+ </petsc>
+ </linear_solver>
+ </linear_solvers>
+</OpenGeoSysProject>
diff --git a/Tests/Data/PhaseField/beam_jupyter_notebook/figures/bar_.png b/Tests/Data/PhaseField/beam_jupyter_notebook/figures/bar_.png
new file mode 100644
index 0000000000000000000000000000000000000000..c4ad4d0b56d620173dc713650f685242dd7de3cd
Binary files /dev/null and b/Tests/Data/PhaseField/beam_jupyter_notebook/figures/bar_.png differ
diff --git a/Tests/Data/PhaseField/beam_jupyter_notebook/figures/beam_final.png b/Tests/Data/PhaseField/beam_jupyter_notebook/figures/beam_final.png
new file mode 100644
index 0000000000000000000000000000000000000000..021254e38dbefa26bc86beeaff3b42aa7a981b8f
Binary files /dev/null and b/Tests/Data/PhaseField/beam_jupyter_notebook/figures/beam_final.png differ
diff --git a/Tests/Data/PhaseField/beam_jupyter_notebook/figures/beam_result_with_analytical.png b/Tests/Data/PhaseField/beam_jupyter_notebook/figures/beam_result_with_analytical.png
new file mode 100644
index 0000000000000000000000000000000000000000..de2e2c62061c1a08d0604624ab32158dbb1b9580
Binary files /dev/null and b/Tests/Data/PhaseField/beam_jupyter_notebook/figures/beam_result_with_analytical.png differ
diff --git a/web/content/docs/benchmarks/phase-field/phasefield/Miao_Biot2017.pdf b/web/content/docs/benchmarks/phase-field/phasefield/Miao_Biot2017.pdf
deleted file mode 100755
index 053152c7ecbae216aa541594d7b28e48e73454cf..0000000000000000000000000000000000000000
Binary files a/web/content/docs/benchmarks/phase-field/phasefield/Miao_Biot2017.pdf and /dev/null differ
diff --git a/web/content/docs/benchmarks/phase-field/phasefield/beam.png b/web/content/docs/benchmarks/phase-field/phasefield/beam.png
deleted file mode 100644
index be65c39911a50838bf73948720cde1a016aad182..0000000000000000000000000000000000000000
Binary files a/web/content/docs/benchmarks/phase-field/phasefield/beam.png and /dev/null differ
diff --git a/web/content/docs/benchmarks/phase-field/phasefield/beam_d.png b/web/content/docs/benchmarks/phase-field/phasefield/beam_d.png
deleted file mode 100644
index eeb41ec784f8200ddcdf0a9fda4180bea6bb599a..0000000000000000000000000000000000000000
Binary files a/web/content/docs/benchmarks/phase-field/phasefield/beam_d.png and /dev/null differ
diff --git a/web/content/docs/benchmarks/phase-field/phasefield/beam_u.png b/web/content/docs/benchmarks/phase-field/phasefield/beam_u.png
deleted file mode 100644
index 0965d63074cd13868724b52bee50549b4057e5ac..0000000000000000000000000000000000000000
Binary files a/web/content/docs/benchmarks/phase-field/phasefield/beam_u.png and /dev/null differ
diff --git a/web/content/docs/benchmarks/phase-field/phasefield/error_u.png b/web/content/docs/benchmarks/phase-field/phasefield/error_u.png
deleted file mode 100644
index 678b9ea8d070a559f7fa9c07930150436549d24c..0000000000000000000000000000000000000000
Binary files a/web/content/docs/benchmarks/phase-field/phasefield/error_u.png and /dev/null differ
diff --git a/web/content/docs/benchmarks/phase-field/phasefield/index.md b/web/content/docs/benchmarks/phase-field/phasefield/index.md
deleted file mode 100644
index 5bfa14755b9d917be648e87350ef56cddad33aa5..0000000000000000000000000000000000000000
--- a/web/content/docs/benchmarks/phase-field/phasefield/index.md
+++ /dev/null
@@ -1,44 +0,0 @@
-+++
-author = "Xing-Yuan Miao"
-date = "2017-05-19T09:10:33+01:00"
-title = "Crack beam under tension"
-image = "beam_u.png"
-+++
-
-{{< data-link >}}
-
-## Problem description
-
-We solve a homogeneous beam model under a given displacement loading. The length of the beam is 2\,mm. Detailed model description can refer [this PDF](Miao_Biot2017.pdf).
-
-## Results and evaluation
-
-Results show crack Phase-Field and displacement field distributions through the length of the beam:
-
-{{< img src="beam.png" >}}
-{{< img src="beam_d.png" >}}
-{{< img src="beam_u.png" >}}
-
-For highlighting the deviation between the analytical and numerical solution, we provide the absolute error of the analytical solution and numerical simulation as follows:
-{{< img src="error_u.png" >}}
-
-The analytical solution is:
-$$
-\begin{equation}
-d (x) = 1 - {\mathrm{e}}^{\frac{- |x|}{2 \varepsilon}}
-\end{equation}
-$$
-$$
-\begin{equation}
-u (x) = \dfrac{\sigma}{E} \int_0^x \dfrac{1}{d (x)^2 + k} \mathrm{d}x
-\end{equation}
-$$
-with
-$$
-\begin{equation}
-\sigma = \dfrac{E u_0}{I (\varepsilon, k)}
-\end{equation}
-$$
-\begin{equation}
-I (\varepsilon, k) = \int_0^1 \dfrac{1}{d (x)^2 + k} \mathrm{d}x
-\end{equation}