b2material_util.H

Go to the documentation of this file.

//------------------------------------------------------------------------
// b2material_util.H --
//
//     Helper classes for elastic material implementations.
//
// written by Thomas Ludwig
//            Neda Ebrahimi Pour <neda.ebrahimipour@dlr.de>
//
// (c) 2015,2016,2017 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 B2_MATERIAL_UTIL_H_
#define B2_MATERIAL_UTIL_H_
 
#include <algorithm>
#include <cmath>
 
#include "b2domain_database.H"
#include "model/b2solver.H"
#include "utils/b2linear_algebra.H"
#include "utils/b2rtable.H"
#include "utils/b2tensor_transformation.H"
 
 
namespace {
 
// Multiply by 1/2 the off-diagonal entries of a strain tensor and
// store it in the gradient_container.
//
template <typename T>
inline void set_field_value_strain(
      b2000::GradientContainer* gradient_container,
      const b2000::GradientContainer::FieldDescription& descr, const T strain[6]) {
    const T strain_storage[6] = {strain[0],         strain[1],         strain[2],
                                 T(.5) * strain[3], T(.5) * strain[4], T(.5) * strain[5]};
    gradient_container->set_field_value(descr, strain_storage);
}
 
template <typename T>
inline void set_field_value_shell_strain(
      b2000::GradientContainer* gradient_container,
      const b2000::GradientContainer::FieldDescription& descr, const T strain[8]) {
    const T strain_storage[8] = {
          strain[0], strain[1], strain[2],         strain[3],
          strain[4], strain[5], T(.5) * strain[6], T(.5) * strain[7],
    };
    gradient_container->set_field_value(descr, strain_storage);
}
 
template <typename T>
inline void set_field_value_strain_2(
      b2000::GradientContainer* gradient_container,
      const b2000::GradientContainer::FieldDescription& descr, const T strain_2d[3],
      const T strain_normal) {
    const T strain_storage[6] = {strain_2d[0],         strain_2d[1], strain_normal,
                                 T(.5) * strain_2d[2], T(0),         T(0)};
    gradient_container->set_field_value(descr, strain_storage);
}
 
template <typename T>
inline void set_field_value_stress_2(
      b2000::GradientContainer* gradient_container,
      const b2000::GradientContainer::FieldDescription& descr, const T stress_2d[3],
      const T stress_normal) {
    const T stress_storage[6] = {stress_2d[0], stress_2d[1], stress_normal,
                                 stress_2d[2], T(0),         T(0)};
    gradient_container->set_field_value(descr, stress_storage);
}
 
}  // namespace
 
namespace b2000 { namespace b2dbv3 {
 
template <typename T>
struct ElasticityHelper {
    bool has_thermal_expansion;
 
    bool has_material_ref;
 
    T M[3][3];
 
    bool has_covariant_base;
 
    T G[3][3];
 
    T G_d[3][3];
 
    T strain_gl_bg[6];
 
    T strain_gl_bg_te[6];
 
    T strain_gl_mat[6];
 
    T strain_gl_mat_te[6];
 
    T stress_pk2_bg[6];
 
    T stress_pk2_mat[6];
 
    T stress_cauchy_bg[6];
 
    T stress_cauchy_mat[6];
 
    T failure_index;
 
    ElasticityHelper(const T material_ref[3][3], const T covariant_base[3][3])
        : has_thermal_expansion(false),
          has_material_ref(material_ref != 0),
          has_covariant_base(covariant_base != 0),
          failure_index(0.) {
        // Material basis (orthogonal reference frame) and its transpose.
        if (has_material_ref) {
            assert(material_ref != 0);
            std::copy(&material_ref[0][0], &material_ref[0][0] + 9, &M[0][0]);
        } else {
            b2000::set_identity_3_3(M);
        }
 
        // Element covariant and contravariant basis.
        if (has_covariant_base) {
            assert(covariant_base != 0);
            b2000::copy_3_3(covariant_base, G);
            b2000::get_base_opposite_variance(G, G_d);
        } else {
            // Everything is aligned with the global coordinate system.
            b2000::set_identity_3_3(G);
            b2000::set_identity_3_3(G_d);
        }
 
        std::fill_n(strain_gl_bg, 6, 0);
        std::fill_n(strain_gl_bg_te, 6, 0);
        std::fill_n(strain_gl_mat, 6, 0);
        std::fill_n(strain_gl_mat_te, 6, 0);
        std::fill_n(stress_pk2_bg, 6, 0);
        std::fill_n(stress_pk2_mat, 6, 0);
        std::fill_n(stress_cauchy_bg, 6, 0);
        std::fill_n(stress_cauchy_mat, 6, 0);
    }
 
