Skip to content

Commit

Permalink
Eq, cargo fmt
Browse files Browse the repository at this point in the history
  • Loading branch information
mscroggs committed Aug 23, 2023
1 parent eacbf9c commit 57a8ef3
Show file tree
Hide file tree
Showing 6 changed files with 82 additions and 38 deletions.
19 changes: 9 additions & 10 deletions element/src/cell.rs
Original file line number Diff line number Diff line change
Expand Up @@ -10,17 +10,16 @@ pub mod cells_3d;
pub use bempp_traits::cell::{InvalidConnectivity, PhysicalCell, ReferenceCell, ReferenceCellType};
pub use cells_3d::*;

pub fn create_cell(cell_type: ReferenceCellType) -> Box<dyn ReferenceCell>
{
pub fn create_cell(cell_type: ReferenceCellType) -> Box<dyn ReferenceCell> {
match cell_type {
ReferenceCellType::Point => Box::new(Point{}),
ReferenceCellType::Interval => Box::new(Interval{}),
ReferenceCellType::Triangle => Box::new(Triangle{}),
ReferenceCellType::Quadrilateral => Box::new(Quadrilateral{}),
ReferenceCellType::Tetrahedron => Box::new(Tetrahedron{}),
ReferenceCellType::Hexahedron => Box::new(Hexahedron{}),
ReferenceCellType::Prism => Box::new(Prism{}),
ReferenceCellType::Pyramid => Box::new(Pyramid{}),
ReferenceCellType::Point => Box::new(Point {}),
ReferenceCellType::Interval => Box::new(Interval {}),
ReferenceCellType::Triangle => Box::new(Triangle {}),
ReferenceCellType::Quadrilateral => Box::new(Quadrilateral {}),
ReferenceCellType::Tetrahedron => Box::new(Tetrahedron {}),
ReferenceCellType::Hexahedron => Box::new(Hexahedron {}),
ReferenceCellType::Prism => Box::new(Prism {}),
ReferenceCellType::Pyramid => Box::new(Pyramid {}),
}
}

Expand Down
67 changes: 50 additions & 17 deletions element/src/element.rs
Original file line number Diff line number Diff line change
@@ -1,14 +1,14 @@
//! Finite Element definitions

use crate::cell::create_cell;
use crate::polynomials::{legendre_shape, polynomial_count, tabulate_legendre_polynomials};
use bempp_tools::arrays::{AdjacencyList, Array2D, Array3D};
use bempp_traits::arrays::{AdjacencyListAccess, Array2DAccess, Array3DAccess, Array4DAccess};
use bempp_traits::cell::ReferenceCellType;
use bempp_traits::element::{ElementFamily, FiniteElement, MapType};
use crate::cell::create_cell;
use crate::polynomials::{polynomial_count, tabulate_legendre_polynomials, legendre_shape};
use rlst_dense::RandomAccessMut;
use rlst_algorithms::linalg::LinAlg;
use rlst_algorithms::traits::inverse::Inverse;
use rlst_dense::RandomAccessMut;
pub mod lagrange;
pub mod raviart_thomas;

Expand Down Expand Up @@ -39,7 +39,6 @@ impl CiarletElement {
discontinuous: bool,
highest_degree: usize,
) -> CiarletElement {

let mut dim = 0;
for epts in &x {
for pts in epts {
Expand Down Expand Up @@ -79,7 +78,9 @@ impl CiarletElement {
}
new_x[tdim].push(all_pts);
new_x
} else { x };
} else {
x
};
let new_m = if discontinuous {
let mut new_m = [vec![], vec![], vec![], vec![]];
let mut pn = 0;
Expand All @@ -91,7 +92,8 @@ impl CiarletElement {
for j in 0..mat.shape().0 {
for k in 0..value_size {
for l in 0..mat.shape().2 {
*all_mat.get_mut(dn + j, k, pn + l).unwrap() = *mat.get(j, k, l).unwrap();
*all_mat.get_mut(dn + j, k, pn + l).unwrap() =
*mat.get(j, k, l).unwrap();
}
}
}
Expand All @@ -101,7 +103,9 @@ impl CiarletElement {
}
new_m[tdim].push(all_mat);
new_m
} else { m };
} else {
m
};

// Compute the dual matrix
let pdim = polynomial_count(cell_type, highest_degree);
Expand All @@ -120,7 +124,8 @@ impl CiarletElement {
let mut value = d_matrix.get_mut(j, l, dof + i).unwrap();
*value = 0.0;
for k in 0..pts.shape().0 {
*value += *mat.get(i, j, k).unwrap() * *table.get(0, l, k).unwrap()
*value +=
*mat.get(i, j, k).unwrap() * *table.get(0, l, k).unwrap()
}
}
}
Expand Down Expand Up @@ -153,21 +158,37 @@ impl CiarletElement {
println!();
}

let mut entity_dofs = [AdjacencyList::<usize>::new(), AdjacencyList::<usize>::new(), AdjacencyList::<usize>::new(), AdjacencyList::<usize>::new()];
let mut entity_dofs = [
AdjacencyList::<usize>::new(),
AdjacencyList::<usize>::new(),
AdjacencyList::<usize>::new(),
AdjacencyList::<usize>::new(),
];
let mut dof = 0;
let coefficients = Array3D::<f64>::new((0,0,0));
let coefficients = Array3D::<f64>::new((0, 0, 0));
for i in 0..4 {
for pts in &new_x[i] {
let dofs: Vec<usize> =(dof..dof + pts.shape().0).collect();
let dofs: Vec<usize> = (dof..dof + pts.shape().0).collect();
entity_dofs[i].add_row(&dofs);
dof += pts.shape().0;
}
}
CiarletElement { cell_type, degree, highest_degree, map_type, value_shape, value_size, family, discontinuous, dim, coefficients, entity_dofs }
CiarletElement {
cell_type,
degree,
highest_degree,
map_type,
value_shape,
value_size,
family,
discontinuous,
dim,
coefficients,
entity_dofs,
}
}
}


