conspire 0.7.1

The Rust interface to conspire.
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94

#[cfg(test)]
mod test;

use crate::{
    geometry::{Coordinate, Coordinates, mesh::Mesh},
    math::{Scalar, TensorArray},
};
use std::collections::HashMap;

#[derive(Clone, Copy)]
pub enum Weighting {
    Uniform,
    Cotangent,
}

fn edge_key(a: usize, b: usize) -> (usize, usize) {
    if a < b { (a, b) } else { (b, a) }
}

impl<const D: usize> Mesh<D> {
    pub fn laplacian(&self, weighting: Weighting) -> Coordinates<D> {
        self.laplacian_over(self.node_node_connectivity(), weighting)
    }
    pub(crate) fn laplacian_over(
        &self,
        adjacency: &[Vec<usize>],
        weighting: Weighting,
    ) -> Coordinates<D> {
        let coordinates = self.coordinates();
        match weighting {
            Weighting::Uniform => adjacency
                .iter()
                .enumerate()
                .map(|(node_a, nodes)| {
                    if nodes.is_empty() {
                        Coordinate::zero()
                    } else {
                        &coordinates[node_a]
                            - nodes
                                .iter()
                                .map(|&node_b| &coordinates[node_b])
                                .sum::<Coordinate<D>>()
                                / (nodes.len() as Scalar)
                    }
                })
                .collect(),
            Weighting::Cotangent => {
                let weights = self.cotangent_weights();
                adjacency
                    .iter()
                    .enumerate()
                    .map(|(node_a, nodes)| {
                        let mut total = 0.0;
                        let displacement = nodes
                            .iter()
                            .map(|&node_b| {
                                let weight = weights[&edge_key(node_a, node_b)];
                                total += weight;
                                (&coordinates[node_a] - &coordinates[node_b]) * weight
                            })
                            .sum::<Coordinate<D>>();
                        if total == 0.0 {
                            Coordinate::zero()
                        } else {
                            displacement / total
                        }
                    })
                    .collect()
            }
        }
    }
    fn cotangent_weights(&self) -> HashMap<(usize, usize), Scalar> {
        let coordinates = self.coordinates();
        let mut weights = HashMap::new();
        for block in self.iter() {
            if block.number_of_nodes_per_element() == Some(3) {
                for element in block.iter() {
                    let triangle = [element[0], element[1], element[2]];
                    for local in 0..3 {
                        let i = triangle[local];
                        let j = triangle[(local + 1) % 3];
                        let k = triangle[(local + 2) % 3];
                        let u = &coordinates[i] - &coordinates[k];
                        let v = &coordinates[j] - &coordinates[k];
                        let dot = &u * &v;
                        let cross = ((&u * &u) * (&v * &v) - dot * dot).sqrt();
                        *weights.entry(edge_key(i, j)).or_insert(0.0) += dot / cross;
                    }
                }
            }
        }
        weights
    }
}