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