import dolfin
import numpy

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

qois_filename     = "Fig8-qois.dat"
res_basename      = "Fig8"
seeds_filename    = "ch8.dat"
mesh_filebasename = "ch8-mesh"


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

def graph_printer(thickness, params, linestyle):
    
    domain_y = 1
    domain_x = domain_y * numpy.sqrt(3)/1.5/2
    # thickness = 0.12

    gen.generate_seeds_semi_regular(
        DoI = 0.,
        row = 1,
        domain_y = domain_y,
        seeds_filename = seeds_filename)
    phi = 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

    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.001, "dt_max":0.005},
        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.0005:
            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(pf_lst, vf_lst, '#084594', label='$\Phi_{f0} = '+str('{:.2f}'.format(phi * 100))+'\%$', linestyle=linestyle)

    plt.xlim(0, 25)
    plt.legend(loc = 'lower right', fontsize=12, shadow=True)

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)

params = [0.08855929243285596, 0.011039510924095856, 0.6281487879627474, 3.409513378002055]
graph_printer(0.18, params, 'dashdot')
graph_printer(0.12, params, None)
graph_printer(0.09, params, 'dashed')
plt.savefig('Porosity_compare.pdf',bbox_inches='tight')
plt.savefig("Porosity_compare.tif", format="tiff", dpi=300, bbox_inches="tight")
plt.show()

Figure 8¶

Imports¶

Loading¶

Printing function¶

Running the printing function for different porosities¶