    void set_strain(
          const bool planestress, const T strain[6], const T temperature, const T alpha[6]) {
        has_thermal_expansion = false;
 
        if (strain != 0) {
            transform_strain_from_base_A_to_I(G, strain, strain_gl_bg);
            transform_strain_from_I_to_base_B(M, strain_gl_bg, strain_gl_mat);
        }
        if (alpha != 0 && temperature != 0) {
            has_thermal_expansion = true;
            strain_gl_mat_te[0] = strain_gl_mat[0] - alpha[0] * temperature;
            strain_gl_mat_te[1] = strain_gl_mat[1] - alpha[1] * temperature;
            if (planestress) {
                strain_gl_mat_te[2] = 0;  // no valid zz strain
            } else {
                strain_gl_mat_te[2] = strain_gl_mat[2] - alpha[2] * temperature;
            }
            strain_gl_mat_te[3] = strain_gl_mat[3] - alpha[3] * temperature;
            strain_gl_mat_te[4] = strain_gl_mat[4] - alpha[4] * temperature;
            strain_gl_mat_te[5] = strain_gl_mat[5] - alpha[5] * temperature;
        } else {
            std::copy(strain_gl_mat, strain_gl_mat + 6, strain_gl_mat_te);
        }
    }
 
    void set_stress_and_or_d_stress_d_strain_mat(
          const T stress_mat[6],  // may be 0
          const T d_stress_d_strain_mat[21], T stress[6], T d_stress_d_strain[21]) {
        if (d_stress_d_strain != 0) {
            assert(d_stress_d_strain_mat != 0);
            transform_constitutive_tensor_from_base_A_ov_to_base_B(
                  M, G_d, d_stress_d_strain_mat, d_stress_d_strain);
        }
 
        //
        // Calculate the stresses.
        //
        if (stress_mat != 0) {
            std::copy(stress_mat, stress_mat + 6, stress_pk2_mat);
        } else {
            // Linear material: Calculate stress from constitutive matrix.
            assert(d_stress_d_strain_mat != 0);
            b2linalg::Matrix<T, b2linalg::Mlower_packed_constref> C_mat(6, d_stress_d_strain_mat);
            b2linalg::Vector<T, b2linalg::Vdense_constref> e_mat(6, strain_gl_mat_te);
            b2linalg::Vector<T, b2linalg::Vdense_ref> s_mat(6, stress_pk2_mat);
            s_mat = C_mat * e_mat;
        }
        transform_stress_from_base_A_ov_to_I(M, stress_pk2_mat, stress_pk2_bg);
        if (stress != 0) { transform_stress_from_I_to_base_B(G_d, stress_pk2_bg, stress); }
    }
 
