Skip to main content

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

1use crate::{
2    fem::{
3        Blocks, ElementModel, ElementModelError, Elements, Model, NodalCoordinates,
4        NodalCoordinatesHistory, NodalVelocities, NodalVelocitiesHistory,
5        block::{finalize_node_neighbors, solver_from_neighbors},
6        solid::{NodalDampingsSolid, NodalForcesSolid, viscoelastic::ViscoelasticElements},
7    },
8    math::{
9        Quantity, Tensor,
10        integrate::{ImplicitDaeSecondOrderMinimize, IntegrationError},
11        optimize::{EqualityConstraint, SecondOrderOptimization},
12    },
13    mechanics::Times,
14    units::{Power, Time},
15};
16
17pub trait ElasticHyperviscousElements<const D: usize>
18where
19    Self: ViscoelasticElements<D>,
20{
21    fn viscous_dissipation(
22        &self,
23        nodal_coordinates: &NodalCoordinates<D>,
24        nodal_velocities: &NodalVelocities<D>,
25    ) -> Result<Quantity<Power>, ElementModelError>;
26    fn dissipation_potential(
27        &self,
28        nodal_coordinates: &NodalCoordinates<D>,
29        nodal_velocities: &NodalVelocities<D>,
30    ) -> Result<Quantity<Power>, ElementModelError>;
31}
32
33impl<B, const D: usize> ElasticHyperviscousElements<D> for Model<B, D>
34where
35    B: ElasticHyperviscousElements<D>,
36{
37    fn viscous_dissipation(
38        &self,
39        nodal_coordinates: &NodalCoordinates<D>,
40        nodal_velocities: &NodalVelocities<D>,
41    ) -> Result<Quantity<Power>, ElementModelError> {
42        self.blocks
43            .viscous_dissipation(nodal_coordinates, nodal_velocities)
44    }
45    fn dissipation_potential(
46        &self,
47        nodal_coordinates: &NodalCoordinates<D>,
48        nodal_velocities: &NodalVelocities<D>,
49    ) -> Result<Quantity<Power>, ElementModelError> {
50        self.blocks
51            .dissipation_potential(nodal_coordinates, nodal_velocities)
52    }
53}
54
55impl<B1, B2, const D: usize> ElasticHyperviscousElements<D> for Blocks<B1, B2>
56where
57    B1: ElasticHyperviscousElements<D>,
58    B2: ElasticHyperviscousElements<D>,
59{
60    fn viscous_dissipation(
61        &self,
62        nodal_coordinates: &NodalCoordinates<D>,
63        nodal_velocities: &NodalVelocities<D>,
64    ) -> Result<Quantity<Power>, ElementModelError> {
65        Ok(self
66            .0
67            .viscous_dissipation(nodal_coordinates, nodal_velocities)?
68            + self
69                .1
70                .viscous_dissipation(nodal_coordinates, nodal_velocities)?)
71    }
72    fn dissipation_potential(
73        &self,
74        nodal_coordinates: &NodalCoordinates<D>,
75        nodal_velocities: &NodalVelocities<D>,
76    ) -> Result<Quantity<Power>, ElementModelError> {
77        Ok(self
78            .0
79            .dissipation_potential(nodal_coordinates, nodal_velocities)?
80            + self
81                .1
82                .dissipation_potential(nodal_coordinates, nodal_velocities)?)
83    }
84}
85
86pub trait SecondOrderMinimize<const D: usize> {
87    fn minimize(
88        &self,
89        equality_constraint: EqualityConstraint,
90        integrator: impl ImplicitDaeSecondOrderMinimize<
91            Quantity<Power>,
92            NodalForcesSolid<D>,
93            NodalDampingsSolid<D>,
94            NodalCoordinates<D>,
95            NodalCoordinatesHistory<D>,
96            NodalVelocitiesHistory<D>,
97        >,
98        time: &[Quantity<Time>],
99        solver: impl SecondOrderOptimization<
100            Quantity<Power>,
101            NodalForcesSolid<D>,
102            NodalDampingsSolid<D>,
103            NodalVelocities<D>,
104        >,
105    ) -> Result<(Times, NodalCoordinatesHistory<D>, NodalVelocitiesHistory<D>), IntegrationError>;
106}
107
108impl<B, const D: usize> SecondOrderMinimize<D> for Model<B, D>
109where
110    B: ElasticHyperviscousElements<D>,
111{
112    fn minimize(
113        &self,
114        equality_constraint: EqualityConstraint,
115        integrator: impl ImplicitDaeSecondOrderMinimize<
116            Quantity<Power>,
117            NodalForcesSolid<D>,
118            NodalDampingsSolid<D>,
119            NodalCoordinates<D>,
120            NodalCoordinatesHistory<D>,
121            NodalVelocitiesHistory<D>,
122        >,
123        time: &[Quantity<Time>],
124        solver: impl SecondOrderOptimization<
125            Quantity<Power>,
126            NodalForcesSolid<D>,
127            NodalDampingsSolid<D>,
128            NodalVelocities<D>,
129        >,
130    ) -> Result<(Times, NodalCoordinatesHistory<D>, NodalVelocitiesHistory<D>), IntegrationError>
131    {
132        let mut neighbors = vec![Vec::new(); self.coordinates().len()];
133        self.node_neighbors(&mut neighbors);
134        finalize_node_neighbors(&mut neighbors);
135        let sparse = solver_from_neighbors(&neighbors, &equality_constraint, D, true);
136        integrator.integrate(
137            |_: Quantity<Time>,
138             nodal_coordinates: &NodalCoordinates<D>,
139             nodal_velocities: &NodalVelocities<D>| {
140                Ok(self.dissipation_potential(nodal_coordinates, nodal_velocities)?)
141            },
142            |_: Quantity<Time>,
143             nodal_coordinates: &NodalCoordinates<D>,
144             nodal_velocities: &NodalVelocities<D>| {
145                Ok(self.nodal_forces(nodal_coordinates, nodal_velocities)?)
146            },
147            |_: Quantity<Time>,
148             nodal_coordinates: &NodalCoordinates<D>,
149             nodal_velocities: &NodalVelocities<D>| {
150                Ok(self.nodal_stiffnesses(nodal_coordinates, nodal_velocities)?)
151            },
152            solver,
153            time,
154            self.coordinates().clone().into(),
155            |_: Quantity<Time>| equality_constraint.clone(),
156            Some(sparse),
157        )
158    }
159}