conspire 0.7.7

//! Elastic-hyperviscous solid constitutive models.
//!
//! ---
//!
#![doc = include_str!("doc.md")]

#[cfg(feature = "doc")]
pub mod doc;

#[cfg(test)]
pub mod test;

mod almansi_hamel;

pub use almansi_hamel::AlmansiHamel;

use super::{
    viscoelastic::{AppliedLoad, Viscoelastic},
    *,
};
use crate::{
    math::{
        ContractWith, Matrix, Quantity, Vector,
        integrate::{ImplicitDaeFirstOrderMinimize, ImplicitDaeSecondOrderMinimize},
        optimize::{EqualityConstraint, FirstOrderOptimization, SecondOrderOptimization},
    },
    units::{Dissipation, Time},
};

/// Required methods for elastic-hyperviscous solid constitutive models.
pub trait ElasticHyperviscous
where
    Self: Viscoelastic,
{
    /// Calculates and returns the dissipation potential.
    ///
    /// ```math
    /// \phi(\mathbf{F},\dot{\mathbf{F}}) = \mathbf{P}^e(\mathbf{F}):\dot{\mathbf{F}} + \psi(\mathbf{F},\dot{\mathbf{F}})
    /// ```
    fn dissipation_potential(
        &self,
        deformation_gradient: &DeformationGradient,
        deformation_gradient_rate: &DeformationGradientRate,
    ) -> Result<Quantity<Dissipation>, ConstitutiveError> {
        Ok(self
            .first_piola_kirchhoff_stress(deformation_gradient, &DeformationGradientRate::zero())?
            .contract_with(deformation_gradient_rate)
            + self.viscous_dissipation(deformation_gradient, deformation_gradient_rate)?)
    }
    /// Calculates and returns the internal dissipation.
    ///
    /// ```math
    /// T\dot{s} = \mathbf{P}(\mathbf{F},\dot{\mathbf{F}}):\dot{\mathbf{F}} - \mathbf{P}^e(\mathbf{F}):\dot{\mathbf{F}}
    /// ```
    fn internal_dissipation(
        &self,
        deformation_gradient: &DeformationGradient,
        deformation_gradient_rate: &DeformationGradientRate,
    ) -> Result<Quantity<Dissipation>, ConstitutiveError> {
        Ok(self
            .first_piola_kirchhoff_stress(deformation_gradient, deformation_gradient_rate)?
            .contract_with(deformation_gradient_rate)
            - self
                .first_piola_kirchhoff_stress(
                    deformation_gradient,
                    &DeformationGradientRate::zero(),
                )?
                .contract_with(deformation_gradient_rate))
    }
    /// Calculates and returns the viscous dissipation.
    ///
    /// ```math
    /// \psi = \psi(\mathbf{F},\dot{\mathbf{F}})
    /// ```
    fn viscous_dissipation(
        &self,
        deformation_gradient: &DeformationGradient,
        deformation_gradient_rate: &DeformationGradientRate,
    ) -> Result<Quantity<Dissipation>, ConstitutiveError>;
}

/// First-order optimization methods for elastic-hyperviscous solid constitutive models.
pub trait FirstOrderMinimize {
    /// Solve for the unknown components of the deformation gradient and rate under an applied load.
    ///
    /// ```math
    /// \Pi(\mathbf{F},\dot{\mathbf{F}},\boldsymbol{\lambda}) = \mathbf{P}^e(\mathbf{F}):\dot{\mathbf{F}} + \psi(\mathbf{F},\dot{\mathbf{F}}) - \boldsymbol{\lambda}:(\dot{\mathbf{F}} - \dot{\mathbf{F}}_0) - \mathbf{P}_0:\dot{\mathbf{F}}
    /// ```
    fn minimize(
        &self,
        applied_load: AppliedLoad,
        integrator: impl ImplicitDaeFirstOrderMinimize<
            Quantity<Dissipation>,
            FirstPiolaKirchhoffStress,
            DeformationGradient,
            DeformationGradients,
            DeformationGradientRates,
        >,
        solver: impl FirstOrderOptimization<
            Quantity<Dissipation>,
            FirstPiolaKirchhoffStress,
            DeformationGradientRate,
        >,
    ) -> Result<(Times, DeformationGradients, DeformationGradientRates), ConstitutiveError>;
}