    void store(
          const bool linear, const T displacement_gradient[3][3],
          const FailureCriterion* failure_criterion, b2000::GradientContainer* gradient_container) {
        // Failure criterion.
        if (failure_criterion != nullptr) {
            if (failure_criterion->get_type() == FailureCriterion::STRESS) {
                failure_index =
                      failure_criterion->process_stress(stress_pk2_mat, gradient_container);
            } else {
                failure_index =
                      failure_criterion->process_strain(strain_gl_mat, false, gradient_container);
            }
        }
 
        if (gradient_container == nullptr) { return; }
 
        // Transform the PK2 stress to the Cauchy stress.
        if (!linear && displacement_gradient != 0) {
            T deformation_gradient[3][3] = {};
            if (has_covariant_base) {
                // T G_inv[3][3];
                // invert_3_3(G, G_inv);
                b2000::inner_product_3_3_NT(displacement_gradient, G_d, deformation_gradient);
            } else {
                std::copy(
                      &displacement_gradient[0][0], &displacement_gradient[0][0] + 9,
                      &deformation_gradient[0][0]);
            }
            deformation_gradient[0][0] += 1;
            deformation_gradient[1][1] += 1;
            deformation_gradient[2][2] += 1;
            const T det = 1 / b2000::determinant_3_3(deformation_gradient);
            if (has_material_ref) {
                T deformation_gradient_aligned_ortho[3][3];
                b2000::inner_product_3_3_NN(
                      deformation_gradient, M, deformation_gradient_aligned_ortho);
                b2000::transform_stress_from_base_A_ov_to_I(
                      deformation_gradient_aligned_ortho, stress_pk2_mat, stress_cauchy_bg);
            } else {
                transform_stress_from_base_A_ov_to_I(
                      deformation_gradient, stress_pk2_mat, stress_cauchy_bg);
            }
            for (int i = 0; i != 6; ++i) { stress_cauchy_bg[i] *= det; }
        } else {
            std::copy(stress_pk2_bg, stress_pk2_bg + 6, stress_cauchy_bg);
        }
 
        //
        // Store MBASE.
        //
        if (has_material_ref) {
            static const b2000::GradientContainer::FieldDescription descr = {
                  "MBASE_IP",
                  "e1x.e1y.e1z.e2x.e2y.e2z.e3x.e3y.e3z",
                  "Material base vectors, undeformed configuration",
                  3,
                  9,
                  1,
                  3,
                  false,
                  typeid(T)};
            if (gradient_container->compute_field_value(descr.name)) {
                gradient_container->set_field_value(descr, &M[0][0]);
            }
        }
 
        //
        // Store the different strain and stress tensors.
        //
        if (1) {
            static const b2000::GradientContainer::FieldDescription descr = {
                  "STRAIN",
                  "Exx.Eyy.Ezz.Exy.Eyz.Exz",
                  "Green-Lagrange/linear strain, global reference frame, no thermal expansion",
                  3,
                  6,
                  2,
                  3,
                  true,
                  typeid(T)};
            if (gradient_container->compute_field_value(descr.name)) {
                set_field_value_strain(gradient_container, descr, strain_gl_bg);
            }
        }
 
        if (0) {
            static const b2000::GradientContainer::FieldDescription descr = {
                  "STRAIN_MATERIAL",
                  "Exx.Eyy.Ezz.Exy.Eyz.Exz",
                  "Green-Lagrange/linear strain, material reference frame, no thermal expansion",
                  3,
                  6,
                  2,
                  3,
                  true,
                  typeid(T)};
            if (gradient_container->compute_field_value(descr.name)) {
                set_field_value_strain(gradient_container, descr, strain_gl_mat);
            }
        }
 
        if (0 && has_thermal_expansion) {
            static const b2000::GradientContainer::FieldDescription descr = {
                  "STRAIN_MECHANICAL",
                  "Exx.Eyy.Ezz.Exy.Eyz.Exz",
                  "Green-Lagrange/linear strain, global reference frame, thermal expansion",
                  3,
                  6,
                  2,
                  3,
                  true,
                  typeid(T)};
            if (gradient_container->compute_field_value(descr.name)) {
                b2000::transform_strain_from_base_A_ov_to_I(M, strain_gl_mat_te, strain_gl_bg_te);
                set_field_value_strain(gradient_container, descr, strain_gl_bg_te);
            }
        }
 
        if (0 && has_thermal_expansion) {
            static const b2000::GradientContainer::FieldDescription descr = {
                  "STRAIN_MECHANICAL_MATERIAL",
                  "Exx.Eyy.Ezz.Exy.Eyz.Exz",
                  "Green-Lagrange/linear strain, material reference frame, thermal expansion",
                  3,
                  6,
                  2,
                  3,
                  true,
                  typeid(T)};
            if (gradient_container->compute_field_value(descr.name)) {
                set_field_value_strain(gradient_container, descr, strain_gl_mat_te);
            }
        }
 
        if (1) {
            static const b2000::GradientContainer::FieldDescription descr = {
                  "STRESS",
                  "Sxx.Syy.Szz.Sxy.Syz.Sxz",
                  "Cauchy/linear stress, global reference frame",
                  3,
                  6,
                  2,
                  3,
                  true,
                  typeid(T)};
            if (gradient_container->compute_field_value(descr.name)) {
                gradient_container->set_field_value(descr, stress_cauchy_bg);
            }
        }
 
        if (0) {
            static const b2000::GradientContainer::FieldDescription descr = {
                  "STRESS_MATERIAL",
                  "Sxx.Syy.Szz.Sxy.Syz.Sxz",
                  "Cauchy/linear stress, material reference frame",
                  3,
                  6,
                  2,
                  3,
                  true,
                  typeid(T)};
            if (gradient_container->compute_field_value(descr.name)) {
                b2000::transform_stress_from_I_to_base_B(M, stress_cauchy_bg, stress_cauchy_mat);
                gradient_container->set_field_value(descr, stress_cauchy_mat);
            }
        }
 
        if (0) {
            static const b2000::GradientContainer::FieldDescription descr = {
                  "STRESS_PK2",
                  "Sxx.Syy.Szz.Sxy.Syz.Sxz",
                  "Second-Piola Kirchhoff stress, global reference frame",
                  3,
                  6,
                  2,
                  3,
                  true,
                  typeid(T)};
            if (gradient_container->compute_field_value(descr.name)) {
                transform_stress_from_base_A_ov_to_I(M, stress_pk2_mat, stress_pk2_bg);
                gradient_container->set_field_value(descr, stress_pk2_bg);
            }
        }
 
        if (0) {
            static const b2000::GradientContainer::FieldDescription descr = {
                  "STRESS_PK2_MATERIAL",
                  "Sxx.Syy.Szz.Sxy.Syz.Sxz",
                  "Second-Piola Kirchhoff stress in the material reference frame",
                  3,
                  6,
                  2,
                  3,
                  true,
                  typeid(T)};
            if (gradient_container->compute_field_value(descr.name)) {
                gradient_container->set_field_value(descr, stress_pk2_mat);
            }
        }
 
        if (0) {
            static const b2000::GradientContainer::FieldDescription descr = {
                  "ENERGY_DENSITY",
                  "Elastic.Kinetic.RateDissipated.Dissipated",
                  "Free elastic, kinetic, rate of dissipated energy and "
                  "dissipated energy, per unit volume, in the undeformed "
                  "configuration",
                  3,
                  4,
                  -1,
                  0,
                  false,
                  typeid(T)};
            if (gradient_container->compute_field_value(descr.name)) {
                T es = 0;
                for (int i = 0; i != 6; ++i) { es += strain_gl_mat[i] * stress_pk2_mat[i]; }
                const T e[4] = {0.5 * es, 0, 0, 0};
                gradient_container->set_field_value(descr, e);
            }
        }
    }
};
 
template <typename T>
struct ElasticityHelperShellResultants {
    bool has_thermal_expansion;
 
