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::{ImplicitDaeSecondOrderMinimize, IntegrationError},
17        optimize::{EqualityConstraint, 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 ImplicitDaeSecondOrderMinimize<
48            Scalar,
49            NodalForcesSolid,
50            NodalStiffnessesSolid,
51            NodalVelocities,
52            NodalVelocitiesHistory,
53        >,
54        time: &[Scalar],
55        solver: impl SecondOrderOptimization<
56            Scalar,
57            NodalForcesSolid,
58            NodalStiffnessesSolid,
59            NodalCoordinates,
60        >,
61    ) -> Result<(Times, NodalCoordinatesHistory, NodalVelocitiesHistory), IntegrationError>;
62}
63
64impl<C, F, const G: usize, const M: usize, const N: usize, const P: usize>
65    ElasticHyperviscousFiniteElementBlock<C, F, G, M, N, P> for Block<C, F, G, M, N, P>
66where
67    C: ElasticHyperviscous,
68    F: ElasticHyperviscousFiniteElement<C, G, M, N, P>,
69    Self: ViscoelasticFiniteElementBlock<C, F, G, M, N, P>,
70{
71    fn viscous_dissipation(
72        &self,
73        nodal_coordinates: &NodalCoordinates,
74        nodal_velocities: &NodalVelocities,
75    ) -> Result<Scalar, FiniteElementBlockError> {
76        match self
77            .elements()
78            .iter()
79            .zip(self.connectivity())
80            .map(|(element, nodes)| {
81                element.viscous_dissipation(
82                    self.constitutive_model(),
83                    &Self::element_coordinates(nodal_coordinates, nodes),
84                    &self.element_velocities(nodal_velocities, nodes),
85                )
86            })
87            .sum()
88        {
89            Ok(viscous_dissipation) => Ok(viscous_dissipation),
90            Err(error) => Err(FiniteElementBlockError::Upstream(
91                format!("{error}"),
92                format!("{self:?}"),
93            )),
94        }
95    }
96    fn dissipation_potential(
97        &self,
98        nodal_coordinates: &NodalCoordinates,
99        nodal_velocities: &NodalVelocities,
100    ) -> Result<Scalar, FiniteElementBlockError> {
101        match self
102            .elements()
103            .iter()
104            .zip(self.connectivity())
105            .map(|(element, nodes)| {
106                element.dissipation_potential(
107                    self.constitutive_model(),
108                    &Self::element_coordinates(nodal_coordinates, nodes),
109                    &self.element_velocities(nodal_velocities, nodes),
110                )
111            })
112            .sum()
113        {
114            Ok(dissipation_potential) => Ok(dissipation_potential),
115            Err(error) => Err(FiniteElementBlockError::Upstream(
116                format!("{error}"),
117                format!("{self:?}"),
118            )),
119        }
120    }
121    fn minimize(
122        &self,
123        equality_constraint: EqualityConstraint,
124        integrator: impl ImplicitDaeSecondOrderMinimize<
125            Scalar,
126            NodalForcesSolid,
127            NodalStiffnessesSolid,
128            NodalVelocities,
129            NodalVelocitiesHistory,
130        >,
131        time: &[Scalar],
132        solver: impl SecondOrderOptimization<
133            Scalar,
134            NodalForcesSolid,
135            NodalStiffnessesSolid,
136            NodalCoordinates,
137        >,
138    ) -> Result<(Times, NodalCoordinatesHistory, NodalVelocitiesHistory), IntegrationError> {
139        let banded = band(
140            self.connectivity(),
141            &equality_constraint,
142            self.coordinates().len(),
143            3,
144        );
145        integrator.integrate(
146            |_: Scalar,
147             nodal_coordinates: &NodalCoordinates,
148             nodal_velocities: &NodalVelocities| {
149                Ok(self.dissipation_potential(nodal_coordinates, nodal_velocities)?)
150            },
151            |_: Scalar,
152             nodal_coordinates: &NodalCoordinates,
153             nodal_velocities: &NodalVelocities| {
154                Ok(self.nodal_forces(nodal_coordinates, nodal_velocities)?)
155            },
156            |_: Scalar,
157             nodal_coordinates: &NodalCoordinates,
158             nodal_velocities: &NodalVelocities| {
159                Ok(self.nodal_stiffnesses(nodal_coordinates, nodal_velocities)?)
160            },
161            solver,
162            time,
163            self.coordinates().clone().into(),
164            |_: Scalar| equality_constraint.clone(),
165            Some(banded),
166        )
167    }
168}