/// Second-order optimization methods for elastic-hyperviscous solid constitutive models.
pub trait SecondOrderMinimize {
    /// Solve for the unknown components of the deformation gradient and rate under an applied load.
    ///
    /// ```math
    /// \Pi(\mathbf{F},\dot{\mathbf{F}},\boldsymbol{\lambda}) = \mathbf{P}^e(\mathbf{F}):\dot{\mathbf{F}} + \psi(\mathbf{F},\dot{\mathbf{F}}) - \boldsymbol{\lambda}:(\dot{\mathbf{F}} - \dot{\mathbf{F}}_0) - \mathbf{P}_0:\dot{\mathbf{F}}
    /// ```
    fn minimize(
        &self,
        applied_load: AppliedLoad,
        integrator: impl ImplicitDaeSecondOrderMinimize<
            Quantity<Dissipation>,
            FirstPiolaKirchhoffStress,
            FirstPiolaKirchhoffRateTangentStiffness,
            DeformationGradient,
            DeformationGradients,
            DeformationGradientRates,
        >,
        solver: impl SecondOrderOptimization<
            Quantity<Dissipation>,
            FirstPiolaKirchhoffStress,
            FirstPiolaKirchhoffRateTangentStiffness,
            DeformationGradientRate,
        >,
    ) -> Result<(Times, DeformationGradients, DeformationGradientRates), ConstitutiveError>;
}

impl<T> FirstOrderMinimize for T
where
    T: ElasticHyperviscous,
{
    fn minimize(
        &self,
        applied_load: AppliedLoad,
        integrator: impl ImplicitDaeFirstOrderMinimize<
            Quantity<Dissipation>,
            FirstPiolaKirchhoffStress,
            DeformationGradient,
            DeformationGradients,
            DeformationGradientRates,
        >,
        solver: impl FirstOrderOptimization<
            Quantity<Dissipation>,
            FirstPiolaKirchhoffStress,
            DeformationGradientRate,
        >,
    ) -> Result<(Times, DeformationGradients, DeformationGradientRates), ConstitutiveError> {
        match applied_load {
            AppliedLoad::UniaxialStress(deformation_gradient_rate_11, time) => {
                let mut matrix = Matrix::zero(4, 9);
                let mut vector = Vector::zero(4);
                matrix[0][0] = 1.0;
                matrix[1][1] = 1.0;
                matrix[2][2] = 1.0;
                matrix[3][5] = 1.0;
                integrator.integrate(
                    |_: Quantity<Time>,
                     deformation_gradient: &DeformationGradient,
                     deformation_gradient_rate: &DeformationGradientRate| {
                        Ok(self.dissipation_potential(
                            deformation_gradient,
                            deformation_gradient_rate,
                        )?)
                    },
                    |_: Quantity<Time>,
                     deformation_gradient: &DeformationGradient,
                     deformation_gradient_rate: &DeformationGradientRate| {
                        Ok(self.first_piola_kirchhoff_stress(
                            deformation_gradient,
                            deformation_gradient_rate,
                        )?)
                    },
                    solver,
                    time,
                    DeformationGradient::identity(),
                    |t: Quantity<Time>| {
                        vector[0] = deformation_gradient_rate_11(t);
                        EqualityConstraint::Linear(matrix.clone(), vector.clone())
                    },
                )
            }
            AppliedLoad::BiaxialStress(
                deformation_gradient_rate_11,
                deformation_gradient_rate_22,
                time,
            ) => {
                let mut matrix = Matrix::zero(5, 9);
                let mut vector = Vector::zero(5);
                matrix[0][0] = 1.0;
                matrix[1][1] = 1.0;
                matrix[2][2] = 1.0;
                matrix[3][5] = 1.0;
                matrix[4][4] = 1.0;
                integrator.integrate(
                    |_: Quantity<Time>,
                     deformation_gradient: &DeformationGradient,
                     deformation_gradient_rate: &DeformationGradientRate| {
                        Ok(self.dissipation_potential(
                            deformation_gradient,
                            deformation_gradient_rate,
                        )?)
                    },
                    |_: Quantity<Time>,
                     deformation_gradient: &DeformationGradient,
                     deformation_gradient_rate: &DeformationGradientRate| {
                        Ok(self.first_piola_kirchhoff_stress(
                            deformation_gradient,
                            deformation_gradient_rate,
                        )?)
                    },
                    solver,
                    time,
                    DeformationGradient::identity(),
                    |t: Quantity<Time>| {
                        vector[0] = deformation_gradient_rate_11(t);
                        vector[4] = deformation_gradient_rate_22(t);
                        EqualityConstraint::Linear(matrix.clone(), vector.clone())
                    },
                )
            }
        }
        .map_err(|error| ConstitutiveError::upstream(error, self))
    }
}

