Figure 6

Imports

import dolfin
import numpy
import scipy.optimize

from scipy.interpolate import interp1d

import dolfin_mech                    as dmech
import matplotlib.pyplot              as plt
import micro_poro_structure_generator as gen

Importing experimental data

smith_PV_inflation_gamma_0   = numpy.load('smith_PV_inflation_gamma_0.npy')
p_smith_PV_inflation_gamma_0 = smith_PV_inflation_gamma_0[:, 0]
v_smith_PV_inflation_gamma_0 = smith_PV_inflation_gamma_0[:, 1]

smith_PV_deflation_gamma_0   = numpy.load('smith_PV_deflation_gamma_0.npy')
p_smith_PV_deflation_gamma_0 = smith_PV_deflation_gamma_0[:, 0]
v_smith_PV_deflation_gamma_0 = smith_PV_deflation_gamma_0[:, 1]

smith_PV_deflation_gamma_0

Defining geometry

seeds_filename    = "Fig6-seeds.dat"
mesh_filebasename = "Fig6-mesh"
qois_filename     = "Fig6-qois.dat"
res_basename      = "Fig6"

domain_y  = 0.1 * 0.8
domain_x  = domain_y * numpy.sqrt(3)/1.5/2
thickness = 0.012 * 0.8

gen.generate_seeds_semi_regular(
    DoI = 0.,
    row = 1,
    domain_y = domain_y,
    seeds_filename = seeds_filename)
gen.generate_mesh_2D_rectangle_w_voronoi_inclusions(
    mesh_filename = mesh_filebasename,
    seeds_filename = seeds_filename,
    h = thickness,
    lcar = thickness/5,
    domain_x = domain_x,
    domain_y = domain_y,
    shift_y = 0.,
    remove_seeds = True)

mesh = dolfin.Mesh()
dolfin.XDMFFile(mesh_filebasename+".xdmf").read(mesh)
dV = dolfin.Measure("dx",domain=mesh)

coord = mesh.coordinates()
xmax = max(coord[:,0]); xmin = min(coord[:,0])
ymax = max(coord[:,1]); ymin = min(coord[:,1])

V   = (xmax - xmin)*(ymax - ymin)
VS0 = dolfin.assemble(dolfin.Constant(1) * dV)
Vf0 = V - VS0

Loading

load_params = {}
load_params["pf"] = 4
load_params["sigma_bar_00"] = 0.0
load_params["sigma_bar_11"] = 0.0
load_params["sigma_bar_01"] = 0.0
load_params["sigma_bar_10"] = 0.0

Identifier function

def parameter_identifier(exp_data, params_initial, bnds, asym_slope):

    p_exp = exp_data[:, 0]
    v_exp = exp_data[:, 1]

    def J_cost(params_dimless):
        params = [params_initial[0]*(1 + params_dimless[0]), params_initial[1]*(1 + params_dimless[1]), params_initial[2]*(1 + params_dimless[2]), params_initial[3]*(1 + params_dimless[3])]
        print(params)


        mat_params = {"model":"exponentialneoHookean", "parameters":{"beta1":params[0], "beta2":params[1], "beta3":params[2], "beta4":100*params[0], "alpha":params[3]}}  


        dmech.run_HollowBox_MicroPoroHyperelasticity(
            dim=2,
            mesh=mesh,
            mat_params=mat_params,
            load_params=load_params,
            step_params={"Deltat":1., "dt_ini":0.1, "dt_min":0.005, "dt_max":0.1},
            res_basename=res_basename,
            write_qois_limited_precision=False,
            verbose=1
        )

        qois_vals = numpy.loadtxt(qois_filename)
        qois_name_list = open(qois_filename).readline().split()
        pf_lst = qois_vals[:, qois_name_list.index("p_f") - 1]*10.20
        vf_lst = qois_vals[:, qois_name_list.index("vf") - 1]

        for i in range(1, len(vf_lst)):
            slope = (vf_lst[i] - vf_lst[i - 1])/(pf_lst[i] - pf_lst[i - 1])
            if slope < asym_slope:
                break

        vf_asym = vf_lst[i]
        vf_lst = [vf_/vf_asym *100 for vf_ in vf_lst]

        model_interpolator = interp1d(pf_lst, vf_lst, kind='cubic')  

        JC = 0
        for i in range(len(exp_data)):
            JC += (v_exp[i] - model_interpolator(p_exp[i]))**2

        print("JC: " +str(JC))
        return JC


    params_ini_dimless = [0, 0,  0, 0]
    params_id= (scipy.optimize.minimize(J_cost, x0=params_ini_dimless, method="TNC", bounds=bnds, tol=1e-3)).x

    params_id = [params_initial[0]*(1 + params_id[0]), params_initial[1]*(1 + params_id[1]), params_initial[2]*(1 + params_id[2]), params_initial[3]*(1 + params_id[3])]

    return params_id

