Skip to main content

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

1use crate::{
2    fem::{
3        Blocks, ElementModel, ElementModelError, Elements, FirstOrderRoot, Model, NodalCoordinates,
4        ZerothOrderRoot,
5        block::{finalize_node_neighbors, solver_from_neighbors},
6        solid::{NodalForcesSolid, NodalStiffnessesSolid},
7    },
8    math::{
9        Tensor,
10        optimize::{
11            EqualityConstraint, FirstOrderRootFinding, OptimizationError, ZerothOrderRootFinding,
12        },
13    },
14};
15
16pub trait ElasticElements<const D: usize>
17where
18    Self: Elements,
19{
20    fn nodal_forces_into(
21        &self,
22        nodal_coordinates: &NodalCoordinates<D>,
23        nodal_forces: &mut NodalForcesSolid<D>,
24    ) -> Result<(), ElementModelError>;
25    fn nodal_forces(
26        &self,
27        nodal_coordinates: &NodalCoordinates<D>,
28    ) -> Result<NodalForcesSolid<D>, ElementModelError> {
29        let mut nodal_forces = NodalForcesSolid::zero(nodal_coordinates.len());
30        self.nodal_forces_into(nodal_coordinates, &mut nodal_forces)?;
31        Ok(nodal_forces)
32    }
33    fn nodal_stiffnesses_into(
34        &self,
35        nodal_coordinates: &NodalCoordinates<D>,
36        nodal_stiffnesses: &mut NodalStiffnessesSolid<D>,
37    ) -> Result<(), ElementModelError>;
38    fn nodal_stiffnesses(
39        &self,
40        nodal_coordinates: &NodalCoordinates<D>,
41    ) -> Result<NodalStiffnessesSolid<D>, ElementModelError> {
42        let mut nodal_stiffnesses = NodalStiffnessesSolid::zero(nodal_coordinates.len());
43        self.nodal_stiffnesses_into(nodal_coordinates, &mut nodal_stiffnesses)?;
44        Ok(nodal_stiffnesses)
45    }
46}
47
48impl<B, const D: usize> ElasticElements<D> for Model<B, D>
49where
50    B: ElasticElements<D>,
51{
52    fn nodal_forces_into(
53        &self,
54        nodal_coordinates: &NodalCoordinates<D>,
55        nodal_forces: &mut NodalForcesSolid<D>,
56    ) -> Result<(), ElementModelError> {
57        self.blocks
58            .nodal_forces_into(nodal_coordinates, nodal_forces)
59    }
60    fn nodal_stiffnesses_into(
61        &self,
62        nodal_coordinates: &NodalCoordinates<D>,
63        nodal_stiffnesses: &mut NodalStiffnessesSolid<D>,
64    ) -> Result<(), ElementModelError> {
65        self.blocks
66            .nodal_stiffnesses_into(nodal_coordinates, nodal_stiffnesses)
67    }
68}
69
70impl<B1, B2, const D: usize> ElasticElements<D> for Blocks<B1, B2>
71where
72    B1: ElasticElements<D>,
73    B2: ElasticElements<D>,
74{
75    fn nodal_forces_into(
76        &self,
77        nodal_coordinates: &NodalCoordinates<D>,
78        nodal_forces: &mut NodalForcesSolid<D>,
79    ) -> Result<(), ElementModelError> {
80        self.0.nodal_forces_into(nodal_coordinates, nodal_forces)?;
81        self.1.nodal_forces_into(nodal_coordinates, nodal_forces)
82    }
83    fn nodal_stiffnesses_into(
84        &self,
85        nodal_coordinates: &NodalCoordinates<D>,
86        nodal_stiffnesses: &mut NodalStiffnessesSolid<D>,
87    ) -> Result<(), ElementModelError> {
88        self.0
89            .nodal_stiffnesses_into(nodal_coordinates, nodal_stiffnesses)?;
90        self.1
91            .nodal_stiffnesses_into(nodal_coordinates, nodal_stiffnesses)
92    }
93}
94
95impl<B, const D: usize> ZerothOrderRoot<NodalCoordinates<D>> for Model<B, D>
96where
97    B: ElasticElements<D>,
98{
99    fn root(
100        &self,
101        equality_constraint: EqualityConstraint,
102        solver: impl ZerothOrderRootFinding<NodalCoordinates<D>>,
103    ) -> Result<NodalCoordinates<D>, OptimizationError> {
104        solver.root(
105            |nodal_coordinates: &NodalCoordinates<D>| Ok(self.nodal_forces(nodal_coordinates)?),
106            self.coordinates().clone().into(),
107            equality_constraint,
108        )
109    }
110}
111
112impl<B, const D: usize>
113    FirstOrderRoot<NodalForcesSolid<D>, NodalStiffnessesSolid<D>, NodalCoordinates<D>>
114    for Model<B, D>
115where
116    B: ElasticElements<D>,
117{
118    fn root(
119        &self,
120        equality_constraint: EqualityConstraint,
121        solver: impl FirstOrderRootFinding<
122            NodalForcesSolid<D>,
123            NodalStiffnessesSolid<D>,
124            NodalCoordinates<D>,
125        >,
126    ) -> Result<NodalCoordinates<D>, OptimizationError> {
127        let mut neighbors = vec![Vec::new(); self.coordinates().len()];
128        self.node_neighbors(&mut neighbors);
129        finalize_node_neighbors(&mut neighbors);
130        let sparse = solver_from_neighbors(&neighbors, &equality_constraint, D, false);
131        solver.root(
132            |nodal_coordinates: &NodalCoordinates<D>| Ok(self.nodal_forces(nodal_coordinates)?),
133            |nodal_coordinates: &NodalCoordinates<D>| {
134                Ok(self.nodal_stiffnesses(nodal_coordinates)?)
135            },
136            self.coordinates().clone().into(),
137            equality_constraint,
138            Some(sparse),
139        )
140    }
141}