fem_2d/
lib.rs

1#![cfg_attr(
2    feature = "max_ortho_basis",
3    feature(const_fn_floating_point_arithmetic)
4)]
5#![doc = include_str!("../README.md")]
6
7/// Structures defining the FEM Domain, Mesh, and Basis Space
8pub mod fem_domain;
9
10/// Structures and Functions to define and solve the FEM Problem
11pub mod fem_problem;
12
13/// Convenient Re-Exports
14pub mod prelude {
15    pub use crate::fem_domain::basis::hierarchical_basis_fns::poly::HierPoly;
16    #[cfg(feature = "max_ortho_basis")]
17    pub use crate::fem_domain::basis::shape_fns::max_ortho::MaxOrthoShapeFn;
18    pub use crate::fem_domain::domain::{
19        dof::{
20            basis_spec::{BSAddress, BasisSpec},
21            DoF,
22        },
23        fields::UniformFieldSpace,
24        mesh::{
25            elem::Elem,
26            h_refinement::{HRef, HRefError},
27            p_refinement::{PRef, Orders, PRefError},
28            Mesh,
29        },
30        ContinuityCondition, Domain,
31    };
32    pub use crate::fem_problem::galerkin::{galerkin_sample_gep_hcurl, GalerkinSamplingError};
33    pub use crate::fem_problem::integration::integrals::{curl_curl::CurlCurl, inner::L2Inner};
34    pub use crate::fem_problem::linalg::{
35        nalgebra_solve::{nalgebra_solve_gep, NalgebraGEPError},
36        slepc_solve::{slepc_solve_gep, SlepcGEPError},
37        EigenPair, GEP,
38    };
39}
40
41#[cfg(test)]
42mod tests {
43    use super::prelude::*;
44
45    #[test]
46    fn nalg_problem() {
47        // Define Mesh
48        let mut mesh = Mesh::from_file("./test_input/test_mesh_a.json").unwrap();
49        mesh.global_p_refinement(PRef::from(2, 2));
50
51        // Construct Domain
52        let domain = Domain::from_mesh(mesh, ContinuityCondition::HCurl);
53        let ndofs = domain.dofs.len();
54        println!("Domain constructed with {} Degrees of Freedom", ndofs);
55
56        // Fill Matrices
57        let eigenproblem =
58            galerkin_sample_gep_hcurl::<HierPoly, CurlCurl, L2Inner>(&domain, Some([8, 8]))
59                .unwrap();
60
61        // Solve Eigenvalue Problem
62        let solution = nalgebra_solve_gep(eigenproblem, 2.64).unwrap();
63        println!("Found eigenvalue: {:.15}", solution.value);
64
65        assert!((solution.value - 2.6479657_f64).abs() < 1e-6);
66        assert_eq!(solution.vector.len(), ndofs);
67
68        let mut field_space = UniformFieldSpace::new(&domain, [8, 8]);
69        let e_field_names = field_space
70            .xy_fields::<HierPoly>("E", solution.normalized_eigenvector())
71            .unwrap();
72        field_space
73            .expression_2arg(e_field_names, "E_mag", |ex, ey| {
74                (ex.powi(2) + ey.powi(2)).sqrt()
75            })
76            .unwrap();
77        field_space
78            .print_all_to_vtk("./test_output/mesh_b_fields.vtk")
79            .unwrap();
80    }
81
82    #[test]
83    fn slepc_problem() {
84        // Define Mesh
85        let mut mesh = Mesh::from_file("./test_input/test_mesh_b.json").unwrap();
86        mesh.global_p_refinement(PRef::from(3, 3));
87        mesh.global_h_refinement(HRef::T);
88        mesh.h_refine_elems(vec![6, 9, 12], HRef::T).unwrap();
89
90        // Construct Domain
91        let domain = Domain::from_mesh(mesh, ContinuityCondition::HCurl);
92        let ndofs = domain.dofs.len();
93        println!("Domain constructed with {} Degrees of Freedom", ndofs);
94
95        // Fill Matrices
96        let eigenproblem =
97            galerkin_sample_gep_hcurl::<HierPoly, CurlCurl, L2Inner>(&domain, Some([8, 8]))
98                .unwrap();
99
100        // Solve Eigenvalue Problem
101        let solution = slepc_solve_gep(eigenproblem, 1.475).unwrap();
102        println!("Found eigenvalue: {:.15}", solution.value);
103
104        assert!((solution.value - 1.4745880937_f64).abs() < 1e-9);
105        assert_eq!(solution.vector.len(), ndofs);
106
107        let mut field_space = UniformFieldSpace::new(&domain, [8, 8]);
108        let e_field_names = field_space
109            .xy_fields::<HierPoly>("E", solution.normalized_eigenvector())
110            .unwrap();
111        field_space
112            .expression_2arg(e_field_names, "E_mag", |ex, ey| {
113                (ex.powi(2) + ey.powi(2)).sqrt()
114            })
115            .unwrap();
116        field_space
117            .print_all_to_vtk("./test_output/mesh_b_fields.vtk")
118            .unwrap();
119    }
120}
fem_2d/lib.rs

fem_2d/
lib.rs
