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

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

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