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::{ImplicitDaeFirstOrderRoot, IntegrationError},
17 optimize::{EqualityConstraint, FirstOrderRootFinding},
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 ImplicitDaeFirstOrderRoot<
57 NodalForcesSolid,
58 NodalStiffnessesSolid,
59 NodalVelocities,
60 NodalVelocitiesHistory,
61 >,
62 time: &[Scalar],
63 solver: impl FirstOrderRootFinding<NodalForcesSolid, NodalStiffnessesSolid, NodalCoordinates>,
64 ) -> Result<(Times, NodalCoordinatesHistory, NodalVelocitiesHistory), IntegrationError>;
65}
66
67impl<C, F, const G: usize, const M: usize, const N: usize, const P: usize>
68 ViscoelasticFiniteElementBlock<C, F, G, M, N, P> for Block<C, F, G, M, N, P>
69where
70 C: Viscoelastic,
71 F: ViscoelasticFiniteElement<C, G, M, N, P>,
72 Self: SolidFiniteElementBlock<C, F, G, M, N, P>,
73{
74 fn deformation_gradient_rates(
75 &self,
76 nodal_coordinates: &NodalCoordinates,
77 nodal_velocities: &NodalVelocities,
78 ) -> Vec<DeformationGradientRateList<G>> {
79 self.elements()
80 .iter()
81 .zip(self.connectivity())
82 .map(|(element, nodes)| {
83 element.deformation_gradient_rates(
84 &Self::element_coordinates(nodal_coordinates, nodes),
85 &self.element_velocities(nodal_velocities, nodes),
86 )
87 })
88 .collect()
89 }
90 fn element_velocities(
91 &self,
92 velocities: &NodalVelocities,
93 nodes: &[usize; N],
94 ) -> ElementNodalVelocities<N> {
95 nodes.iter().map(|&node| velocities[node].clone()).collect()
96 }
97 fn nodal_forces(
98 &self,
99 nodal_coordinates: &NodalCoordinates,
100 nodal_velocities: &NodalVelocities,
101 ) -> Result<NodalForcesSolid, FiniteElementBlockError> {
102 let mut nodal_forces = NodalForcesSolid::zero(nodal_coordinates.len());
103 match self
104 .elements()
105 .iter()
106 .zip(self.connectivity())
107 .try_for_each(|(element, nodes)| {
108 element
109 .nodal_forces(
110 self.constitutive_model(),
111 &Self::element_coordinates(nodal_coordinates, nodes),
112 &self.element_velocities(nodal_velocities, nodes),
113 )?
114 .iter()
115 .zip(nodes)
116 .for_each(|(nodal_force, &node)| nodal_forces[node] += nodal_force);
117 Ok::<(), FiniteElementError>(())
118 }) {
119 Ok(()) => Ok(nodal_forces),
120 Err(error) => Err(FiniteElementBlockError::Upstream(
121 format!("{error}"),
122 format!("{self:?}"),
123 )),
124 }
125 }
126 fn nodal_stiffnesses(
127 &self,
128 nodal_coordinates: &NodalCoordinates,
129 nodal_velocities: &NodalVelocities,
130 ) -> Result<NodalStiffnessesSolid, FiniteElementBlockError> {
131 let mut nodal_stiffnesses = NodalStiffnessesSolid::zero(nodal_coordinates.len());
132 match self
133 .elements()
134 .iter()
135 .zip(self.connectivity())
136 .try_for_each(|(element, nodes)| {
137 element
138 .nodal_stiffnesses(
139 self.constitutive_model(),
140 &Self::element_coordinates(nodal_coordinates, nodes),
141 &self.element_velocities(nodal_velocities, nodes),
142 )?
143 .iter()
144 .zip(nodes)
145 .for_each(|(object, &node_a)| {
146 object
147 .iter()
148 .zip(nodes)
149 .for_each(|(nodal_stiffness, &node_b)| {
150 nodal_stiffnesses[node_a][node_b] += nodal_stiffness
151 })
152 });
153 Ok::<(), FiniteElementError>(())
154 }) {
155 Ok(()) => Ok(nodal_stiffnesses),
156 Err(error) => Err(FiniteElementBlockError::Upstream(
157 format!("{error}"),
158 format!("{self:?}"),
159 )),
160 }
161 }
162 fn root(
163 &self,
164 equality_constraint: EqualityConstraint,
165 integrator: impl ImplicitDaeFirstOrderRoot<
166 NodalForcesSolid,
167 NodalStiffnessesSolid,
168 NodalVelocities,
169 NodalVelocitiesHistory,
170 >,
171 time: &[Scalar],
172 solver: impl FirstOrderRootFinding<NodalForcesSolid, NodalStiffnessesSolid, NodalCoordinates>,
173 ) -> Result<(Times, NodalCoordinatesHistory, NodalVelocitiesHistory), IntegrationError> {
174 integrator.integrate(
175 |_: Scalar,
176 nodal_coordinates: &NodalCoordinates,
177 nodal_velocities: &NodalVelocities| {
178 Ok(self.nodal_forces(nodal_coordinates, nodal_velocities)?)
179 },
180 |_: Scalar,
181 nodal_coordinates: &NodalCoordinates,
182 nodal_velocities: &NodalVelocities| {
183 Ok(self.nodal_stiffnesses(nodal_coordinates, nodal_velocities)?)
184 },
185 solver,
186 time,
187 self.coordinates().clone().into(),
188 |_: Scalar| equality_constraint.clone(),
189 )
190 }
191}