fenris 0.0.32

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
134
135
136
137
138
139
140
141
142
143

use crate::element::{Tet4Element, Tri3d2Element, Tri3d3Element};
use crate::geometry::Orientation::Counterclockwise;
use crate::mesh::procedural::create_rectangular_uniform_quad_mesh_2d;
use crate::mesh::QuadMesh2d;
use ::proptest::prelude::*;
use fenris_geometry::proptest::Triangle2dParams;
use fenris_geometry::Triangle2d;
use nalgebra::{Point2, Point3, Vector2};
use std::cmp::max;

pub fn point2() -> impl Strategy<Value = Point2<f64>> {
    // Pick a reasonably small range to pick coordinates from,
    // otherwise we can easily get floating point numbers that are
    // so ridiculously large as to break anything we might want to do with them
    let range = -10.0..10.0;
    [range.clone(), range.clone()].prop_map(|[x, y]| Point2::new(x, y))
}

pub fn point3() -> impl Strategy<Value = Point3<f64>> {
    // Pick a reasonably small range to pick coordinates from,
    // otherwise we can easily get floating point numbers that are
    // so ridiculously large as to break anything we might want to do with them
    let range = -10.0..10.0;
    [range.clone(), range.clone(), range.clone()].prop_map(|[x, y, z]| Point3::new(x, y, z))
}

impl Arbitrary for Tri3d2Element<f64> {
    type Parameters = ();
    type Strategy = BoxedStrategy<Self>;

    fn arbitrary_with(_args: Self::Parameters) -> Self::Strategy {
        any_with::<Triangle2d<f64>>(Triangle2dParams::default().with_orientation(Counterclockwise))
            .prop_map(|triangle| Self::from(triangle))
            .boxed()
    }
}

impl Arbitrary for Tri3d3Element<f64> {
    type Parameters = ();
    type Strategy = BoxedStrategy<Self>;

    fn arbitrary_with(_args: Self::Parameters) -> Self::Strategy {
        let vertices: [_; 3] = std::array::from_fn(|_| point3());
        vertices
            .prop_map(|vertices| Tri3d3Element::from_vertices(vertices))
            .boxed()
    }
}

impl Arbitrary for Tet4Element<f64> {
    // TODO: Reasonable parameters?
    type Parameters = ();
    type Strategy = BoxedStrategy<Self>;

    fn arbitrary_with(_args: Self::Parameters) -> Self::Strategy {
        let l = 5.0;
        (Tri3d3Element::arbitrary(), -l..l, -l..l, 0.0..l)
            .prop_map(|(tri_element, t1_scale, t2_scale, n_scale)| {
                let [a, b, c] = tri_element.vertices().clone();
                let t1 = b - a;
                let t2 = c - a;
                let n = t1.cross(&t2);
                let d = a + t1_scale * t1 + t2_scale * t2 + n_scale * n;
                Tet4Element::from_vertices([a, b, c, d])
            })
            .boxed()
    }
}

// Returns a strategy in which each value is a triplet (cells_per_unit, units_x, units_y)
// such that cells_per_unit^2 * units_x * units_y <= max_cells
fn rectangular_uniform_mesh_cell_distribution_strategy(
    max_cells: usize,
) -> impl Strategy<Value = (usize, usize, usize)> {
    let max_cells_per_unit = f64::floor(f64::sqrt(max_cells as f64)) as usize;
    (1..=max(1, max_cells_per_unit))
        .prop_flat_map(move |cells_per_unit| (Just(cells_per_unit), 0..=max_cells / (cells_per_unit * cells_per_unit)))
        .prop_flat_map(move |(cells_per_unit, units_x)| {
            let units_y_strategy = 0..=max_cells / (cells_per_unit * cells_per_unit * max(1, units_x));
            (Just(cells_per_unit), Just(units_x), units_y_strategy)
        })
}

pub fn rectangular_uniform_mesh_strategy(unit_length: f64, max_cells: usize) -> impl Strategy<Value = QuadMesh2d<f64>> {
    rectangular_uniform_mesh_cell_distribution_strategy(max_cells).prop_map(
        move |(cells_per_unit, units_x, units_y)| {
            create_rectangular_uniform_quad_mesh_2d(
                unit_length,
                units_x,
                units_y,
                cells_per_unit,
                &Vector2::new(0.0, 0.0),
            )
        },
    )
}

#[cfg(test)]
mod tests {
    use super::rectangular_uniform_mesh_cell_distribution_strategy;
    use crate::geometry::proptest::{
        convex_quad2d_strategy_f64, nondegenerate_convex_quad2d_strategy_f64, nondegenerate_triangle2d_strategy_f64,
    };
    use crate::geometry::Orientation;
    use proptest::prelude::*;

    proptest! {
        #[test]
        fn rectangular_uniform_mesh_cell_distribution_strategy_respects_max_cells(
            (max_cells, cells_per_unit, units_x, units_y)
             in (0..20usize).prop_flat_map(|max_cells| {
                rectangular_uniform_mesh_cell_distribution_strategy(max_cells)
                    .prop_map(move |(cells_per_unit, units_x, units_y)| {
                    (max_cells, cells_per_unit, units_x, units_y)
                })
             })
        ) {
            // Test that the distribution strategy for rectangular meshes
            // respects the maximum number of cells given
            prop_assert!(cells_per_unit * cells_per_unit * units_x * units_y <= max_cells);
        }

        #[test]
        fn convex_quads_are_convex(quad in convex_quad2d_strategy_f64()) {
            prop_assert!(quad.concave_corner().is_none());
        }

        #[test]
        fn nondegenerate_triangles_have_positive_area(
            triangle in nondegenerate_triangle2d_strategy_f64()
        ){
            prop_assert!(triangle.area() > 0.0);
            prop_assert!(triangle.orientation() == Orientation::Counterclockwise);
        }

        #[test]
        fn nondegenerate_quads_have_positive_area(
            quad in nondegenerate_convex_quad2d_strategy_f64()
        ){
            prop_assert!(quad.area() > 0.0);
        }
    }
}