LargeDeformations: Torsion of a hyperelastic bar

import os
from pathlib import Path

import matplotlib.pyplot as plt
import numpy as np
import ogstools as ot
import pandas as pd
import pyvista as pv
from ogstools.logparser import fill_ogs_context, parse_file

pv.set_plot_theme("document")
pv.set_jupyter_backend("static")

# Some plot settings
plt.rcParams["lines.linewidth"] = 2.0
plt.rcParams["lines.color"] = "black"
plt.rcParams["legend.frameon"] = True
plt.rcParams["font.family"] = "serif"
plt.rcParams["legend.fontsize"] = 14
plt.rcParams["font.size"] = 14
plt.rcParams["axes.spines.right"] = False
plt.rcParams["axes.spines.top"] = False
plt.rcParams["axes.spines.left"] = True
plt.rcParams["axes.spines.bottom"] = True
plt.rcParams["axes.axisbelow"] = True
plt.rcParams["figure.figsize"] = (8, 6)

def R(x, y):
    return np.sqrt(x**2 + y**2)


def phi(x, y):
    theta = np.arctan2(y, x)
    return np.where(theta < 0, theta + 2 * np.pi, theta)


def u(x, y):
    ux = R(x, y) * np.cos(phi(x, y) + np.pi / 20) - x
    uy = R(x, y) * np.sin(phi(x, y) + np.pi / 20) - y
    return [ux, uy]

x = np.linspace(-0.5, 0.5, 10)
y = x.copy()
X, Y = np.meshgrid(x, y)
fig, ax = plt.subplots()
ax.plot(X, Y, marker="s", color="red", ls="", alpha=0.5)
ax.quiver(X, Y, u(X, Y)[0], u(X, Y)[1], pivot="mid")
ax.set_aspect("equal")
ax.set_xlabel("$x$ / m")
ax.set_ylabel("$y$ / m")
ax.set_title("Displacement-controlled torsion on face of prismatic bar")
fig.tight_layout()

out_dir = Path(os.environ.get("OGS_TESTRUNNER_OUT_DIR", "_out"))
if not out_dir.exists():
    out_dir.mkdir(parents=True)

model = ot.Project(input_file="bar1to6_torsion.prj", output_file="bar1to6_torsion.prj")

model.run_model(logfile=f"{out_dir}/out.txt", args=f"-o {out_dir}")

reader = pv.get_reader(f"{out_dir}/bar1to6_torsion.pvd")
print(reader.time_values)
reader.set_active_time_value(0.05)
mesh = reader.read()[0]

plotter = pv.Plotter(shape=(1, 2), window_size=[1000, 500])

warped = mesh.warp_by_vector("displacement")
plotter.subplot(0, 0)
plotter.add_mesh(mesh, show_edges=True, show_scalar_bar=False, color=None, scalars=None)
plotter.show_bounds(
    ticks="outside", xtitle="x / m", ytitle="y / m", ztitle="z / m", font_size=10
)
plotter.add_axes()
plotter.view_isometric()
plotter.add_text("undeformed", font_size=10)

plotter.subplot(0, 1)
plotter.add_mesh(
    warped, show_edges=True, show_scalar_bar=None, color=None, scalars="displacement"
)
plotter.view_isometric()
plotter.add_text("deformed", font_size=10)
plotter.show()

log_file_raw = parse_file(f"{out_dir}/out.txt")
log_file_df = pd.DataFrame(log_file_raw)
ogs_log_context = fill_ogs_context(log_file_df)

last_ts = ogs_log_context["time_step"].to_numpy()[-1]
fig, ax = plt.subplots()
for i in range(1, last_ts + 1):
    timestep = ogs_log_context[ogs_log_context["time_step"] == i]
    ax.plot(
        timestep["iteration_number"],
        timestep["dx"],
        label=f"timestep {i}",
        color="darkblue",
        alpha=i / (last_ts + 1),
    )

ax.axhline(1e-13, ls="--", label="convergence criterion")
ax.plot([7, 8, 7, 7], [1e-6, 1e-12, 1e-12, 1e-6], color="black")
ax.text(7.3, 2e-12, "1")
ax.text(6.7, 1e-9, "2")
ax.set_yscale("log")
ax.set_xlabel("iteration number")
ax.set_ylabel("$|| \\Delta \\vec{u}$ || / m")
ax.legend()
fig.tight_layout()