b2viscoelasticity.H

//------------------------------------------------------------------------
// b2visoelastity.H --
//
// written by Mathias Doreille
//            Neda Ebrahimi Pour <neda.ebrahimipour@dlr.de>
//
// (c) 2011, 2018 SMR Engineering & Development SA
//               2502 Bienne, Switzerland
//
// (c) 2023 Deutsches Zentrum für Luft- und Raumfahrt (DLR) e.V.
//          Linder Höhe, 51147 Köln
//
// All Rights Reserved.  Proprietary source code.  The contents of
// this file may not be disclosed to third parties, copied or
// duplicated in any form, in whole or in part, without the prior
// written permission of SMR.
//------------------------------------------------------------------------
 
#ifndef B2VISOELASTITY_H_
#define B2VISOELASTITY_H_
 
#include "b2solid_mechanics.H"
#include "utils/b2exception.H"
#include "utils/b2linear_algebra.H"
#include "utils/b2object.H"
#include "utils/b2tensor_calculus.H"
#include "utils/b2type.H"
 
namespace {
 
template <typename T>
inline void mult_mlp_v6(const T C[21], const T strain[6], T stress[6]) {
    const b2000::b2linalg::Matrix<T, b2000::b2linalg::Mlower_packed_constref> C_tmp(6, C);
    const b2000::b2linalg::Vector<T, b2000::b2linalg::Vdense_constref> strain_tmp(6, strain);
    b2000::b2linalg::Vector<T, b2000::b2linalg::Vdense_ref> stress_tmp(6, stress);
    stress_tmp = C_tmp * strain_tmp;
}
 
template <typename T>
inline void mult_mlp_v8(const T C[36], const T strain[8], T stress[8]) {
    const b2000::b2linalg::Matrix<T, b2000::b2linalg::Mlower_packed_constref> C_tmp(8, C);
    const b2000::b2linalg::Vector<T, b2000::b2linalg::Vdense_constref> strain_tmp(8, strain);
    b2000::b2linalg::Vector<T, b2000::b2linalg::Vdense_ref> stress_tmp(8, stress);
    stress_tmp = C_tmp * strain_tmp;
}
 
}  // namespace
 
