conspire/domain/fem/solid/elastic_hyperviscous/
mod.rs1use 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}