conspire/domain/fem/solid/elastic/internal_variables/
mod.rs

1use crate::{
2    fem::{
3        ElementModel, ElementModelError, Elements, Model, NodalCoordinates,
4        block::{
5            element::solid::elastic::internal_variables::InternalVariables,
6            finalize_node_neighbors, solver_from_neighbors,
7        },
8        solid::{NodalForcesSolid, NodalStiffnessesSolid},
9    },
10    math::{
11        Scalar, Tensor, TensorVector, Vector,
12        optimize::{
13            EqualityConstraint, FirstOrderRootFinding, FirstOrderRootFindingIncremental,
14            NewtonRaphson, OptimizationError, SolveStrategy,
15        },
16    },
17};
18use std::cell::{Ref, RefCell};
19
20/// The internal variables held at every integration point of every element.
21pub type InternalVariablesField<const G: usize, V> = TensorVector<InternalVariables<G, V>>;
22
23pub trait ElasticIVElements<const G: usize, V, const D: usize>
24where
25    Self: Elements,
26    V: Tensor,
27{
28    /// The internal variables every integration point starts from.
29    fn internal_variables_initial(&self) -> InternalVariablesField<G, V>;
30    /// Steps the internal variables everywhere alongside a decrement of the
31    /// nodal coordinates, rather than solving them where the coordinates now
32    /// are.
33    ///
34    /// This is the second row of the Newton system the solver never assembles,
35    /// each integration point taking its own share of the increment.
36    fn internal_variables_increment(
37        &self,
38        nodal_coordinates: &NodalCoordinates<D>,
39        internal_variables: &InternalVariablesField<G, V>,
40        nodal_decrement: &NodalCoordinates<D>,
41        step: Scalar,
42    ) -> Result<InternalVariablesField<G, V>, ElementModelError>;
43    /// Solves the internal variables everywhere, holding the deformation fixed.
44    ///
45    /// Every integration point is independent of every other, so this is a
46    /// gather of local solves rather than a system.
47    fn internal_variables_root(
48        &self,
49        local_solver: &NewtonRaphson,
50        nodal_coordinates: &NodalCoordinates<D>,
51        internal_variables: &InternalVariablesField<G, V>,
52    ) -> Result<InternalVariablesField<G, V>, ElementModelError>;
53    fn nodal_forces_into(
54        &self,
55        nodal_coordinates: &NodalCoordinates<D>,
56        internal_variables: &InternalVariablesField<G, V>,
57        nodal_forces: &mut NodalForcesSolid<D>,
58    ) -> Result<(), ElementModelError>;
59    fn nodal_forces(
60        &self,
61        nodal_coordinates: &NodalCoordinates<D>,
62        internal_variables: &InternalVariablesField<G, V>,
63    ) -> Result<NodalForcesSolid<D>, ElementModelError> {
64        let mut nodal_forces = NodalForcesSolid::zero(nodal_coordinates.len());
65        self.nodal_forces_into(nodal_coordinates, internal_variables, &mut nodal_forces)?;
66        Ok(nodal_forces)
67    }
68    /// The nodal forces with the residual of the internal variables eliminated
69    /// into them, for when they are carried rather than solved.
70    ///
71    /// ```math
72    /// \mathbf{r}_u - \mathcal{K}_{uv}\mathcal{K}_{vv}^{-1}\mathbf{r}_v
73    /// ```
74    fn nodal_forces_eliminated_into(
75        &self,
76        nodal_coordinates: &NodalCoordinates<D>,
77        internal_variables: &InternalVariablesField<G, V>,
78        nodal_forces: &mut NodalForcesSolid<D>,
79    ) -> Result<(), ElementModelError>;
80    fn nodal_forces_eliminated(
81        &self,
82        nodal_coordinates: &NodalCoordinates<D>,
83        internal_variables: &InternalVariablesField<G, V>,
84    ) -> Result<NodalForcesSolid<D>, ElementModelError> {
85        let mut nodal_forces = NodalForcesSolid::zero(nodal_coordinates.len());
86        self.nodal_forces_eliminated_into(
87            nodal_coordinates,
88            internal_variables,
89            &mut nodal_forces,
90        )?;
91        Ok(nodal_forces)
92    }
93    fn nodal_stiffnesses_into(
94        &self,
95        nodal_coordinates: &NodalCoordinates<D>,
96        internal_variables: &InternalVariablesField<G, V>,
97        nodal_stiffnesses: &mut NodalStiffnessesSolid<D>,
98    ) -> Result<(), ElementModelError>;
99    fn nodal_stiffnesses(
100        &self,
101        nodal_coordinates: &NodalCoordinates<D>,
102        internal_variables: &InternalVariablesField<G, V>,
103    ) -> Result<NodalStiffnessesSolid<D>, ElementModelError> {
104        let mut nodal_stiffnesses = NodalStiffnessesSolid::zero(nodal_coordinates.len());
105        self.nodal_stiffnesses_into(
106            nodal_coordinates,
107            internal_variables,
108            &mut nodal_stiffnesses,
109        )?;
110        Ok(nodal_stiffnesses)
111    }
112}
113
114impl<B, const G: usize, V, const D: usize> ElasticIVElements<G, V, D> for Model<B, D>
115where
116    B: ElasticIVElements<G, V, D>,
117    V: Tensor,
118{
119    fn internal_variables_initial(&self) -> InternalVariablesField<G, V> {
120        self.blocks.internal_variables_initial()
121    }
122    fn internal_variables_increment(
123        &self,
124        nodal_coordinates: &NodalCoordinates<D>,
125        internal_variables: &InternalVariablesField<G, V>,
126        nodal_decrement: &NodalCoordinates<D>,
127        step: Scalar,
128    ) -> Result<InternalVariablesField<G, V>, ElementModelError> {
129        self.blocks.internal_variables_increment(
130            nodal_coordinates,
131            internal_variables,
132            nodal_decrement,
133            step,
134        )
135    }
136    fn internal_variables_root(
137        &self,
138        local_solver: &NewtonRaphson,
139        nodal_coordinates: &NodalCoordinates<D>,
140        internal_variables: &InternalVariablesField<G, V>,
141    ) -> Result<InternalVariablesField<G, V>, ElementModelError> {
142        self.blocks
143            .internal_variables_root(local_solver, nodal_coordinates, internal_variables)
144    }
145    fn nodal_forces_eliminated_into(
146        &self,
147        nodal_coordinates: &NodalCoordinates<D>,
148        internal_variables: &InternalVariablesField<G, V>,
149        nodal_forces: &mut NodalForcesSolid<D>,
150    ) -> Result<(), ElementModelError> {
151        self.blocks.nodal_forces_eliminated_into(
152            nodal_coordinates,
153            internal_variables,
154            nodal_forces,
155        )
156    }
157    fn nodal_forces_into(
158        &self,
159        nodal_coordinates: &NodalCoordinates<D>,
160        internal_variables: &InternalVariablesField<G, V>,
161        nodal_forces: &mut NodalForcesSolid<D>,
162    ) -> Result<(), ElementModelError> {
163        self.blocks
164            .nodal_forces_into(nodal_coordinates, internal_variables, nodal_forces)
165    }
166    fn nodal_stiffnesses_into(
167        &self,
168        nodal_coordinates: &NodalCoordinates<D>,
169        internal_variables: &InternalVariablesField<G, V>,
170        nodal_stiffnesses: &mut NodalStiffnessesSolid<D>,
171    ) -> Result<(), ElementModelError> {
172        self.blocks
173            .nodal_stiffnesses_into(nodal_coordinates, internal_variables, nodal_stiffnesses)
174    }
175}
176
177/// The internal variables as a function of the nodal coordinates, solved
178/// wherever they are asked for.
179///
180/// The residual and the tangent are asked for at the same coordinates, so the
181/// solve is remembered rather than repeated, and what it converged to last is
182/// where the next one starts.
183pub struct SolvedInternalVariables<'a, B, const G: usize, V, const D: usize>
184where
185    V: Tensor,
186{
187    initial: InternalVariablesField<G, V>,
188    local_solver: &'a NewtonRaphson,
189    model: &'a Model<B, D>,
190    solved: RefCell<Option<(NodalCoordinates<D>, InternalVariablesField<G, V>)>>,
191}
192
193impl<'a, B, const G: usize, V, const D: usize> SolvedInternalVariables<'a, B, G, V, D>
194where
195    B: ElasticIVElements<G, V, D>,
196    V: Tensor,
197{
198    pub fn new(
199        model: &'a Model<B, D>,
200        local_solver: &'a NewtonRaphson,
201        initial: InternalVariablesField<G, V>,
202    ) -> Self {
203        Self {
204            initial,
205            local_solver,
206            model,
207            solved: RefCell::new(None),
208        }
209    }
210    pub fn at(
211        &self,
212        nodal_coordinates: &NodalCoordinates<D>,
213    ) -> Result<InternalVariablesField<G, V>, ElementModelError> {
214        if let Some((ref at, ref variables)) = *self.solved.borrow()
215            && at == nodal_coordinates
216        {
217            return Ok(variables.clone());
218        }
219        let warm = match *self.solved.borrow() {
220            Some((_, ref variables)) => variables.clone(),
221            None => self.initial.clone(),
222        };
223        let variables =
224            self.model
225                .internal_variables_root(self.local_solver, nodal_coordinates, &warm)?;
226        *self.solved.borrow_mut() = Some((nodal_coordinates.clone(), variables.clone()));
227        Ok(variables)
228    }
229}
230
231/// The internal variables carried alongside the nodal coordinates, stepped by
232/// their share of whatever increment the solver lends out.
233///
234/// A step is offered before it is taken, so what the increment is measured from
235/// and what the residual is evaluated at come apart while one is being
236/// considered. Increments are always taken from the state last kept, which is
237/// what makes offering several of them in a row harmless.
238pub struct CarriedInternalVariables<'a, B, const G: usize, V, const D: usize>
239where
240    V: Tensor,
241{
242    committed: RefCell<InternalVariablesField<G, V>>,
243    model: &'a Model<B, D>,
244    stepped: RefCell<InternalVariablesField<G, V>>,
245}
246
247impl<'a, B, const G: usize, V, const D: usize> CarriedInternalVariables<'a, B, G, V, D>
248where
249    B: ElasticIVElements<G, V, D>,
250    V: Tensor,
251{
252    pub fn new(model: &'a Model<B, D>, initial: InternalVariablesField<G, V>) -> Self {
253        Self {
254            committed: RefCell::new(initial.clone()),
255            model,
256            stepped: RefCell::new(initial),
257        }
258    }
259    pub fn stepped(&self) -> Ref<'_, InternalVariablesField<G, V>> {
260        self.stepped.borrow()
261    }
262    pub fn step(
263        &self,
264        nodal_coordinates: &NodalCoordinates<D>,
265        nodal_decrement: &Vector,
266        step: Scalar,
267        commit: bool,
268    ) -> Result<(), ElementModelError> {
269        let stepped = self.model.internal_variables_increment(
270            nodal_coordinates,
271            &self.committed.borrow(),
272            &NodalCoordinates::from(nodal_decrement.clone()),
273            step,
274        )?;
275        if commit {
276            *self.committed.borrow_mut() = stepped.clone()
277        }
278        *self.stepped.borrow_mut() = stepped;
279        Ok(())
280    }
281}
282
283/// First-order root-finding for elastic models whose internal variables are
284/// condensed out at every integration point.
285pub trait FirstOrderRootIV<const G: usize, V, const D: usize>
286where
287    V: Tensor,
288{
289    fn root(
290        &self,
291        equality_constraint: EqualityConstraint,
292        solver: impl FirstOrderRootFinding<
293            NodalForcesSolid<D>,
294            NodalStiffnessesSolid<D>,
295            NodalCoordinates<D>,
296        > + FirstOrderRootFindingIncremental<
297            NodalForcesSolid<D>,
298            NodalStiffnessesSolid<D>,
299            NodalCoordinates<D>,
300        >,
301        strategy: SolveStrategy,
302    ) -> Result<NodalCoordinates<D>, OptimizationError>;
303}
304
305impl<B, const G: usize, V, const D: usize> FirstOrderRootIV<G, V, D> for Model<B, D>
306where
307    B: ElasticIVElements<G, V, D>,
308    V: Tensor,
309{
310    fn root(
311        &self,
312        equality_constraint: EqualityConstraint,
313        solver: impl FirstOrderRootFinding<
314            NodalForcesSolid<D>,
315            NodalStiffnessesSolid<D>,
316            NodalCoordinates<D>,
317        > + FirstOrderRootFindingIncremental<
318            NodalForcesSolid<D>,
319            NodalStiffnessesSolid<D>,
320            NodalCoordinates<D>,
321        >,
322        strategy: SolveStrategy,
323    ) -> Result<NodalCoordinates<D>, OptimizationError> {
324        let mut neighbors = vec![Vec::new(); self.coordinates().len()];
325        self.node_neighbors(&mut neighbors);
326        finalize_node_neighbors(&mut neighbors);
327        //
328        // Either way the solver only ever sees the nodal coordinates, so the
329        // sparsity is that of an ordinary mesh. The internal variables are
330        // element-local, so what they contribute to the tangent lands where the
331        // nodes of their own element already meet.
332        //
333        let sparse = solver_from_neighbors(&neighbors, &equality_constraint, D, false);
334        let initial = self.internal_variables_initial();
335        match strategy {
336            //
337            // The internal variables are solved before they are used, so the
338            // residual is one of the nodal coordinates alone, the solver being
339            // free to evaluate wherever it likes.
340            //
341            SolveStrategy::Condensed(ref local_solver) => {
342                let solved = SolvedInternalVariables::new(self, local_solver, initial);
343                solver.root(
344                    |nodal_coordinates: &NodalCoordinates<D>| {
345                        Ok(self.nodal_forces(nodal_coordinates, &solved.at(nodal_coordinates)?)?)
346                    },
347                    |nodal_coordinates: &NodalCoordinates<D>| {
348                        Ok(self
349                            .nodal_stiffnesses(nodal_coordinates, &solved.at(nodal_coordinates)?)?)
350                    },
351                    self.coordinates().clone().into(),
352                    equality_constraint,
353                    Some(sparse),
354                )
355            }
356            //
357            // The internal variables are carried instead, stepped once per
358            // iteration by the increment the solver lends out. They are never
359            // at their own root along the way, so their residual has to be
360            // eliminated into the one the solver does see.
361            //
362            SolveStrategy::Monolithic { elimination: true } => {
363                let carried = CarriedInternalVariables::new(self, initial);
364                solver.root_incremental(
365                    |nodal_coordinates: &NodalCoordinates<D>| {
366                        Ok(self.nodal_forces_eliminated(nodal_coordinates, &carried.stepped())?)
367                    },
368                    |nodal_coordinates: &NodalCoordinates<D>| {
369                        Ok(self.nodal_stiffnesses(nodal_coordinates, &carried.stepped())?)
370                    },
371                    |nodal_coordinates: &NodalCoordinates<D>,
372                     decrement: &Vector,
373                     step: Scalar,
374                     commit: bool| {
375                        Ok(carried.step(nodal_coordinates, decrement, step, commit)?)
376                    },
377                    self.coordinates().clone().into(),
378                    equality_constraint,
379                    Some(sparse),
380                )
381            }
382            SolveStrategy::Monolithic { elimination: false } => unimplemented!(
383                "The internal variables must be unknowns of the solver to be solved with it."
384            ),
385        }
386    }
387}
conspire/domain/fem/solid/elastic/internal_variables/mod.rs

conspire/domain/fem/solid/elastic/internal_variables/
mod.rs
