# fem_2d A Rust library for 2D Finite Element Method computations, featuring: - Highly flexible *hp*-Refinement - Isotropic & Anisotropic *h*-refinements (with support for n-irregularity) - Isotropic & Anisotropic *p*-refinements - Generic shape function evaluation - You can use one of the two built in sets of H(curl) conforming Shape Functions - Or you can define your own by implementing the `ShapeFn` Trait - Two Eigensolvers - Sparse: Using an external Slepc Solver (code and installation instructions found [here](https://github.com/jeremiah-corrado/slepc_gep_solver)) - Dense: Using [Nalgebra](https://nalgebra.org/docs/user_guide/decompositions_and_lapack#eigendecomposition-of-a-hermitian-matrix)'s Eigen-Decomposition (not recommended for large or ill-conditioned problems) - Expressive Solution Evaluation - Field solutions can easily be generated from an eigenvector - Arbitrary functions of solutions can also be evaluated (ex: magnitude of a field) - Solutions and expressions are easily printed to `.vtk` files for plotting (using [VISIT](https://visit-dav.github.io/visit-website/index.html) or similar tools) ## Usage Include this line in your `Cargo.toml` file under **\[dependencies\]**: ```toml fem_2d = "0.1.0" ``` Please include one or more of the following citations in any academic or commercial work based on this repository: * Corrado, Jeremiah; Harmon, Jake; Notaros, Branislav; Ilic, Milan M. (2022): FEM_2D: A Rust Package for 2D Finite Element Method Computations with Extensive Support for hp-refinement. TechRxiv. Preprint. [https://doi.org/10.36227/techrxiv.19166339.v1](https://doi.org/10.36227/techrxiv.19166339.v1) * Corrado, Jeremiah; Harmon, Jake; Notaros, Branislav (2021): A Refinement-by-Superposition Approach to Fully Anisotropic hp-Refinement for Improved Efficiency in CEM. TechRxiv. Preprint. [https://doi.org/10.36227/techrxiv.16695163.v1](https://doi.org/10.36227/techrxiv.16695163.v1) * Harmon, Jake; Corrado, Jeremiah; Notaros, Branislav (2021): A Refinement-by-Superposition hp-Method for H(curl)- and H(div)-Conforming Discretizations. TechRxiv. Preprint. [https://doi.org/10.36227/techrxiv.14807895.v1](https://doi.org/10.36227/techrxiv.14807895.v1) ## Documentation The latest Documentation can be found [here](https://docs.rs/fem_2d/0.2.1/fem_2d/) ## Example Solve the Maxwell Eigenvalue Problem on a standard Waveguide and print the Electric Fields to a VTK file. This example encompasses most of the functionality of the library: ```rust use fem_2d::prelude::*; fn solve_basic_problem() -> Result<(), Box> { // Load a standard air-filled waveguide mesh from a JSON file let mut mesh = Mesh::from_file("./test_input/test_mesh_a.json")?; // Set the polynomial expansion order to 4 in both directions on all Elems mesh.set_global_expansion_orders(Orders::new(4, 4)); // Isotropically refine all Elems mesh.global_h_refinement(HRef::t()); // Then anisotropically refine the resultant Elems in the center of the mesh let cenral_node_id = mesh.elems[0].nodes[3]; mesh.h_refine_with_filter(|elem| { if elem.nodes.contains(&cenral_node_id) { Some(HRef::u()) } else { None } }); // Construct a Domain with H(curl) continuity conditions let domain = Domain::from_mesh(mesh, ContinuityCondition::HCurl); println!("Constructed Domain with {} DoFs", domain.dofs.len()); // Construct a generalized eigenvalue problem for the Electric Field // (in parallel using the Rayon Global ThreadPool) let gep = galerkin_sample_gep_hcurl::< HierPoly, CurlCurl, L2Inner >(&domain, Some([8, 8]))?; // Solve the generalized eigenvalue problem using Nalgebra's Eigen-Decomposition // look for an eigenvalue close to 10.0 let solution = nalgebra_solve_gep(gep, 10.0)?; println!("Found Eigenvalue: {:.15}", solution.value); // Construct a solution-field-space over the Domain with 64 samples on each "leaf" Elem let mut field_space = UniformFieldSpace::new(&domain, [8, 8]); // Compute the Electric Field in the X- and Y-directions (using the same ShapeFns as above) let e_field_names = field_space.xy_fields:: ("E", solution.normalized_eigenvector())?; // Compute the magnitude of the Electric Field field_space.expression_2arg(e_field_names, "E_mag", |ex, ey| { (ex.powi(2) + ey.powi(2)).sqrt() })?; // Print "E_x", "E_y" and "E_mag" to a VTK file field_space.print_all_to_vtk("./test_output/electric_field_solution.vtk")?; Ok(()) } ``` ## Mesh Refinement A `Mesh` structure keeps track of the geometric layout of the finite elements (designated as `Elem`s in the library), as well as the polynomial expansion orders on each element. These can be updated using h- and p-refinements respectively ### *h*-Refinement: *h*-Refinements* are implemented using the Refinement by Superposition (RBS) method > Technical details can be found in this paper: [A Refinement-by-Superposition Approach to Fully Anisotropichp-Refinement for Improved Efficiency in CEM](https://www.techrxiv.org/articles/preprint/A_Refinement-by-Superposition_Approach_to_Fully_Anisotropic_hp-Refinement_for_Improved_Efficiency_in_CEM/16695163) Three types of h-refinement are supported: * **T**: Elements are superimposed with 4 equally sized child elements * **U**: Elements are superimposed with 2 child elements, such that the resolution is improved in the x-direction * **V**: Elements are superimposed with 2 child elements, such that the resolution is improved in the y-direction These are designated as an Enum: `HRef`, located in the `h_refinements` module. They can be executed by constructing a refinement as follows: ```Rust let h_iso = HRef::T; let h_aniso_u = HRef::U(None); let h_aniso_v = HRef::V(None); ``` ...and applying it to an element or group of elements using one of the many *h*-refinement methods on `Mesh`. Multi-step anisotropic *h*-refinements can be executed by constructing the U or V variant with `Some(0)` or `Some(1)`. This will cause the 0th or 1st resultant child element to be anisotropically refined in the opposite direction. Mesh coarsening is not currently supported ### *p*-Refinement: *p*-Refinements allow elements to support a range of expansion orders in the X and Y directions. These can be modified separately for greater control over resource usage and solution accuracy. As a `Domain` is constructed from a `Mesh`, Basis Functions are constructed based on the elements expansion orders. Expansion orders can be increased or decreased by constructing a `PRef`, located in the `p_refinements` module: ```Rust let uv_plus_2 = PRef::from(2, 2); let u_plus_1_v_minus_3 = PRef::from(1, -3); ``` ...and applying it to an element or group of elements using one of the many *p*-refinement methods on `Mesh`. ## JSON Mesh Files fem_2d uses `.json` files to import and export Mesh layouts. * The input format is simplified and describes the geometry of the problem only. * The output format describes the Mesh in it's refined state which is usefull for debugging and observing the refinement state. ### Input Mesh Files A `Mesh` can be constructed from a JSON file with the following format: ```JSON { "Elements": [ { "materials": [eps_rel_re, eps_rel_im, mu_rel_re, mu_rel_im], "node_ids": [node_0_id, node_1_id, node_2_id, node_3_id], }, { "materials": [1.0, 0.0, 1.0, 0.0], "node_ids": [1, 2, 4, 5], }, { "materials": [1.2, 0.0, 0.9999, 0.0], "node_ids": [2, 3, 5, 6], } ], "Nodes": [ [x_coordinate, y_coordinate], [0.0, 0.0], [1.0, 0.0], [2.0, 0.0], [0.0, 0.5], [1.0, 0.5], [2.0, 0.5], ] } ``` (The first Element and Node are simply there to explain what variables mean. Those should not be included in an actual mesh file!) The above file corresponds to this 2 Element mesh (with Node indices labeled): ```text 3 4 5 0.5 *---------------*---------------* | | | | air | teflon | | | | 0.0 *---------------*---------------* y 0 1 2 x 0.0 1.0 2.0 Both Elements start with a polynomial expansion order of 1 upon construction. ``` This library does not yet support curvilinear elements. When that feature is added, this file format will also be extended to describe higher-order geometry. ### Output Mesh Files A refined `Mesh` can also be exported and visualized using [this](https://github.com/jeremiah-corrado/fem_2d_mesh_plot) tool: ## Solution Plotting Solution figures like the ones below can be generated by exporting solutions as `.vtk` files and plotting using a compatible tool (such as [Visit](https://visit-dav.github.io/visit-website/index.html)): ## Community Guidelines / Code of Conduct Contributions, questions, and bug-reports are welcome! Any pull requests, issues, etc. must adhere to Rust's [Code of Conduct](https://www.rust-lang.org/policies/code-of-conduct). More details about how to contribute to this project can be found [here](https://github.com/jeremiah-corrado/fem_2d/blob/main/rm_figs/contributing.md). Questions and bug-reports should be directed to the Issues tab.