params_initial = [0.12414156917794349, 0.011853212744289299, 0.6282749906141952, 3.124581670582977]
bnds = [(-0.9, 10), (-0.8, 10), (-0.2, 0.2), (-0.3, 4)]
asym_slope = 0.1 * Vf0
params_id_simth =  parameter_identifier(numpy.vstack((smith_PV_deflation_gamma_0, smith_PV_inflation_gamma_0)), params_initial, bnds, asym_slope)

Model response

def graph_printer_smith(params):
    mat_params = {"model":"exponentialneoHookean", "parameters":{"beta1":params[0], "beta2":params[1], "beta3":params[2], "beta4":100*params[0], "alpha":params[3]}}    
    dmech.run_HollowBox_MicroPoroHyperelasticity(
        dim=2,
        mesh=mesh,
        mat_params=mat_params,
        load_params=load_params,
        step_params={"Deltat":1., "dt_ini":0.1, "dt_min":0.005, "dt_max":0.1},
        res_basename=res_basename,
        write_qois_limited_precision=False,
        verbose=1
    )

    qois_vals = numpy.loadtxt(qois_filename)
    qois_name_list = open(qois_filename).readline().split()
    pf_lst = qois_vals[:, qois_name_list.index("p_f") - 1]*10.20
    vf_lst = qois_vals[:, qois_name_list.index("vf") - 1]

    for i in range(1, len(vf_lst)):
        slope = (vf_lst[i] - vf_lst[i - 1])/(pf_lst[i] - pf_lst[i - 1])
        if slope < 0.1 * Vf0:
            break

    vf_asym = vf_lst[i]
    vf_lst = [vf_/vf_asym *100 for vf_ in vf_lst]

    print("pf_asym: {:.2f} cmH2O".format(pf_lst[i]))

    plt.plot()
    plt.rc('xtick' , labelsize=14)
    plt.rc('ytick' , labelsize=14)
    plt.rc('legend', fontsize=12)
    plt.xlabel(r'$p_f~(cm H_2O)$', fontsize=16)
    plt.ylabel(r'$Volume~(\% TLC)$', fontsize=16)


    plt.plot(p_smith_PV_deflation_gamma_0, v_smith_PV_deflation_gamma_0, '#D94801', label='[Smith, 1986], deflation, $\gamma = 0~dyn/cm$', linestyle='dashed')
    plt.plot(p_smith_PV_inflation_gamma_0, v_smith_PV_inflation_gamma_0, '#D94801', label='[Smith, 1986], inflation, $\gamma = 0~dyn/cm$')

    plt.plot(pf_lst[0:80], vf_lst[0:80], '#084594', label='Model, $\gamma = 0~dyn/cm$')

    plt.xlim(0, 25)
    plt.ylim(0, 110)
    plt.legend(loc = 'lower right', fontsize=12, shadow=True)
    plt.savefig('Smith_iden.pdf',bbox_inches='tight')
    plt.savefig("Smith_iden.tif", format="tiff", dpi=300, bbox_inches="tight")
    plt.show()

Plot

graph_printer_smith(params_id_simth)