conspire/domain/fem/block/solid/hyperelastic/
mod.rs1use 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}