conspire/domain/fem/block/solid/elastic_hyperviscous/
mod.rs

1use crate::{
2    constitutive::solid::elastic_hyperviscous::ElasticHyperviscous,
3    fem::{
4        NodalCoordinates, NodalCoordinatesHistory, NodalVelocities, NodalVelocitiesHistory,
5        block::{
6            Block, FiniteElementBlockError, band,
7            element::solid::elastic_hyperviscous::ElasticHyperviscousFiniteElement,
8            solid::{
9                NodalForcesSolid, NodalStiffnessesSolid,
10                viscoelastic::ViscoelasticFiniteElementBlock,
11            },
12        },
13    },
14    math::{
15        Scalar, Tensor,
16        integrate::{Explicit, IntegrationError},
17        optimize::{EqualityConstraint, OptimizationError, SecondOrderOptimization},
18    },
19    mechanics::Times,
20};
21
22pub trait ElasticHyperviscousFiniteElementBlock<
23    C,
24    F,
25    const G: usize,
26    const M: usize,
27    const N: usize,
28    const P: usize,
29> where
30    C: ElasticHyperviscous,
31    F: ElasticHyperviscousFiniteElement<C, G, M, N, P>,
32    Self: ViscoelasticFiniteElementBlock<C, F, G, M, N, P>,
33{
34    fn viscous_dissipation(
35        &self,
36        nodal_coordinates: &NodalCoordinates,
37        nodal_velocities: &NodalVelocities,
38    ) -> Result<Scalar, FiniteElementBlockError>;
39    fn dissipation_potential(
40        &self,
41        nodal_coordinates: &NodalCoordinates,
42        nodal_velocities: &NodalVelocities,
43    ) -> Result<Scalar, FiniteElementBlockError>;
44    fn minimize(
45        &self,
46        equality_constraint: EqualityConstraint,
47        integrator: impl Explicit<NodalVelocities, NodalVelocitiesHistory>,
48        time: &[Scalar],
49        solver: impl SecondOrderOptimization<
50            Scalar,
51            NodalForcesSolid,
52            NodalStiffnessesSolid,
53            NodalCoordinates,
54        >,
55    ) -> Result<(Times, NodalCoordinatesHistory, NodalVelocitiesHistory), IntegrationError>;
56    #[doc(hidden)]
57    fn minimize_inner(
58        &self,
59        equality_constraint: &EqualityConstraint,
60        nodal_coordinates: &NodalCoordinates,
61        solver: &impl SecondOrderOptimization<
62            Scalar,
63            NodalForcesSolid,
64            NodalStiffnessesSolid,
65            NodalCoordinates,
66        >,
67        initial_guess: &NodalVelocities,
68    ) -> Result<NodalVelocities, OptimizationError>;
69}
70
71impl<C, F, const G: usize, const M: usize, const N: usize, const P: usize>
72    ElasticHyperviscousFiniteElementBlock<C, F, G, M, N, P> for Block<C, F, G, M, N, P>
73where
74    C: ElasticHyperviscous,
75    F: ElasticHyperviscousFiniteElement<C, G, M, N, P>,
76    Self: ViscoelasticFiniteElementBlock<C, F, G, M, N, P>,
77{
78    fn viscous_dissipation(
79        &self,
80        nodal_coordinates: &NodalCoordinates,
81        nodal_velocities: &NodalVelocities,
82    ) -> Result<Scalar, FiniteElementBlockError> {
83        match self
84            .elements()
85            .iter()
86            .zip(self.connectivity())
87            .map(|(element, nodes)| {
88                element.viscous_dissipation(
89                    self.constitutive_model(),
90                    &Self::element_coordinates(nodal_coordinates, nodes),
91                    &self.element_velocities(nodal_velocities, nodes),
92                )
93            })
94            .sum()
95        {
96            Ok(viscous_dissipation) => Ok(viscous_dissipation),
97            Err(error) => Err(FiniteElementBlockError::Upstream(
98                format!("{error}"),
99                format!("{self:?}"),
100            )),
101        }
102    }
103    fn dissipation_potential(
104        &self,
105        nodal_coordinates: &NodalCoordinates,
106        nodal_velocities: &NodalVelocities,
107    ) -> Result<Scalar, FiniteElementBlockError> {
108        match self
109            .elements()
110            .iter()
111            .zip(self.connectivity())
112            .map(|(element, nodes)| {
113                element.dissipation_potential(
114                    self.constitutive_model(),
115                    &Self::element_coordinates(nodal_coordinates, nodes),
116                    &self.element_velocities(nodal_velocities, nodes),
117                )
118            })
119            .sum()
120        {
121            Ok(dissipation_potential) => Ok(dissipation_potential),
122            Err(error) => Err(FiniteElementBlockError::Upstream(
123                format!("{error}"),
124                format!("{self:?}"),
125            )),
126        }
127    }
128    fn minimize(
129        &self,
130        equality_constraint: EqualityConstraint,
131        integrator: impl Explicit<NodalVelocities, NodalVelocitiesHistory>,
132        time: &[Scalar],
133        solver: impl SecondOrderOptimization<
134            Scalar,
135            NodalForcesSolid,
136            NodalStiffnessesSolid,
137            NodalCoordinates,
138        >,
139    ) -> Result<(Times, NodalCoordinatesHistory, NodalVelocitiesHistory), IntegrationError> {
140        let mut solution = NodalVelocities::zero(self.coordinates().len());
141        integrator.integrate(
142            |_: Scalar, nodal_coordinates: &NodalCoordinates| {
143                solution = self.minimize_inner(
144                    &equality_constraint,
145                    nodal_coordinates,
146                    &solver,
147                    &solution,
148                )?;
149                Ok(solution.clone())
150            },
151            time,
152            self.coordinates().clone().into(),
153        )
154    }
155    fn minimize_inner(
156        &self,
157        equality_constraint: &EqualityConstraint,
158        nodal_coordinates: &NodalCoordinates,
159        solver: &impl SecondOrderOptimization<
160            Scalar,
161            NodalForcesSolid,
162            NodalStiffnessesSolid,
163            NodalVelocities,
164        >,
165        initial_guess: &NodalVelocities,
166    ) -> Result<NodalVelocities, OptimizationError> {
167        let banded = band(
168            self.connectivity(),
169            equality_constraint,
170            nodal_coordinates.len(),
171            3,
172        );
173        solver.minimize(
174            |nodal_velocities: &NodalVelocities| {
175                Ok(self.dissipation_potential(nodal_coordinates, nodal_velocities)?)
176            },
177            |nodal_velocities: &NodalVelocities| {
178                Ok(self.nodal_forces(nodal_coordinates, nodal_velocities)?)
179            },
180            |nodal_velocities: &NodalVelocities| {
181                Ok(self.nodal_stiffnesses(nodal_coordinates, nodal_velocities)?)
182            },
183            initial_guess.clone(),
184            equality_constraint.clone(),
185            Some(banded),
186        )
187    }
188}