conspire/domain/cbm/block/node/
mod.rs

1pub mod solid;
2#[cfg(test)]
3mod test;
4mod tetrahedron;
5
6use crate::{
7    domain::NodalReferenceCoordinates,
8    geometry::mesh::PrimitiveConnectivity,
9    math::{
10        CrossProduct, Quantity, Reference, Scalar, Tensor, TensorRank1, TensorRank1List,
11        TensorRank1Vec,
12    },
13    units::{ReciprocalLength, UnitMul, Volume},
14};
15use std::collections::HashMap;
16use tetrahedron::{ElementNodalReferenceCoordinates, GradientVectors, Tetrahedron};
17
18#[derive(Clone, Copy, Debug, Default)]
19pub enum Weighting {
20    SolidAngle,
21    #[default]
22    Uniform,
23}
24
25impl Weighting {
26    fn weights(&self, coordinates: &ElementNodalReferenceCoordinates) -> [Scalar; 4] {
27        match self {
28            Self::Uniform => [0.25; 4],
29            Self::SolidAngle => solid_angle_weights(coordinates),
30        }
31    }
32}
33
34fn solid_angle_weights(coordinates: &ElementNodalReferenceCoordinates) -> [Scalar; 4] {
35    let solid_angle = |apex: usize, others: [usize; 3]| -> Scalar {
36        let a = &coordinates[others[0]] - &coordinates[apex];
37        let b = &coordinates[others[1]] - &coordinates[apex];
38        let c = &coordinates[others[2]] - &coordinates[apex];
39        let numerator = (&a * &b.cross(&c)).value().abs();
40        let (norm_a, norm_b, norm_c) = (a.norm().value(), b.norm().value(), c.norm().value());
41        let denominator = norm_a * norm_b * norm_c
42            + (&a * &b).value() * norm_c
43            + (&a * &c).value() * norm_b
44            + (&b * &c).value() * norm_a;
45        2.0 * numerator.atan2(denominator)
46    };
47    let angles = [
48        solid_angle(0, [1, 2, 3]),
49        solid_angle(1, [0, 2, 3]),
50        solid_angle(2, [0, 1, 3]),
51        solid_angle(3, [0, 1, 2]),
52    ];
53    let sum: Scalar = angles.iter().sum();
54    angles.map(|angle| angle / sum)
55}
56
57pub(crate) type BondGradientVector = TensorRank1<3, Reference, ReciprocalLength>;
58type UnnormalizedBondGradientVector =
59    TensorRank1<3, Reference, <ReciprocalLength as UnitMul<Volume>>::Output>;
60
61pub(crate) struct Node {
62    volume: Quantity<Volume>,
63    neighbors: Vec<usize>,
64    gradient_vectors: Vec<BondGradientVector>,
65}
66
67impl Node {
68    pub(crate) fn neighbors(&self) -> &[usize] {
69        &self.neighbors
70    }
71    pub(crate) fn gradient_vectors(&self) -> &[BondGradientVector] {
72        &self.gradient_vectors
73    }
74    fn element_coordinates<const D: usize, I, U>(
75        coordinates: &TensorRank1Vec<D, I, U>,
76        nodes: &[usize; 4],
77    ) -> TensorRank1List<D, I, 4, U> {
78        nodes
79            .iter()
80            .map(|&node| coordinates[node].clone())
81            .collect()
82    }
83    fn accumulate_bonds(
84        connectivity: &PrimitiveConnectivity<3, 4>,
85        elements: &[Tetrahedron],
86        weightings: &[[Scalar; 4]],
87    ) -> HashMap<(usize, usize), UnnormalizedBondGradientVector> {
88        let mut bonds: HashMap<(usize, usize), UnnormalizedBondGradientVector> = HashMap::new();
89        connectivity
90            .iter()
91            .zip(elements.iter())
92            .zip(weightings.iter())
93            .for_each(|((nodes, element), weights)| {
94                let gradient_vectors: &GradientVectors = element.gradient_vectors();
95                nodes.iter().enumerate().for_each(|(index_a, &node_a)| {
96                    let weight = element.volume() * weights[index_a];
97                    nodes.iter().zip(gradient_vectors.iter()).for_each(
98                        |(&node_b, gradient_vector_b)| {
99                            if node_a != node_b {
100                                let contribution = gradient_vector_b * weight;
101                                bonds
102                                    .entry((node_a, node_b))
103                                    .and_modify(|gradient_vector| *gradient_vector += &contribution)
104                                    .or_insert(contribution);
105                            }
106                        },
107                    )
108                })
109            });
110        bonds
111    }
112    pub(crate) fn vec_from(
113        connectivity: &PrimitiveConnectivity<3, 4>,
114        reference_coordinates: &NodalReferenceCoordinates<3>,
115        weighting: Weighting,
116    ) -> Vec<Self> {
117        let (elements, weightings): (Vec<Tetrahedron>, Vec<[Scalar; 4]>) = connectivity
118            .iter()
119            .map(|nodes| {
120                let coordinates: ElementNodalReferenceCoordinates =
121                    Self::element_coordinates(reference_coordinates, nodes);
122                let weights = weighting.weights(&coordinates);
123                (Tetrahedron::from(coordinates), weights)
124            })
125            .unzip();
126        let mut volumes = vec![Quantity::<Volume>::new(0.0); reference_coordinates.len()];
127        connectivity
128            .iter()
129            .zip(elements.iter())
130            .zip(weightings.iter())
131            .for_each(|((nodes, element), weights)| {
132                nodes.iter().enumerate().for_each(|(index, &node)| {
133                    volumes[node] += &(element.volume() * weights[index])
134                })
135            });
136        let bonds = Self::accumulate_bonds(connectivity, &elements, &weightings);
137        let mut unnormalized: Vec<Vec<(usize, UnnormalizedBondGradientVector)>> =
138            vec![Vec::new(); volumes.len()];
139        bonds
140            .into_iter()
141            .for_each(|((node_a, node_b), bond)| unnormalized[node_a].push((node_b, bond)));
142        unnormalized
143            .iter_mut()
144            .enumerate()
145            .for_each(|(node, bonds)| {
146                let self_gradient_vector = -bonds
147                    .iter()
148                    .map(|(_, gradient_vector)| gradient_vector)
149                    .sum::<UnnormalizedBondGradientVector>();
150                bonds.push((node, self_gradient_vector));
151            });
152        unnormalized
153            .into_iter()
154            .zip(volumes)
155            .map(|(bonds, volume)| {
156                let (neighbors, gradient_vectors) = bonds
157                    .into_iter()
158                    .map(|(node, bond)| (node, bond / volume))
159                    .unzip();
160                Node {
161                    volume,
162                    neighbors,
163                    gradient_vectors,
164                }
165            })
166            .collect()
167    }
168}
conspire/domain/cbm/block/node/mod.rs

conspire/domain/cbm/block/node/
mod.rs
