Skip to main content

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

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