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

1use crate::{
2    constitutive::solid::elastic::Elastic,
3    fem::{
4        NodalCoordinates,
5        block::{
6            Block, FiniteElementBlockError, FirstOrderRoot, ZerothOrderRoot,
7            element::{FiniteElementError, solid::elastic::ElasticFiniteElement},
8            solid::{NodalForcesSolid, NodalStiffnessesSolid, SolidFiniteElementBlock},
9        },
10    },
11    math::{
12        Tensor,
13        optimize::{
14            EqualityConstraint, FirstOrderRootFinding, OptimizationError, ZerothOrderRootFinding,
15        },
16    },
17};
18
19pub trait ElasticFiniteElementBlock<
20    C,
21    F,
22    const G: usize,
23    const M: usize,
24    const N: usize,
25    const P: usize,
26> where
27    C: Elastic,
28    F: ElasticFiniteElement<C, G, M, N, P>,
29{
30    fn nodal_forces(
31        &self,
32        nodal_coordinates: &NodalCoordinates,
33    ) -> Result<NodalForcesSolid, FiniteElementBlockError>;
34    fn nodal_stiffnesses(
35        &self,
36        nodal_coordinates: &NodalCoordinates,
37    ) -> Result<NodalStiffnessesSolid, FiniteElementBlockError>;
38}
39
40impl<C, F, const G: usize, const M: usize, const N: usize, const P: usize>
41    ElasticFiniteElementBlock<C, F, G, M, N, P> for Block<C, F, G, M, N, P>
42where
43    C: Elastic,
44    F: ElasticFiniteElement<C, G, M, N, P>,
45    Self: SolidFiniteElementBlock<C, F, G, M, N, P>,
46{
47    fn nodal_forces(
48        &self,
49        nodal_coordinates: &NodalCoordinates,
50    ) -> Result<NodalForcesSolid, FiniteElementBlockError> {
51        let mut nodal_forces = NodalForcesSolid::zero(nodal_coordinates.len());
52        match self
53            .elements()
54            .iter()
55            .zip(self.connectivity())
56            .try_for_each(|(element, nodes)| {
57                element
58                    .nodal_forces(
59                        self.constitutive_model(),
60                        &Self::element_coordinates(nodal_coordinates, nodes),
61                    )?
62                    .iter()
63                    .zip(nodes)
64                    .for_each(|(nodal_force, &node)| nodal_forces[node] += nodal_force);
65                Ok::<(), FiniteElementError>(())
66            }) {
67            Ok(()) => Ok(nodal_forces),
68            Err(error) => Err(FiniteElementBlockError::Upstream(
69                format!("{error}"),
70                format!("{self:?}"),
71            )),
72        }
73    }
74    fn nodal_stiffnesses(
75        &self,
76        nodal_coordinates: &NodalCoordinates,
77    ) -> Result<NodalStiffnessesSolid, FiniteElementBlockError> {
78        let mut nodal_stiffnesses = NodalStiffnessesSolid::zero(nodal_coordinates.len());
79        match self
80            .elements()
81            .iter()
82            .zip(self.connectivity())
83            .try_for_each(|(element, nodes)| {
84                element
85                    .nodal_stiffnesses(
86                        self.constitutive_model(),
87                        &Self::element_coordinates(nodal_coordinates, nodes),
88                    )?
89                    .iter()
90                    .zip(nodes)
91                    .for_each(|(object, &node_a)| {
92                        object
93                            .iter()
94                            .zip(nodes)
95                            .for_each(|(nodal_stiffness, &node_b)| {
96                                nodal_stiffnesses[node_a][node_b] += nodal_stiffness
97                            })
98                    });
99                Ok::<(), FiniteElementError>(())
100            }) {
101            Ok(()) => Ok(nodal_stiffnesses),
102            Err(error) => Err(FiniteElementBlockError::Upstream(
103                format!("{error}"),
104                format!("{self:?}"),
105            )),
106        }
107    }
108}
109
110impl<C, F, const G: usize, const M: usize, const N: usize, const P: usize>
111    ZerothOrderRoot<C, F, G, M, N, NodalCoordinates> for Block<C, F, G, M, N, P>
112where
113    C: Elastic,
114    F: ElasticFiniteElement<C, G, M, N, P>,
115{
116    fn root(
117        &self,
118        equality_constraint: EqualityConstraint,
119        solver: impl ZerothOrderRootFinding<NodalCoordinates>,
120    ) -> Result<NodalCoordinates, OptimizationError> {
121        solver.root(
122            |nodal_coordinates: &NodalCoordinates| Ok(self.nodal_forces(nodal_coordinates)?),
123            self.coordinates().clone().into(),
124            equality_constraint,
125        )
126    }
127}
128
129impl<C, F, const G: usize, const M: usize, const N: usize, const P: usize>
130    FirstOrderRoot<C, F, G, M, N, NodalForcesSolid, NodalStiffnessesSolid, NodalCoordinates>
131    for Block<C, F, G, M, N, P>
132where
133    C: Elastic,
134    F: ElasticFiniteElement<C, G, M, N, P>,
135{
136    fn root(
137        &self,
138        equality_constraint: EqualityConstraint,
139        solver: impl FirstOrderRootFinding<NodalForcesSolid, NodalStiffnessesSolid, NodalCoordinates>,
140    ) -> Result<NodalCoordinates, OptimizationError> {
141        solver.root(
142            |nodal_coordinates: &NodalCoordinates| Ok(self.nodal_forces(nodal_coordinates)?),
143            |nodal_coordinates: &NodalCoordinates| Ok(self.nodal_stiffnesses(nodal_coordinates)?),
144            self.coordinates().clone().into(),
145            equality_constraint,
146        )
147    }
148}