Skip to content

Commit

Permalink
Fix arbitrary degree Lagrange elements
Browse files Browse the repository at this point in the history
  • Loading branch information
@mscroggs
mscroggs committed Aug 30, 2023
1 parent eee6ffa commit 93e2aca
Show file tree
Hide file tree
Showing 5 changed files with 166 additions and 29 deletions.
26 changes: 16 additions & 10 deletions element/src/element.rs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8,7 +8,7 @@ use bempp_traits::cell::ReferenceCellType;
use bempp_traits::element::{ElementFamily, FiniteElement, MapType};
use rlst_algorithms::linalg::LinAlg;
use rlst_algorithms::traits::inverse::Inverse;
use rlst_dense::{RandomAccessMut, RandomAccessByRef};
use rlst_dense::{RandomAccessByRef, RandomAccessMut};
pub mod lagrange;
pub mod raviart_thomas;

Expand Down Expand Up @@ -41,6 +41,7 @@ impl CiarletElement {
) -> CiarletElement {
let mut dim = 0;
let mut npts = 0;

for emats in &m {
for mat in emats {
dim += mat.shape().0;
Expand All @@ -67,9 +68,13 @@ impl CiarletElement {
let mut new_x = [vec![], vec![], vec![], vec![]];
let mut pn = 0;
let mut all_pts = Array2D::<f64>::new((npts, tdim));
for i in 0..tdim {
for _pts in &x[i] {
new_x[i].push(Array2D::<f64>::new((0, tdim)));
}
}
for i in 0..tdim + 1 {
for pts in &x[i] {
new_x[i].push(Array2D::<f64>::new((0, tdim)));
for j in 0..pts.shape().0 {
for k in 0..tdim {
*all_pts.get_mut(pn + j, k).unwrap() = *pts.get(j, k).unwrap();
Expand All @@ -88,9 +93,13 @@ impl CiarletElement {
let mut pn = 0;
let mut dn = 0;
let mut all_mat = Array3D::<f64>::new((dim, value_size, npts));
for i in 0..tdim {
for _mat in &m[i] {
new_m[i].push(Array3D::<f64>::new((0, value_size, 0)));
}
}
for i in 0..tdim + 1 {
for mat in &m[i] {
new_m[i].push(Array3D::<f64>::new((0, value_size, 0)));
for j in 0..mat.shape().0 {
for k in 0..value_size {
for l in 0..mat.shape().2 {
Expand Down Expand Up @@ -156,7 +165,8 @@ impl CiarletElement {
for i in 0..dim {
for j in 0..value_size {
for k in 0..dim / value_size {
*coefficients.get_mut(i, j, k).unwrap() = *inverse.get(i, j * dim / value_size + k).unwrap()
*coefficients.get_mut(i, j, k).unwrap() =
*inverse.get(i, j * dim / value_size + k).unwrap()
}
}
}
Expand Down Expand Up @@ -194,7 +204,6 @@ impl CiarletElement {
pub fn value_shape(&self) -> &Vec<usize> {
&self.value_shape
}

}

impl FiniteElement for CiarletElement {
Expand Down Expand Up @@ -249,11 +258,8 @@ impl FiniteElement for CiarletElement {
let value = data.get_mut(d, p, b, j).unwrap();
*value = 0.0;
for i in 0..table.shape().1 {
*value +=
*self.coefficients.get(b, j, i).unwrap()
* *table.get_mut(d, i, p).unwrap();
println!("data({d}, {p}, {b}, {j}) += {} * {}", *self.coefficients.get(b, j, i).unwrap(), *table.get_mut(d, i, p).unwrap());
println!("value = {}", *value);
*value += *self.coefficients.get(b, j, i).unwrap()
* *table.get_mut(d, i, p).unwrap();
}
}
}
Expand Down
155 changes: 147 additions & 8 deletions element/src/element/lagrange.rs
View file
Original file line number Diff line number Diff line change
@@ -1,9 +1,9 @@
//! Lagrange elements

use crate::element::OldCiarletElement;
use crate::element::{CiarletElement, create_cell};
use crate::element::{create_cell, CiarletElement};
use crate::polynomials::polynomial_count;
use bempp_tools::arrays::{AdjacencyList, Array3D, Array2D};
use bempp_tools::arrays::{AdjacencyList, Array2D, Array3D};
use bempp_traits::arrays::Array3DAccess;
use bempp_traits::cell::ReferenceCellType;
use bempp_traits::element::{ElementFamily, MapType};
Expand Down Expand Up @@ -198,13 +198,115 @@ pub fn create_new(
x[d].push(Array2D::<f64>::new((0, tdim)));
m[d].push(Array3D::<f64>::new((0, 1, 0)));
}
x[tdim].push(Array2D::<f64>::from_data(cell.midpoint(), (1, tdim)));
m[tdim].push(Array3D::<f64>::from_data(vec![1.0], (1, 1, 1)));
}
x[tdim].push(Array2D::<f64>::from_data(cell.midpoint(), (1, tdim)));
m[tdim].push(Array3D::<f64>::from_data(vec![1.0], (1, 1, 1)));
} else {

// TODO: GLL points
for e in 0..cell.entity_count(0) {
let mut pts = vec![0.0; tdim];
let vertex = &cell.vertices()[e * tdim..(e + 1) * tdim];
for i in 0..tdim {
pts[i] = vertex[i];
}
x[0].push(Array2D::<f64>::from_data(pts, (1, tdim)));
m[0].push(Array3D::<f64>::from_data(vec![1.0], (1, 1, 1)));
}
for e in 0..cell.entity_count(1) {
let mut pts = vec![0.0; tdim * (degree - 1)];
let mut ident = vec![0.0; (degree - 1).pow(2)];
let vn0 = cell.edges()[2 * e];
let vn1 = cell.edges()[2 * e + 1];
let v0 = &cell.vertices()[vn0 * tdim..(vn0 + 1) * tdim];
let v1 = &cell.vertices()[vn1 * tdim..(vn1 + 1) * tdim];
let mut n = 0;
for i in 1..degree {
ident[n * degree] = 1.0;
for j in 0..tdim {
pts[n * tdim + j] = v0[j] + i as f64 / degree as f64 * (v1[j] - v0[j]);
}
n += 1;
}
x[1].push(Array2D::<f64>::from_data(pts, (degree - 1, tdim)));
m[1].push(Array3D::<f64>::from_data(
ident,
(degree - 1, 1, degree - 1),
));
}
let mut start = 0;
for e in 0..cell.entity_count(2) {
let nvertices = cell.faces_nvertices()[e];
let npts = if nvertices == 3 {
if degree > 2 {
(degree - 1) * (degree - 2) / 2
} else {
0
}
} else if nvertices == 4 {
(degree - 1).pow(2)
} else {
panic!("Unsupported face type");
};
let mut pts = vec![0.0; tdim * npts];
let mut ident = vec![0.0; npts.pow(2)];

let vn0 = cell.faces()[start];
let vn1 = cell.faces()[start + 1];
let vn2 = cell.faces()[start + 2];
let v0 = &cell.vertices()[vn0 * tdim..(vn0 + 1) * tdim];
let v1 = &cell.vertices()[vn1 * tdim..(vn1 + 1) * tdim];
let v2 = &cell.vertices()[vn2 * tdim..(vn2 + 1) * tdim];

if nvertices == 3 {
// Triangle
let mut n = 0;
for i0 in 1..degree {
for i1 in 1..degree - i0 {
for j in 0..tdim {
pts[n * tdim + j] = v0[j]
+ i0 as f64 / degree as f64 * (v1[j] - v0[j])
+ i1 as f64 / degree as f64 * (v2[j] - v0[j]);
}
n += 1;
}
}
} else if nvertices == 4 {
// Quadrilateral
let mut n = 0;
for i0 in 1..degree {
for i1 in 1..degree {
for j in 0..tdim {
pts[n * tdim + j] = v0[j]
+ i0 as f64 / degree as f64 * (v1[j] - v0[j])
+ i1 as f64 / degree as f64 * (v2[j] - v0[j]);
}
n += 1;
}
}
} else {
panic!("Unsupported face type.");
}

for i in 0..npts {
ident[i * npts + i] = 1.0;
}
x[2].push(Array2D::<f64>::from_data(pts, (npts, tdim)));
m[2].push(Array3D::<f64>::from_data(ident, (npts, 1, npts)));
start += nvertices;
}
}
CiarletElement::create(ElementFamily::Lagrange, cell_type, degree, vec![], wcoeffs, x, m, MapType::Identity, discontinuous, degree)
CiarletElement::create(
ElementFamily::Lagrange,
cell_type,
degree,
vec![],
wcoeffs,
x,
m,
MapType::Identity,
discontinuous,
degree,
)
}

#[cfg(test)]
Expand Down Expand Up @@ -420,8 +522,6 @@ mod test {
check_dofs(e);
}



#[test]
fn test_lagrange_0_interval() {
let e = create_new(ReferenceCellType::Interval, 0, true);
Expand Down Expand Up @@ -493,6 +593,45 @@ mod test {
check_dofs(e);
}

#[test]
fn test_lagrange_higher_degree_triangle() {
create_new(ReferenceCellType::Triangle, 2, false);
create_new(ReferenceCellType::Triangle, 3, false);
create_new(ReferenceCellType::Triangle, 4, false);
create_new(ReferenceCellType::Triangle, 5, false);

create_new(ReferenceCellType::Triangle, 2, true);
create_new(ReferenceCellType::Triangle, 3, true);
create_new(ReferenceCellType::Triangle, 4, true);
create_new(ReferenceCellType::Triangle, 5, true);
}

#[test]
fn test_lagrange_higher_degree_interval() {
create_new(ReferenceCellType::Interval, 2, false);
create_new(ReferenceCellType::Interval, 3, false);
create_new(ReferenceCellType::Interval, 4, false);
create_new(ReferenceCellType::Interval, 5, false);

create_new(ReferenceCellType::Interval, 2, true);
create_new(ReferenceCellType::Interval, 3, true);
create_new(ReferenceCellType::Interval, 4, true);
create_new(ReferenceCellType::Interval, 5, true);
}

#[test]
fn test_lagrange_higher_degree_quadrilateral() {
create_new(ReferenceCellType::Quadrilateral, 2, false);
create_new(ReferenceCellType::Quadrilateral, 3, false);
create_new(ReferenceCellType::Quadrilateral, 4, false);
create_new(ReferenceCellType::Quadrilateral, 5, false);

create_new(ReferenceCellType::Quadrilateral, 2, true);
create_new(ReferenceCellType::Quadrilateral, 3, true);
create_new(ReferenceCellType::Quadrilateral, 4, true);
create_new(ReferenceCellType::Quadrilateral, 5, true);
}

#[test]
fn test_lagrange_0_quadrilateral() {
let e = create_new(ReferenceCellType::Quadrilateral, 0, true);
Expand Down
8 changes: 0 additions & 8 deletions element/src/polynomials.rs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -68,12 +68,6 @@ fn tabulate_legendre_polynomials_quadrilateral<'a>(
assert_eq!(data.shape().2, points.shape().0);
assert_eq!(points.shape().1, 2);

for j in 0..derivatives + 1 {
for i in 0..derivatives + 1 - j {
println!("{i} {j} {}", tri_index(i, j));
}
}

for i in 0..data.shape().2 {
*data
.get_mut(tri_index(0, 0), quad_index(0, 0, degree), i)
Expand Down Expand Up @@ -471,7 +465,6 @@ mod test {
product +=
data.get(0, i, k).unwrap() * data.get(0, j, k).unwrap() * rule.weights[k];
}
println!("{i} {j} {product}");
}
}
for i in 0..data.shape().1 {
Expand Down Expand Up @@ -567,7 +560,6 @@ mod test {
let mut index = 0;
for i in 0..10 {
for j in 0..10 - i {
println!("{index}");
p[6 * index] = i as f64 / 10.0;
p[6 * index + 1] = j as f64 / 10.0;
p[6 * index + 2] = p[6 * index] + epsilon;
Expand Down
2 changes: 1 addition & 1 deletion tools/src/arrays.rs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,7 +6,7 @@ use std::clone::Clone;
/// A two-dimensional rectangular array
pub struct Array2D<T: Num> {
/// The data in the array, in row-major order
data: Vec<T>,
pub data: Vec<T>,
/// The shape of the array
shape: (usize, usize),
}
Expand Down
4 changes: 2 additions & 2 deletions traits/src/cell.rs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -41,7 +41,7 @@ pub trait ReferenceCell {
let vertices = self.vertices();
let mut m = vec![0.0; dim];
for (i, v) in vertices.iter().enumerate() {
m[i%dim] += v;
m[i % dim] += v;
}
for i in 0..dim {
m[i] /= (vertices.len() / dim) as f64;
Expand All @@ -56,7 +56,7 @@ pub trait ReferenceCell {

/// The faces of the cell
///
/// The first `self.faces_nvertices()[0]` components are the vertex numbers of vertices of the first face, the next `self.faces_nvertices()[1]` the second edge, and so on.
/// The first `self.faces_nvertices()[0]` components are the vertex numbers of vertices of the first face, the next `self.faces_nvertices()[1]` the second face, and so on.
fn faces(&self) -> &[usize];

/// The number of vertices adjacent to each face
Expand Down

0 comments on commit 93e2aca

Please sign in to comment.