From d1ad13c2e01f7edb7e745dbf7a83c70ea1cb51d4 Mon Sep 17 00:00:00 2001 From: Timo Betcke Date: Wed, 1 Jan 2025 19:21:53 +0100 Subject: [PATCH] WIP: space_to_point_map --- examples/test_green_evaluators.rs | 10 +- src/evaluator_tools.rs | 158 ++++++++++++++++++++++++++++++ src/lib.rs | 1 + 3 files changed, 165 insertions(+), 4 deletions(-) create mode 100644 src/evaluator_tools.rs diff --git a/examples/test_green_evaluators.rs b/examples/test_green_evaluators.rs index 6b62e0cd..fa6f04f3 100644 --- a/examples/test_green_evaluators.rs +++ b/examples/test_green_evaluators.rs @@ -45,10 +45,12 @@ fn main() { let mut charges = space.zero(); - charges - .view_mut() - .local_mut() - .fill_from_equally_distributed(&mut rng); + // charges + // .view_mut() + // .local_mut() + // .fill_from_equally_distributed(&mut rng); + + charges.view_mut().local_mut().set_one(); // Create the dense evaluator. diff --git a/src/evaluator_tools.rs b/src/evaluator_tools.rs new file mode 100644 index 00000000..75c9ddd1 --- /dev/null +++ b/src/evaluator_tools.rs @@ -0,0 +1,158 @@ +//! Various helper functions to support evaluators. + +use itertools::izip; +use mpi::traits::{Communicator, Equivalence}; +use ndelement::traits::FiniteElement; +use ndgrid::{ + traits::{Entity, GeometryMap, Grid}, + types::Ownership, +}; +use rlst::{ + rlst_array_from_slice2, rlst_dynamic_array1, rlst_dynamic_array2, rlst_dynamic_array3, + rlst_dynamic_array4, DefaultIterator, DistributedCsrMatrix, IndexLayout, RawAccess, + RawAccessMut, RlstScalar, +}; + +use crate::function::FunctionSpaceTrait; + +/// Create a linear operator from the map of a basis to points. +pub fn basis_to_point_map< + 'a, + C: Communicator + 'a, + DomainLayout: IndexLayout, + RangeLayout: IndexLayout, + T: RlstScalar + Equivalence, + Space: FunctionSpaceTrait, +>( + function_space: &Space, + domain_layout: &'a DomainLayout, + range_layout: &'a RangeLayout, + quadrature_points: &[T::Real], + quadrature_weights: &[T::Real], + return_transpose: bool, +) -> DistributedCsrMatrix<'a, DomainLayout, RangeLayout, T, C> +where + T::Real: Equivalence, +{ + // Get the grid. + let grid = function_space.grid(); + + // Topological dimension of the grid. + let tdim = grid.topology_dim(); + + // The method is currently restricted to single element type grids. + // So let's make sure that this is the case. + + assert_eq!(grid.entity_types(tdim).len(), 1); + + let reference_cell = grid.entity_types(tdim)[0]; + + // Number of cells. We are only interested in owned cells. + let n_cells = grid + .entity_iter(tdim) + .filter(|e| matches!(e.ownership(), Ownership::Owned)) + .count(); + + // Get the number of quadrature points and check that weights + // and points have compatible dimensions. + + let n_quadrature_points = quadrature_weights.len(); + assert_eq!(quadrature_points.len() % tdim, 0); + assert_eq!(quadrature_points.len() / tdim, n_quadrature_points); + + // Check that domain space and function space are compatible. + + let n_domain_dofs = domain_layout.number_of_global_indices(); + let n_range_dofs = range_layout.number_of_global_indices(); + + assert_eq!(function_space.local_size(), n_domain_dofs,); + + assert_eq!(n_cells * n_quadrature_points, n_range_dofs); + + // All the dimensions are OK. Let's get to work. We need to iterate through the elements, + // get the attached global dofs and the corresponding jacobian map. + + // Let's first tabulate the basis function values at the quadrature points on the reference element. + + // Quadrature point is needed here as a RLST matrix. + + let quadrature_points = rlst_array_from_slice2!(quadrature_points, [tdim, n_quadrature_points]); + + let mut table = rlst_dynamic_array4!( + T, + function_space + .element(reference_cell) + .tabulate_array_shape(0, n_quadrature_points) + ); + function_space + .element(reference_cell) + .tabulate(&quadrature_points, 0, &mut table); + + // We have tabulated the basis functions on the reference element. Now need + // the map to physical elements. + + let geometry_evaluator = grid.geometry_map(reference_cell, quadrature_points.data()); + + // The following arrays hold the jacobians, their determinants and the normals. + + let mut jacobians = + rlst_dynamic_array3![T::Real, [grid.geometry_dim(), tdim, n_quadrature_points]]; + let mut jdets = rlst_dynamic_array1![T::Real, [n_quadrature_points]]; + let mut normals = rlst_dynamic_array2![T::Real, [grid.geometry_dim(), n_quadrature_points]]; + + // Now iterate through the cells of the grid, get the attached dofs and evaluate the geometry map. + + // These arrays store the data of the transformation matrix. + let mut rows = Vec::::default(); + let mut cols = Vec::::default(); + let mut data = Vec::::default(); + + for cell in grid + .entity_iter(tdim) + .filter(|e| matches!(e.ownership(), Ownership::Owned)) + { + // Get the Jacobians + geometry_evaluator.jacobians_dets_normals( + cell.local_index(), + jacobians.data_mut(), + jdets.data_mut(), + normals.data_mut(), + ); + // Get the global dofs of the cell. + let global_dofs = function_space + .cell_dofs(cell.local_index()) + .unwrap() + .iter() + .map(|local_dof_index| function_space.global_dof_index(*local_dof_index)) + .collect::>(); + + for (qindex, (jdet, qweight)) in izip!(jdets.iter(), quadrature_weights.iter()).enumerate() + { + for global_dof in global_dofs.iter() { + rows.push(n_quadrature_points * cell.global_index() + qindex); + cols.push(*global_dof); + data.push(T::from_real(jdet * *qweight)); + } + } + } + + if return_transpose { + DistributedCsrMatrix::from_aij( + &domain_layout, + &range_layout, + &cols, + &rows, + &data, + domain_layout.comm(), + ) + } else { + DistributedCsrMatrix::from_aij( + &domain_layout, + &range_layout, + &rows, + &cols, + &data, + domain_layout.comm(), + ) + } +} diff --git a/src/lib.rs b/src/lib.rs index 74cc3080..ffcb8094 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -5,6 +5,7 @@ //pub mod bindings; pub mod boundary_assemblers; pub mod boundary_evaluators; +pub mod evaluator_tools; pub mod function; pub mod greens_function_evaluators; pub mod helmholtz;