1use crate::{
2 constitutive::solid::elastic_viscoplastic::ElasticViscoplastic,
3 fem::{
4 NodalCoordinates, NodalCoordinatesHistory,
5 block::{
6 Block, FiniteElementBlockError,
7 element::{
8 FiniteElementError, solid::elastic_viscoplastic::ElasticViscoplasticFiniteElement,
9 },
10 solid::{NodalForcesSolid, NodalStiffnessesSolid, SolidFiniteElementBlock},
11 },
12 },
13 math::{
14 Scalar, Tensor, TensorTupleListVec, TensorTupleListVec2D,
15 integrate::{ExplicitDaeFirstOrderRoot, IntegrationError},
16 optimize::{EqualityConstraint, FirstOrderRootFinding},
17 },
18 mechanics::{DeformationGradientPlastic, Times},
19};
20use std::array::from_fn;
21
22pub type ViscoplasticStateVariables<const G: usize, Y> =
23 TensorTupleListVec<DeformationGradientPlastic, Y, G>;
24
25pub type ViscoplasticStateVariablesHistory<const G: usize, Y> =
26 TensorTupleListVec2D<DeformationGradientPlastic, Y, G>;
27
28pub type ElasticViscoplasticBCs = fn(Scalar) -> EqualityConstraint;
29
30pub trait ElasticViscoplasticFiniteElementBlock<
31 C,
32 F,
33 const G: usize,
34 const M: usize,
35 const N: usize,
36 const P: usize,
37 Y,
38> where
39 C: ElasticViscoplastic<Y>,
40 F: ElasticViscoplasticFiniteElement<C, G, M, N, P, Y>,
41 Self: SolidFiniteElementBlock<C, F, G, M, N, P>,
42 Y: Tensor,
43{
44 fn nodal_forces(
45 &self,
46 nodal_coordinates: &NodalCoordinates,
47 state_variables: &ViscoplasticStateVariables<G, Y>,
48 ) -> Result<NodalForcesSolid, FiniteElementBlockError>;
49 fn nodal_stiffnesses(
50 &self,
51 nodal_coordinates: &NodalCoordinates,
52 state_variables: &ViscoplasticStateVariables<G, Y>,
53 ) -> Result<NodalStiffnessesSolid, FiniteElementBlockError>;
54 fn state_variables_evolution(
55 &self,
56 nodal_coordinates: &NodalCoordinates,
57 state_variables: &ViscoplasticStateVariables<G, Y>,
58 ) -> Result<ViscoplasticStateVariables<G, Y>, FiniteElementBlockError>;
59 fn root(
60 &self,
61 integrator: impl ExplicitDaeFirstOrderRoot<
62 NodalForcesSolid,
63 NodalStiffnessesSolid,
64 ViscoplasticStateVariables<G, Y>,
65 NodalCoordinates,
66 ViscoplasticStateVariablesHistory<G, Y>,
67 NodalCoordinatesHistory,
68 >,
69 solver: impl FirstOrderRootFinding<NodalForcesSolid, NodalStiffnessesSolid, NodalCoordinates>,
70 time: &[Scalar],
71 bcs: ElasticViscoplasticBCs,
72 ) -> Result<
73 (
74 Times,
75 NodalCoordinatesHistory,
76 ViscoplasticStateVariablesHistory<G, Y>,
77 ),
78 IntegrationError,
79 >;
80}
81
82impl<C, F, const G: usize, const M: usize, const N: usize, const P: usize, Y>
83 ElasticViscoplasticFiniteElementBlock<C, F, G, M, N, P, Y> for Block<C, F, G, M, N, P>
84where
85 C: ElasticViscoplastic<Y>,
86 F: ElasticViscoplasticFiniteElement<C, G, M, N, P, Y>,
87 Self: SolidFiniteElementBlock<C, F, G, M, N, P>,
88 Y: Tensor,
89{
90 fn nodal_forces(
91 &self,
92 nodal_coordinates: &NodalCoordinates,
93 state_variables: &ViscoplasticStateVariables<G, Y>,
94 ) -> Result<NodalForcesSolid, FiniteElementBlockError> {
95 let mut nodal_forces = NodalForcesSolid::zero(nodal_coordinates.len());
96 match self
97 .elements()
98 .iter()
99 .zip(self.connectivity())
100 .zip(state_variables)
101 .try_for_each(|((element, nodes), state_variables_element)| {
102 element
103 .nodal_forces(
104 self.constitutive_model(),
105 &Self::element_coordinates(nodal_coordinates, nodes),
106 state_variables_element,
107 )?
108 .iter()
109 .zip(nodes)
110 .for_each(|(nodal_force, &node)| nodal_forces[node] += nodal_force);
111 Ok::<(), FiniteElementError>(())
112 }) {
113 Ok(()) => Ok(nodal_forces),
114 Err(error) => Err(FiniteElementBlockError::Upstream(
115 format!("{error}"),
116 format!("{self:?}"),
117 )),
118 }
119 }
120 fn nodal_stiffnesses(
121 &self,
122 nodal_coordinates: &NodalCoordinates,
123 state_variables: &ViscoplasticStateVariables<G, Y>,
124 ) -> Result<NodalStiffnessesSolid, FiniteElementBlockError> {
125 let mut nodal_stiffnesses = NodalStiffnessesSolid::zero(nodal_coordinates.len());
126 match self
127 .elements()
128 .iter()
129 .zip(self.connectivity())
130 .zip(state_variables)
131 .try_for_each(|((element, nodes), state_variables_element)| {
132 element
133 .nodal_stiffnesses(
134 self.constitutive_model(),
135 &Self::element_coordinates(nodal_coordinates, nodes),
136 state_variables_element,
137 )?
138 .iter()
139 .zip(nodes)
140 .for_each(|(object, &node_a)| {
141 object
142 .iter()
143 .zip(nodes)
144 .for_each(|(nodal_stiffness, &node_b)| {
145 nodal_stiffnesses[node_a][node_b] += nodal_stiffness
146 })
147 });
148 Ok::<(), FiniteElementError>(())
149 }) {
150 Ok(()) => Ok(nodal_stiffnesses),
151 Err(error) => Err(FiniteElementBlockError::Upstream(
152 format!("{error}"),
153 format!("{self:?}"),
154 )),
155 }
156 }
157 fn state_variables_evolution(
158 &self,
159 nodal_coordinates: &NodalCoordinates,
160 state_variables: &ViscoplasticStateVariables<G, Y>,
161 ) -> Result<ViscoplasticStateVariables<G, Y>, FiniteElementBlockError> {
162 match self
163 .elements()
164 .iter()
165 .zip(self.connectivity())
166 .zip(state_variables)
167 .map(|((element, nodes), element_state_variables)| {
168 element.state_variables_evolution(
169 self.constitutive_model(),
170 &Self::element_coordinates(nodal_coordinates, nodes),
171 element_state_variables,
172 )
173 })
174 .collect()
175 {
176 Ok(state_variables_evolution) => Ok(state_variables_evolution),
177 Err(error) => Err(FiniteElementBlockError::Upstream(
178 format!("{error}"),
179 format!("{self:?}"),
180 )),
181 }
182 }
183 fn root(
184 &self,
185 integrator: impl ExplicitDaeFirstOrderRoot<
186 NodalForcesSolid,
187 NodalStiffnessesSolid,
188 ViscoplasticStateVariables<G, Y>,
189 NodalCoordinates,
190 ViscoplasticStateVariablesHistory<G, Y>,
191 NodalCoordinatesHistory,
192 >,
193 solver: impl FirstOrderRootFinding<NodalForcesSolid, NodalStiffnessesSolid, NodalCoordinates>,
194 time: &[Scalar],
195 bcs: ElasticViscoplasticBCs,
196 ) -> Result<
197 (
198 Times,
199 NodalCoordinatesHistory,
200 ViscoplasticStateVariablesHistory<G, Y>,
201 ),
202 IntegrationError,
203 > {
204 let (time_history, state_variables_history, _, nodal_coordinates_history) = integrator
205 .integrate(
206 |_: Scalar,
207 state_variables: &ViscoplasticStateVariables<G, Y>,
208 nodal_coordinates: &NodalCoordinates| {
209 Ok(self.state_variables_evolution(nodal_coordinates, state_variables)?)
210 },
211 |_t: Scalar,
212 state_variables: &ViscoplasticStateVariables<G, Y>,
213 nodal_coordinates: &NodalCoordinates| {
214 Ok(self.nodal_forces(nodal_coordinates, state_variables)?)
215 },
216 |_t: Scalar,
217 state_variables: &ViscoplasticStateVariables<G, Y>,
218 nodal_coordinates: &NodalCoordinates| {
219 Ok(self.nodal_stiffnesses(nodal_coordinates, state_variables)?)
220 },
221 solver,
222 time,
223 (
224 self.elements()
225 .iter()
226 .map(|_| from_fn(|_| self.constitutive_model().initial_state()).into())
227 .collect(),
228 self.coordinates().clone().into(),
229 ),
230 bcs,
231 )?;
232 Ok((
233 time_history,
234 nodal_coordinates_history,
235 state_variables_history,
236 ))
237 }
238}