    bool has_material_ref;
 
    T M[3][3];
 
    bool has_covariant_base;
 
    T G[3][3];
 
    T G_d[3][3];
 
    T strain_gl_mat[8];
 
    T strain_gl_mat_te[8];
 
    T stress_pk2_mat[8];
 
    T failure_index;
 
    ElasticityHelperShellResultants(const T material_ref[3][3], const T covariant_base[3][3])
        : has_thermal_expansion(false),
          has_material_ref(material_ref != 0),
          has_covariant_base(covariant_base != 0),
          failure_index(0.) {
        // Material basis (orthogonal reference frame) and its transpose.
        if (has_material_ref) {
            assert(material_ref != 0);
            std::copy(&material_ref[0][0], &material_ref[0][0] + 9, &M[0][0]);
        } else {
            set_identity_3_3(M);
        }
 
        // Element covariant and contravariant basis.
        if (has_covariant_base) {
            assert(covariant_base != 0);
            copy_3_3(covariant_base, G);
            get_base_opposite_variance(G, G_d);
        } else {
            // Everything is aligned with the global coordinate system.
            set_identity_3_3(G);
            set_identity_3_3(G_d);
        }
 
        std::fill_n(strain_gl_mat, 8, 0);
        std::fill_n(strain_gl_mat_te, 8, 0);
        std::fill_n(stress_pk2_mat, 8, 0);
    }
 
    void set_strain(const T strain[8], const T temperature, const T alpha[8]) {
        has_thermal_expansion = false;
 
        if (strain != 0) {
            b2000::transform_shell_strain_from_base_A_to_base_B(G, M, strain, strain_gl_mat);
        }
        if (alpha != 0 && temperature != 0) {
            has_thermal_expansion = true;
            strain_gl_mat_te[0] = strain_gl_mat[0] - alpha[0] * temperature;
            strain_gl_mat_te[1] = strain_gl_mat[1] - alpha[1] * temperature;
            strain_gl_mat_te[2] = strain_gl_mat[2] - alpha[2] * temperature;
            strain_gl_mat_te[3] = strain_gl_mat[3] - alpha[3] * temperature;
            strain_gl_mat_te[4] = strain_gl_mat[4] - alpha[4] * temperature;
            strain_gl_mat_te[5] = strain_gl_mat[5] - alpha[5] * temperature;
            strain_gl_mat_te[6] = strain_gl_mat[6] - alpha[6] * temperature;
            strain_gl_mat_te[7] = strain_gl_mat[7] - alpha[7] * temperature;
        } else {
            std::copy(strain_gl_mat, strain_gl_mat + 8, strain_gl_mat_te);
        }
    }
 
 
    void set_stress_and_or_d_stress_d_strain_mat(
          const T stress_mat[8],  // may be 0
          const T d_stress_d_strain_mat[36], T stress[8], T d_stress_d_strain[36]) {
        if (d_stress_d_strain != 0) {
            assert(d_stress_d_strain_mat != 0);
            transform_shell_constitutive_tensor_from_base_A_ov_to_base_B(
                  M, G_d, d_stress_d_strain_mat, d_stress_d_strain);
        }
 
        //
        // Calculate the stresses.
        //
        if (stress_mat != 0) {
            std::copy(stress_mat, stress_mat + 8, stress_pk2_mat);
        } else {
            // Linear material: Calculate stress from constitutive matrix.
            assert(d_stress_d_strain_mat != 0);
            b2linalg::Matrix<T, b2linalg::Mlower_packed_constref> C_mat(8, d_stress_d_strain_mat);
            b2linalg::Vector<T, b2linalg::Vdense_constref> e_mat(8, strain_gl_mat_te);
            b2linalg::Vector<T, b2linalg::Vdense_ref> s_mat(8, stress_pk2_mat);
            s_mat = C_mat * e_mat;
        }
        if (stress != 0) {
            b2000::transform_shell_stress_from_base_A_to_base_B(M, G_d, stress_pk2_mat, stress);
        }
    }
 
