Figure 8

Contents

Figure 8#

Imports#

import dolfin
import matplotlib.pyplot as plt
import numpy
import pandas

import dolfin_mech as dmech
import get_invariants

Parameters#

Mesh#

cube_params = {"path_and_mesh_name":"Data/generic_mesh.xdmf"}

Material#

params = {
    "alpha":0.16,     # kPa
    "gamma":0.5,      # [-]
    "c1":0.8,         # kPa
    "c2":0.2,         # kPa
    "kappa":1e2,      # kPa
    "eta":1e-5,       # kPa
    "rho_solid":1e-6} # g/mm3
mat_params = {"scaling":"linear", "parameters":params}

Loading#

pe = -0.5    # kPa
pi = -3.0    # kPa
g  = +9.81e3 # mm/s2

h_lst = [0., -0.005, +0.005, -0.0065, +0.0065]

Computing strain fields#

invariants = {}

for h in h_lst:

    ### computing the unloaded configuration from a generic end-exhalation configuration
    Uref, phisref_computation, dV_expi_ini = dmech.run_RivlinCube_PoroHyperelasticity(
        inverse=1,
        porosity_params={"type":"mesh_function_random_xml", "val":0.3},
        cube_params=cube_params,
        mat_params=mat_params,
        inertia_params={"applied":True, "rho_val":1e-8},
        step_params={"dt_min":1e-4, "dt_ini":1},
        load_params={"type":"p_boundary_condition0", "f":g, "P0":float(pe), "H":h},
        res_basename="Fig8-unloaded"+str(h),
        get_results=1,
        verbose=1)

    ### computing the end-exhalation configuration
    U_exhal, phisexhal, dV_unloaded = dmech.run_RivlinCube_PoroHyperelasticity(
        cube_params=cube_params,
        move_params = {"move":True, "U":Uref},
        porosity_params={"type":"mesh_function_random_xml", "val":0.5},
        mat_params=mat_params,
        inertia_params={"applied":True, "rho_val":1e-8},
        step_params={"dt_min":1e-4, "dt_ini":0.125},
        load_params={"type":"p_boundary_condition", "f":g, "P0":float(pe), "H":h},
        res_basename="Fig8-exhalation"+str(h),
        get_results=1,
        verbose=1)

    ### computing the end-inhalation configuration
    U_inhal, phisinhal, dV_unloaded = dmech.run_RivlinCube_PoroHyperelasticity(
        cube_params=cube_params,
        move_params = {"move":True, "U":Uref},
        porosity_params={"type":"mesh_function_random_xml", "val":0.5},
        mat_params=mat_params,
        inertia_params={"applied":True, "rho_val":1e-8},
        step_params={"dt_min":1e-4, "dt_ini":0.125},
        load_params={"type":"p_boundary_condition", "f":g, "P0":float(pi), "H":h},
        res_basename="Fig8-inhalation"+str(h),
        get_results=1,
        verbose=1)

    mesh = dolfin.Mesh()
    mesh_name = str(cube_params["path_and_mesh_name"])
    dolfin.XDMFFile(mesh_name).read(mesh)
    dolfin.ALE.move(mesh, Uref)

    invariants_mesh = get_invariants.get_invariants(U_exhal=U_exhal, U_inhal=U_inhal, mesh=mesh)
    invariants["h="+str(h)] = invariants_mesh

Generating plots#

[Hurtado et al.]’s data#

J  = pandas.read_table("Data/hurtado_J.dat" , sep="\s+", usecols=["zone", "mean", "mean+std", "mean-std"]).to_dict("list")
I1 = pandas.read_table("Data/hurtado_I1.dat", sep="\s+", usecols=["zone", "mean", "mean+std", "mean-std"]).to_dict("list")
I2 = pandas.read_table("Data/hurtado_I2.dat", sep="\s+", usecols=["zone", "mean", "mean+std", "mean-std"]).to_dict("list")

I1#

fig, ax = plt.subplots()

### plotting parameters
plt.rc("xtick", labelsize=20)
plt.rc("ytick", labelsize=20)
ax.set_xticks(numpy.linspace(1, 10, 10))
plt.xlim([0.9, 10.1])
fig.set_size_inches(6, 6)
plt.grid()

### axis labels
ax.set_xlabel("Zone", fontsize=20, labelpad=10)
ax.set_ylabel("Log Normalized First Invariant", fontsize=20, labelpad=10)

### error bars
down_up_lst = [[abs(J["mean-std"][i]-J["mean"][i])for i in range(len(J["zone"]))], [abs(J["mean+std"][i]-J["mean"][i]) for i in range(len(J["zone"]))]]
mean_lst=[J["mean"][i] for i in range(len(J["zone"]))]
plt.errorbar(invariants["h=0.0"]["zone"], mean_lst, yerr=down_up_lst, capsize=4, fmt="s", ecolor = "black", markeredgecolor='black', markerfacecolor='black')