namespace b2000 {
 
/* Reference:
 * Computational Inelasticity, Juan C. Simó, Thomas J. R. Hughes, 1998.
 * Section 10.3 (for the deviatoric viscosity only). The "full" elasticity is
 * easily deduced.
 *
 * This is for linear viscoleasticity with the small strain
 * assumption. The extension to finite strain (section 10.5 of the
 * same book) is not implemented here.
 */
template <typename T>
inline void visco_elastic_get_stress(
      const T E[3], const T P[3], const T G[3], const bool shell, const T sh,
      const std::vector<T>& relative_moduli, const std::vector<T>& relaxation_time,
      const bool only_deviatoric_viscosity, b2linalg::Matrix<T>& internal_viscous_stress,
      const T strain[6], const T delta_time, const EquilibriumSolution equilibrium_solution,
      T stress[6], T d_stress_d_strain[21], T& stored_energy_per_unit_volume,
      T& dissipated_energy_per_unit_volume) {
    if (internal_viscous_stress.size2() == 0) {
        internal_viscous_stress.resize(6, relative_moduli.size());
    }
 
    T d_stress_d_strain_tmp[21];
    T* C = d_stress_d_strain ? d_stress_d_strain : d_stress_d_strain_tmp;
 
    T ematrix[3][3];
    T cmatrix[3][3];
    std::fill_n(C, 21, T(0));
    if (shell) {
        cmatrix[0][0] = T(1) / E[0];
        cmatrix[1][1] = T(1) / E[1];
        cmatrix[2][2] = T(0);
        cmatrix[0][1] = cmatrix[1][0] = -P[0] * cmatrix[0][0];
        cmatrix[0][2] = cmatrix[2][0] = T(0);
        cmatrix[1][2] = cmatrix[2][1] = T(0);
 
        const T inv_det = T(1) / (cmatrix[0][0] * cmatrix[1][1] - cmatrix[0][1] * cmatrix[0][1]);
        C[0] = cmatrix[1][1] * inv_det;
        C[1] = -cmatrix[0][1] * inv_det;
        C[6] = cmatrix[0][0] * inv_det;
 
        // Set the diagonal values for the last 3 rows and columns.
        C[15] = G[0];
        C[18] = G[1] * sh;
        C[20] = G[2] * sh;
    } else {
        cmatrix[0][0] = T(1) / E[0];
        cmatrix[1][1] = T(1) / E[1];
        cmatrix[2][2] = T(1) / E[2];
        cmatrix[0][1] = cmatrix[1][0] = -P[0] * cmatrix[0][0];
        cmatrix[0][2] = cmatrix[2][0] = -P[1] * cmatrix[0][0];
        cmatrix[1][2] = cmatrix[2][1] = -P[2] * cmatrix[1][1];
 
        // Invert the compliance matrix to obtain the first three
        // columns and first three rows of the elasticity matrix.
        invert_3_3(cmatrix, ematrix);
 
        // Copy the 3x3 part of the elasticity matrix to the array.
        C[0] = ematrix[0][0];
        C[1] = ematrix[0][1];
        C[2] = ematrix[0][2];
        C[6] = ematrix[1][1];
        C[7] = ematrix[1][2];
        C[11] = ematrix[2][2];
 
        // Set the diagonal values for the last 3 rows and columns.
        C[15] = G[0];
        C[18] = G[1];
        C[20] = G[2];
    }
 
    if (stress || equilibrium_solution) {
        if (stress) { std::fill_n(stress, 6, T(0)); }
        if (equilibrium_solution) {
            stored_energy_per_unit_volume = T(0);
            dissipated_energy_per_unit_volume = T(0);
        }
        double stress_0[6] = {};
        if (only_deviatoric_viscosity) {
            const T strain_vol = (strain[0] + strain[1] + strain[2]) / 3;
            if (stress) {
                stress[0] = strain_vol * (C[0] + C[1] + C[2]);
                stress[1] = strain_vol * (C[1] + C[6] + C[7]);
                stress[2] = strain_vol * (C[2] + C[7] + C[11]);
            }
            if (equilibrium_solution) {
                stored_energy_per_unit_volume =
                      0.5 * strain_vol * strain_vol
                      * (C[0] + C[1] + C[2] + C[1] + C[6] + C[7] + C[2] + C[7] + C[11]);
            }
            T strain_dev[6];
            for (int i = 0; i != 3; ++i) { strain_dev[i] = strain[i] - strain_vol; }
            for (int i = 3; i != 6; ++i) { strain_dev[i] = strain[i]; }
            mult_mlp_v6(C, strain_dev, stress_0);
            const T stress_0_vol = (stress_0[0] + stress_0[1] + stress_0[2]) / 3;
            for (int i = 0; i != 3; ++i) { stress_0[i] -= stress_0_vol; }
        } else {
            mult_mlp_v6(C, strain, stress_0);
        }
 
        T dstress_0[6];
        const size_t n = relaxation_time.size();
        for (int i = 0; i != 6; ++i) { dstress_0[i] = stress_0[i] - internal_viscous_stress(i, n); }
        if (stress) {
            for (int i = 0; i != 6; ++i) { stress[i] += relative_moduli[n] * stress_0[i]; }
        }
        if (equilibrium_solution) {
            T etmp = 0;
            for (int i = 0; i != 6; ++i) {
                internal_viscous_stress(i, n) = stress_0[i];
                etmp += stress_0[i] * strain[i];
            }
            stored_energy_per_unit_volume += 0.5 * etmp * relative_moduli[n];
        }
 
        T* internal_viscous_stress_ptr = &internal_viscous_stress(0, 0);
        for (size_t j = 0; j != n; ++j) {
            const T dr = delta_time / relaxation_time[j];
            const T a = exp(-dr);
            const T b = (1 - a) / dr;
            for (int i = 0; i != 6; ++i, ++internal_viscous_stress_ptr) {
                const T h = a * *internal_viscous_stress_ptr + b * dstress_0[i];
                if (equilibrium_solution) { *internal_viscous_stress_ptr = h; }
                if (stress) { stress[i] += relative_moduli[j] * h; }
            }
            if (equilibrium_solution) {
                T strain_tmp[3];
                inner_product_3_1_NN(cmatrix, internal_viscous_stress_ptr - 6, strain_tmp);
                if (only_deviatoric_viscosity) {
                    const T strain_tmp_vol = (strain_tmp[0] + strain_tmp[1] + strain_tmp[2]) / 3;
                    for (int i = 0; i != 3; ++i) { strain_tmp[i] -= strain_tmp_vol; }
                }
                T etmp = dot_3(strain_tmp, internal_viscous_stress_ptr - 6);
                for (int i = 0; i != 3; ++i) {
                    etmp += internal_viscous_stress_ptr[i - 3] * internal_viscous_stress_ptr[i - 3]
                            / G[i];
                }
                stored_energy_per_unit_volume += relative_moduli[j] * 0.5 * etmp;
                assert(etmp >= 0);
                dissipated_energy_per_unit_volume += relative_moduli[j] * etmp / relaxation_time[j];
            }
        }
    }
 
    if (d_stress_d_strain) {
        T s = relative_moduli.back();
        for (size_t j = 0; j != relaxation_time.size(); ++j) {
            const T dr = delta_time / relaxation_time[j];
            s += relative_moduli[j] * (1 - exp(-dr)) / dr;
        }
        if (only_deviatoric_viscosity) {
            const T C_dd_I[6] = {
                  (C[0] + C[1] + C[2]) / 3,
                  (C[1] + C[6] + C[7]) / 3,
                  (C[2] + C[7] + C[11]) / 3,
                  C[15] / 3,
                  C[18] / 3,
                  C[20] / 3};
            T I_dd_C_dd_I = 0;
            for (int i = 0; i != 6; ++i) { I_dd_C_dd_I += C_dd_I[i]; }
            I_dd_C_dd_I /= 3;
            for (int j = 0; j != 6; ++j) {
                int i = j;
                for (; i < 3; ++i, ++C) {
                    *C = (C_dd_I[i] + C_dd_I[j]) / 2
                         + s * (*C - C_dd_I[i] - C_dd_I[j] + I_dd_C_dd_I);
                }
                for (; i < 6; ++i, ++C) { *C = s * (*C - C_dd_I[i] - C_dd_I[j] + I_dd_C_dd_I); }
            }
        } else {
            for (int i = 0; i != 21; ++i) { C[i] *= s; }
        }
    }
}
 
/* Reference:
 * Computational Inelasticity, Juan C. Simó, Thomas J. R. Hughes, 1998.
 * Section 10.3 (for the deviatoric viscosity only). The "full" elasticity is
 * easily deduced.
 *
 * This is for linear orthotropic viscoleasticity with the small strain
 * assumption. The extension to finite strain (section 10.5 of the
 * same book) is not implemented here.
 */
template <typename T>
inline void visco_elastic_get_stress(
      const T E[3], const T P[3], const T G[3], const bool shell, const T sh,
      const b2linalg::Matrix<T>& relative_moduli, const b2linalg::Matrix<T>& relaxation_time,
      const bool only_deviatoric_viscosity, b2linalg::Matrix<T>& internal_viscous_stress,
      const T strain[6], const T delta_time, const EquilibriumSolution equilibrium_solution,
      T stress[6], T d_stress_d_strain[21], T& stored_energy_per_unit_volume,
      T& dissipated_energy_per_unit_volume) {
    if (internal_viscous_stress.size2() == 0) {
        internal_viscous_stress.resize(6, relative_moduli.size2());
    }
 
    T d_stress_d_strain_tmp[21];
    T* C = d_stress_d_strain ? d_stress_d_strain : d_stress_d_strain_tmp;
 
    T ematrix[3][3];
    T cmatrix[3][3];
    std::fill_n(C, 21, T(0));
    if (shell) {
        cmatrix[0][0] = T(1) / E[0];
        cmatrix[1][1] = T(1) / E[1];
        cmatrix[2][2] = T(0);
        cmatrix[0][1] = cmatrix[1][0] = -P[0] * cmatrix[0][0];
        cmatrix[0][2] = cmatrix[2][0] = T(0);
        cmatrix[1][2] = cmatrix[2][1] = T(0);
 
        const T inv_det = 1 / (cmatrix[0][0] * cmatrix[1][1] - cmatrix[0][1] * cmatrix[0][1]);
        C[0] = cmatrix[1][1] * inv_det;
        C[1] = -cmatrix[0][1] * inv_det;
        C[6] = cmatrix[0][0] * inv_det;
 
        // Set the diagonal values for the last 3 rows and columns.
        C[15] = G[0];
        C[18] = G[1] * sh;
        C[20] = G[2] * sh;
    } else {
        cmatrix[0][0] = T(1) / E[0];
        cmatrix[1][1] = T(1) / E[1];
        cmatrix[2][2] = T(1) / E[2];
        cmatrix[0][1] = cmatrix[1][0] = -P[0] * cmatrix[0][0];
        cmatrix[0][2] = cmatrix[2][0] = -P[1] * cmatrix[0][0];
        cmatrix[1][2] = cmatrix[2][1] = -P[2] * cmatrix[1][1];
 
        // Invert the compliance matrix to obtain the first three
        // columns and first three rows of the elasticity matrix.
        invert_3_3(cmatrix, ematrix);
 
        // Copy the 3x3 part of the elasticity matrix to the array.
        C[0] = ematrix[0][0];
        C[1] = ematrix[0][1];
        C[2] = ematrix[0][2];
        C[6] = ematrix[1][1];
        C[7] = ematrix[1][2];
        C[11] = ematrix[2][2];
 
        // Set the diagonal values for the last 3 rows and columns.
        C[15] = G[0];
        C[18] = G[1];
        C[20] = G[2];
    }
 
    if (stress || equilibrium_solution) {
        if (stress) { std::fill_n(stress, 6, T(0)); }
        if (equilibrium_solution) {
            stored_energy_per_unit_volume = T(0);
            dissipated_energy_per_unit_volume = T(0);
        }
        T stress_0[6];
        if (only_deviatoric_viscosity) {
            const T strain_vol = (strain[0] + strain[1] + strain[2]) / 3;
            if (stress) {
                stress[0] = strain_vol * (C[0] + C[1] + C[2]);
                stress[1] = strain_vol * (C[1] + C[6] + C[7]);
                stress[2] = strain_vol * (C[2] + C[7] + C[11]);
            }
            if (equilibrium_solution) {
                stored_energy_per_unit_volume =
                      0.5 * strain_vol * strain_vol
                      * (C[0] + C[1] + C[2] + C[1] + C[6] + C[7] + C[2] + C[7] + C[11]);
            }
            T strain_dev[6];
            for (int i = 0; i != 3; ++i) { strain_dev[i] = strain[i] - strain_vol; }
            for (int i = 3; i != 6; ++i) { strain_dev[i] = strain[i]; }
            mult_mlp_v6(C, strain_dev, stress_0);
            const T stress_0_vol = (stress_0[0] + stress_0[1] + stress_0[2]) / 3;
            for (int i = 0; i != 3; ++i) { stress_0[i] -= stress_0_vol; }
        } else {
            mult_mlp_v6(C, strain, stress_0);
        }
 
        T dstress_0[6];
        const size_t n = relaxation_time.size2();
        for (int i = 0; i != 6; ++i) { dstress_0[i] = stress_0[i] - internal_viscous_stress(i, n); }
        if (stress) {
            for (int i = 0; i != 6; ++i) { stress[i] += relative_moduli(i, n) * stress_0[i]; }
        }
        if (equilibrium_solution) {
            for (int i = 0; i != 6; ++i) {
                internal_viscous_stress(i, n) = stress_0[i];
                const T etmp = stress_0[i] * strain[i];
                stored_energy_per_unit_volume += 0.5 * etmp * relative_moduli(i, n);
            }
        }
 
        T* internal_viscous_stress_ptr = &internal_viscous_stress(0, 0);
        for (size_t j = 0; j != n; ++j) {
            for (int i = 0; i != 6; ++i, ++internal_viscous_stress_ptr) {
                const T dr = delta_time / relaxation_time(i, j);
                const T a = exp(-dr);
                const T b = (1 - a) / dr;
                const T h = a * *internal_viscous_stress_ptr + b * dstress_0[i];
                if (equilibrium_solution) { *internal_viscous_stress_ptr = h; }
                if (stress) { stress[i] += relative_moduli(i, j) * h; }
            }
            if (equilibrium_solution) {
                T strain_tmp[6];
                inner_product_3_1_NN(cmatrix, internal_viscous_stress_ptr - 6, strain_tmp);
                if (only_deviatoric_viscosity) {
                    const T strain_tmp_vol = (strain_tmp[0] + strain_tmp[1] + strain_tmp[2]) / 3;
                    for (int i = 0; i != 3; ++i) { strain_tmp[i] -= strain_tmp_vol; }
                }
                for (int i = 0; i != 3; ++i) {
                    strain_tmp[i + 3] = internal_viscous_stress_ptr[i - 3] / G[i];
                }
                for (int i = 0; i != 6; ++i) {
                    stored_energy_per_unit_volume += relative_moduli(i, j) * 0.5 * strain_tmp[i]
                                                     * internal_viscous_stress_ptr[i - 6];
                    dissipated_energy_per_unit_volume += relative_moduli(i, j)
                                                         / relaxation_time(i, j) * strain_tmp[i]
                                                         * internal_viscous_stress_ptr[i - 6];
                }
            }
        }
    }
 
    if (d_stress_d_strain) {
        T s[6];
        for (int i = 0; i != 6; ++i) {
            s[i] = relative_moduli(i, relative_moduli.size2() - 1);
            for (size_t j = 0; j != relaxation_time.size2(); ++j) {
                const T dr = delta_time / relaxation_time(i, j);
                s[i] += relative_moduli(i, j) * (1 - exp(-dr)) / dr;
            }
        }
        if (only_deviatoric_viscosity) {
            const T C_dd_I[6] = {
                  (C[0] + C[1] + C[2]) / 3,
                  (C[1] + C[6] + C[7]) / 3,
                  (C[2] + C[7] + C[11]) / 3,
                  C[15] / 3,
                  C[18] / 3,
                  C[20] / 3};
            T I_dd_C_dd_I = 0;
            for (int i = 0; i != 6; ++i) { I_dd_C_dd_I += C_dd_I[i]; }
            I_dd_C_dd_I /= 3;
            for (int j = 0; j != 6; ++j) {
                int i = j;
                for (; i < 3; ++i, ++C) {
                    *C = (C_dd_I[i] + C_dd_I[j]) / 2
                         + 0.5 * (s[i] + s[j]) * (*C - C_dd_I[i] - C_dd_I[j] + I_dd_C_dd_I);
                }
                for (; i < 6; ++i, ++C) {
                    *C = 0.5 * (s[i] + s[j]) * (*C - C_dd_I[i] - C_dd_I[j] + I_dd_C_dd_I);
                }
            }
        } else {
            for (int j = 0; j != 6; ++j) {
                for (int i = j; i != 6; ++i, ++C) { *C *= 0.5 * (s[i] + s[j]); }
            }
        }
    }
}
 
/*
 * Same implementation for shell strain and shell stress
 */
template <typename T>
inline void visco_elastic_get_stress(
      const T d_stress_d_strain_0[36], const T inv_d_stress_d_strain_0[36],
      const std::vector<T>& relative_moduli, const std::vector<T>& relaxation_time,
      b2linalg::Matrix<T>& internal_viscous_stress, const T strain[8], const T delta_time,
      const EquilibriumSolution equilibrium_solution, T stress[8], T d_stress_d_strain[36],
      T& stored_energy_per_unit_volume, T& dissipated_energy_per_unit_volume) {
    if (internal_viscous_stress.size2() == 0) {
        internal_viscous_stress.resize(8, relative_moduli.size());
    }
 
    if (stress || equilibrium_solution) {
        if (stress) { std::fill_n(stress, 8, T(0)); }
        if (equilibrium_solution) {
            stored_energy_per_unit_volume = T(0);
            dissipated_energy_per_unit_volume = T(0);
        }
        T stress_0[8];
        mult_mlp_v8(d_stress_d_strain_0, strain, stress_0);
 
        T dstress_0[8];
        const size_t n = relaxation_time.size();
        for (int i = 0; i != 8; ++i) { dstress_0[i] = stress_0[i] - internal_viscous_stress(i, n); }
        if (stress) {
            for (int i = 0; i != 8; ++i) { stress[i] += relative_moduli[n] * stress_0[i]; }
        }
        if (equilibrium_solution) {
            T etmp = 0;
            for (int i = 0; i != 8; ++i) {
                internal_viscous_stress(i, n) = stress_0[i];
                etmp += stress_0[i] * strain[i];
            }
            stored_energy_per_unit_volume += 0.5 * etmp * relative_moduli[n];
        }
 
        T* internal_viscous_stress_ptr = &internal_viscous_stress(0, 0);
        for (size_t j = 0; j != n; ++j) {
            const T dr = delta_time / relaxation_time[j];
            const T a = exp(-dr);
            const T b = (1 - a) / dr;
            for (int i = 0; i != 8; ++i, ++internal_viscous_stress_ptr) {
                const T h = a * *internal_viscous_stress_ptr + b * dstress_0[i];
                if (equilibrium_solution) { *internal_viscous_stress_ptr = h; }
                if (stress) { stress[i] += relative_moduli[j] * h; }
            }
            if (equilibrium_solution) {
                T strain_tmp[8];
                mult_mlp_v8(inv_d_stress_d_strain_0, internal_viscous_stress_ptr - 8, strain_tmp);
                T etmp = T(0);
                for (int i = 0; i != 8; ++i) {
                    etmp += strain_tmp[i] * internal_viscous_stress_ptr[i - 8];
                }
                stored_energy_per_unit_volume += relative_moduli[j] * 0.5 * etmp;
                assert(etmp >= 0);
                dissipated_energy_per_unit_volume += relative_moduli[j] * etmp / relaxation_time[j];
            }
        }
    }
 
    if (d_stress_d_strain) {
        T s = relative_moduli.back();
        for (size_t j = 0; j != relaxation_time.size(); ++j) {
            const T dr = delta_time / relaxation_time[j];
            s += relative_moduli[j] * (1 - exp(-dr)) / dr;
        }
        for (int i = 0; i != 36; ++i) { d_stress_d_strain[i] = d_stress_d_strain_0[i] * s; }
    }
}
 
// same implementation for rod
template <typename T>
inline void visco_elastic_get_stress(
      const T E, const std::vector<T>& relative_moduli, const std::vector<T>& relaxation_time,
      b2linalg::Vector<T>& internal_viscous_stress, const T strain, const T delta_time,
      const EquilibriumSolution equilibrium_solution, T& stress, T& d_stress_d_strain,
      T& stored_energy_per_unit_volume, T& dissipated_energy_per_unit_volume) {
    if (internal_viscous_stress.size() == 0) {
        internal_viscous_stress.resize(relative_moduli.size());
    }
 
    {
        if (equilibrium_solution) {
            stored_energy_per_unit_volume = T(0);
            dissipated_energy_per_unit_volume = T(0);
        }
        T stress_0 = strain * E;
        T dstress_0;
        const size_t n = relaxation_time.size();
        dstress_0 = stress_0 - internal_viscous_stress[n];
        stress = relative_moduli[n] * stress_0;
        if (equilibrium_solution) {
            internal_viscous_stress[n] = stress_0;
            stored_energy_per_unit_volume += 0.5 * stress_0 * strain * relative_moduli[n];
        }
 
        for (size_t j = 0; j != n; ++j) {
            const T dr = delta_time / relaxation_time[j];
            const T a = exp(-dr);
            const T b = (1 - a) / dr;
            const T h = a * internal_viscous_stress[j] + b * dstress_0;
            if (equilibrium_solution) { internal_viscous_stress[j] = h; }
            stress += relative_moduli[j] * h;
            if (equilibrium_solution) {
                const T etmp = internal_viscous_stress[j] * internal_viscous_stress[j] / E;
                stored_energy_per_unit_volume += relative_moduli[j] * 0.5 * etmp;
                assert(etmp >= 0);
                dissipated_energy_per_unit_volume += relative_moduli[j] * etmp / relaxation_time[j];
            }
        }
    }
 
    {
        T s = relative_moduli.back();
        for (size_t j = 0; j != relaxation_time.size(); ++j) {
            const T dr = delta_time / relaxation_time[j];
            s += relative_moduli[j] * (1 - exp(-dr)) / dr;
        }
        d_stress_d_strain = E * s;
    }
}
 
class ViscoElasticParameterIdentification {
public:
    ViscoElasticParameterIdentification(
          int nb_maxwell_element_, double freq_start_, double freq_end_, bool logarithm_precision_,
          double epsilon_integration_ = 1e-5, double tol_opt_ = 1e-5)
        : scale_real(0), scale_imag(0), x_nlopt(nullptr), modulus_at_relaxation(0) {
        nb_maxwell_element = nb_maxwell_element_;
        freq_start = freq_start_;
        logarithm_precision = logarithm_precision_;
        epsilon_integration = epsilon_integration_;
        tol_opt = tol_opt_;
        freq_start = freq_start_;
        freq_end = freq_end_;
        x.resize(nb_maxwell_element_ * 2);
    }
 