    void store(
          const bool linear, const b2dbv3::FailureCriterion* failure_criterion,
          b2000::GradientContainer* gradient_container) {
        // // Failure criterion.
        // if (failure_criterion != 0) {
        //     if (failure_criterion->get_type() == FailureCriterion::STRESS)
        //         failure_index = failure_criterion->process_shell_stress(
        //             stress_pk2_mat, gradient_container);
        //     else
        //         failure_index = failure_criterion->process_shell_strain(
        //             strain_gl_mat, false, gradient_container);
        // }
 
        if (gradient_container == nullptr) { return; }
 
        //
        // Store MBASE.
        //
        if (has_material_ref) {
            static const b2000::GradientContainer::FieldDescription descr = {
                  "MBASE_IP",
                  "e1x.e1y.e1z.e2x.e2y.e2z.e3x.e3y.e3z",
                  "Material base vectors, undeformed configuration",
                  3,
                  9,
                  1,
                  3,
                  false,
                  typeid(T)};
            if (gradient_container->compute_field_value(descr.name)) {
                gradient_container->set_field_value(descr, &M[0][0]);
            }
        }
 
        //
        // Store the different strain and stress tensors.
        //
        {
            static const b2000::GradientContainer::FieldDescription descr = {
                  "STRAIN_SHELL_MATERIAL",
                  "Mxx.Myy.Mxy.Bxx.Byy.Bxy.Tyz.Txz",
                  "Green-Lagrange/linear strain, material reference "
                  "frame, no thermal expansion",
                  3,
                  8,
                  -1,
                  -1,
                  false,
                  typeid(T)};
            if (gradient_container->compute_field_value(descr.name)) {
                set_field_value_shell_strain(gradient_container, descr, strain_gl_mat);
            }
        }
        {
            static const b2000::GradientContainer::FieldDescription descr = {
                  "STRESS_SHELL_MATERIAL",
                  "Mxx.Myy.Mxy.Bxx.Byy.Bxy.Tyz.Txz",
                  "Second Piola-Kirchhoff/linear stress, material "
                  "reference frame, no thermal expansion",
                  3,
                  8,
                  -1,
                  -1,
                  false,
                  typeid(T)};
            if (gradient_container->compute_field_value(descr.name)) {
                gradient_container->set_field_value(descr, stress_pk2_mat);
            }
        }
    }
};
 
template <typename T>
struct ElasticityHelper2D {
    bool has_thermal_expansion;
 
    bool has_material_ref;
 
    T M[2][2];
 
    bool has_covariant_base;
 
    T G[2][2];
 
    T G_d[2][2];
 
    T strain_gl_bg[3];
 
    T strain_gl_bg_te[3];
 
    T strain_gl_mat[3];
 
    T strain_gl_mat_te[3];
 
    T strain_normal;
 
    T stress_pk2_bg[3];
 
    T stress_pk2_mat[3];
 
    T stress_cauchy_bg[3];
 
    T stress_cauchy_mat[3];
 
    T stress_normal;
 
    T failure_index;
 
    ElasticityHelper2D(const T material_ref[2][2], const T covariant_base[2][2])
        : has_thermal_expansion(false),
          has_material_ref(material_ref != 0),
          has_covariant_base(covariant_base != 0),
          strain_normal(0),
          stress_normal(0),
          failure_index(0) {
        if (has_material_ref) {
            assert(material_ref != 0);
            std::copy(&material_ref[0][0], &material_ref[0][0] + 4, &M[0][0]);
        } else {
            set_identity_2_2(M);
        }
 
        if (has_covariant_base) {
            assert(covariant_base != 0);
            b2000::copy_2_2(covariant_base, G);
            b2000::get_base_opposite_variance(G, G_d);
        } else {
            // Everything is aligned with the global coordinate system.
            b2000::set_identity_2_2(G);
            b2000::set_identity_2_2(G_d);
        }
 
        std::fill_n(strain_gl_bg, 3, 0);
        std::fill_n(strain_gl_bg_te, 3, 0);
        std::fill_n(strain_gl_mat, 3, 0);
        std::fill_n(strain_gl_mat_te, 3, 0);
        std::fill_n(stress_pk2_bg, 3, 0);
        std::fill_n(stress_pk2_mat, 3, 0);
        std::fill_n(stress_cauchy_bg, 3, 0);
        std::fill_n(stress_cauchy_mat, 3, 0);
    }
 
    void set_strain(const T strain[3], const T temperature, const T alpha[3]) {
        // Compute the Linear/Green-Lagrange strain in the branch-global
        // and material reference frames.
        has_thermal_expansion = false;
        if (strain != 0) {
            transform_strain_from_base_A_to_I(G, strain, strain_gl_bg);
            transform_strain_from_I_to_base_B(M, strain_gl_bg, strain_gl_mat);
        }
        if (alpha != 0) {
            has_thermal_expansion = true;
            strain_gl_mat_te[0] = strain_gl_mat[0] - alpha[0] * temperature;
            strain_gl_mat_te[1] = strain_gl_mat[1] - alpha[1] * temperature;
            strain_gl_mat_te[2] = strain_gl_mat[2] - alpha[2] * temperature;
        } else {
            std::copy(strain_gl_mat, strain_gl_mat + 3, strain_gl_mat_te);
        }
    }
 
    void set_stress_and_or_d_stress_d_strain_mat(
          const T stress_mat[3],  // may be 0
          const T d_stress_d_strain_mat[6], T stress[3], T d_stress_d_strain[6]) {
        if (d_stress_d_strain) {
            assert(d_stress_d_strain_mat != 0);
            transform_constitutive_tensor_from_base_A_ov_to_base_B(
                  M, G_d, d_stress_d_strain_mat, d_stress_d_strain);
        }
 
        //
        // Calculate the stresses.
        //
        if (stress_mat != 0) {
            std::copy(stress_mat, stress_mat + 3, stress_pk2_mat);
        } else {
            // Linear material: Calculate stress from constitutive matrix.
            assert(d_stress_d_strain_mat != 0);
            b2linalg::Matrix<T, b2linalg::Mlower_packed_constref> C_mat(3, d_stress_d_strain_mat);
            b2linalg::Vector<T, b2linalg::Vdense_constref> e_mat(3, strain_gl_mat_te);
            b2linalg::Vector<T, b2linalg::Vdense_ref> s_mat(3, stress_pk2_mat);
            s_mat = C_mat * e_mat;
        }
        transform_stress_from_base_A_ov_to_I(M, stress_pk2_mat, stress_pk2_bg);
        if (stress != 0) { transform_stress_from_I_to_base_B(G_d, stress_pk2_bg, stress); }
    }
 
