fenris 0.0.33

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
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178

//! Tests for estimation of L2/H1 seminorm errors approximated by using a higher-resolution
//! reference solution.
use fenris::assembly::local::UniformQuadratureTable;
use fenris::error::{estimate_H1_seminorm_error, estimate_L2_error, SpaceInterpolationFn};
use fenris::mesh::procedural::create_unit_square_uniform_tri_mesh_2d;
use fenris::quadrature;
use fenris::space::SpatiallyIndexed;
use fenris::util::global_vector_from_point_fn;
use itertools::izip;
use nalgebra::coordinates::XY;
use nalgebra::{vector, Point2, Vector1, Vector2, U1};
use serde::{Deserialize, Serialize};
use std::f64::consts::PI;
use std::fs::{create_dir_all, File};
use std::ops::Deref;
use std::path::PathBuf;

fn sin(x: f64) -> f64 {
    x.sin()
}
fn cos(x: f64) -> f64 {
    x.cos()
}

fn u_2d(x: &Point2<f64>) -> Vector1<f64> {
    let &XY { x, y } = x.coords.deref();
    vector![sin(PI * x) * sin(PI * y)]
}

fn u_grad_2d(x: &Point2<f64>) -> Vector2<f64> {
    let &XY { x, y } = x.coords.deref();
    let u_x = PI * cos(PI * x) * sin(PI * y);
    let u_y = PI * sin(PI * x) * cos(PI * y);
    Vector2::new(u_x, u_y)
}

#[derive(Debug, Default, Clone, PartialEq, Serialize, Deserialize)]
#[allow(non_snake_case)]
pub struct ErrorSample {
    pub coarse_res: usize,
    /// The resolution of the fine mesh used for error estimation, or `0` if analytic
    pub fine_res: usize,
    pub L2_error: f64,
    pub H1_semi_error: f64,
}

/// For serializing to JSON for subsequent analysis/plots
#[derive(Debug, Default, PartialEq, Serialize, Deserialize)]
#[allow(non_snake_case)]
pub struct ErrorSummary {
    pub element_name: String,
    pub samples: Vec<ErrorSample>,
}

fn export_summary(summary: &ErrorSummary) {
    let base_path = PathBuf::from("data/convergence_tests/error_estimation/");
    let path = base_path.join(format!("{}_summary.json", summary.element_name.to_ascii_lowercase()));
    create_dir_all(path.parent().unwrap()).unwrap();
    let mut file = File::create(path).unwrap();
    serde_json::to_writer_pretty(&mut file, summary).unwrap();
}

fn load_reference_summary(element_name: &str) -> ErrorSummary {
    let path = format!(
        "tests/convergence_tests/reference_values/error_estimation_{}_summary.json",
        element_name.to_ascii_lowercase()
    );
    let json_str = std::fs::read_to_string(&path).expect(&format!(
        "I/O error: Failed to load reference summary for \
                          element {} at location \"{}\"",
        element_name, &path
    ));
    match serde_json::from_str(&json_str) {
        Ok(summary) => summary,
        Err(error) => panic!(
            "Failed to deserialize reference summary for element {}. Error: {}",
            element_name, error
        ),
    }
}

fn rel_error(value: f64, reference: f64) -> f64 {
    (value - reference).abs() / reference.abs()
}

#[allow(non_snake_case)]
fn assert_summary_close_to_reference(summary: &ErrorSummary) {
    let reference = load_reference_summary(&summary.element_name);

    assert_eq!(summary.element_name, reference.element_name, "Element name mismatch");
    assert_eq!(summary.samples.len(), reference.samples.len(), "Sample count mismatch");

    for (summary_sample, reference_sample) in izip!(&summary.samples, &reference.samples) {
        assert_eq!(
            summary_sample.coarse_res, reference_sample.coarse_res,
            "Coarse resolution mismatch"
        );
        assert_eq!(
            summary_sample.fine_res, reference_sample.fine_res,
            "Fine resolution mismatch"
        );

        let L2_rel_error = rel_error(summary_sample.L2_error, reference_sample.L2_error);
        let H1_semi_rel_error = rel_error(summary_sample.H1_semi_error, reference_sample.H1_semi_error);
        assert!(
            L2_rel_error <= 0.01,
            "L2 error deviates by more than 1% (summary: {}, reference: {}, rel error: {}",
            summary_sample.L2_error,
            reference_sample.L2_error,
            L2_rel_error
        );
        assert!(
            H1_semi_rel_error <= 0.01,
            "H1 seminorm error deviates by more than 1% (summary: {}, reference: {}, rel error: {}",
            summary_sample.H1_semi_error,
            reference_sample.H1_semi_error,
            H1_semi_rel_error
        );
    }
}

#[allow(non_snake_case)]
#[test]
fn tri3_error_estimation() {
    let fine_resolutions = vec![64, 71, 91, 128, 131, 512, 1024];
    let coarse_resolutions = vec![1, 2, 3, 4, 8, 16, 32];

    let base_error_quadrature = quadrature::total_order::triangle(20).unwrap();
    let qtable = UniformQuadratureTable::from_quadrature(base_error_quadrature);

    let mut summary = ErrorSummary::default();
    summary.element_name = "Tri3".into();

    for &fine_res in &fine_resolutions {
        let mesh_fine = create_unit_square_uniform_tri_mesh_2d(fine_res);
        let u_weights_fine = global_vector_from_point_fn(mesh_fine.vertices(), u_2d);
        let fine_interpolator = SpatiallyIndexed::from_space(mesh_fine);
        let u_fine = SpaceInterpolationFn(&fine_interpolator, &u_weights_fine);
        for &coarse_res in &coarse_resolutions {
            let mesh_coarse = create_unit_square_uniform_tri_mesh_2d(coarse_res);
            let u_weights_coarse = global_vector_from_point_fn(mesh_coarse.vertices(), u_2d);
            // TODO: Can we make SolutionDim unambiguous, so that we don't have to specify generic types?
            // I tried to make SolutionDim an associated type of SolutionFunction/Gradient,
            // but unfortunately then we run into some compiler bugs where rustc gets
            // the types wrong :-(
            let L2_error_approx =
                estimate_L2_error::<_, U1, _, _>(&mesh_coarse, &u_fine, &u_weights_coarse, &qtable).unwrap();
            let H1_semi_error_approx =
                estimate_H1_seminorm_error::<_, U1, _, _>(&mesh_coarse, &u_fine, &u_weights_coarse, &qtable).unwrap();
            let sample = ErrorSample {
                coarse_res,
                fine_res,
                L2_error: L2_error_approx,
                H1_semi_error: H1_semi_error_approx,
            };

            summary.samples.push(sample);
        }
    }

    for &coarse_res in &coarse_resolutions {
        let mesh_coarse = create_unit_square_uniform_tri_mesh_2d(coarse_res);
        let u_weights_coarse = global_vector_from_point_fn(mesh_coarse.vertices(), u_2d);
        let L2_error = estimate_L2_error::<_, U1, _, _>(&mesh_coarse, &u_2d, &u_weights_coarse, &qtable).unwrap();
        let H1_semi_error =
            estimate_H1_seminorm_error::<_, U1, _, _>(&mesh_coarse, &u_grad_2d, &u_weights_coarse, &qtable).unwrap();
        let sample = ErrorSample {
            coarse_res,
            fine_res: 0,
            L2_error,
            H1_semi_error,
        };
        summary.samples.push(sample);
    }

    export_summary(&summary);
    assert_summary_close_to_reference(&summary);
}