    void set_expected_modulus_as_expresion(
          const double modulus, const std::string& expresion_storage_modulus,
          const std::string& expresion_loss_moduluis, bool loss_modulus_as_tan_modulus = false);
 
    void set_expected_modulus_as_tabulated_values(
          const double modulus, const std::vector<double>& storage_modulus,
          const std::vector<double>& lose_modulus, bool loss_modulus_as_tan_modulus = false,
          int interpolation_order = 1);
 
    void get_relativ_modulus(
          std::vector<double>& relative_moduli, std::vector<double>& relaxation_time,
          const bool modulus_at_relaxation_ = true);
 
private:
    double f(const double w, const unsigned int w_level);
 
    double adaptiveSimpsonsAux(
          double a, double b, unsigned int a_level, unsigned int b_level, double epsilon, double S,
          double fa, double fb, double fc);
 
    static double myfunc(unsigned n, const double* x, double* grad, void* my_func_data);
 
    double adaptiveSimpsonsAux_log(
          double a, double b, unsigned int a_level, unsigned int b_level, double epsilon, double S,
          double fa, double fb, double fc);
 
    static double myfunc_log(unsigned n, const double* x, double* grad, void* my_func_data);
 
    static const int max_level = 12;
    static const unsigned int pow_max_level = 4096;  // = 2 * max_level
    int nb_maxwell_element;
    double freq_start;
    double freq_end;
    bool logarithm_precision;
    double epsilon_integration;
    double tol_opt;
    double scale_real;
    double scale_imag;
    std::pair<double, double> expected_value[pow_max_level + 1];
    std::vector<double> x;
    const double* x_nlopt;
    bool modulus_at_relaxation;
};
 
}  // namespace b2000
 
#endif  // B2VISOELASTITY_H_