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

1use crate::{
2    fem::{
3        ElementModel, ElementModelError, Elements, Model, NodalCoordinates,
4        block::{finalize_node_neighbors, solver_from_neighbors},
5        solid::{
6            NodalForcesSolid, NodalStiffnessesSolid,
7            elastic::internal_variables::{
8                CarriedInternalVariables, ElasticIVElements, InternalVariablesField,
9                SolvedInternalVariables,
10            },
11        },
12    },
13    math::{
14        Quantity, Scalar, Tensor, Vector,
15        optimize::{
16            EqualityConstraint, OptimizationError, SecondOrderOptimization,
17            SecondOrderOptimizationIncremental, SolveStrategy,
18        },
19    },
20    units::Energy,
21};
22
23pub trait HyperelasticIVElements<const G: usize, V, const D: usize>
24where
25    Self: ElasticIVElements<G, V, D>,
26    V: Tensor,
27{
28    fn helmholtz_free_energy(
29        &self,
30        nodal_coordinates: &NodalCoordinates<D>,
31        internal_variables: &InternalVariablesField<G, V>,
32    ) -> Result<Quantity<Energy>, ElementModelError>;
33}
34
35impl<B, const G: usize, V, const D: usize> HyperelasticIVElements<G, V, D> for Model<B, D>
36where
37    B: HyperelasticIVElements<G, V, D>,
38    V: Tensor,
39{
40    fn helmholtz_free_energy(
41        &self,
42        nodal_coordinates: &NodalCoordinates<D>,
43        internal_variables: &InternalVariablesField<G, V>,
44    ) -> Result<Quantity<Energy>, ElementModelError> {
45        self.blocks
46            .helmholtz_free_energy(nodal_coordinates, internal_variables)
47    }
48}
49
50/// Second-order minimization for hyperelastic models whose internal variables
51/// are condensed out at every integration point.
52pub trait SecondOrderMinimizeIV<const G: usize, V, const D: usize>
53where
54    V: Tensor,
55{
56    fn minimize(
57        &self,
58        equality_constraint: EqualityConstraint,
59        solver: impl SecondOrderOptimization<
60            Quantity<Energy>,
61            NodalForcesSolid<D>,
62            NodalStiffnessesSolid<D>,
63            NodalCoordinates<D>,
64        > + SecondOrderOptimizationIncremental<
65            Quantity<Energy>,
66            NodalForcesSolid<D>,
67            NodalStiffnessesSolid<D>,
68            NodalCoordinates<D>,
69        >,
70        strategy: SolveStrategy,
71    ) -> Result<NodalCoordinates<D>, OptimizationError>;
72}
73
74impl<B, const G: usize, V, const D: usize> SecondOrderMinimizeIV<G, V, D> for Model<B, D>
75where
76    B: HyperelasticIVElements<G, V, D>,
77    V: Tensor,
78{
79    fn minimize(
80        &self,
81        equality_constraint: EqualityConstraint,
82        solver: impl SecondOrderOptimization<
83            Quantity<Energy>,
84            NodalForcesSolid<D>,
85            NodalStiffnessesSolid<D>,
86            NodalCoordinates<D>,
87        > + SecondOrderOptimizationIncremental<
88            Quantity<Energy>,
89            NodalForcesSolid<D>,
90            NodalStiffnessesSolid<D>,
91            NodalCoordinates<D>,
92        >,
93        strategy: SolveStrategy,
94    ) -> Result<NodalCoordinates<D>, OptimizationError> {
95        let initial = self.internal_variables_initial();
96        let mut neighbors = vec![Vec::new(); self.coordinates().len()];
97        self.node_neighbors(&mut neighbors);
98        finalize_node_neighbors(&mut neighbors);
99        //
100        // Condensation preserves symmetry: the local block and the couplings
101        // are those of one energy Hessian, so the Schur complement is symmetric
102        // wherever the unreduced tangent was.
103        //
104        let sparse = solver_from_neighbors(&neighbors, &equality_constraint, D, true);
105        match strategy {
106            //
107            // The internal variables are solved before they are used, so the
108            // energy is one of the nodal coordinates alone, the solver being
109            // free to evaluate wherever it likes.
110            //
111            SolveStrategy::Condensed(ref local_solver) => {
112                let solved = SolvedInternalVariables::new(self, local_solver, initial);
113                solver.minimize(
114                    |nodal_coordinates: &NodalCoordinates<D>| {
115                        Ok(self.helmholtz_free_energy(
116                            nodal_coordinates,
117                            &solved.at(nodal_coordinates)?,
118                        )?)
119                    },
120                    |nodal_coordinates: &NodalCoordinates<D>| {
121                        Ok(self.nodal_forces(nodal_coordinates, &solved.at(nodal_coordinates)?)?)
122                    },
123                    |nodal_coordinates: &NodalCoordinates<D>| {
124                        Ok(self
125                            .nodal_stiffnesses(nodal_coordinates, &solved.at(nodal_coordinates)?)?)
126                    },
127                    self.coordinates().clone().into(),
128                    equality_constraint,
129                    Some(sparse),
130                )
131            }
132            //
133            // The internal variables are carried instead, stepped by their
134            // share of whatever step the nodal coordinates take. What the line
135            // search weighs is the energy of the whole state, so they are moved
136            // to where a step would put them before the energy there is asked
137            // for, and left there only if that step is the one taken.
138            //
139            SolveStrategy::Monolithic { elimination: true } => {
140                let carried = CarriedInternalVariables::new(self, initial);
141                solver.minimize_incremental(
142                    |nodal_coordinates: &NodalCoordinates<D>| {
143                        Ok(self.helmholtz_free_energy(nodal_coordinates, &carried.stepped())?)
144                    },
145                    |nodal_coordinates: &NodalCoordinates<D>| {
146                        Ok(self.nodal_forces_eliminated(nodal_coordinates, &carried.stepped())?)
147                    },
148                    |nodal_coordinates: &NodalCoordinates<D>| {
149                        Ok(self.nodal_stiffnesses(nodal_coordinates, &carried.stepped())?)
150                    },
151                    |nodal_coordinates: &NodalCoordinates<D>,
152                     decrement: &Vector,
153                     step: Scalar,
154                     commit: bool| {
155                        Ok(carried.step(nodal_coordinates, decrement, step, commit)?)
156                    },
157                    self.coordinates().clone().into(),
158                    equality_constraint,
159                    Some(sparse),
160                )
161            }
162            SolveStrategy::Monolithic { elimination: false } => unimplemented!(
163                "The internal variables must be unknowns of the solver to be solved with it."
164            ),
165        }
166    }
167}
conspire/domain/fem/solid/hyperelastic/internal_variables/mod.rs

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