### curves
plt.plot(invariants["h=0.0"]["zone"]    , invariants["h=0.0"]["I1^"]    , color="green"        , linestyle="-", marker="s", label="Model (w gravity), H=0"  )
plt.plot(invariants["h=0.005"]["zone"]  , invariants["h=0.005"]["I1^"]  , color="red"          , linestyle="-", marker="s", label="Model (w gravity), H=+H1")
plt.plot(invariants["h=-0.005"]["zone"] , invariants["h=-0.005"]["I1^"] , color="darkturquoise", linestyle="-", marker="s", label="Model (w gravity), H=–H1")
plt.plot(invariants["h=0.0065"]["zone"] , invariants["h=0.0065"]["I1^"] , color="purple"       , linestyle="-", marker="s", label="Model (w gravity), H=+H2")
plt.plot(invariants["h=-0.0065"]["zone"], invariants["h=-0.0065"]["I1^"], color="orange"       , linestyle="-", marker="s", label="Model (w gravity), H=–H2")
plt.legend(loc="upper right", fontsize=14)
plt.show()

I2#

fig, ax = plt.subplots()

### plotting parameters
plt.rc("xtick", labelsize=20)
plt.rc("ytick", labelsize=20)
plt.xlim([0.9, 10.1])
fig.set_size_inches(6, 6)
ax.set_xticks(numpy.linspace(1, 10, 10))
plt.grid()

### axis labels
ax.set_xlabel("Zone", fontsize=20, labelpad=10)
ax.set_ylabel("Log Normalized Second Invariant", fontsize=20, labelpad=10)

### error bars
down_up_lst = [[abs(J["mean-std"][i]-J["mean"][i])for i in range(len(J["zone"]))], [abs(J["mean+std"][i]-J["mean"][i]) for i in range(len(J["zone"]))]]
mean_lst=[J["mean"][i] for i in range(len(J["zone"]))]
plt.errorbar(invariants["h=0.0"]["zone"], mean_lst, yerr=down_up_lst, capsize=4, fmt="s", ecolor = "black", markeredgecolor='black', markerfacecolor='black')

### curves
plt.plot(invariants["h=0.0"]["zone"]    , invariants["h=0.0"]["I2^"]    , color="green"        , linestyle="-", marker="s", label="Model (w gravity), H=0"  )
plt.plot(invariants["h=0.005"]["zone"]  , invariants["h=0.005"]["I2^"]  , color="red"          , linestyle="-", marker="s", label="Model (w gravity), H=+H1")
plt.plot(invariants["h=-0.005"]["zone"] , invariants["h=-0.005"]["I2^"] , color="darkturquoise", linestyle="-", marker="s", label="Model (w gravity), H=–H1")
plt.plot(invariants["h=0.0065"]["zone"] , invariants["h=0.0065"]["I2^"] , color="purple"       , linestyle="-", marker="s", label="Model (w gravity), H=+H2")
plt.plot(invariants["h=-0.0065"]["zone"], invariants["h=-0.0065"]["I2^"], color="orange"       , linestyle="-", marker="s", label="Model (w gravity), H=–H2")
plt.legend(loc="upper right", fontsize=14)
plt.show()

J#

fig, ax = plt.subplots()

### plotting parameters
plt.rc("xtick", labelsize=20)
plt.rc("ytick", labelsize=20)
ax.set_xticks(numpy.linspace(1, 10, 10))
plt.xlim([0.9, 10.1])
fig.set_size_inches(6, 6)
plt.grid()

### axis labels
ax.set_xlabel("Zone", fontsize=20, labelpad=10)
ax.set_ylabel("Log Normalized Third Invariant", fontsize=20, labelpad=10)

### error bars
down_up_lst = [[abs(J["mean-std"][i]-J["mean"][i])for i in range(len(J["zone"]))], [abs(J["mean+std"][i]-J["mean"][i]) for i in range(len(J["zone"]))]]
mean_lst=[J["mean"][i] for i in range(len(J["zone"]))]
plt.errorbar(invariants["h=0.0"]["zone"], mean_lst, yerr=down_up_lst, capsize=4, fmt="s", ecolor = "black", markeredgecolor='black', markerfacecolor='black')

### curves
plt.plot(invariants["h=0.0"]["zone"]    , invariants["h=0.0"]["J^"]    , color="green"        , linestyle="-", marker="s", label="Model (w gravity), H=0"  )
plt.plot(invariants["h=0.005"]["zone"]  , invariants["h=0.005"]["J^"]  , color="red"          , linestyle="-", marker="s", label="Model (w gravity), H=+H1")
plt.plot(invariants["h=-0.005"]["zone"] , invariants["h=-0.005"]["J^"] , color="darkturquoise", linestyle="-", marker="s", label="Model (w gravity), H=–H1")
plt.plot(invariants["h=0.0065"]["zone"] , invariants["h=0.0065"]["J^"] , color="purple"       , linestyle="-", marker="s", label="Model (w gravity), H=+H2")
plt.plot(invariants["h=-0.0065"]["zone"], invariants["h=-0.0065"]["J^"], color="orange"       , linestyle="-", marker="s", label="Model (w gravity), H=–H2")
plt.legend(loc="upper right", fontsize=14)
plt.show()