conspire/geometry/mesh/differential/laplace/
mod.rs

1#[cfg(test)]
2mod test;
3
4use crate::{
5    geometry::{
6        Coordinate, Coordinates,
7        mesh::{Connectivity, Mesh},
8    },
9    math::{Scalar, TensorArray},
10};
11use std::collections::HashMap;
12
13#[derive(Clone, Copy)]
14pub enum Weighting {
15    Uniform,
16    Cotangent,
17}
18
19fn edge_key(a: usize, b: usize) -> (usize, usize) {
20    if a < b { (a, b) } else { (b, a) }
21}
22
23impl<const D: usize> Mesh<D> {
24    pub fn laplacian(&self, weighting: Weighting) -> Result<Coordinates<D>, &'static str> {
25        self.laplacian_over(self.node_node_connectivity(), weighting)
26    }
27    pub(crate) fn laplacian_over(
28        &self,
29        adjacency: &[Vec<usize>],
30        weighting: Weighting,
31    ) -> Result<Coordinates<D>, &'static str> {
32        let coordinates = self.coordinates();
33        Ok(match weighting {
34            Weighting::Uniform => adjacency
35                .iter()
36                .enumerate()
37                .map(|(node_a, nodes)| {
38                    if nodes.is_empty() {
39                        Coordinate::zero()
40                    } else {
41                        &coordinates[node_a]
42                            - nodes
43                                .iter()
44                                .map(|&node_b| &coordinates[node_b])
45                                .sum::<Coordinate<D>>()
46                                / (nodes.len() as Scalar)
47                    }
48                })
49                .collect(),
50            Weighting::Cotangent => {
51                if !self
52                    .iter()
53                    .all(|block| matches!(block, Connectivity::Triangular(_)))
54                {
55                    return Err("cotangent weighting requires an all-triangular mesh");
56                }
57                let weights = self.cotangent_weights();
58                adjacency
59                    .iter()
60                    .enumerate()
61                    .map(|(node_a, nodes)| {
62                        let mut total = 0.0;
63                        let displacement = nodes
64                            .iter()
65                            .map(|&node_b| {
66                                let weight = weights[&edge_key(node_a, node_b)];
67                                total += weight;
68                                (&coordinates[node_a] - &coordinates[node_b]) * weight
69                            })
70                            .sum::<Coordinate<D>>();
71                        if total == 0.0 {
72                            Coordinate::zero()
73                        } else {
74                            displacement / total
75                        }
76                    })
77                    .collect()
78            }
79        })
80    }
81    fn cotangent_weights(&self) -> HashMap<(usize, usize), Scalar> {
82        let coordinates = self.coordinates();
83        let mut weights = HashMap::new();
84        for block in self.iter() {
85            if block.number_of_nodes_per_element() == Some(3) {
86                for element in block.iter() {
87                    let triangle = [element[0], element[1], element[2]];
88                    for local in 0..3 {
89                        let i = triangle[local];
90                        let j = triangle[(local + 1) % 3];
91                        let k = triangle[(local + 2) % 3];
92                        let u = &coordinates[i] - &coordinates[k];
93                        let v = &coordinates[j] - &coordinates[k];
94                        let dot = (&u * &v).value();
95                        let cross = ((&u * &u).value() * (&v * &v).value() - dot * dot).sqrt();
96                        *weights.entry(edge_key(i, j)).or_insert(0.0) += dot / cross;
97                    }
98                }
99            }
100        }
101        weights
102    }
103}
conspire/geometry/mesh/differential/laplace/mod.rs

conspire/geometry/mesh/differential/laplace/
mod.rs
