Figure 10

Figure 10#

Imports#

import matplotlib
import matplotlib.pyplot as plt
import numpy

import dolfin_estim as destim
import dolfin_mech as dmech

### disable deprecation warning to avoid heavy output
import warnings
from ffc.quadrature.deprecation import QuadratureRepresentationDeprecationWarning
warnings.simplefilter("ignore", QuadratureRepresentationDeprecationWarning)

Parameters#

Geometry and varying parameters#

mesh_size = 0.1 ### mesh size for the resolution
l = 0.6 ### characteristic length of the cube
delta_lst = [1000, 10, 5, 2.5, 1.7, 1.1, 0.9, 0.099] ### list of wave lengths chosen 
noise_lst = [0.0, 0.01, 0.05, 0.1] ### noise levels chosen

Material behavior#

params = {
    "E":1.,        # kPa
    "nu":0.3}      # [-]
mat_params = {"model":"CGNH", "parameters":params} ### hyperelastic law

Loading#

load_params_boundary = {"type":"pres", "f":0.3} ### studied in the case of a boundary force, f = 3 kPa

Generating plots#

Helpers#

cube_params = {"X0":0.2, "Y0":0.2, "X1":0.8, "Y1":0.8, "l":0.01} # exact solution --without noise-- computed on a very fine mesh for better resolution, computed here and not in dolfin_estim in order to avoid heavy computations

u, v = dmech.run_RivlinCube_Hyperelasticity(
    dim          = 2,
    cube_params  = cube_params,
    mat_params   = mat_params,
    step_params  = {"dt_ini":1/20},
    const_params = {"type":"blox"},
    load_params  = load_params_boundary,
    get_results  = 1,
    res_basename = "./")

SNR_lst = [] ### defining SNRs -- Signal-to-Noise ratios--
for noise in noise_lst:
    if noise == 0.:
        SNR_lst.append(40.) ### setting the SNR arbitrarily when should be +∞
    else:
        SNR_lst.append(1/noise)

def writing_results_to_pdf(results_all = [], delta_lst = [], SNR_lst = []): ### function to write results --boxplots--
    
    fig, ax = plt.subplots()
    ### plotting parameters
    plt.rc("xtick", labelsize = 20)
    plt.rc("ytick", labelsize = 20)
    plt.rc("legend", fontsize = 12)
    plt.xlabel("Signal to Noise Ratio (SNR)", fontsize = 12)
    plt.ylabel("Estimation error (%)", fontsize = 12)
    fig.set_size_inches(10, 6)
    plt.grid()
    plt.xlim([5, 2000.])
    plt.ylim([-100, 100.])

    nb_deltas = len(delta_lst)
    spacing = numpy.linspace(-0.0305*(nb_deltas//2), 0.0305*(nb_deltas//2), nb_deltas) ### spacing between the different boxplots
    color_lst = matplotlib.cm.viridis(numpy.linspace(0, 1, len(delta_lst))) #### colors in the viridis scale

    count_spacing = 0
    for delta in delta_lst:
        results_delta = results_all[str(delta_lst[count_spacing])] ### getting results for each value of delta, stored into the dictionary "results_delta"
        down_up_lst = [[abs(results_delta["E_-"][i]-results_delta["E_average"][i])for i in range(len(results_delta["noise"]))], [abs(results_delta["E_+"][i]-results_delta["E_average"][i]) for i in range(len(results_delta["noise"]))]] #### to create y-error bars, defining up and down whisker
        mean_lst = [results_delta["E_average"][i] for i in range(len(results_delta["noise"]))]
        SNR_lst_scaled = [SNR_lst[i]*numpy.exp((spacing[count_spacing]))**numpy.log(10) for i in range(len(SNR_lst))] ### to include spacing, that will remain even with the log scale
        plt.errorbar(SNR_lst_scaled, mean_lst, yerr = down_up_lst, capsize = 4, fmt = "s", ecolor = color_lst[-1], markeredgecolor = color_lst[-1], markerfacecolor = color_lst[-1], label = r'$\Delta$'+"  = " +'{0:.1f}'.format(delta)+"l")
        plt.gca().set_xscale('log') ### plotting y-error bars
        color_lst = color_lst[:-1] ### different color for each error bar
        count_spacing+= 1
    ax.set_xticks([])
    ax.set_xticks([], minor = True)
    plt.xticks(SNR_lst, ['$\infty$', 100, 20, 10]) 
    plt.legend(loc = "upper right", fontsize = 13, ncol = 2)

    plt.savefig("./different_values_beta_VFM_plane_waves.pdf", bbox_inches = 'tight')
    plt.show()

Plots#

results_all = {}
for delta in delta_lst: ### for each value of the period
    ### creating arrays and lists to store results
    results_std = {}
    results = {}
    noise_results = []
    E_all = []
    results_std["noise"] = noise_lst
    E_average, E_plus, E_minus = [], [], []
    for noise in noise_lst:
        E_results = []
        for i in range(1, 50): #### many iterations for convergence of the distribution
            noise_results.append(noise)
            E = destim.identifying_parameter(method = "VFM", delta = delta*l, load_type = "boundary_force", load_params = load_params_boundary, mesh_size = mesh_size, cube_params = cube_params, refine = False, noise_from_images = False, noise = noise,  u_params = {"u":u, "v":v}) #  this computation was only conducted for a boundary force test, without refining the mesh, and with synthetic data generated by adding noise to displacements
            E_error = (E-1)*100 ### computing the error compared to the ground truth, i.e., E = 1 kPa
            E_all.append(E_error)
            E_results.append(E_error)
        E_average.append(numpy.average(E_results))
        E_plus.append(numpy.average(E_results)+numpy.std(E_results))
        E_minus.append(numpy.average(E_results)-numpy.std(E_results))
    results_std["E_+"] = E_plus
    results_std["E_-"] = E_minus
    results_std["E_average"] = E_average
    results_all[str(delta)] = results_std
    results["noise"] = noise_results
    results["E"] = E_all
writing_results_to_pdf(results_all = results_all, delta_lst = delta_lst, SNR_lst = SNR_lst)
Contents