    void store(
          const bool linear, const T displacement_gradient[2][2],
          const FailureCriterion* failure_criterion, b2000::GradientContainer* gradient_container) {
        // Failure criterion.
        if (failure_criterion) {
            if (failure_criterion->get_type() == FailureCriterion::STRESS) {
                const T stress_pk2_mat_6[6] = {
                      stress_pk2_mat[0], stress_pk2_mat[1], stress_normal, stress_pk2_mat[2], 0, 0};
                failure_index =
                      failure_criterion->process_stress(stress_pk2_mat_6, gradient_container);
            } else {
                const T strain_gl_mat_6[6] = {
                      strain_gl_mat[0], strain_gl_mat[1], strain_normal, strain_gl_mat[2], 0, 0};
                failure_index =
                      failure_criterion->process_strain(strain_gl_mat_6, false, gradient_container);
            }
        }
 
        if (!gradient_container) { return; }
 
        // Transform the PK2 stress to the Cauchy stress.
        if (!linear && displacement_gradient) {
            T deformation_gradient[2][2] = {};
            if (has_covariant_base) {
                b2000::inner_product_2_2_NT(displacement_gradient, G_d, deformation_gradient);
            } else {
                std::copy(
                      &displacement_gradient[0][0], &displacement_gradient[0][0] + 4,
                      &deformation_gradient[0][0]);
            }
            deformation_gradient[0][0] += 1;
            deformation_gradient[1][1] += 1;
            const T det = 1 / b2000::determinant_2_2(deformation_gradient);
            if (has_material_ref) {
                T deformation_gradient_aligned_ortho[2][2];
                b2000::inner_product_2_2_NN(
                      deformation_gradient, M, deformation_gradient_aligned_ortho);
                b2000::transform_stress_from_base_A_ov_to_I(
                      deformation_gradient_aligned_ortho, stress_pk2_mat, stress_cauchy_bg);
            } else {
                b2000::transform_stress_from_base_A_ov_to_I(
                      deformation_gradient, stress_pk2_mat, stress_cauchy_bg);
            }
            for (int i = 0; i != 3; ++i) { stress_cauchy_bg[i] *= det; }
        } else {
            std::copy(stress_pk2_bg, stress_pk2_bg + 3, stress_cauchy_bg);
        }
 
        //
        // Store MBASE.
        //
        if (has_material_ref) {
            static const b2000::GradientContainer::FieldDescription descr = {
                  "MBASE_IP",
                  "e1x.e1y.e1z.e2x.e2y.e2z.e3x.e3y.e3z",
                  "Material base vectors, undeformed configuration",
                  3,
                  9,
                  1,
                  3,
                  false,
                  typeid(T)};
            if (gradient_container->compute_field_value(descr.name)) {
                T mbase_3_3[3][3] = {
                      {M[0][0], M[0][1], 0.},
                      {M[1][0], M[1][1], 0.},
                      {0., 0., 1.},
                };
                gradient_container->set_field_value(descr, &mbase_3_3[0][0]);
            }
        }
 
        //
        // Store the different strain and stress tensors.
        //
        if (1) {
            static const b2000::GradientContainer::FieldDescription descr = {
                  "STRAIN",
                  "Exx.Eyy.Ezz.Exy.Eyz.Exz",
                  "Green-Lagrange/linear strain, global reference "
                  "frame, no thermal expansion",
                  3,
                  6,
                  2,
                  3,
                  true,
                  typeid(T)};
            if (gradient_container->compute_field_value(descr.name)) {
                set_field_value_strain_2(gradient_container, descr, strain_gl_bg, strain_normal);
            }
        }
 
        if (0) {
            static const b2000::GradientContainer::FieldDescription descr = {
                  "STRAIN_MATERIAL",
                  "Exx.Eyy.Ezz.Exy.Eyz.Exz",
                  "Green-Lagrange/linear strain, material reference "
                  "frame, no thermal expansion",
                  3,
                  6,
                  2,
                  3,
                  true,
                  typeid(T)};
            if (gradient_container->compute_field_value(descr.name)) {
                set_field_value_strain_2(gradient_container, descr, strain_gl_mat, strain_normal);
            }
        }
 
        if (0 && has_thermal_expansion) {
            static const b2000::GradientContainer::FieldDescription descr = {
                  "STRAIN_MECHANICAL",
                  "Exx.Eyy.Ezz.Exy.Eyz.Exz",
                  "Green-Lagrange/linear strain, global reference "
                  "frame, thermal expansion",
                  3,
                  6,
                  2,
                  3,
                  true,
                  typeid(T)};
            if (gradient_container->compute_field_value(descr.name)) {
                b2000::transform_strain_from_base_A_ov_to_I(M, strain_gl_mat_te, strain_gl_bg_te);
                set_field_value_strain_2(gradient_container, descr, strain_gl_bg_te, strain_normal);
            }
        }
 
        if (0 && has_thermal_expansion) {
            static const b2000::GradientContainer::FieldDescription descr = {
                  "STRAIN_MECHANICAL_MATERIAL",
                  "Exx.Eyy.Ezz.Exy.Eyz.Exz",
                  "Green-Lagrange/linear strain, material reference "
                  "frame, thermal expansion",
                  3,
                  6,
                  2,
                  3,
                  true,
                  typeid(T)};
            if (gradient_container->compute_field_value(descr.name)) {
                set_field_value_strain_2(
                      gradient_container, descr, strain_gl_mat_te, strain_normal);
            }
        }
 
        if (1) {
            static const b2000::GradientContainer::FieldDescription descr = {
                  "STRESS",
                  "Sxx.Syy.Szz.Sxy.Syz.Sxz",
                  "Cauchy/linear stress, global reference frame",
                  3,
                  6,
                  2,
                  3,
                  true,
                  typeid(T)};
            if (gradient_container->compute_field_value(descr.name)) {
                set_field_value_stress_2(
                      gradient_container, descr, stress_cauchy_bg, stress_normal);
            }
        }
 
        if (0) {
            static const b2000::GradientContainer::FieldDescription descr = {
                  "STRESS_MATERIAL",
                  "Sxx.Syy.Szz.Sxy.Syz.Sxz",
                  "Cauchy/linear stress, material reference frame",
                  3,
                  6,
                  2,
                  3,
                  true,
                  typeid(T)};
            if (gradient_container->compute_field_value(descr.name)) {
                b2000::transform_stress_from_I_to_base_B(M, stress_cauchy_bg, stress_cauchy_mat);
                set_field_value_stress_2(
                      gradient_container, descr, stress_cauchy_mat, stress_normal);
            }
        }
 
        if (0) {
            static const b2000::GradientContainer::FieldDescription descr = {
                  "STRESS_PK2",
                  "Sxx.Syy.Szz.Sxy.Syz.Sxz",
                  "Second-Piola Kirchhoff stress, global reference frame",
                  3,
                  6,
                  2,
                  3,
                  true,
                  typeid(T)};
            if (gradient_container->compute_field_value(descr.name)) {
                b2000::transform_stress_from_base_A_ov_to_I(M, stress_pk2_mat, stress_pk2_bg);
                set_field_value_stress_2(gradient_container, descr, stress_pk2_bg, stress_normal);
            }
        }
 
        if (0) {
            static const b2000::GradientContainer::FieldDescription descr = {
                  "STRESS_PK2_MATERIAL",
                  "Sxx.Syy.Szz.Sxy.Syz.Sxz",
                  "Second-Piola Kirchhoff stress, material reference frame",
                  3,
                  6,
                  2,
                  3,
                  true,
                  typeid(T)};
            if (gradient_container->compute_field_value(descr.name)) {
                set_field_value_stress_2(gradient_container, descr, stress_pk2_mat, stress_normal);
            }
        }
 
        if (0) {
            static const b2000::GradientContainer::FieldDescription descr = {
                  "ENERGY_DENSITY",
                  "Elastic.Kinetic.RateDissipated.Dissipated",
                  "Free elastic, kinetic, rate of dissipated energy and "
                  "dissipated energy, per unit volume, in the undeformed "
                  "configuration",
                  3,
                  4,
                  -1,
                  0,
                  false,
                  typeid(T)};
            if (gradient_container->compute_field_value(descr.name)) {
                const T e[4] = {0.5 * dot_3(strain_gl_mat, stress_pk2_mat), 0, 0, 0};
                gradient_container->set_field_value(descr, e);
            }
        }
    }
};
 
template <typename T>
struct LayerData {
    std::string not_supported;
    bool isotropic;
    bool orthotropic;
    bool anisotropic;
    T d_stress_d_strain_mat[21];
    T alpha[6];
    T density;
 
