fenris 0.0.21

A library for advanced finite element computations in Rust
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
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133

use fenris::model::FiniteElementInterpolator;
use nalgebra::{DVector, Vector2};

use proptest::prelude::*;

use matrixcompare::assert_scalar_eq;

#[test]
fn interpolate_into() {
    let values = vec![1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0];
    let supported_nodes = vec![
        0, 1, // Interpolation point 1
        0, 3, 4, // Interpolation point 2
        1, 2, 4, // Interpolation point 3
    ]; // Interpolation point 4
    let supported_node_offsets = vec![0, 2, 5, 8, 8];

    let interpolator = FiniteElementInterpolator::from_compressed_values(
        values.into_iter().zip(supported_nodes).collect(),
        supported_node_offsets,
    );

    let u = DVector::from_column_slice(&[1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0]);

    let mut interpolated_values = vec![Vector2::zeros(); 4];

    interpolator.interpolate_into(&mut interpolated_values, &u);

    let v = interpolated_values;
    assert_scalar_eq!(v[0].x, 7.0);
    assert_scalar_eq!(v[0].y, 10.0);
    assert_scalar_eq!(v[1].x, 76.0);
    assert_scalar_eq!(v[1].y, 88.0);
    assert_scalar_eq!(v[2].x, 125.0);
    assert_scalar_eq!(v[2].y, 146.0);
    assert_scalar_eq!(v[3].x, 0.0);
    assert_scalar_eq!(v[3].y, 0.0);
}

// TODO: Use this in rewritten tests?
// fn find_interior_point<T>(polygon: &ConvexPolygon<T>) -> Point2<T>
// where
//     T: Real,
// {
//     // Find an interior point by averaging all vertices (note: this is not in general the centroid)
//     let num_vertices = polygon.vertices().len();
//     let vertex_sum = polygon
//         .vertices()
//         .iter()
//         .map(|p| p.coords)
//         .fold(Vector2::zeros(), |a, b| a + b);
//     Point2::from(vertex_sum / T::from_usize(num_vertices).unwrap())
// }

proptest! {

    // TODO: Rewrite this test without NodalModel/Quad4Model
    // #[test]
    // fn finite_element_interpolator_is_identity_for_node_vertices(
    //     mesh in rectangular_uniform_mesh_strategy(1.0, 8))
    // {
    //     prop_assume!(mesh.connectivity().len() > 0);
    //
    //     let quadrature = quad_quadrature_strength_5_f64();
    //     let model = Quad4Model::from_mesh_and_quadrature(mesh.clone(), quadrature);
    //     let interpolation_vertices = mesh.vertices();
    //     let interpolator = model.make_interpolator(interpolation_vertices).unwrap();
    //
    //     // Assume that the solution variable is x, with x corresponding to the position
    //     // of vertices. So x_h(X) = sum_i N_i(X) x_i, where x_h(X) is the deformed
    //     // position at reference position X. Then, clearly, if x_i == X_i, with X_i
    //     // being a node in the Lagrange finite element basis, then
    //     // x_h(X_i) == X_i.
    //     let vertices: Vec<_> = mesh.vertices()
    //         .iter()
    //         .map(|x| x.coords)
    //         .collect();
    //
    //     let solution_variables = flatten_vertically(&vertices).unwrap();
    //     let mut result = vec![Vector2::zeros(); interpolation_vertices.len()];
    //
    //     interpolator.interpolate_into(&mut result, &solution_variables);
    //
    //     assert_eq!(result.len(), interpolation_vertices.len());
    //     for (v_result, v_expected) in result.iter().zip(interpolation_vertices) {
    //         // TODO: Use matrixcompare
    //         let diff = v_result - v_expected.coords;
    //         prop_assert!(diff.norm() <= 1e-6);
    //     }
    // }

    // TODO: Rewrite this test without NodalModel/Quad4Model
    // #[test]
    // fn finite_element_interpolation_at_interior_points_is_bounded_by_nodal_values(
    //     (mesh, u) in rectangular_uniform_mesh_strategy(1.0, 8).prop_flat_map(|mesh| {
    //         let ndof = mesh.vertices().len();
    //         let u_strategy = vec(-10.0..10.0, ndof).prop_map(move |v| DVector::from_iterator(ndof, v));
    //         (Just(mesh), u_strategy)
    //     }))
    // {
    //     assert_eq!(mesh.vertices().len(), u.len());
    //     prop_assume!(mesh.connectivity().len() > 0);
    //
    //     // It is easy to show that |u_h(x)| = |sum_i N_i(x) u_i| <= max_i |u_i|,
    //     // where u_i are the nodal values of the nodes for which x is in the support of the
    //     // nodal basis function.
    //
    //     let quadrature = quad_quadrature_strength_5_f64();
    //     let model = Quad4Model::from_mesh_and_quadrature(mesh.clone(), quadrature);
    //
    //     let interpolation_vertices: Vec<_> = mesh.connectivity().iter()
    //         .map(|connectivity| connectivity.cell(mesh.vertices()).unwrap())
    //         .map(|cell| ConvexPolygon::try_from(cell).expect("Meshes should only have convex cells"))
    //         .map(|polygon| find_interior_point(&polygon))
    //         .collect();
    //
    //     let interpolator = model.make_interpolator(&interpolation_vertices).unwrap();
    //     let mut result = vec![OVector::<f64, U1>::zeros(); interpolation_vertices.len()];
    //
    //     interpolator.interpolate_into(&mut result, &u);
    //
    //     prop_assert_eq!(result.len(), interpolation_vertices.len());
    //
    //     for (i, u_result) in result.iter().map(|p| p[0]).enumerate() {
    //         let cell_connectivity = mesh.connectivity()[i];
    //         let smaller_than_neighbor_nodes = cell_connectivity.vertex_indices()
    //             .iter()
    //             .map(|idx| u[*idx])
    //             .any(|u_j| u_result.abs() <= u_j.abs());
    //         prop_assert!(smaller_than_neighbor_nodes);
    //     }
    // }
}