diff --git a/.github/workflows/run-tests.yml b/.github/workflows/run-tests.yml index 36ec4561..b03bd914 100644 --- a/.github/workflows/run-tests.yml +++ b/.github/workflows/run-tests.yml @@ -14,7 +14,7 @@ jobs: runs-on: ubuntu-latest strategy: matrix: - rust-version: ["stable"] + rust-version: ["stable", "nightly"] mpi: ['mpich', 'openmpi'] steps: - name: Set up Rust @@ -34,28 +34,28 @@ jobs: sudo apt-get install libopenblas-dev liblapack-dev - name: Style checks run: | - cargo fmt -- --check - cargo clippy -- -D warnings - cargo clippy --tests -- -D warnings - cargo clippy --examples -- -D warnings + cargo +${{ matrix.rust-version }} fmt -- --check + cargo +${{ matrix.rust-version }} clippy -- -D warnings + cargo +${{ matrix.rust-version }} clippy --tests -- -D warnings + cargo +${{ matrix.rust-version }} clippy --examples -- -D warnings - name: Build rust library - run: cargo build --features "strict" + run: cargo +${{ matrix.rust-version }} build --features "strict" - name: Run unit tests - run: cargo test --lib --features "strict" + run: cargo +${{ matrix.rust-version }} test --lib --features "strict" - name: Run unit tests (with mpi enabled) - run: cargo test --lib --features "mpi,strict" + run: cargo +${{ matrix.rust-version }} test --lib --features "mpi,strict" - name: Run tests - run: cargo test --examples --release --features "mpi,strict" + run: cargo +${{ matrix.rust-version }} test --examples --release --features "mpi,strict" - name: Run examples run: | - python3 find_examples.py + python3 find_examples.py --cargo-version ${{ matrix.rust-version }} chmod +x examples.sh ./examples.sh - name: Build docs - run: cargo doc --features "strict" --no-deps + run: cargo +${{ matrix.rust-version }} doc --features "strict" --no-deps run-tests-python: name: Run tests (Python) diff --git a/.github/workflows/run-weekly-tests.yml b/.github/workflows/run-weekly-tests.yml index 302e2ca7..b2c3a973 100644 --- a/.github/workflows/run-weekly-tests.yml +++ b/.github/workflows/run-weekly-tests.yml @@ -30,23 +30,20 @@ jobs: sudo apt-get install libopenblas-dev liblapack-dev - name: Build rust library (debug) - run: cargo build --features "strict" + run: cargo +${{ matrix.rust-version }} build --features "strict" - name: Build rust library (release) - run: cargo build --release --features "strict" + run: cargo +${{ matrix.rust-version }} build --release --features "strict" - name: Build rust library (release with mpi) - run: cargo build --release --features "strict,mpi" + run: cargo +${{ matrix.rust-version }} build --release --features "strict,mpi" - name: Run unit tests - run: cargo test --lib --features "strict" + run: cargo +${{ matrix.rust-version }} test --lib --features "strict" - name: Run unit tests (with mpi enabled) - run: cargo test --lib --features "mpi,strict" + run: cargo +${{ matrix.rust-version }} test --lib --features "mpi,strict" - name: Run tests - run: cargo test --examples --release --features "mpi,strict" + run: cargo +${{ matrix.rust-version }} test --examples --release --features "mpi,strict" - name: Run examples run: | - python3 find_examples.py + python3 find_examples.py --cargo-version ${{ matrix.rust-version }} chmod +x examples.sh ./examples.sh - - - name: Build docs - run: cargo doc --features "strict" --no-deps diff --git a/Cargo.toml b/Cargo.toml index 7bfcdc7b..200d5330 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -4,6 +4,7 @@ members = [ "bem", "element", "fmm", + "generation", "grid", "hyksort", "quadrature", @@ -17,6 +18,7 @@ default-members = [ "bem", "element", "fmm", + "generation", "grid", "hyksort", "quadrature", diff --git a/find_examples.py b/find_examples.py index c5f78f31..23d54a82 100644 --- a/find_examples.py +++ b/find_examples.py @@ -18,8 +18,12 @@ parser = argparse.ArgumentParser(description='Process some integers.') parser.add_argument('--features', default=None, help='feature flags to pass to the examples') +parser.add_argument('--cargo-version', default="stable", + help='cargo version to use (stable/beta/nightly') -raw_features = parser.parse_args().features +a = parser.parse_args() +raw_features = a.features +cargo_version = a.cargo_version features = [] if raw_features is not None: @@ -70,7 +74,7 @@ else: options += f" --features \"{','.join(features)}\"" - command = f"cargo {cmd} --example {example_name} --release" + command = f"cargo +{cargo_version} {cmd} --example {example_name} --release" if options is not None: command += f" {options}" if "{{NPROCESSES}}" in command: diff --git a/generation/Cargo.toml b/generation/Cargo.toml new file mode 100644 index 00000000..b6060d08 --- /dev/null +++ b/generation/Cargo.toml @@ -0,0 +1,28 @@ +[features] +strict = [] +slice-static-tables = [] + +[lib] +proc-macro = true + +[package] +name = "bempp-generation" +version = "0.1.0" +edition = "2021" + +# See more keys and their definitions at https://doc.rust-lang.org/cargo/reference/manifest.html + +[dependencies] +bempp-tools = { path = "../tools" } +bempp-traits = { path = "../traits" } +bempp-element = { path = "../element" } +bempp-bem = { path = "../bem" } +bempp-grid = { path = "../grid" } +bempp-quadrature = { path = "../quadrature" } +quote = { version = "1.0" } +syn = { version = "2.0.22" } +num = "0.4" +approx = "0.5" + +[build-dependencies] +rustversion = "1.0" diff --git a/generation/build.rs b/generation/build.rs new file mode 100644 index 00000000..5ee3c555 --- /dev/null +++ b/generation/build.rs @@ -0,0 +1,7 @@ +#[rustversion::not(nightly)] +fn main() {} + +#[rustversion::nightly] +fn main() { + println!("cargo:rustc-cfg=nightly"); +} diff --git a/generation/examples/kernels.rs b/generation/examples/kernels.rs new file mode 100644 index 00000000..4b17b1b8 --- /dev/null +++ b/generation/examples/kernels.rs @@ -0,0 +1,102 @@ +#![cfg_attr(nightly, feature(core_intrinsics))] + +use bempp_generation::generate_kernels; +use bempp_traits::element::FiniteElement; + +fn main() { + generate_kernels!( + dp0_dp0_kernel, + "Lagrange", + "Triangle", + 0, + true, + "Lagrange", + "Triangle", + 0, + true, + "Lagrange", + "Triangle", + 1, + false, + "Lagrange", + "Triangle", + 1, + false + ); + generate_kernels!( + p1_p1_kernel, + "Lagrange", + "Triangle", + 1, + false, + "Lagrange", + "Triangle", + 1, + false, + "Lagrange", + "Triangle", + 1, + false, + "Lagrange", + "Triangle", + 1, + false + ); + generate_kernels!( + p1_p2_kernel, + "Lagrange", + "Triangle", + 1, + false, + "Lagrange", + "Triangle", + 2, + false, + "Lagrange", + "Triangle", + 1, + false, + "Lagrange", + "Triangle", + 1, + false + ); + generate_kernels!( + high_order_geometry_kernel, + "Lagrange", + "Triangle", + 1, + false, + "Lagrange", + "Triangle", + 1, + false, + "Lagrange", + "Triangle", + 2, + false, + "Lagrange", + "Triangle", + 2, + false + ); + generate_kernels!( + mixed_order_geometry_kernel, + "Lagrange", + "Triangle", + 1, + false, + "Lagrange", + "Triangle", + 1, + false, + "Lagrange", + "Triangle", + 2, + false, + "Lagrange", + "Triangle", + 1, + false + ); +} diff --git a/generation/examples/single_layer_assembly.rs b/generation/examples/single_layer_assembly.rs new file mode 100644 index 00000000..8c8d2bfa --- /dev/null +++ b/generation/examples/single_layer_assembly.rs @@ -0,0 +1,213 @@ +#![cfg_attr(nightly, feature(core_intrinsics))] + +use approx::*; +use bempp_bem::function_space::SerialFunctionSpace; +use bempp_generation::generate_kernels; +use bempp_grid::shapes::regular_sphere; +use bempp_tools::arrays::Array2D; +use bempp_traits::arrays::Array2DAccess; +use bempp_traits::bem::{DofMap, FunctionSpace, Kernel}; +use bempp_traits::element::FiniteElement; +use bempp_traits::grid::{Geometry, Grid, Topology}; +use std::time::Instant; + +fn main() { + generate_kernels!( + bem_kernel, "Lagrange", "Triangle", 0, true, "Lagrange", "Triangle", 0, true, "Lagrange", + "Triangle", 1, false, "Lagrange", "Triangle", 1, false + ); + + let grid = regular_sphere(0); + let space = SerialFunctionSpace::new(&grid, &bem_kernel.test_element); + + let ndofs = space.dofmap().global_size(); + + let mut matrix = Array2D::::new((ndofs, ndofs)); + + assemble(&mut matrix, &space, &space, &bem_kernel); + + // Compare to result from bempp-cl + #[rustfmt::skip] + let from_cl = vec![vec![0.1854538822982487, 0.08755414595678074, 0.05963897421514472, 0.08755414595678074, 0.08755414595678074, 0.05963897421514473, 0.04670742127454548, 0.05963897421514472], vec![0.08755414595678074, 0.1854538822982487, 0.08755414595678074, 0.05963897421514472, 0.05963897421514472, 0.08755414595678074, 0.05963897421514473, 0.04670742127454548], vec![0.05963897421514472, 0.08755414595678074, 0.1854538822982487, 0.08755414595678074, 0.04670742127454548, 0.05963897421514472, 0.08755414595678074, 0.05963897421514473], vec![0.08755414595678074, 0.05963897421514472, 0.08755414595678074, 0.1854538822982487, 0.05963897421514473, 0.04670742127454548, 0.05963897421514472, 0.08755414595678074], vec![0.08755414595678074, 0.05963897421514472, 0.046707421274545476, 0.05963897421514473, 0.1854538822982487, 0.08755414595678074, 0.05963897421514472, 0.08755414595678074], vec![0.05963897421514473, 0.08755414595678074, 0.05963897421514472, 0.046707421274545476, 0.08755414595678074, 0.1854538822982487, 0.08755414595678074, 0.05963897421514472], vec![0.046707421274545476, 0.05963897421514473, 0.08755414595678074, 0.05963897421514472, 0.05963897421514472, 0.08755414595678074, 0.1854538822982487, 0.08755414595678074], vec![0.05963897421514472, 0.046707421274545476, 0.05963897421514473, 0.08755414595678074, 0.08755414595678074, 0.05963897421514472, 0.08755414595678074, 0.1854538822982487]]; + + for (i, row) in from_cl.iter().enumerate() { + for (j, entry) in row.iter().enumerate() { + assert_relative_eq!(*matrix.get(i, j).unwrap(), entry, epsilon = 1e-4); + } + } + + for i in 2..6 { + let grid = regular_sphere(i); + let space = SerialFunctionSpace::new(&grid, &bem_kernel.test_element); + + let ndofs = space.dofmap().global_size(); + + let mut matrix = Array2D::::new((ndofs, ndofs)); + let start = Instant::now(); + + assemble(&mut matrix, &space, &space, &bem_kernel); + println!( + "Time taken to assemble dense {} by {} matrix: {}ms", + grid.topology().entity_count(2), + grid.topology().entity_count(2), + start.elapsed().as_millis() + ); + } +} + +fn assemble<'a, E: FiniteElement>( + matrix: &mut impl Array2DAccess<'a, f64>, + test_space: &SerialFunctionSpace, + trial_space: &SerialFunctionSpace, + bem_kernel: &impl Kernel, +) { + let mut local_result = Array2D::::new(( + bem_kernel.test_element_dim(), + bem_kernel.trial_element_dim(), + )); + let grid = test_space.grid(); + + let mut test_vertices = Array2D::::new(( + bem_kernel.test_geometry_element_dim(), + grid.geometry().dim(), + )); + let mut trial_vertices = Array2D::::new(( + bem_kernel.trial_geometry_element_dim(), + grid.geometry().dim(), + )); + + // Test and trial cells are equal + for test_cell in 0..grid.geometry().cell_count() { + let test_cell_tindex = grid.topology().index_map()[test_cell]; + let test_dofs = test_space.dofmap().cell_dofs(test_cell_tindex).unwrap(); + let trial_dofs = trial_space.dofmap().cell_dofs(test_cell_tindex).unwrap(); + let test_cell_dofs = grid.geometry().cell_vertices(test_cell).unwrap(); + + for (i, dof) in test_cell_dofs.iter().enumerate() { + for (j, coord) in grid.geometry().point(*dof).unwrap().iter().enumerate() { + *test_vertices.get_mut(i, j).unwrap() = *coord; + } + } + bem_kernel.same_cell_kernel( + &mut local_result.data, + &test_vertices.data, + &test_vertices.data, + ); + for (test_i, test_dof) in test_dofs.iter().enumerate() { + for (trial_i, trial_dof) in trial_dofs.iter().enumerate() { + *matrix.get_mut(*test_dof, *trial_dof).unwrap() = + *local_result.get(test_i, trial_i).unwrap() + } + } + } + + for (test_cell, trial_cell, edge_info) in grid.topology().facet_adjacent_cells().iter() { + let test_cell_tindex = grid.topology().index_map()[*test_cell]; + let test_dofs = test_space.dofmap().cell_dofs(test_cell_tindex).unwrap(); + let test_cell_dofs = grid.geometry().cell_vertices(*test_cell).unwrap(); + + for i in 0..3 { + for j in 0..3 { + *test_vertices.get_mut(i, j).unwrap() = + grid.geometry().point(test_cell_dofs[i]).unwrap()[j]; + } + } + + let trial_cell_tindex = grid.topology().index_map()[*trial_cell]; + let trial_dofs = trial_space.dofmap().cell_dofs(trial_cell_tindex).unwrap(); + let trial_cell_dofs = grid.geometry().cell_vertices(*trial_cell).unwrap(); + + for i in 0..3 { + for j in 0..3 { + *trial_vertices.get_mut(i, j).unwrap() = + grid.geometry().point(trial_cell_dofs[i]).unwrap()[j]; + } + } + + bem_kernel.shared_edge_kernel( + &mut local_result.data, + &test_vertices.data, + &trial_vertices.data, + *edge_info, + ); + for (test_i, test_dof) in test_dofs.iter().enumerate() { + for (trial_i, trial_dof) in trial_dofs.iter().enumerate() { + *matrix.get_mut(*test_dof, *trial_dof).unwrap() = + *local_result.get(test_i, trial_i).unwrap() + } + } + } + + for (test_cell, trial_cell, vertex_info) in grid.topology().ridge_adjacent_cells().iter() { + let test_cell_tindex = grid.topology().index_map()[*test_cell]; + let test_dofs = test_space.dofmap().cell_dofs(test_cell_tindex).unwrap(); + let test_cell_dofs = grid.geometry().cell_vertices(*test_cell).unwrap(); + + for i in 0..3 { + for j in 0..3 { + *test_vertices.get_mut(i, j).unwrap() = + grid.geometry().point(test_cell_dofs[i]).unwrap()[j]; + } + } + + let trial_cell_tindex = grid.topology().index_map()[*trial_cell]; + let trial_dofs = trial_space.dofmap().cell_dofs(trial_cell_tindex).unwrap(); + let trial_cell_dofs = grid.geometry().cell_vertices(*trial_cell).unwrap(); + + for i in 0..3 { + for j in 0..3 { + *trial_vertices.get_mut(i, j).unwrap() = + grid.geometry().point(trial_cell_dofs[i]).unwrap()[j]; + } + } + + bem_kernel.shared_vertex_kernel( + &mut local_result.data, + &test_vertices.data, + &trial_vertices.data, + *vertex_info, + ); + for (test_i, test_dof) in test_dofs.iter().enumerate() { + for (trial_i, trial_dof) in trial_dofs.iter().enumerate() { + *matrix.get_mut(*test_dof, *trial_dof).unwrap() = + *local_result.get(test_i, trial_i).unwrap() + } + } + } + + for (test_cell, trial_cell) in grid.topology().nonadjacent_cells().iter() { + let test_cell_tindex = grid.topology().index_map()[*test_cell]; + let test_dofs = test_space.dofmap().cell_dofs(test_cell_tindex).unwrap(); + let test_cell_dofs = grid.geometry().cell_vertices(*test_cell).unwrap(); + + for i in 0..3 { + for j in 0..3 { + *test_vertices.get_mut(i, j).unwrap() = + grid.geometry().point(test_cell_dofs[i]).unwrap()[j]; + } + } + + let trial_cell_tindex = grid.topology().index_map()[*trial_cell]; + let trial_dofs = trial_space.dofmap().cell_dofs(trial_cell_tindex).unwrap(); + let trial_cell_dofs = grid.geometry().cell_vertices(*trial_cell).unwrap(); + + for i in 0..3 { + for j in 0..3 { + *trial_vertices.get_mut(i, j).unwrap() = + grid.geometry().point(trial_cell_dofs[i]).unwrap()[j]; + } + } + + bem_kernel.nonneighbour_kernel( + &mut local_result.data, + &test_vertices.data, + &trial_vertices.data, + ); + for (test_i, test_dof) in test_dofs.iter().enumerate() { + for (trial_i, trial_dof) in trial_dofs.iter().enumerate() { + *matrix.get_mut(*test_dof, *trial_dof).unwrap() = + *local_result.get(test_i, trial_i).unwrap() + } + } + } +} diff --git a/generation/src/lib.rs b/generation/src/lib.rs new file mode 100644 index 00000000..d1f098a5 --- /dev/null +++ b/generation/src/lib.rs @@ -0,0 +1,1683 @@ +//! Generation of boundary element kernels +#![cfg_attr(feature = "strict", deny(warnings))] + +// TODO: make a struct that stores all information about the kernel, then look up (eg) tables names from it + +extern crate proc_macro; +use bempp_element::element::create_element; +use bempp_quadrature::duffy::triangle::triangle_duffy; +use bempp_quadrature::simplex_rules::simplex_rule; +use bempp_quadrature::types::{ + CellToCellConnectivity, NumericalQuadratureDefinition, TestTrialNumericalQuadratureDefinition, +}; +use bempp_tools::arrays::{Array2D, Array4D}; +use bempp_traits::arrays::{Array2DAccess, Array4DAccess}; +use bempp_traits::cell::ReferenceCellType; +use bempp_traits::element::ElementFamily; +use bempp_traits::element::FiniteElement; +use num::traits::real::Real; +use num::Num; +use proc_macro::TokenStream; +use quote::quote; +use std::collections::HashMap; +use std::fmt::Debug; + +use syn::{parse::Parse, parse_macro_input, Expr, Token}; + +#[derive(Debug, PartialEq, Eq, Clone, Copy, Hash)] +#[repr(u8)] +enum FunctionType { + Test = 0, + Trial = 1, +} + +#[derive(Debug, PartialEq, Eq, Clone, Copy, Hash)] +#[repr(u8)] +enum ValueType { + Value = 0, + DX = 1, + DY = 2, +} + +struct GenerationInput { + kernel_name: Expr, + _comma0: Token![,], + test_element_family: Expr, + _comma1: Token![,], + test_element_cell: Expr, + _comma2: Token![,], + test_element_degree: Expr, + _comma3: Token![,], + test_element_discontinuous: Expr, + _comma4: Token![,], + trial_element_family: Expr, + _comma5: Token![,], + trial_element_cell: Expr, + _comma6: Token![,], + trial_element_degree: Expr, + _comma7: Token![,], + trial_element_discontinuous: Expr, + _comma8: Token![,], + test_geometry_element_family: Expr, + _comma9: Token![,], + test_geometry_element_cell: Expr, + _comma10: Token![,], + test_geometry_element_degree: Expr, + _comma11: Token![,], + test_geometry_element_discontinuous: Expr, + _comma12: Token![,], + trial_geometry_element_family: Expr, + _comma13: Token![,], + trial_geometry_element_cell: Expr, + _comma14: Token![,], + trial_geometry_element_degree: Expr, + _comma15: Token![,], + trial_geometry_element_discontinuous: Expr, +} + +impl Parse for GenerationInput { + fn parse(input: syn::parse::ParseStream) -> syn::Result { + Ok(Self { + kernel_name: input.parse()?, + _comma0: input.parse()?, + test_element_family: input.parse()?, + _comma1: input.parse()?, + test_element_cell: input.parse()?, + _comma2: input.parse()?, + test_element_degree: input.parse()?, + _comma3: input.parse()?, + test_element_discontinuous: input.parse()?, + _comma4: input.parse()?, + trial_element_family: input.parse()?, + _comma5: input.parse()?, + trial_element_cell: input.parse()?, + _comma6: input.parse()?, + trial_element_degree: input.parse()?, + _comma7: input.parse()?, + trial_element_discontinuous: input.parse()?, + _comma8: input.parse()?, + test_geometry_element_family: input.parse()?, + _comma9: input.parse()?, + test_geometry_element_cell: input.parse()?, + _comma10: input.parse()?, + test_geometry_element_degree: input.parse()?, + _comma11: input.parse()?, + test_geometry_element_discontinuous: input.parse()?, + _comma12: input.parse()?, + trial_geometry_element_family: input.parse()?, + _comma13: input.parse()?, + trial_geometry_element_cell: input.parse()?, + _comma14: input.parse()?, + trial_geometry_element_degree: input.parse()?, + _comma15: input.parse()?, + trial_geometry_element_discontinuous: input.parse()?, + }) + } +} + +fn parse_family(family: &Expr) -> ElementFamily { + let family_str = format!("{}", quote! { #family }); + if family_str == "\"Lagrange\"" { + ElementFamily::Lagrange + } else if family_str == "\"Raviart-Thomas\"" { + ElementFamily::RaviartThomas + } else { + panic!("Unsupported element family: {}", family_str); + } +} + +fn parse_cell(cell: &Expr) -> ReferenceCellType { + let cell_str = format!("{}", quote! { #cell }); + if cell_str == "\"Triangle\"" { + ReferenceCellType::Triangle + } else if cell_str == "\"Quadrilateral\"" { + ReferenceCellType::Quadrilateral + } else { + panic!("Unsupported cell type: {}", cell_str); + } +} + +fn parse_int(int: &Expr) -> usize { + format!("{}", quote! { #int }).parse().unwrap() +} + +fn parse_bool(bool: &Expr) -> bool { + format!("{}", quote! { #bool }).parse().unwrap() +} +fn parse_string(string: &Expr) -> String { + format!("{}", quote! { #string }) +} + +fn assert_sizes( + test_element: &impl FiniteElement, + trial_element: &impl FiniteElement, + test_geometry_element: &impl FiniteElement, + trial_geometry_element: &impl FiniteElement, + gdim: usize, + level: usize, +) -> String { + let mut code = String::new(); + code += &indent(level); + code += &format!( + "assert_eq!(result.len(), {});\n", + trial_element.dim() * test_element.dim() + ); + code += &indent(level); + code += &format!( + "assert_eq!(test_vertices.len(), {});\n", + test_geometry_element.dim() * gdim + ); + code += &indent(level); + code += &format!( + "assert_eq!(trial_vertices.len(), {});\n", + trial_geometry_element.dim() * gdim + ); + code +} + +fn table(table: &[T]) -> String { + let mut t = String::new(); + t += "["; + for (i, w) in table.iter().enumerate() { + if i > 0 { + t += ", "; + } + t += &format!("{w:?}"); + } + t += "]"; + t +} + +fn eval_table_transpose( + table: &Array4D, + deriv: usize, + component: usize, +) -> String { + let mut t = String::new(); + t += "["; + for j in 0..table.shape().2 { + for i in 0..table.shape().1 { + if i > 0 || j > 0 { + t += ", "; + } + t += &format!("{:?}", table.get(deriv, i, j, component).unwrap()); + } + } + t += "]"; + t +} + +fn tables( + test_element: &impl FiniteElement, + trial_element: &impl FiniteElement, + test_geometry_element: &impl FiniteElement, + trial_geometry_element: &impl FiniteElement, + test_points: &HashMap>, + trial_points: &HashMap>, + level: usize, +) -> String { + type T = f64; + let typename = std::any::type_name::().to_string(); + + let test_npts = if let Some((_qid, pts)) = test_points.iter().next() { + pts.shape().0 + } else { + panic!("No points."); + }; + let trial_npts = if let Some((_qid, pts)) = trial_points.iter().next() { + pts.shape().0 + } else { + panic!("No points."); + }; + + let test_derivs = if test_geometry_element.degree() == 1 { + 0 + } else { + 1 + }; + let trial_derivs = if trial_geometry_element.degree() == 1 { + 0 + } else { + 1 + }; + + let mut test_eval_table = Array4D::::new(test_element.tabulate_array_shape(0, test_npts)); + let mut test_geometry_eval_table = + Array4D::::new(test_geometry_element.tabulate_array_shape(test_derivs, test_npts)); + let mut trial_eval_table = Array4D::::new(trial_element.tabulate_array_shape(0, trial_npts)); + let mut trial_geometry_eval_table = + Array4D::::new(trial_geometry_element.tabulate_array_shape(trial_derivs, trial_npts)); + let mut out = String::new(); + if test_points.len() == 1 { + if let Some((_qid, pts)) = test_points.iter().next() { + test_element.tabulate(pts, 0, &mut test_eval_table); + test_geometry_element.tabulate(pts, test_derivs, &mut test_geometry_eval_table); + } + if let Some((_qid, pts)) = trial_points.iter().next() { + trial_element.tabulate(pts, 0, &mut trial_eval_table); + trial_geometry_element.tabulate(pts, trial_derivs, &mut trial_geometry_eval_table); + } + out += &indent(level); + out += &format!( + "let TEST_EVALS: [{typename}; {}] = ", + test_npts * test_element.dim() + ); + out += &eval_table_transpose(&test_eval_table, 0, 0); + out += ";\n"; + out += &indent(level); + out += &format!( + "let TEST_GEOMETRY_EVALS: [{typename}; {}] = ", + test_npts * test_geometry_element.dim() + ); + out += &eval_table_transpose(&test_geometry_eval_table, 0, 0); + out += ";\n"; + if test_derivs > 0 { + out += &indent(level); + out += &format!( + "let TEST_GEOMETRY_EVALS_DX: [{typename}; {}] = ", + test_npts * test_geometry_element.dim() + ); + out += &eval_table_transpose(&test_geometry_eval_table, 1, 0); + out += ";\n"; + out += &indent(level); + out += &format!( + "let TEST_GEOMETRY_EVALS_DY: [{typename}; {}] = ", + test_npts * test_geometry_element.dim() + ); + out += &eval_table_transpose(&test_geometry_eval_table, 2, 0); + out += ";\n"; + } + out += &indent(level); + out += &format!( + "let TRIAL_EVALS: [{typename}; {}] = ", + trial_npts * trial_element.dim() + ); + out += &eval_table_transpose(&trial_eval_table, 0, 0); + out += ";\n"; + out += &indent(level); + out += &format!( + "let TRIAL_GEOMETRY_EVALS: [{typename}; {}] = ", + trial_npts * trial_geometry_element.dim() + ); + out += &eval_table_transpose(&trial_geometry_eval_table, 0, 0); + out += ";\n"; + if trial_derivs > 0 { + out += &indent(level); + out += &format!( + "let TRIAL_GEOMETRY_EVALS_DX: [{typename}; {}] = ", + trial_npts * trial_geometry_element.dim() + ); + out += &eval_table_transpose(&trial_geometry_eval_table, 1, 0); + out += ";\n"; + out += &indent(level); + out += &format!( + "let TRIAL_GEOMETRY_EVALS_DY: [{typename}; {}] = ", + trial_npts * trial_geometry_element.dim() + ); + out += &eval_table_transpose(&trial_geometry_eval_table, 2, 0); + out += ";\n"; + } + } else { + out += &indent(level); + out += &format!( + "let TEST_EVALS: [{typename}; {}];\n", + test_npts * test_element.dim() + ); + out += &indent(level); + out += &format!( + "let TEST_GEOMETRY_EVALS: [{typename}; {}];\n", + test_npts * test_geometry_element.dim() + ); + if test_derivs > 0 { + out += &indent(level); + out += &format!( + "let TEST_GEOMETRY_EVALS_DX: [{typename}; {}];\n", + test_npts * test_geometry_element.dim() + ); + out += &indent(level); + out += &format!( + "let TEST_GEOMETRY_EVALS_DY: [{typename}; {}];\n", + test_npts * test_geometry_element.dim() + ); + } + out += &indent(level); + out += &format!( + "let TRIAL_EVALS: [{typename}; {}];\n", + trial_npts * trial_element.dim() + ); + out += &indent(level); + out += &format!( + "let TRIAL_GEOMETRY_EVALS: [{typename}; {}];\n", + trial_npts * trial_geometry_element.dim() + ); + if trial_derivs > 0 { + out += &indent(level); + out += &format!( + "let TRIAL_GEOMETRY_EVALS_DX: [{typename}; {}];\n", + trial_npts * trial_geometry_element.dim() + ); + out += &indent(level); + out += &format!( + "let TRIAL_GEOMETRY_EVALS_DY: [{typename}; {}];\n", + trial_npts * trial_geometry_element.dim() + ); + } + out += &indent(level); + out += "match qid {\n"; + for (qid, test_pts) in test_points.iter() { + let trial_pts = trial_points.get(qid).unwrap(); + test_element.tabulate(test_pts, test_derivs, &mut test_eval_table); + test_geometry_element.tabulate(test_pts, test_derivs, &mut test_geometry_eval_table); + trial_element.tabulate(trial_pts, trial_derivs, &mut trial_eval_table); + trial_geometry_element.tabulate( + trial_pts, + trial_derivs, + &mut trial_geometry_eval_table, + ); + out += &indent(level + 1); + out += &format!("{qid} => {{\n"); + out += &indent(level + 2); + out += "TEST_EVALS = "; + out += &eval_table_transpose(&test_eval_table, 0, 0); + out += ";\n"; + out += &indent(level + 2); + out += "TEST_GEOMETRY_EVALS = "; + out += &eval_table_transpose(&test_geometry_eval_table, 0, 0); + out += ";\n"; + if test_derivs > 0 { + out += &indent(level + 2); + out += "TEST_GEOMETRY_EVALS_DX = "; + out += &eval_table_transpose(&test_geometry_eval_table, 1, 0); + out += ";\n"; + out += &indent(level + 2); + out += "TEST_GEOMETRY_EVALS_DY = "; + out += &eval_table_transpose(&test_geometry_eval_table, 2, 0); + out += ";\n"; + } + out += &indent(level + 2); + out += "TRIAL_EVALS = "; + out += &eval_table_transpose(&trial_eval_table, 0, 0); + out += ";\n"; + out += &indent(level + 2); + out += "TRIAL_GEOMETRY_EVALS = "; + out += &eval_table_transpose(&trial_geometry_eval_table, 0, 0); + out += ";\n"; + if trial_derivs > 0 { + out += &indent(level + 2); + out += "TRIAL_GEOMETRY_EVALS_DX = "; + out += &eval_table_transpose(&trial_geometry_eval_table, 1, 0); + out += ";\n"; + out += &indent(level + 2); + out += "TRIAL_GEOMETRY_EVALS_DY = "; + out += &eval_table_transpose(&trial_geometry_eval_table, 2, 0); + out += ";\n"; + } + out += &indent(level + 1); + out += "},\n"; + } + out += &indent(level + 1); + out += "_ => { panic!(\"Unsupported quadrature rule.\"); },\n"; + out += &indent(level); + out += "}\n"; + } + + out +} + +fn set_zero(len: usize, level: usize) -> String { + let mut out = String::new(); + + if len == 1 { + out += &indent(level); + out += "result[0] = 0.0;\n"; + } else { + out += &indent(level); + out += &format!("for i in 0..{len} {{\n"); + out += &indent(level + 1); + out += "result[i] = 0.0;\n"; + out += &indent(level); + out += "}\n"; + } + out +} + +#[allow(clippy::too_many_arguments)] +#[cfg(feature = "slice-static-tables")] +fn take_slices( + test_element: &impl FiniteElement, + trial_element: &impl FiniteElement, + test_geometry_element: &impl FiniteElement, + trial_geometry_element: &impl FiniteElement, + test_npts: usize, + trial_npts: usize, + gdim: usize, + level: usize, +) -> String { + let mut out = String::new(); + for b in 0..test_geometry_element.dim() { + out += &indent(level); + out += &format!( + "let tsv{b} = &test_vertices[{}..{}];\n", + b * gdim, + (b + 1) * gdim + ); + } + for b in 0..test_geometry_element.dim() { + out += &indent(level); + out += &format!( + "let tsg{b} = &TEST_GEOMETRY_EVALS[{}..{}];\n", + b * test_npts, + (b + 1) * test_npts, + ); + } + if test_geometry_element.degree() > 1 { + for b in 0..test_geometry_element.dim() { + out += &indent(level); + out += &format!( + "let tsgx{b} = &TEST_GEOMETRY_EVALS_DX[{}..{}];\n", + b * test_npts, + (b + 1) * test_npts, + ); + out += &indent(level); + out += &format!( + "let tsgy{b} = &TEST_GEOMETRY_EVALS_DY[{}..{}];\n", + b * test_npts, + (b + 1) * test_npts, + ); + } + } + for b in 0..trial_geometry_element.dim() { + out += &indent(level); + out += &format!( + "let trv{b} = &trial_vertices[{}..{}];\n", + b * gdim, + (b + 1) * gdim + ); + } + for b in 0..trial_geometry_element.dim() { + out += &indent(level); + out += &format!( + "let trg{b} = &TRIAL_GEOMETRY_EVALS[{}..{}];\n", + b * trial_npts, + (b + 1) * trial_npts, + ); + } + if trial_geometry_element.degree() > 1 { + for b in 0..trial_geometry_element.dim() { + out += &indent(level); + out += &format!( + "let trgx{b} = &TRIAL_GEOMETRY_EVALS_DX[{}..{}];\n", + b * trial_npts, + (b + 1) * trial_npts, + ); + out += &indent(level); + out += &format!( + "let trgy{b} = &TRIAL_GEOMETRY_EVALS_DY[{}..{}];\n", + b * trial_npts, + (b + 1) * trial_npts, + ); + } + } + if test_element.degree() > 0 { + for b in 0..test_element.dim() { + out += &indent(level); + out += &format!( + "let ts{b} = &TEST_EVALS[{}..{}];\n", + b * test_npts, + (b + 1) * test_npts + ); + } + } + if trial_element.degree() > 0 { + for b in 0..trial_element.dim() { + out += &indent(level); + out += &format!( + "let tr{b} = &TRIAL_EVALS[{}..{}];\n", + b * trial_npts, + (b + 1) * trial_npts + ); + } + } + out +} + +#[allow(clippy::too_many_arguments)] +#[cfg(not(feature = "slice-static-tables"))] +fn take_slices( + _test_element: &impl FiniteElement, + _trial_element: &impl FiniteElement, + test_geometry_element: &impl FiniteElement, + trial_geometry_element: &impl FiniteElement, + _test_npts: usize, + _trial_npts: usize, + gdim: usize, + level: usize, +) -> String { + let mut out = String::new(); + for b in 0..test_geometry_element.dim() { + out += &indent(level); + out += &format!( + "let tsv{b} = &test_vertices[{}..{}];\n", + b * gdim, + (b + 1) * gdim + ); + } + for b in 0..trial_geometry_element.dim() { + out += &indent(level); + out += &format!( + "let trv{b} = &trial_vertices[{}..{}];\n", + b * gdim, + (b + 1) * gdim + ); + } + out +} + +trait Index { + fn is_int(&self) -> bool; + fn str(&self) -> String; + fn int(&self) -> usize; +} +struct IntIndex { + i: usize, +} +impl Index for IntIndex { + fn is_int(&self) -> bool { + true + } + fn str(&self) -> String { + self.i.to_string() + } + fn int(&self) -> usize { + self.i + } +} +struct StrIndex { + i: String, +} +impl Index for StrIndex { + fn is_int(&self) -> bool { + false + } + fn str(&self) -> String { + self.i.to_string() + } + fn int(&self) -> usize { + panic!("Not an integer."); + } +} +fn vertex(t: FunctionType, v: &impl Index, d: &impl Index, gdim: usize) -> String { + if v.is_int() { + format!( + "{}v{}[{}]", + match t { + FunctionType::Test => "ts", + FunctionType::Trial => "tr", + }, + v.str(), + d.str() + ) + } else { + format!( + "{}_vertices[{} * {gdim} + {}]", + match t { + FunctionType::Test => "test", + FunctionType::Trial => "trial", + }, + v.str(), + d.str() + ) + } +} + +#[cfg(features = "slice-static-tables")] +fn function_value( + t: FunctionType, + v: ValueType, + b: &impl Index, + p: &impl Index, + edim: usize, +) -> String { + if v.is_int() { + format!( + "{}{}{}[{}]", + match t { + FunctionType::Test => "ts", + FunctionType::Trial => "tr", + }, + match v { + ValueType::Value => "", + ValueType::DX => "x", + ValueType::DY => "y", + }, + b.str(), + p.str() + ) + } else { + if b.is_int() { + format!( + "{}_EVALS{}[{} + {}]", + match t { + FunctionType::Test => "TEST", + FunctionType::Trial => "TRIAL", + }, + match v { + ValueType::Value => "", + ValueType::DX => "_DX", + ValueType::DY => "_DY", + }, + b.int() * edim, + p.str() + ) + } else { + format!( + "{}_EVALS{}[{} * {edim} + {}]", + match t { + FunctionType::Test => "TEST", + FunctionType::Trial => "TRIAL", + }, + match v { + ValueType::Value => "", + ValueType::DX => "_DX", + ValueType::DY => "_DY", + }, + b.str(), + p.str() + ) + } + } +} +#[cfg(not(features = "slice-static-tables"))] +fn function_value( + t: FunctionType, + v: ValueType, + b: &impl Index, + p: &impl Index, + edim: usize, +) -> String { + if b.is_int() { + format!( + "{}_EVALS{}[{} + {}]", + match t { + FunctionType::Test => "TEST", + FunctionType::Trial => "TRIAL", + }, + match v { + ValueType::Value => "", + ValueType::DX => "_DX", + ValueType::DY => "_DY", + }, + b.int() * edim, + p.str() + ) + } else { + format!( + "{}_EVALS{}[{} * {edim} + {}]", + match t { + FunctionType::Test => "TEST", + FunctionType::Trial => "TRIAL", + }, + match v { + ValueType::Value => "", + ValueType::DX => "_DX", + ValueType::DY => "_DY", + }, + b.str(), + p.str() + ) + } +} + +#[cfg(features = "slice-static-tables")] +fn geometry(t: FunctionType, v: ValueType, b: &impl Index, p: &impl Index, npts: usize) -> String { + if v.is_int() { + format!( + "{}g{}{}[{}]", + match t { + FunctionType::Test => "ts", + FunctionType::Trial => "tr", + }, + match v { + ValueType::Value => "", + ValueType::DX => "x", + ValueType::DY => "y", + }, + b.str(), + p.str() + ) + } else { + if b.is_int() { + format!( + "{}_GEOMETRY_EVALS{}[{} + {}]", + match t { + FunctionType::Test => "TEST", + FunctionType::Trial => "TRIAL", + }, + match v { + ValueType::Value => "", + ValueType::DX => "_DX", + ValueType::DY => "_DY", + }, + b.int() * npts, + p.str() + ) + } else { + format!( + "{}_GEOMETRY_EVALS{}[{} * {npts} + {}]", + match t { + FunctionType::Test => "TEST", + FunctionType::Trial => "TRIAL", + }, + match v { + ValueType::Value => "", + ValueType::DX => "_DX", + ValueType::DY => "_DY", + }, + b.str(), + p.str() + ) + } + } +} +#[cfg(not(features = "slice-static-tables"))] +fn geometry(t: FunctionType, v: ValueType, b: &impl Index, p: &impl Index, npts: usize) -> String { + if b.is_int() { + format!( + "{}_GEOMETRY_EVALS{}[{} + {}]", + match t { + FunctionType::Test => "TEST", + FunctionType::Trial => "TRIAL", + }, + match v { + ValueType::Value => "", + ValueType::DX => "_DX", + ValueType::DY => "_DY", + }, + b.int() * npts, + p.str() + ) + } else { + format!( + "{}_GEOMETRY_EVALS{}[{} * {npts} + {}]", + match t { + FunctionType::Test => "TEST", + FunctionType::Trial => "TRIAL", + }, + match v { + ValueType::Value => "", + ValueType::DX => "_DX", + ValueType::DY => "_DY", + }, + b.str(), + p.str() + ) + } +} + +fn linear_derivative( + name: String, + vertices: &String, + gdim: usize, + table: &Array4D, + deriv: usize, +) -> String { + let mut code = String::new(); + code += &format!("let {name} = ["); + for d in 0..gdim { + if d > 0 { + code += ", "; + } + let mut started = false; + for f in 0..table.shape().2 { + let value = table.get(deriv, 0, f, 0).unwrap(); + if value.abs() > (T::one() + T::one()).powi(-15) { + if (*value - T::one()).abs() > (T::one() + T::one()).powi(-15) { + // value is 1 + if started { + code += " + "; + } + } else if (*value + T::one()).abs() > (T::one() + T::one()).powi(-15) { + // value is -1 + code += " - "; + } else { + if started { + code += " + "; + } + code += &format!("{:?} * ", value); + } + code += &format!("{vertices}{f}[{d}]"); + started = true; + } + } + } + code += "];\n"; + code +} + +fn linear_jacobian( + name: String, + vertices: String, + tdim: usize, + gdim: usize, + table: &Array4D, + level: usize, +) -> String { + let mut code = String::new(); + if tdim == 2 && gdim == 3 { + code += &indent(level); + code += &linear_derivative("dx".to_string(), &vertices, gdim, table, 1); + code += &indent(level); + code += &linear_derivative("dy".to_string(), &vertices, gdim, table, 2); + code += &indent(level); + code += "let j = [dx[1] * dy[2] - dx[2] * dy[1], dx[2] * dy[0] - dx[0] * dy[2], dx[0] * dy[1] - dx[1] * dy[0]];\n"; + code += &indent(level); + code += &format!("let {name} = (j[0].powi(2) + j[1].powi(2) + j[2].powi(2)).sqrt();\n"); + } else { + panic!( + "Jacobian computation not implemented for tdim {} and gdim {}.", + tdim, gdim + ); + } + code +} + +fn geometry_physical( + t: FunctionType, + v: ValueType, + point: &String, + gdim: usize, + basis_count: usize, +) -> String { + let mut code = String::new(); + code += "["; + for d in 0..gdim { + if d > 0 { + code += ", "; + } + for f in 0..basis_count { + if f > 0 { + code += " + "; + } + code += &geometry( + t, + v, + &IntIndex { i: f }, + &StrIndex { + i: point.to_string(), + }, + basis_count, + ); + code += " * "; + code += &vertex(t, &IntIndex { i: f }, &IntIndex { i: d }, gdim); + } + } + code += "]"; + code +} + +fn jacobian( + t: FunctionType, + geometry_element: &impl FiniteElement, + point: String, + gdim: usize, + tdim: usize, + level: usize, +) -> String { + let mut out = String::new(); + if geometry_element.degree() > 1 { + if tdim == 2 && gdim == 3 { + out += &indent(level); + out += "let dx = "; + out += &geometry_physical(t, ValueType::DX, &point, gdim, geometry_element.dim()); + out += ";\n"; + out += &indent(level); + out += "let dy = "; + out += &geometry_physical(t, ValueType::DY, &point, gdim, geometry_element.dim()); + out += ";\n"; + out += &indent(level); + out += "let j = [dx[1] * dy[2] - dx[2] * dy[1], dx[2] * dy[0] - dx[0] * dy[2], dx[0] * dy[1] - dx[1] * dy[0]];\n"; + out += &indent(level); + out += &format!( + "let {}_jdet = (j[0].powi(2) + j[1].powi(2) + j[2].powi(2)).sqrt();\n", + match t { + FunctionType::Test => "test", + FunctionType::Trial => "trial", + } + ); + } else { + panic!( + "Jacobian computation not implemented for tdim {} and gdim {}.", + tdim, gdim + ); + } + } + out +} + +#[cfg(not(nightly))] +fn fadd(a: String, b: String) -> String { + format!("{a} += {b};") +} + +#[cfg(nightly)] +fn fadd(a: String, b: String) -> String { + format!("{a} = unsafe {{ std::intrinsics::fadd_fast({a}, {b}) }};") +} + +#[allow(clippy::too_many_arguments)] +fn singular_kernel( + name: String, + quadrules: HashMap, + test_element: &impl FiniteElement, + trial_element: &impl FiniteElement, + test_geometry_element: &impl FiniteElement, + trial_geometry_element: &impl FiniteElement, + tdim: usize, + gdim: usize, +) -> String { + // TODO: make this into an input + type T = f64; + + let typename = std::any::type_name::().to_string(); + let npts = if let Some((_qid, qrule)) = quadrules.iter().next() { + qrule.npoints + } else { + panic!("No quadrature rule given."); + }; + + let mut code = String::new(); + // Function name and inputs + code += &format!("fn {name}(&self, result: &mut [{typename}], test_vertices: &[{typename}], trial_vertices: &[{typename}]"); + if quadrules.len() > 1 { + code += ", qid: u8"; + } + code += ") {\n"; + + // Assert that inputs have correct size + code += &assert_sizes( + test_element, + trial_element, + test_geometry_element, + trial_geometry_element, + gdim, + 1, + ); + + // Write quadrature weights + code += &indent(1); + code += &format!("let WTS: [{typename}; {npts}] = "); + if let Some((_qid, qrule)) = quadrules.iter().next() { + code += &table(&qrule.weights); + } else { + panic!("No quadrature rule given."); + } + code += ";\n"; + + // Write tabulated elements + let mut test_points = HashMap::new(); + let mut trial_points = HashMap::new(); + for (qid, qrule) in quadrules.iter() { + let mut ts = Array2D::::new((npts, tdim)); + let mut tr = Array2D::::new((npts, tdim)); + for p in 0..npts { + for d in 0..tdim { + *ts.get_mut(p, d).unwrap() = qrule.test_points[p * tdim + d]; + *tr.get_mut(p, d).unwrap() = qrule.trial_points[p * tdim + d]; + } + } + test_points.insert(*qid, ts); + trial_points.insert(*qid, tr); + } + code += &tables( + test_element, + trial_element, + test_geometry_element, + trial_geometry_element, + &test_points, + &trial_points, + 1, + ); + code += "\n"; + + // Set result to zero + code += &set_zero(test_element.dim() * trial_element.dim(), 1); + + // Write slices + code += &take_slices( + test_element, + trial_element, + test_geometry_element, + trial_geometry_element, + npts, + npts, + gdim, + 1, + ); + + // Quadrature loop + code += &indent(1); + code += &format!("for q in 0..{npts} {{\n"); + + // Compute the distance between points + code += &indent(2); + code += &format!("let mut sum_squares: {typename} = 0.0;\n"); + code += &indent(2); + code += &format!("for d in 0..{gdim} {{\n"); + code += &indent(3); + code += "let x = "; + for b in 0..test_geometry_element.dim() { + if b > 0 { + code += " + "; + } + code += &vertex( + FunctionType::Test, + &IntIndex { i: b }, + &StrIndex { i: "d".to_string() }, + gdim, + ); + code += " * "; + code += &geometry( + FunctionType::Test, + ValueType::Value, + &IntIndex { i: b }, + &StrIndex { i: "q".to_string() }, + npts, + ); + } + code += ";\n"; + code += &indent(3); + code += "let y = "; + for b in 0..trial_geometry_element.dim() { + if b > 0 { + code += " + "; + } + code += &vertex( + FunctionType::Trial, + &IntIndex { i: b }, + &StrIndex { i: "d".to_string() }, + gdim, + ); + code += " * "; + code += &geometry( + FunctionType::Trial, + ValueType::Value, + &IntIndex { i: b }, + &StrIndex { i: "q".to_string() }, + npts, + ); + } + code += ";\n"; + code += &indent(3); + code += &fadd("sum_squares".to_string(), "(x - y).powi(2)".to_string()); + code += "\n"; + code += &indent(2); + code += "}\n"; + code += &indent(2); + code += "let distance = sum_squares.sqrt();\n"; + + // Compute jacobians + code += &jacobian( + FunctionType::Test, + test_geometry_element, + "q".to_string(), + gdim, + tdim, + 2, + ); + code += &jacobian( + FunctionType::Trial, + trial_geometry_element, + "q".to_string(), + gdim, + tdim, + 2, + ); + + code += &indent(2); + code += "let c = WTS[q]"; + if test_geometry_element.degree() > 1 { + code += " * test_jdet"; + } + if trial_geometry_element.degree() > 1 { + code += " * trial_jdet"; + } + code += " / distance;\n"; + for i in 0..test_element.dim() { + for j in 0..trial_element.dim() { + let mut term = String::new(); + term += "c"; + if test_element.degree() > 0 { + term += " * "; + term += &function_value( + FunctionType::Test, + ValueType::Value, + &IntIndex { i }, + &StrIndex { i: "q".to_string() }, + test_element.dim(), + ); + } + if trial_element.degree() > 0 { + term += " * "; + term += &function_value( + FunctionType::Trial, + ValueType::Value, + &IntIndex { i }, + &StrIndex { i: "q".to_string() }, + trial_element.dim(), + ); + } + code += &indent(2); + code += &fadd(format!("result[{}]", i * trial_element.dim() + j), term); + code += "\n"; + } + } + + // End quadrature loop + code += &indent(1); + code += "}\n"; + + if test_geometry_element.degree() == 1 { + let pts = Array2D::from_data(vec![0.0; tdim], (1, tdim)); + let mut test_geometry_eval_table = + Array4D::::new(test_geometry_element.tabulate_array_shape(1, 1)); + test_geometry_element.tabulate(&pts, 1, &mut test_geometry_eval_table); + code += &linear_jacobian( + "test_jdet".to_string(), + "tsv".to_string(), + tdim, + gdim, + &test_geometry_eval_table, + 1, + ); + } + if trial_geometry_element.degree() == 1 { + let pts = Array2D::from_data(vec![0.0; tdim], (1, tdim)); + let mut trial_geometry_eval_table = + Array4D::::new(trial_geometry_element.tabulate_array_shape(1, 1)); + trial_geometry_element.tabulate(&pts, 1, &mut trial_geometry_eval_table); + code += &linear_jacobian( + "trial_jdet".to_string(), + "trv".to_string(), + tdim, + gdim, + &trial_geometry_eval_table, + 1, + ); + } + code += &indent(1); + code += "let multiplier: f64 = 0.25 * std::f64::consts::FRAC_1_PI"; + if test_geometry_element.degree() == 1 { + code += " * test_jdet"; + } + if trial_geometry_element.degree() == 1 { + code += " * trial_jdet"; + } + code += ";\n"; + + code += &indent(1); + code += &format!( + "for i in 0..{} {{\n", + test_element.dim() * trial_element.dim() + ); + code += &indent(2); + code += "result[i] *= multiplier;\n"; + code += &indent(1); + code += "}\n"; + + code += "}\n\n"; + // If you want to see the code that this macro generates, uncomment this line + // println!("{code}"); + code +} + +#[allow(clippy::too_many_arguments)] +fn nonsingular_kernel( + name: String, + test_rule: NumericalQuadratureDefinition, + trial_rule: NumericalQuadratureDefinition, + test_element: &impl FiniteElement, + trial_element: &impl FiniteElement, + test_geometry_element: &impl FiniteElement, + trial_geometry_element: &impl FiniteElement, + tdim: usize, + gdim: usize, +) -> String { + // TODO: make this into an input + type T = f64; + + let typename = std::any::type_name::().to_string(); + let test_npts = test_rule.npoints; + let trial_npts = trial_rule.npoints; + + let mut code = String::new(); + // Function name and inputs + code += &format!("fn {name}(&self, result: &mut [{typename}], test_vertices: &[{typename}], trial_vertices: &[{typename}]) {{\n"); + + // Assert that inputs have correct size + code += &assert_sizes( + test_element, + trial_element, + test_geometry_element, + trial_geometry_element, + gdim, + 1, + ); + + // Write quadrature weights + code += &indent(1); + code += &format!("let TEST_WTS: [{typename}; {test_npts}] = "); + code += &table(&test_rule.weights); + code += ";\n"; + code += &indent(1); + code += &format!("let TRIAL_WTS: [{typename}; {trial_npts}] = "); + code += &table(&trial_rule.weights); + code += ";\n"; + + // Write tabulated elements + let mut test_points = HashMap::new(); + let mut trial_points = HashMap::new(); + test_points.insert( + 0, + Array2D::::from_data(test_rule.points, (test_npts, tdim)), + ); + trial_points.insert( + 0, + Array2D::::from_data(trial_rule.points, (trial_npts, tdim)), + ); + + code += &tables( + test_element, + trial_element, + test_geometry_element, + trial_geometry_element, + &test_points, + &trial_points, + 1, + ); + + code += "\n"; + + // Set result to zero + code += &set_zero(test_element.dim() * trial_element.dim(), 1); + + // Write slices + code += &take_slices( + test_element, + trial_element, + test_geometry_element, + trial_geometry_element, + test_npts, + trial_npts, + gdim, + 1, + ); + + code += &indent(1); + code += &format!("let mut x = [0.0; {gdim}];\n"); + code += &indent(1); + code += &format!("for test_q in 0..{test_npts} {{\n"); + + if test_geometry_element.degree() > 1 { + code += &jacobian( + FunctionType::Test, + test_geometry_element, + "test_q".to_string(), + gdim, + tdim, + 2, + ); + } + code += &indent(2); + code += &format!("for d in 0..{gdim} {{\n"); + code += &indent(3); + code += "x[d] = "; + for b in 0..test_geometry_element.dim() { + if b > 0 { + code += " + "; + } + code += &vertex( + FunctionType::Test, + &IntIndex { i: b }, + &StrIndex { i: "d".to_string() }, + gdim, + ); + code += " * "; + code += &geometry( + FunctionType::Test, + ValueType::Value, + &IntIndex { i: b }, + &StrIndex { + i: "test_q".to_string(), + }, + test_npts, + ); + } + code += ";\n"; + code += &indent(2); + code += "}\n"; + code += &indent(2); + code += &format!("for trial_q in 0..{trial_npts} {{\n"); + code += &indent(3); + code += "let mut sum_squares = 0.0;\n"; + code += &indent(3); + code += &format!("for d in 0..{gdim} {{\n"); + code += &indent(4); + code += "let y = "; + for b in 0..trial_geometry_element.dim() { + if b > 0 { + code += " + "; + } + code += &vertex( + FunctionType::Trial, + &IntIndex { i: b }, + &StrIndex { i: "d".to_string() }, + gdim, + ); + code += " * "; + code += &geometry( + FunctionType::Trial, + ValueType::Value, + &IntIndex { i: b }, + &StrIndex { + i: "trial_q".to_string(), + }, + trial_npts, + ); + } + code += ";\n"; + code += &indent(4); + code += &fadd("sum_squares".to_string(), "(x[d] - y).powi(2)".to_string()); + code += "\n"; + code += &indent(3); + code += "}\n"; + code += &indent(3); + code += "let distance = sum_squares.sqrt();\n"; + if trial_geometry_element.degree() > 1 { + code += &jacobian( + FunctionType::Trial, + trial_geometry_element, + "trial_q".to_string(), + gdim, + tdim, + 2, + ); + } + code += &indent(3); + code += "let c = TEST_WTS[test_q] * TRIAL_WTS[trial_q]"; + if test_geometry_element.degree() > 1 { + code += " * test_jdet"; + } + if trial_geometry_element.degree() > 1 { + code += " * trial_jdet"; + } + code += " / distance;\n"; + for i in 0..test_element.dim() { + for j in 0..trial_element.dim() { + let mut term = String::new(); + term += "c"; + if test_element.degree() > 0 { + term += " * "; + term += &function_value( + FunctionType::Test, + ValueType::Value, + &IntIndex { i }, + &StrIndex { + i: "test_q".to_string(), + }, + test_element.dim(), + ); + } + if trial_element.degree() > 0 { + term += " * "; + term += &function_value( + FunctionType::Trial, + ValueType::Value, + &IntIndex { i }, + &StrIndex { + i: "trial_q".to_string(), + }, + trial_element.dim(), + ); + } + code += &indent(3); + code += &fadd(format!("result[{}]", i * trial_element.dim() + j), term); + code += "\n"; + } + } + code += &indent(2); + code += "}\n"; + code += &indent(1); + code += "}\n"; + + if test_geometry_element.degree() == 1 { + let pts = Array2D::from_data(vec![0.0; tdim], (1, tdim)); + let mut test_geometry_eval_table = + Array4D::::new(test_geometry_element.tabulate_array_shape(1, 1)); + test_geometry_element.tabulate(&pts, 1, &mut test_geometry_eval_table); + code += &linear_jacobian( + "test_jdet".to_string(), + "tsv".to_string(), + tdim, + gdim, + &test_geometry_eval_table, + 1, + ); + } + if trial_geometry_element.degree() == 1 { + let pts = Array2D::from_data(vec![0.0; tdim], (1, tdim)); + let mut trial_geometry_eval_table = + Array4D::::new(trial_geometry_element.tabulate_array_shape(1, 1)); + trial_geometry_element.tabulate(&pts, 1, &mut trial_geometry_eval_table); + code += &linear_jacobian( + "trial_jdet".to_string(), + "trv".to_string(), + tdim, + gdim, + &trial_geometry_eval_table, + 1, + ); + } + + code += &indent(1); + code += "let multiplier: f64 = 0.25 * std::f64::consts::FRAC_1_PI"; + if test_geometry_element.degree() == 1 { + code += " * test_jdet"; + } + if trial_geometry_element.degree() == 1 { + code += " * trial_jdet"; + } + code += ";\n"; + + code += &indent(1); + code += &format!( + "for i in 0..{} {{\n", + test_element.dim() * trial_element.dim() + ); + code += &indent(2); + code += "result[i] *= multiplier;\n"; + code += &indent(1); + code += "}\n"; + + code += "}\n\n"; + // If you want to see the code that this macro generates, uncomment this line + // println!("{code}"); + code +} + +fn indent(level: usize) -> String { + let mut t = String::new(); + for _ in 0..level * 4 { + t += " "; + } + t +} + +#[proc_macro] +pub fn generate_kernels(input: TokenStream) -> TokenStream { + // TODO: make these into inputs + type T = f64; + let tdim = 2; + let gdim = 3; + let singular_quad_degree = 3; + let nonsingular_quad_degree = 7; + + let typename = std::any::type_name::().to_string(); + let es = parse_macro_input!(input as GenerationInput); + + let kernel_name = parse_string(&es.kernel_name); + + let test_element = create_element( + parse_family(&es.test_element_family), + parse_cell(&es.test_element_cell), + parse_int(&es.test_element_degree), + parse_bool(&es.test_element_discontinuous), + ); + let trial_element = create_element( + parse_family(&es.trial_element_family), + parse_cell(&es.trial_element_cell), + parse_int(&es.trial_element_degree), + parse_bool(&es.trial_element_discontinuous), + ); + let test_geometry_element = create_element( + parse_family(&es.test_geometry_element_family), + parse_cell(&es.test_geometry_element_cell), + parse_int(&es.test_geometry_element_degree), + parse_bool(&es.test_geometry_element_discontinuous), + ); + let trial_geometry_element = create_element( + parse_family(&es.trial_geometry_element_family), + parse_cell(&es.trial_geometry_element_cell), + parse_int(&es.trial_geometry_element_degree), + parse_bool(&es.trial_geometry_element_discontinuous), + ); + + let mut code = String::new(); + + code += &format!("struct _BemppKernel_{kernel_name} {{\n"); + code += " test_element: bempp_element::element::CiarletElement,\n"; + code += " trial_element: bempp_element::element::CiarletElement,\n"; + code += " test_geometry_element: bempp_element::element::CiarletElement,\n"; + code += " trial_geometry_element: bempp_element::element::CiarletElement,\n"; + code += "};\n\n"; + code += + &format!("impl bempp_traits::bem::Kernel<{typename}> for _BemppKernel_{kernel_name} {{"); + code += "\n\n"; + code += "fn test_element_dim(&self) -> usize { self.test_element.dim() }\n"; + code += "fn trial_element_dim(&self) -> usize { self.trial_element.dim() }\n"; + code += "fn test_geometry_element_dim(&self) -> usize { self.test_geometry_element.dim() }\n"; + code += "fn trial_geometry_element_dim(&self) -> usize { self.trial_geometry_element.dim() }\n"; + + // TODO: quads + let mut quadrules = HashMap::new(); + quadrules.insert( + 0, + triangle_duffy( + &CellToCellConnectivity { + connectivity_dimension: 2, + local_indices: vec![(0, 0), (1, 1), (2, 2)], + }, + singular_quad_degree, + ) + .unwrap(), + ); + code += &singular_kernel( + "same_cell_kernel".to_string(), + quadrules, + &test_element, + &trial_element, + &test_geometry_element, + &trial_geometry_element, + tdim, + gdim, + ); + + let mut quadrules = HashMap::new(); + for i in 0..3 { + for j in i + 1..3 { + for k in 0..3 { + for l in 0..3 { + if k != l { + quadrules.insert( + i * 64 + k * 16 + j * 4 + l, + triangle_duffy( + &CellToCellConnectivity { + connectivity_dimension: 1, + local_indices: vec![(i, k), (j, l)], + }, + singular_quad_degree, + ) + .unwrap(), + ); + } + } + } + } + } + code += &singular_kernel( + "shared_edge_kernel".to_string(), + quadrules, + &test_element, + &trial_element, + &test_geometry_element, + &trial_geometry_element, + tdim, + gdim, + ); + + let mut quadrules = HashMap::new(); + for i in 0..3 { + for j in 0..3 { + quadrules.insert( + 4 * j + i, + triangle_duffy( + &CellToCellConnectivity { + connectivity_dimension: 0, + local_indices: vec![(i, j)], + }, + singular_quad_degree, + ) + .unwrap(), + ); + } + } + code += &singular_kernel( + "shared_vertex_kernel".to_string(), + quadrules, + &test_element, + &trial_element, + &test_geometry_element, + &trial_geometry_element, + tdim, + gdim, + ); + + let test_rule = simplex_rule(test_element.cell_type(), nonsingular_quad_degree).unwrap(); + let trial_rule = simplex_rule(trial_element.cell_type(), nonsingular_quad_degree).unwrap(); + + code += &nonsingular_kernel( + "nonneighbour_kernel".to_string(), + test_rule, + trial_rule, + &test_element, + &trial_element, + &test_geometry_element, + &trial_geometry_element, + tdim, + gdim, + ); + + code += "}\n\n"; + code += &format!("let {kernel_name} = _BemppKernel_{kernel_name} {{\n"); + code += &format!(" test_element: bempp_element::element::create_element(bempp_traits::element::ElementFamily::{:?}, bempp_traits::cell::ReferenceCellType::{:?}, {}, {}),\n", + test_element.family(), test_element.cell_type(), test_element.degree(), test_element.discontinuous()); + code += &format!(" trial_element: bempp_element::element::create_element(bempp_traits::element::ElementFamily::{:?}, bempp_traits::cell::ReferenceCellType::{:?}, {}, {}),\n", + trial_element.family(), trial_element.cell_type(), trial_element.degree(), trial_element.discontinuous()); + code += &format!(" test_geometry_element: bempp_element::element::create_element(bempp_traits::element::ElementFamily::{:?}, bempp_traits::cell::ReferenceCellType::{:?}, {}, {}),\n", + test_geometry_element.family(), test_geometry_element.cell_type(), test_geometry_element.degree(), test_geometry_element.discontinuous()); + code += &format!(" trial_geometry_element: bempp_element::element::create_element(bempp_traits::element::ElementFamily::{:?}, bempp_traits::cell::ReferenceCellType::{:?}, {}, {}),\n", + trial_geometry_element.family(), trial_geometry_element.cell_type(), trial_geometry_element.degree(), trial_geometry_element.discontinuous()); + code += "};"; + // If you want to print the generated code, uncomment this line + // println!("{code}"); + code.parse().unwrap() +} diff --git a/tools/src/arrays.rs b/tools/src/arrays.rs index 3e13efe0..23f422bc 100644 --- a/tools/src/arrays.rs +++ b/tools/src/arrays.rs @@ -6,7 +6,7 @@ use std::clone::Clone; /// A two-dimensional rectangular array pub struct Array2D { /// The data in the array, in row-major order - data: Vec, + pub data: Vec, /// The shape of the array shape: (usize, usize), } diff --git a/traits/src/bem.rs b/traits/src/bem.rs index 985d7c4d..88f384d9 100644 --- a/traits/src/bem.rs +++ b/traits/src/bem.rs @@ -1,5 +1,6 @@ use crate::element::FiniteElement; use crate::grid::Grid; +use num::Num; pub trait DofMap { /// Get the DOF numbers on the local process associated with the given entity @@ -40,3 +41,26 @@ pub trait FunctionSpace<'a> { self.dofmap().is_serial() && self.grid().is_serial() } } + +pub trait Kernel { + fn test_element_dim(&self) -> usize; + fn trial_element_dim(&self) -> usize; + fn test_geometry_element_dim(&self) -> usize; + fn trial_geometry_element_dim(&self) -> usize; + fn same_cell_kernel(&self, result: &mut [T], test_vertices: &[T], trial_vertices: &[T]); + fn shared_edge_kernel( + &self, + result: &mut [T], + test_vertices: &[T], + trial_vertices: &[T], + qid: u8, + ); + fn shared_vertex_kernel( + &self, + result: &mut [T], + test_vertices: &[T], + trial_vertices: &[T], + qid: u8, + ); + fn nonneighbour_kernel(&self, result: &mut [T], test_vertices: &[T], trial_vertices: &[T]); +}