1use crate::{
2 constitutive::solid::elastic_hyperviscous::ElasticHyperviscous,
3 fem::{
4 NodalCoordinates, NodalCoordinatesHistory, NodalVelocities, NodalVelocitiesHistory,
5 block::{
6 Block, FiniteElementBlockError, band,
7 element::solid::elastic_hyperviscous::ElasticHyperviscousFiniteElement,
8 solid::{
9 NodalForcesSolid, NodalStiffnessesSolid,
10 viscoelastic::ViscoelasticFiniteElementBlock,
11 },
12 },
13 },
14 math::{
15 Scalar, Tensor,
16 integrate::{Explicit, IntegrationError},
17 optimize::{EqualityConstraint, OptimizationError, SecondOrderOptimization},
18 },
19 mechanics::Times,
20};
21
22pub trait ElasticHyperviscousFiniteElementBlock<
23 C,
24 F,
25 const G: usize,
26 const M: usize,
27 const N: usize,
28 const P: usize,
29> where
30 C: ElasticHyperviscous,
31 F: ElasticHyperviscousFiniteElement<C, G, M, N, P>,
32 Self: ViscoelasticFiniteElementBlock<C, F, G, M, N, P>,
33{
34 fn viscous_dissipation(
35 &self,
36 nodal_coordinates: &NodalCoordinates,
37 nodal_velocities: &NodalVelocities,
38 ) -> Result<Scalar, FiniteElementBlockError>;
39 fn dissipation_potential(
40 &self,
41 nodal_coordinates: &NodalCoordinates,
42 nodal_velocities: &NodalVelocities,
43 ) -> Result<Scalar, FiniteElementBlockError>;
44 fn minimize(
45 &self,
46 equality_constraint: EqualityConstraint,
47 integrator: impl Explicit<NodalVelocities, NodalVelocitiesHistory>,
48 time: &[Scalar],
49 solver: impl SecondOrderOptimization<
50 Scalar,
51 NodalForcesSolid,
52 NodalStiffnessesSolid,
53 NodalCoordinates,
54 >,
55 ) -> Result<(Times, NodalCoordinatesHistory, NodalVelocitiesHistory), IntegrationError>;
56 #[doc(hidden)]
57 fn minimize_inner(
58 &self,
59 equality_constraint: &EqualityConstraint,
60 nodal_coordinates: &NodalCoordinates,
61 solver: &impl SecondOrderOptimization<
62 Scalar,
63 NodalForcesSolid,
64 NodalStiffnessesSolid,
65 NodalCoordinates,
66 >,
67 initial_guess: &NodalVelocities,
68 ) -> Result<NodalVelocities, OptimizationError>;
69}
70
71impl<C, F, const G: usize, const M: usize, const N: usize, const P: usize>
72 ElasticHyperviscousFiniteElementBlock<C, F, G, M, N, P> for Block<C, F, G, M, N, P>
73where
74 C: ElasticHyperviscous,
75 F: ElasticHyperviscousFiniteElement<C, G, M, N, P>,
76 Self: ViscoelasticFiniteElementBlock<C, F, G, M, N, P>,
77{
78 fn viscous_dissipation(
79 &self,
80 nodal_coordinates: &NodalCoordinates,
81 nodal_velocities: &NodalVelocities,
82 ) -> Result<Scalar, FiniteElementBlockError> {
83 match self
84 .elements()
85 .iter()
86 .zip(self.connectivity())
87 .map(|(element, nodes)| {
88 element.viscous_dissipation(
89 self.constitutive_model(),
90 &Self::element_coordinates(nodal_coordinates, nodes),
91 &self.element_velocities(nodal_velocities, nodes),
92 )
93 })
94 .sum()
95 {
96 Ok(viscous_dissipation) => Ok(viscous_dissipation),
97 Err(error) => Err(FiniteElementBlockError::Upstream(
98 format!("{error}"),
99 format!("{self:?}"),
100 )),
101 }
102 }
103 fn dissipation_potential(
104 &self,
105 nodal_coordinates: &NodalCoordinates,
106 nodal_velocities: &NodalVelocities,
107 ) -> Result<Scalar, FiniteElementBlockError> {
108 match self
109 .elements()
110 .iter()
111 .zip(self.connectivity())
112 .map(|(element, nodes)| {
113 element.dissipation_potential(
114 self.constitutive_model(),
115 &Self::element_coordinates(nodal_coordinates, nodes),
116 &self.element_velocities(nodal_velocities, nodes),
117 )
118 })
119 .sum()
120 {
121 Ok(dissipation_potential) => Ok(dissipation_potential),
122 Err(error) => Err(FiniteElementBlockError::Upstream(
123 format!("{error}"),
124 format!("{self:?}"),
125 )),
126 }
127 }
128 fn minimize(
129 &self,
130 equality_constraint: EqualityConstraint,
131 integrator: impl Explicit<NodalVelocities, NodalVelocitiesHistory>,
132 time: &[Scalar],
133 solver: impl SecondOrderOptimization<
134 Scalar,
135 NodalForcesSolid,
136 NodalStiffnessesSolid,
137 NodalCoordinates,
138 >,
139 ) -> Result<(Times, NodalCoordinatesHistory, NodalVelocitiesHistory), IntegrationError> {
140 let mut solution = NodalVelocities::zero(self.coordinates().len());
141 integrator.integrate(
142 |_: Scalar, nodal_coordinates: &NodalCoordinates| {
143 solution = self.minimize_inner(
144 &equality_constraint,
145 nodal_coordinates,
146 &solver,
147 &solution,
148 )?;
149 Ok(solution.clone())
150 },
151 time,
152 self.coordinates().clone().into(),
153 )
154 }
155 fn minimize_inner(
156 &self,
157 equality_constraint: &EqualityConstraint,
158 nodal_coordinates: &NodalCoordinates,
159 solver: &impl SecondOrderOptimization<
160 Scalar,
161 NodalForcesSolid,
162 NodalStiffnessesSolid,
163 NodalVelocities,
164 >,
165 initial_guess: &NodalVelocities,
166 ) -> Result<NodalVelocities, OptimizationError> {
167 let banded = band(
168 self.connectivity(),
169 equality_constraint,
170 nodal_coordinates.len(),
171 3,
172 );
173 solver.minimize(
174 |nodal_velocities: &NodalVelocities| {
175 Ok(self.dissipation_potential(nodal_coordinates, nodal_velocities)?)
176 },
177 |nodal_velocities: &NodalVelocities| {
178 Ok(self.nodal_forces(nodal_coordinates, nodal_velocities)?)
179 },
180 |nodal_velocities: &NodalVelocities| {
181 Ok(self.nodal_stiffnesses(nodal_coordinates, nodal_velocities)?)
182 },
183 initial_guess.clone(),
184 equality_constraint.clone(),
185 Some(banded),
186 )
187 }
188}