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