Skip to main content

conspire/domain/solid/viscoelastic/
mod.rs

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