    LayerData() : isotropic(true), orthotropic(false), anisotropic(false), density(0) {
        std::fill_n(d_stress_d_strain_mat, 21, 0.);
        std::fill_n(alpha, 6, 0.);
    }
 
    void init_25d(const RTable& rtable, const std::string& prefix = "") {
        const std::string sub_type = rtable.get_string(prefix + "SUBTYPE");
        isotropic = (sub_type == "ISOTROPIC");
        orthotropic = (sub_type == "ORTHOTROPIC" || sub_type == "ORTHOTROPIC25D");
        anisotropic = (sub_type == "ANISOTROPIC");
 
        const T sh = rtable.get<T>(prefix + "SHEAR_CORRECTION_FACTOR", T(5) / T(6));
        if (anisotropic) {
            // Anisotropic material.  Material matrix is defined in
            // database and is stored in row-major format, upper-triangle.
            // This is equal to column-major format, lower-triangle.
            d_stress_d_strain_mat[0] = rtable.get<T>(prefix + "C1111", 0.);
            d_stress_d_strain_mat[1] = rtable.get<T>(prefix + "C1122", 0.);
            d_stress_d_strain_mat[2] = rtable.get<T>(prefix + "C1133", 0.);
            d_stress_d_strain_mat[3] = rtable.get<T>(prefix + "C1112", 0.);
            d_stress_d_strain_mat[4] = rtable.get<T>(prefix + "C1123", 0.);
            d_stress_d_strain_mat[5] = rtable.get<T>(prefix + "C1113", 0.);
            d_stress_d_strain_mat[6] = rtable.get<T>(prefix + "C2222", 0.);
            d_stress_d_strain_mat[7] = rtable.get<T>(prefix + "C2233", 0.);
            d_stress_d_strain_mat[8] = rtable.get<T>(prefix + "C2212", 0.);
            d_stress_d_strain_mat[9] = rtable.get<T>(prefix + "C2223", 0.);
            d_stress_d_strain_mat[10] = rtable.get<T>(prefix + "C2213", 0.);
            d_stress_d_strain_mat[11] = rtable.get<T>(prefix + "C3333", 0.);
            d_stress_d_strain_mat[12] = rtable.get<T>(prefix + "C3312", 0.);
            d_stress_d_strain_mat[13] = rtable.get<T>(prefix + "C3323", 0.);
            d_stress_d_strain_mat[14] = rtable.get<T>(prefix + "C3313", 0.);
            d_stress_d_strain_mat[15] = rtable.get<T>(prefix + "C1212", 0.);
            d_stress_d_strain_mat[16] = rtable.get<T>(prefix + "C1223", 0.);
            d_stress_d_strain_mat[17] = rtable.get<T>(prefix + "C1213", 0.);
            d_stress_d_strain_mat[18] = rtable.get<T>(prefix + "C2323", 0.);
            d_stress_d_strain_mat[19] = rtable.get<T>(prefix + "C2313", 0.);
            d_stress_d_strain_mat[20] = rtable.get<T>(prefix + "C1313", 0.);
        } else if (isotropic || orthotropic) {
            // Extract parameters.
            T e1, e2, e3;
            T nu12, nu23, nu13;
            T g12, g23, g13;
            if (orthotropic) {
                e1 = rtable.get<T>(prefix + "E1");
                e2 = rtable.get<T>(prefix + "E2");
                e3 = rtable.get<T>(prefix + "E3", e2);
                nu12 = rtable.get<T>(prefix + "NU12");
                nu13 = rtable.get<T>(prefix + "NU13", nu12);
                nu23 = rtable.get<T>(prefix + "NU23", nu13);
                g12 = rtable.get<T>(prefix + "G12");
                g23 = rtable.get<T>(prefix + "G23", g12);
                g13 = rtable.get<T>(prefix + "G13", g12);
            } else {
                e1 = e2 = e3 = rtable.get<T>(prefix + "E");
                nu12 = nu23 = nu13 = rtable.get<T>(prefix + "NU");
                g12 = g23 = g13 = rtable.get<T>(prefix + "G", e1 / (T(2) * (T(1) + nu12)));
            }
 
            // Initialize constitutive matrix.
            std::fill_n(d_stress_d_strain_mat, 21, T(0));
            {
                const T c11 = T(1) / e1;
                const T c22 = T(1) / e2;
                const T c12 = -nu12 * c11;
                const T inv_det = T(1) / (c11 * c22 - c12 * c12);
                d_stress_d_strain_mat[0] = c22 * inv_det;
                d_stress_d_strain_mat[1] = -c12 * inv_det;
                d_stress_d_strain_mat[6] = c11 * inv_det;
            }
            d_stress_d_strain_mat[15] = g12;
            d_stress_d_strain_mat[18] = g23 * sh;
            d_stress_d_strain_mat[20] = g13 * sh;
        } else {
            std::ostringstream o;
            o << rtable.errmsg_context << "The material type " << sub_type
              << " cannot be used with shell elements.";
            not_supported = o.str();
            return;
        }
 
        const T alpha_r = rtable.get<T>(prefix + "ALPHA", T(0));
        alpha[0] = rtable.get<T>(prefix + "ALPHA11", alpha_r);
        alpha[1] = rtable.get<T>(prefix + "ALPHA22", alpha_r);
        alpha[2] = T(0);  // no valid transverse normal strain
        const T alpha12 = rtable.get<T>(prefix + "ALPHA12", T(0));
        alpha[3] = T(2) * alpha12;
        alpha[4] = T(2) * rtable.get<T>(prefix + "ALPHA23", T(0));
        alpha[5] = T(2) * rtable.get<T>(prefix + "ALPHA13", T(0));
 
        density = rtable.get<T>(prefix + "DENSITY", T(0));
    }
};
 
}}  // namespace b2000::b2dbv3
 
#endif  // B2_MATERIAL_UTIL_H_