fenris 0.0.28

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

use eyre::eyre;
use fenris::assembly::global::{
    apply_homogeneous_dirichlet_bc_csr, apply_homogeneous_dirichlet_bc_rhs, CsrAssembler, VectorAssembler,
};
use fenris::assembly::local::{
    ElementEllipticAssemblerBuilder, ElementSourceAssemblerBuilder, SourceFunction, UniformQuadratureTable,
};
use fenris::assembly::operators::{LaplaceOperator, Operator};
use fenris::io::vtk::FiniteElementMeshDataSetBuilder;
use fenris::mesh::procedural::create_unit_square_uniform_quad_mesh_2d;
use fenris::mesh::QuadMesh2d;
use fenris::nalgebra::{DMatrix, DVector, Point2, U1, U2};
use fenris::nalgebra_sparse::CsrMatrix;
use fenris::quadrature;
use nalgebra::Vector1;

fn main() -> eyre::Result<()> {
    // TODO: Make it easy to construct triangle meshes as well.
    // Need to make it easy to convert between different meshes, such as Quad2d -> Tri2d
    let mesh: QuadMesh2d<f64> = create_unit_square_uniform_quad_mesh_2d(4);

    let (a, b) = assemble_linear_system(&mesh)?;
    let u = solve_linear_system(&a, &b)?;

    FiniteElementMeshDataSetBuilder::from_mesh(&mesh)
        .with_title("Poisson 2D")
        .with_point_scalar_attributes("u", 1, u.as_slice())
        .try_export("poisson2d.vtu")?;

    Ok(())
}

fn assemble_linear_system(mesh: &QuadMesh2d<f64>) -> eyre::Result<(CsrMatrix<f64>, DVector<f64>)> {
    // A quadrature table is responsible for providing each element with a quadrature rule.
    // Since we want the same quadrature rule per element, we use a uniform quadrature table.
    let (weights, points) = quadrature::tensor::quadrilateral_gauss(2);
    let quadrature = UniformQuadratureTable::from_points_and_weights(points, weights);

    // TODO: This isn't actually needed. Get rid of it by introducing a separate trait
    // for linear contractions
    let u = DVector::<f64>::zeros(mesh.vertices().len());

    // Set up global assemblers. These are responsible for adding element matrix/vector entries from
    // each element into the global matrix/vector
    let vector_assembler = VectorAssembler::<f64>::default();
    let matrix_assembler = CsrAssembler::default();

    // Set up local matrix assembler for the Laplace operator. This assembler is responsible
    // for building the local element matrices corresponding to each element
    let laplace_assembler = ElementEllipticAssemblerBuilder::new()
        .with_finite_element_space(mesh)
        .with_operator(&LaplaceOperator)
        .with_quadrature_table(&quadrature)
        // TODO: If the operator is linear, u is not actually needed...
        // How to reflect this in the API?
        .with_u(&u)
        .build();

    let mut a_global = matrix_assembler.assemble(&laplace_assembler)?;

    // Set up a local vector assembler for the "source term", which corresponds to the weak form
    // term (f, v), where v is a test function and f is the source function `f` in the
    // Poisson equation delta u = f.
    let source_assembler = ElementSourceAssemblerBuilder::new()
        .with_finite_element_space(mesh)
        .with_quadrature_table(&quadrature)
        .with_source(&PoissonProblemSourceFunction)
        .build();

    let mut b_global = vector_assembler.assemble_vector(&source_assembler)?;

    // We want to have a Dirichlet boundary for x == 0. To account for slight numerical errors,
    // we determine the indices of the Dirichlet nodes by extracting those node indices
    // which satisfy x < eps, for some small epsilon.
    let dirichlet_nodes: Vec<_> = mesh
        .vertices()
        .iter()
        .enumerate()
        .filter_map(|(idx, v)| (v.x < 1e-12).then(|| idx))
        .collect();

    apply_homogeneous_dirichlet_bc_csr(&mut a_global, &dirichlet_nodes, 1);
    apply_homogeneous_dirichlet_bc_rhs(&mut b_global, &dirichlet_nodes, 1);

    Ok((a_global, b_global))
}

fn solve_linear_system(matrix: &CsrMatrix<f64>, rhs: &DVector<f64>) -> eyre::Result<DVector<f64>> {
    // TODO: Use sparse solver
    let matrix = DMatrix::from(matrix);
    // The discrete Laplace operator is positive definite (given appropriate boundary conditions),
    // so we can use a Cholesky factorization
    let cholesky = matrix
        .cholesky()
        .ok_or_else(|| eyre!("Failed to solve linear system"))?;
    Ok(cholesky.solve(rhs))
}

/// Represents the source function `f` in the poisson equation - Delta u = f.
struct PoissonProblemSourceFunction;

impl Operator<f64, U2> for PoissonProblemSourceFunction {
    /// f maps R^2 to R^1 (U1)
    type SolutionDim = U1;
    type Parameters = ();
}

impl SourceFunction<f64, U2> for PoissonProblemSourceFunction {
    fn evaluate(&self, _coords: &Point2<f64>, _data: &Self::Parameters) -> Vector1<f64> {
        // TODO: Use a more interesting function than a constant function
        Vector1::new(1.0)
    }
}