Skip to main content

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

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