impl FiniteElement for CiarletElement {
fn value_size(&self) -> usize {
self.value_size
Expand Down Expand Up @@ -200,14 +221,27 @@ impl FiniteElement for CiarletElement {
nderivs: usize,
data: &mut impl Array4DAccess<f64>,
) {
let mut table = Array3D::<f64>::new(legendre_shape(self.cell_type, points, self.highest_degree, nderivs));
tabulate_legendre_polynomials(self.cell_type, points, self.highest_degree, nderivs, &mut table);
let mut table = Array3D::<f64>::new(legendre_shape(
self.cell_type,
points,
self.highest_degree,
nderivs,
));
tabulate_legendre_polynomials(
self.cell_type,
points,
self.highest_degree,
nderivs,
&mut table,
);
for d in 0..table.shape().0 {
for i in 0..table.shape().1 {
for p in 0..points.shape().0 {
for j in 0..self.value_size {
for b in 0..self.dim {
*data.get_mut(d, p, b, j).unwrap() += *self.coefficients.get(b, j, p).unwrap() * *table.get_mut(d, i, p).unwrap();
*data.get_mut(d, p, b, j).unwrap() +=
*self.coefficients.get(b, j, p).unwrap()
* *table.get_mut(d, i, p).unwrap();
}
}
}
Expand All @@ -219,7 +253,6 @@ impl FiniteElement for CiarletElement {
}
}


pub struct OldCiarletElement {
cell_type: ReferenceCellType,
degree: usize,
Expand Down
6 changes: 5 additions & 1 deletion element/src/element/lagrange.rs
Original file line number Diff line number Diff line change
Expand Up @@ -7,7 +7,11 @@ use bempp_traits::cell::ReferenceCellType;
use bempp_traits::element::{ElementFamily, MapType};

/// Create a Lagrange element
pub fn create(cell_type: ReferenceCellType, degree: usize, discontinuous: bool) -> OldCiarletElement {
pub fn create(
cell_type: ReferenceCellType,
degree: usize,
discontinuous: bool,
) -> OldCiarletElement {
if degree == 0 && !discontinuous {
panic!("Cannot create continuous degree 0 element");
}
Expand Down
6 changes: 5 additions & 1 deletion element/src/element/raviart_thomas.rs
Original file line number Diff line number Diff line change
Expand Up @@ -7,7 +7,11 @@ use bempp_traits::cell::ReferenceCellType;
use bempp_traits::element::{ElementFamily, MapType};

/// Create a Raviart-Thomas element
pub fn create(cell_type: ReferenceCellType, degree: usize, discontinuous: bool) -> OldCiarletElement {
pub fn create(
cell_type: ReferenceCellType,
degree: usize,
discontinuous: bool,
) -> OldCiarletElement {
let coefficients = match cell_type {
ReferenceCellType::Triangle => match degree {
// Basis = {(-x, -y), (x-1,y), (-x, 1-y)}
Expand Down
16 changes: 10 additions & 6 deletions element/src/polynomials.rs
Original file line number Diff line number Diff line change
Expand Up @@ -352,10 +352,10 @@ fn tabulate_legendre_polynomials_triangle<'a>(
}
}

pub fn polynomial_count(cell_type: ReferenceCellType,degree: usize) -> usize {
pub fn polynomial_count(cell_type: ReferenceCellType, degree: usize) -> usize {
match cell_type {
ReferenceCellType::Interval => degree + 1 ,
ReferenceCellType::Triangle => (degree + 1) * (degree + 2) / 2,
ReferenceCellType::Interval => degree + 1,
ReferenceCellType::Triangle => (degree + 1) * (degree + 2) / 2,
ReferenceCellType::Quadrilateral => (degree + 1) * (degree + 1),
_ => {
panic!("Unsupported cell type");
Expand All @@ -365,8 +365,8 @@ pub fn polynomial_count(cell_type: ReferenceCellType,degree: usize) -> usize {

pub fn derivative_count(cell_type: ReferenceCellType, derivatives: usize) -> usize {
match cell_type {
ReferenceCellType::Interval => derivatives + 1 ,
ReferenceCellType::Triangle => (derivatives + 1) * (derivatives + 2) / 2,
ReferenceCellType::Interval => derivatives + 1,
ReferenceCellType::Triangle => (derivatives + 1) * (derivatives + 2) / 2,
ReferenceCellType::Quadrilateral => (derivatives + 1) * (derivatives + 2) / 2,
_ => {
panic!("Unsupported cell type");
Expand All @@ -380,7 +380,11 @@ pub fn legendre_shape<'a>(
degree: usize,
derivatives: usize,
) -> (usize, usize, usize) {
(derivative_count(cell_type, derivatives), polynomial_count(cell_type, degree), points.shape().0)
(
derivative_count(cell_type, derivatives),
polynomial_count(cell_type, degree),
points.shape().0,
)
}

/// Tabulate orthonormal polynomials
Expand Down
6 changes: 3 additions & 3 deletions tools/src/types.rs
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
//! Representation of basic types

// Specification of data types
#[derive(PartialEq, Debug, Clone, Copy)]
#[derive(Eq, PartialEq, Debug, Clone, Copy)]
#[repr(u8)]
pub enum DTYPE {
/// 32 bit float
Expand Down Expand Up @@ -40,15 +40,15 @@ pub enum DTYPE {
// pub(crate) use iterate_over_type;

// Mutability Property
#[derive(PartialEq, Debug, Clone, Copy)]
#[derive(Eq, PartialEq, Debug, Clone, Copy)]
#[repr(u8)]
pub enum MUTABILITY {
NotMutable = 0,
Mutable = 1,
}

// Ownership Property
#[derive(PartialEq, Debug, Clone, Copy)]
#[derive(Eq, PartialEq, Debug, Clone, Copy)]
#[repr(u8)]
pub enum OWNERSHIP {
NotOwner = 0,
Expand Down

0 comments on commit 57a8ef3

Please sign in to comment.