Figure 9

Imports

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

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.4,         # 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 = -2.0    # kPa
g  = +9.81e3 # mm/s2

Computing strain fields in prone and supine positions

invariants = {}

gravity_lst = [-1.,+1.]
for gravity_cste in gravity_lst:

    ### computing the unloaded configuration from a generic end-exhalation configuration
    Uref, phisref_computation, dV_expi_ini = dmech.run_RivlinCube_PoroHyperelasticity(
        inverse=1,
        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":gravity_cste*g, "P0":float(pe)},
        res_basename="Fig9-unloaded",
        get_results=1,
        verbose=1)

    phisref_imposed = [numpy.random.uniform(low=0.4, high=0.6) for i in range(len(phisref_computation))]

    ### 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":"function_xml_from_array", "val":phisref_imposed},
        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":gravity_cste*g, "P0":float(pe)},
        res_basename="Fig9-exhalation",
        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":"function_xml_from_array", "val":phisref_imposed},
        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":gravity_cste*g, "P0":float(pi)},
        res_basename="Fig9-inhalation",
        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, lognorm=False)
    invariants["g="+str(gravity_cste)] = invariants_mesh

Post-processing

Reading displacement field from [Patte et al.] and computing associated invariants

mesh = dolfin.Mesh()
dolfin.XDMFFile("Data/patte_mesh.xdmf").read(mesh)
fe = dolfin.VectorElement(
    family="CG",
    cell=mesh.ufl_cell(),
    degree=1)
U_fs= dolfin.FunctionSpace(
    mesh,
    fe)

U_inhal_patte = dolfin.Function(U_fs, name="U")
U_exhal_patte = dolfin.Function(U_fs, name="U")

mesh_meshio = meshio.read("Data/patte_mesh.xdmf")
u_meshio = mesh_meshio.point_data["U"]
u_meshio = u_meshio.tolist()
u_meshio = [item for sublist in u_meshio for item in sublist] 

c = 0
for dof in dolfin.vertex_to_dof_map(U_fs):
    U_inhal_patte.vector()[dof] = u_meshio[c]
    U_exhal_patte.vector()[dof] = 0.
    c += 1

invariants_mesh_patte = get_invariants.get_invariants(U_exhal=U_exhal_patte, U_inhal=U_inhal_patte, mesh=mesh, lognorm=False)

Generating plots

I1

fig, ax = plt.subplots()

### plotting parameters
plt.rc('xtick', labelsize=16) 
plt.rc('ytick', labelsize=16) 
plt.rc('axes', labelsize=16)
plt.ylim([1, 1.9])
plt.xlim([1, 10])
fig.set_size_inches(6, 6)
plt.legend(loc="upper right", fontsize=16)
plt.grid()

### curves
plt.plot(invariants_mesh_patte["zone"], invariants_mesh_patte["I1"], color="green", linestyle="-", marker="s", label="Patte et al., 2022")
plt.plot(invariants["g=1.0"]["zone"]  , invariants["g=1.0"]["I1"]  , color="blue" , linestyle="-", marker="s", label="Supine position"   )
plt.plot(invariants["g=-1.0"]["zone"] , invariants["g=-1.0"]["I1"] , color="red"  , linestyle="-", marker="s", label="Prone position"    )

### axis name
ax.set_ylabel("Average of the first invariant per zone")
ax.set_xlabel("Zone")

plt.savefig("fig9-I1.pdf", bbox_inches='tight')
plt.show()

I2

fig, ax = plt.subplots()

### plotting parameters
plt.rc('xtick', labelsize=16) 
plt.rc('ytick', labelsize=16) 
plt.rc('axes', labelsize=16)
plt.ylim([1, 1.9])
plt.xlim([1, 10])
fig.set_size_inches(6, 6)
plt.legend(loc="upper right", fontsize=16)
plt.grid()

### curves
plt.plot(invariants_mesh_patte["zone"], invariants_mesh_patte["I2"], color="green", linestyle="-", marker="s", label="Patte et al., 2022")
plt.plot(invariants["g=1.0"]["zone"]  , invariants["g=1.0"]["I2"]  , color="blue" , linestyle="-", marker="s", label="Supine position"   )
plt.plot(invariants["g=-1.0"]["zone"] , invariants["g=-1.0"]["I2"] , color="red"  , linestyle="-", marker="s", label="Prone position"    )

### axis name
ax.set_ylabel("Average of the first invariant per zone")
ax.set_xlabel("Zone")

plt.savefig("fig9-I1.pdf", bbox_inches='tight')
plt.show()

J

fig, ax = plt.subplots()

### plotting parameters
plt.rc('xtick', labelsize=16) 
plt.rc('ytick', labelsize=16) 
plt.rc('axes', labelsize=16)
plt.ylim([1, 1.9])
plt.xlim([1, 10])
fig.set_size_inches(6, 6)
plt.legend(loc="upper right", fontsize=16)
plt.grid()

### curves
plt.plot(invariants_mesh_patte["zone"], invariants_mesh_patte["J"], color="green", linestyle="-", marker="s", label="Patte et al., 2022")
plt.plot(invariants["g=1.0"]["zone"]  , invariants["g=1.0"]["J"]  , color="blue" , linestyle="-", marker="s", label="Supine position"   )
plt.plot(invariants["g=-1.0"]["zone"] , invariants["g=-1.0"]["J"] , color="red"  , linestyle="-", marker="s", label="Prone position"    )

### axis name
ax.set_ylabel("Average of the first invariant per zone")
ax.set_xlabel("Zone")

plt.savefig("fig9-I1.pdf", bbox_inches='tight')
plt.show()