impl<T> SecondOrderMinimize for T
where
    T: ElasticHyperviscous,
{
    fn minimize(
        &self,
        applied_load: AppliedLoad,
        integrator: impl ImplicitDaeSecondOrderMinimize<
            Quantity<Dissipation>,
            FirstPiolaKirchhoffStress,
            FirstPiolaKirchhoffRateTangentStiffness,
            DeformationGradient,
            DeformationGradients,
            DeformationGradientRates,
        >,
        solver: impl SecondOrderOptimization<
            Quantity<Dissipation>,
            FirstPiolaKirchhoffStress,
            FirstPiolaKirchhoffRateTangentStiffness,
            DeformationGradientRate,
        >,
    ) -> Result<(Times, DeformationGradients, DeformationGradientRates), ConstitutiveError> {
        match applied_load {
            AppliedLoad::UniaxialStress(deformation_gradient_rate_11, time) => {
                let mut matrix = Matrix::zero(4, 9);
                let mut vector = Vector::zero(4);
                matrix[0][0] = 1.0;
                matrix[1][1] = 1.0;
                matrix[2][2] = 1.0;
                matrix[3][5] = 1.0;
                integrator.integrate(
                    |_: Quantity<Time>,
                     deformation_gradient: &DeformationGradient,
                     deformation_gradient_rate: &DeformationGradientRate| {
                        Ok(self.dissipation_potential(
                            deformation_gradient,
                            deformation_gradient_rate,
                        )?)
                    },
                    |_: Quantity<Time>,
                     deformation_gradient: &DeformationGradient,
                     deformation_gradient_rate: &DeformationGradientRate| {
                        Ok(self.first_piola_kirchhoff_stress(
                            deformation_gradient,
                            deformation_gradient_rate,
                        )?)
                    },
                    |_: Quantity<Time>,
                     deformation_gradient: &DeformationGradient,
                     deformation_gradient_rate: &DeformationGradientRate| {
                        Ok(self.first_piola_kirchhoff_rate_tangent_stiffness(
                            deformation_gradient,
                            deformation_gradient_rate,
                        )?)
                    },
                    solver,
                    time,
                    DeformationGradient::identity(),
                    |t: Quantity<Time>| {
                        vector[0] = deformation_gradient_rate_11(t);
                        EqualityConstraint::Linear(matrix.clone(), vector.clone())
                    },
                    None,
                )
            }
            AppliedLoad::BiaxialStress(
                deformation_gradient_rate_11,
                deformation_gradient_rate_22,
                time,
            ) => {
                let mut matrix = Matrix::zero(5, 9);
                let mut vector = Vector::zero(5);
                matrix[0][0] = 1.0;
                matrix[1][1] = 1.0;
                matrix[2][2] = 1.0;
                matrix[3][5] = 1.0;
                matrix[4][4] = 1.0;
                integrator.integrate(
                    |_: Quantity<Time>,
                     deformation_gradient: &DeformationGradient,
                     deformation_gradient_rate: &DeformationGradientRate| {
                        Ok(self.dissipation_potential(
                            deformation_gradient,
                            deformation_gradient_rate,
                        )?)
                    },
                    |_: Quantity<Time>,
                     deformation_gradient: &DeformationGradient,
                     deformation_gradient_rate: &DeformationGradientRate| {
                        Ok(self.first_piola_kirchhoff_stress(
                            deformation_gradient,
                            deformation_gradient_rate,
                        )?)
                    },
                    |_: Quantity<Time>,
                     deformation_gradient: &DeformationGradient,
                     deformation_gradient_rate: &DeformationGradientRate| {
                        Ok(self.first_piola_kirchhoff_rate_tangent_stiffness(
                            deformation_gradient,
                            deformation_gradient_rate,
                        )?)
                    },
                    solver,
                    time,
                    DeformationGradient::identity(),
                    |t: Quantity<Time>| {
                        vector[0] = deformation_gradient_rate_11(t);
                        vector[4] = deformation_gradient_rate_22(t);
                        EqualityConstraint::Linear(matrix.clone(), vector.clone())
                    },
                    None,
                )
            }
        }
        .map_err(|error| ConstitutiveError::upstream(error, self))
    }
}