Skip to main content

conspire/domain/solid/hyperelastic/
mod.rs

1use crate::{
2    domain::{
3        Blocks, ElementModel, ElementModelError, FirstOrderMinimize, Model, NodalCoordinates,
4        SecondOrderMinimize,
5        block::{element::Elements, finalize_node_neighbors, solver_from_neighbors},
6        solid::{NodalForcesSolid, NodalStiffnessesSolidSymmetric, elastic::ElasticElements},
7    },
8    math::{
9        Quantity, Tensor,
10        optimize::{
11            EqualityConstraint, FirstOrderOptimization, OptimizationError, SecondOrderOptimization,
12        },
13    },
14    units::Energy,
15};
16
17pub trait HyperelasticElements<const D: usize>
18where
19    Self: ElasticElements<D>,
20{
21    fn helmholtz_free_energy(
22        &self,
23        nodal_coordinates: &NodalCoordinates<D>,
24    ) -> Result<Quantity<Energy>, ElementModelError>;
25    fn nodal_stiffnesses_symmetric_into(
26        &self,
27        nodal_coordinates: &NodalCoordinates<D>,
28        nodal_stiffnesses: &mut NodalStiffnessesSolidSymmetric<D>,
29    ) -> Result<(), ElementModelError>;
30    fn nodal_stiffnesses_symmetric(
31        &self,
32        nodal_coordinates: &NodalCoordinates<D>,
33    ) -> Result<NodalStiffnessesSolidSymmetric<D>, ElementModelError> {
34        let mut nodal_stiffnesses = NodalStiffnessesSolidSymmetric::zero(nodal_coordinates.len());
35        self.nodal_stiffnesses_symmetric_into(nodal_coordinates, &mut nodal_stiffnesses)?;
36        Ok(nodal_stiffnesses)
37    }
38}
39
40impl<B, const D: usize> HyperelasticElements<D> for Model<B, D>
41where
42    B: HyperelasticElements<D>,
43{
44    fn helmholtz_free_energy(
45        &self,
46        nodal_coordinates: &NodalCoordinates<D>,
47    ) -> Result<Quantity<Energy>, ElementModelError> {
48        self.blocks.helmholtz_free_energy(nodal_coordinates)
49    }
50    fn nodal_stiffnesses_symmetric_into(
51        &self,
52        nodal_coordinates: &NodalCoordinates<D>,
53        nodal_stiffnesses: &mut NodalStiffnessesSolidSymmetric<D>,
54    ) -> Result<(), ElementModelError> {
55        self.blocks
56            .nodal_stiffnesses_symmetric_into(nodal_coordinates, nodal_stiffnesses)
57    }
58}
59
60impl<B1, B2, const D: usize> HyperelasticElements<D> for Blocks<B1, B2>
61where
62    B1: HyperelasticElements<D>,
63    B2: HyperelasticElements<D>,
64{
65    fn helmholtz_free_energy(
66        &self,
67        nodal_coordinates: &NodalCoordinates<D>,
68    ) -> Result<Quantity<Energy>, ElementModelError> {
69        Ok(self.0.helmholtz_free_energy(nodal_coordinates)?
70            + self.1.helmholtz_free_energy(nodal_coordinates)?)
71    }
72    fn nodal_stiffnesses_symmetric_into(
73        &self,
74        nodal_coordinates: &NodalCoordinates<D>,
75        nodal_stiffnesses: &mut NodalStiffnessesSolidSymmetric<D>,
76    ) -> Result<(), ElementModelError> {
77        self.0
78            .nodal_stiffnesses_symmetric_into(nodal_coordinates, nodal_stiffnesses)?;
79        self.1
80            .nodal_stiffnesses_symmetric_into(nodal_coordinates, nodal_stiffnesses)
81    }
82}
83
84impl<B, const D: usize>
85    FirstOrderMinimize<Quantity<Energy>, NodalForcesSolid<D>, NodalCoordinates<D>> for Model<B, D>
86where
87    B: HyperelasticElements<D>,
88{
89    fn minimize(
90        &self,
91        equality_constraint: EqualityConstraint,
92        solver: impl FirstOrderOptimization<Quantity<Energy>, NodalForcesSolid<D>, NodalCoordinates<D>>,
93    ) -> Result<NodalCoordinates<D>, OptimizationError> {
94        solver.minimize(
95            |nodal_coordinates: &NodalCoordinates<D>| {
96                Ok(self.helmholtz_free_energy(nodal_coordinates)?)
97            },
98            |nodal_coordinates: &NodalCoordinates<D>| Ok(self.nodal_forces(nodal_coordinates)?),
99            self.coordinates().clone().into(),
100            equality_constraint,
101        )
102    }
103}
104
105impl<B, const D: usize>
106    SecondOrderMinimize<
107        Quantity<Energy>,
108        NodalForcesSolid<D>,
109        NodalStiffnessesSolidSymmetric<D>,
110        NodalCoordinates<D>,
111    > for Model<B, D>
112where
113    B: HyperelasticElements<D>,
114{
115    fn minimize(
116        &self,
117        equality_constraint: EqualityConstraint,
118        solver: impl SecondOrderOptimization<
119            Quantity<Energy>,
120            NodalForcesSolid<D>,
121            NodalStiffnessesSolidSymmetric<D>,
122            NodalCoordinates<D>,
123        >,
124    ) -> Result<NodalCoordinates<D>, OptimizationError> {
125        let mut neighbors = vec![Vec::new(); self.coordinates().len()];
126        self.node_neighbors(&mut neighbors);
127        finalize_node_neighbors(&mut neighbors);
128        let sparse = solver_from_neighbors(&neighbors, &equality_constraint, D, true);
129        solver.minimize(
130            |nodal_coordinates: &NodalCoordinates<D>| {
131                Ok(self.helmholtz_free_energy(nodal_coordinates)?)
132            },
133            |nodal_coordinates: &NodalCoordinates<D>| Ok(self.nodal_forces(nodal_coordinates)?),
134            |nodal_coordinates: &NodalCoordinates<D>| {
135                Ok(self.nodal_stiffnesses_symmetric(nodal_coordinates)?)
136            },
137            self.coordinates().clone().into(),
138            equality_constraint,
139            Some(sparse),
140        )
141    }
142}