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