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