Figure 4

Figure 4#

Imports#

import copy
import matplotlib.pyplot as plt
import numpy

import dolfin_mech as dmech

Parameters#

Mesh#

cube_params = {"path_and_mesh_name":"Data/generic_mesh.xdmf", "refine":True} # the mesh is refined only for plotting purposes

Material#

params = {
    "alpha":0.16,     # kPa
    "gamma":0.5,      # [-]
    "c1":0.2,         # kPa
    "c2":0.4,         # 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 porosity distributions#

gravity_lst = [0,1]
for gravity_ in gravity_lst:

    ### computing the unloaded configuration from a generic end-exhalation configuration
    Uref, phisref_computation, dVexpiini = dmech.run_RivlinCube_PoroHyperelasticity(
        inverse=1,
        cube_params=cube_params,
        mat_params=mat_params,
        porosity_params={"type":"mesh_function_random_xml", "val":0.3},
        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_*g, "P0":float(pe)},
        res_basename="Fig4-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
    Uexhal, phisexhal, dVunloaded = 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_*g, "P0":float(pe)},
        res_basename="Fig4-exhalation",
        get_results=1,
        verbose=1)

    ### computing the end-inhalation configuration
    Uinhal, phisinhal, dVunloaded = 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_*g, "P0":float(pi)},
        res_basename="Fig4-inhalation",
        get_results=1,
        verbose=1)

    if (gravity_ == 0):
        phisexhal_g0 = copy.deepcopy(phisexhal)
        phisinhal_g0 = copy.deepcopy(phisinhal)

porosity_lst = numpy.linspace(0, 1, 300)

porosity_exhal_g0 = []
porosity_inhal_g0 = []
porosity_exhal_g  = []
porosity_inhal_g  = []
porosity_ref      = []
porosity_plot     = []
for c in range(0, len(porosity_lst)-1):
    min = porosity_lst[c]
    max = porosity_lst[c+1]
    porosity_exhal_g0.append(numpy.sum([min<=i<max for i in phisexhal_g0]))
    porosity_inhal_g0.append(numpy.sum([min<=i<max for i in phisinhal_g0]))
    porosity_exhal_g.append(numpy.sum([min<=i<max for i in phisexhal]))
    porosity_inhal_g.append(numpy.sum([min<=i<max for i in phisinhal]))
    porosity_ref.append(numpy.sum([min<=i<max for i in phisref_imposed]))
    porosity_plot.append(1-(min+(max-min)/2))

Generating plots#

End-exhalation#

fig, ax = plt.subplots()

### plotting parameters
plt.rc('xtick', labelsize=18) 
plt.rc('ytick', labelsize=18) 
plt.rc('axes', labelsize=18)
fig.set_size_inches(8, 6)
plt.xlim([0.3, 0.9])
plt.ylim([0., 14000])
ax.set_ylabel("Frequency", fontsize=22, labelpad=10)
ax.set_xlabel("Porosity", fontsize=22, labelpad=10)

##### plot results
width = max-min
ax.bar(porosity_plot, porosity_ref    , color="green" , label='$\Phi_{f0}$'  , alpha=0.5, edgecolor="black", width=width)
ax.bar(porosity_plot, porosity_exhal_g, color="tomato", label="$\Phi_{fe}$ w/gravity"  , alpha=0.5, edgecolor="black", width=width)
ax.bar(porosity_plot, porosity_exhal_g0, color="blue" , label="$\Phi_{fe}$ wo/gravity", alpha=0.4, edgecolor="black", width=width)
ax.legend()

plt.legend(loc="upper left", fontsize=18)
plt.show()

End-inhalation#

fig, ax = plt.subplots()

### plotting parameters
plt.rc('xtick', labelsize=18) 
plt.rc('ytick', labelsize=18) 
plt.rc('axes', labelsize=18)
fig.set_size_inches(8, 6)
plt.xlim([0.3, 0.9])
plt.ylim([0., 14000])
ax.set_ylabel("Frequency", fontsize=22, labelpad=10)
ax.set_xlabel("Porosity", fontsize=22, labelpad=10)

width = max-min
ax.bar(porosity_plot, porosity_ref     , color="green" , label='$\Phi_{f0}$'  , alpha=0.5, edgecolor="black", width=width)
ax.bar(porosity_plot, porosity_inhal_g , color="tomato", label='$\Phi_{fi}$ w/gravity'  , alpha=0.5, edgecolor="black", width=width)
ax.bar(porosity_plot, porosity_inhal_g0, color="blue"  , label='$\Phi_{fi}$ wo/gravity', alpha=0.4, edgecolor="black", width=width)
ax.legend()

plt.legend(loc="upper left", fontsize=18)
plt.show()
Contents