conspire/domain/fem/block/solid/hyperelastic/
mod.rs

1use crate::{
2    constitutive::solid::hyperelastic::Hyperelastic,
3    fem::{
4        NodalCoordinates,
5        block::{
6            Block, FiniteElementBlockError, FirstOrderMinimize, SecondOrderMinimize, band,
7            element::solid::hyperelastic::HyperelasticFiniteElement,
8            solid::{NodalForcesSolid, NodalStiffnessesSolid, elastic::ElasticFiniteElementBlock},
9        },
10    },
11    math::{
12        Scalar, Tensor,
13        optimize::{
14            EqualityConstraint, FirstOrderOptimization, OptimizationError, SecondOrderOptimization,
15        },
16    },
17};
18
19pub trait HyperelasticFiniteElementBlock<
20    C,
21    F,
22    const G: usize,
23    const M: usize,
24    const N: usize,
25    const P: usize,
26> where
27    C: Hyperelastic,
28    F: HyperelasticFiniteElement<C, G, M, N, P>,
29    Self: ElasticFiniteElementBlock<C, F, G, M, N, P>,
30{
31    fn helmholtz_free_energy(
32        &self,
33        nodal_coordinates: &NodalCoordinates,
34    ) -> Result<Scalar, FiniteElementBlockError>;
35}
36
37impl<C, F, const G: usize, const M: usize, const N: usize, const P: usize>
38    HyperelasticFiniteElementBlock<C, F, G, M, N, P> for Block<C, F, G, M, N, P>
39where
40    C: Hyperelastic,
41    F: HyperelasticFiniteElement<C, G, M, N, P>,
42    Self: ElasticFiniteElementBlock<C, F, G, M, N, P>,
43{
44    fn helmholtz_free_energy(
45        &self,
46        nodal_coordinates: &NodalCoordinates,
47    ) -> Result<Scalar, FiniteElementBlockError> {
48        match self
49            .elements()
50            .iter()
51            .zip(self.connectivity())
52            .map(|(element, nodes)| {
53                element.helmholtz_free_energy(
54                    self.constitutive_model(),
55                    &Self::element_coordinates(nodal_coordinates, nodes),
56                )
57            })
58            .sum()
59        {
60            Ok(helmholtz_free_energy) => Ok(helmholtz_free_energy),
61            Err(error) => Err(FiniteElementBlockError::Upstream(
62                format!("{error}"),
63                format!("{self:?}"),
64            )),
65        }
66    }
67}
68
69impl<C, F, const G: usize, const M: usize, const N: usize, const P: usize>
70    FirstOrderMinimize<C, F, G, M, N, NodalCoordinates> for Block<C, F, G, M, N, P>
71where
72    C: Hyperelastic,
73    F: HyperelasticFiniteElement<C, G, M, N, P>,
74{
75    fn minimize(
76        &self,
77        equality_constraint: EqualityConstraint,
78        solver: impl FirstOrderOptimization<Scalar, NodalForcesSolid>,
79    ) -> Result<NodalCoordinates, OptimizationError> {
80        solver.minimize(
81            |nodal_coordinates: &NodalCoordinates| {
82                Ok(self.helmholtz_free_energy(nodal_coordinates)?)
83            },
84            |nodal_coordinates: &NodalCoordinates| Ok(self.nodal_forces(nodal_coordinates)?),
85            self.coordinates().clone().into(),
86            equality_constraint,
87        )
88    }
89}
90
91impl<C, F, const G: usize, const M: usize, const N: usize, const P: usize>
92    SecondOrderMinimize<C, F, G, M, N, NodalForcesSolid, NodalStiffnessesSolid, NodalCoordinates>
93    for Block<C, F, G, M, N, P>
94where
95    C: Hyperelastic,
96    F: HyperelasticFiniteElement<C, G, M, N, P>,
97{
98    fn minimize(
99        &self,
100        equality_constraint: EqualityConstraint,
101        solver: impl SecondOrderOptimization<
102            Scalar,
103            NodalForcesSolid,
104            NodalStiffnessesSolid,
105            NodalCoordinates,
106        >,
107    ) -> Result<NodalCoordinates, OptimizationError> {
108        let banded = band(
109            self.connectivity(),
110            &equality_constraint,
111            self.coordinates().len(),
112            3,
113        );
114        solver.minimize(
115            |nodal_coordinates: &NodalCoordinates| {
116                Ok(self.helmholtz_free_energy(nodal_coordinates)?)
117            },
118            |nodal_coordinates: &NodalCoordinates| Ok(self.nodal_forces(nodal_coordinates)?),
119            |nodal_coordinates: &NodalCoordinates| Ok(self.nodal_stiffnesses(nodal_coordinates)?),
120            self.coordinates().clone().into(),
121            equality_constraint,
122            Some(banded),
123        )
124    }
125}