formoniq 0.2.0

Finite Element Exterior Calculus (FEEC): PDEs formulated with differential forms and solved on simplicial pseudo-Riemannian manifolds of arbitrary dimension, intrinsically and coordinate-free.
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

extern crate nalgebra as na;

use formoniq::operators::{self, ElMatProvider};
use simplicial::linalg::Matrix;
use simplicial::{atlas::refsimp_vol, geometry::metric::simplex::SimplexLengthsSq, Dim};

use approx::assert_relative_eq;

fn check_ref_elmat<F, G, E>(elmat: G, ref_elmat: F)
where
  E: ElMatProvider,
  F: Fn(Dim) -> Option<Matrix>,
  G: Fn(Dim) -> E,
{
  for dim in 1..=10 {
    let Some(expected_elmat) = ref_elmat(dim) else {
      continue;
    };
    let elmat = elmat(dim);

    let refcell = SimplexLengthsSq::standard(dim);
    let computed_elmat = elmat.eval(&refcell.metric());

    assert_relative_eq!(&computed_elmat, &expected_elmat);
  }
}

#[test]
fn laplacian_refcell() {
  check_ref_elmat(|dim| operators::CodifDifElmat::new(dim, 0), ref_laplacian);
}
fn ref_laplacian(dim: Dim) -> Option<Matrix> {
  let ndofs = dim + 1;
  let mut expected_elmat = Matrix::zeros(ndofs, ndofs);
  expected_elmat[(0, 0)] = dim as i32;
  for i in 1..ndofs {
    expected_elmat[(i, 0)] = -1;
    expected_elmat[(0, i)] = -1;
    expected_elmat[(i, i)] = 1;
  }

  Some(expected_elmat.cast::<f64>() * refsimp_vol(dim))
}

#[test]
fn mass_refcell() {
  check_ref_elmat(|dim| operators::HodgeMassElmat::new(dim, 0), ref_mass);
}
fn ref_mass(dim: Dim) -> Option<Matrix> {
  #[rustfmt::skip]
  let mats = [
    na::dmatrix![1.0],
    na::dmatrix![
      1.0/3.0, 1.0/6.0;
      1.0/6.0, 1.0/3.0;
    ],
    na::dmatrix![
      1.0/12.0, 1.0/24.0, 1.0/24.0;
      1.0/24.0, 1.0/12.0, 1.0/24.0;
      1.0/24.0, 1.0/24.0, 1.0/12.0;
    ],
    na::dmatrix![
      1.0/60.0, 1.0/120.0, 1.0/120.0, 1.0/120.0;
      1.0/120.0, 1.0/60.0, 1.0/120.0, 1.0/120.0;
      1.0/120.0, 1.0/120.0, 1.0/60.0, 1.0/120.0;
      1.0/120.0, 1.0/120.0, 1.0/120.0, 1.0/60.0;
    ],
  ];
  mats.get(dim).cloned()
}

#[test]
fn lumped_mass_refcell() {
  check_ref_elmat(|_| operators::ScalarLumpedMassElmat, ref_lumped_mass);
}
fn ref_lumped_mass(dim: Dim) -> Option<Matrix> {
  let nvertices = dim + 1;
  let ndofs = nvertices;
  Some(refsimp_vol(dim) / ndofs as f64 * Matrix::identity(ndofs, ndofs))
}