Skip to content

Commit

Permalink
formatting
Browse files Browse the repository at this point in the history
  • Loading branch information
@Andlon
Andlon committed Jul 18, 2023
1 parent 8b1f21c commit 9060e95
Show file tree
Hide file tree
Showing 3 changed files with 59 additions and 30 deletions.
44 changes: 29 additions & 15 deletions tests/integration_tests/interpolation.rs
View file
Original file line number Diff line number Diff line change
@@ -1,18 +1,24 @@
use crate::integration_tests::data_output_path;
use fenris::assembly::buffers::{BufferUpdate, InterpolationBuffer};
use fenris::connectivity::{Tri3d2Connectivity};
use fenris::connectivity::Tri3d2Connectivity;
use fenris::io::vtk::FiniteElementMeshDataSetBuilder;
use fenris::mesh::procedural::{create_unit_box_uniform_tet_mesh_3d, create_unit_square_uniform_tri_mesh_2d};
use fenris::mesh::refinement::refine_uniformly_repeat;
use fenris::mesh::{Mesh, Tet4Mesh, TriangleMesh2d};
use fenris::space::{FindClosestElement, FiniteElementConnectivity, FiniteElementSpace, InterpolateGradientInSpace, InterpolateInSpace, SpatiallyIndexed};
use fenris::quadrature::Quadrature;
use fenris::space::{
FindClosestElement, FiniteElementConnectivity, FiniteElementSpace, InterpolateGradientInSpace, InterpolateInSpace,
SpatiallyIndexed,
};
use fenris::util::global_vector_from_point_fn;
use fenris::{quadrature, SmallDim};
use fenris_traits::allocators::{BiDimAllocator, TriDimAllocator};
use itertools::{izip, Itertools};
use matrixcompare::assert_matrix_eq;
use nalgebra::{vector, DVectorView, DefaultAllocator, OMatrix, OVector, Point2, Vector1, Vector2, U1, U2, Point3, Vector3, OPoint, U3};
use fenris::quadrature::Quadrature;
use nalgebra::{
vector, DVectorView, DefaultAllocator, OMatrix, OPoint, OVector, Point2, Point3, Vector1, Vector2, Vector3, U1, U2,
U3,
};
use util::flatten_vertically;

fn u_scalar_2d(p: &Point2<f64>) -> Vector1<f64> {
Expand Down Expand Up @@ -145,7 +151,8 @@ fn spatially_indexed_interpolation_trimesh() {
// space), so that we already know the correct answer.
{
// For interior quadrature points, we check both values and gradients of the scalar function
let values = compute_expected_interpolation_test_values::<_, U1, _>(&space, &interior_points, &u_weights_scalar);
let values =
compute_expected_interpolation_test_values::<_, U1, _>(&space, &interior_points, &u_weights_scalar);
let iter = izip!(
values.u_interpolated,
values.grad_u_interpolated,
Expand All @@ -161,7 +168,8 @@ fn spatially_indexed_interpolation_trimesh() {
{
// For boundary points, we only check values since gradients are discontinuous
// at element interfaces
let values = compute_expected_interpolation_test_values::<_, U1, _>(&space, &interface_points, &u_weights_scalar);
let values =
compute_expected_interpolation_test_values::<_, U1, _>(&space, &interface_points, &u_weights_scalar);
let iter = izip!(values.u_interpolated, values.u_expected);
for (u, u_expected) in iter {
assert_matrix_eq!(u, u_expected, comp = abs, tol = 1e-12);
Expand All @@ -170,7 +178,8 @@ fn spatially_indexed_interpolation_trimesh() {

{
// Repeat interior quadrature points for vector function
let values = compute_expected_interpolation_test_values::<_, U2, _>(&space, &interior_points, &u_weights_vector);
let values =
compute_expected_interpolation_test_values::<_, U2, _>(&space, &interior_points, &u_weights_vector);
let iter = izip!(
values.u_interpolated,
values.grad_u_interpolated,
Expand All @@ -185,7 +194,8 @@ fn spatially_indexed_interpolation_trimesh() {

{
// Repeat interface quadrature points for vector function
let values = compute_expected_interpolation_test_values::<_, U2, _>(&space, &interface_points, &u_weights_vector);
let values =
compute_expected_interpolation_test_values::<_, U2, _>(&space, &interface_points, &u_weights_vector);
let iter = izip!(values.u_interpolated, values.u_expected);
for (u, u_expected) in iter {
assert_matrix_eq!(u, u_expected, comp = abs, tol = 1e-12);
Expand Down Expand Up @@ -301,7 +311,7 @@ fn spatially_indexed_interpolation_tet4() {
[-1.0, -1.0, 1.0],
[-(1.0 / 3.0), -(1.0 / 3.0), -(1.0 / 3.0)],
]
.map(Point3::from);
.map(Point3::from);

// For debugging
FiniteElementMeshDataSetBuilder::from_mesh(space.space())
Expand All @@ -313,7 +323,8 @@ fn spatially_indexed_interpolation_tet4() {
// space), so that we already know the correct answer.
{
// For interior quadrature points, we check both values and gradients of the scalar function
let values = compute_expected_interpolation_test_values::<_, U1, _>(&space, &interior_points, &u_weights_scalar);
let values =
compute_expected_interpolation_test_values::<_, U1, _>(&space, &interior_points, &u_weights_scalar);
let u_interpolated = flatten_vertically(&values.u_interpolated).unwrap();
let u_expected = flatten_vertically(&values.u_expected).unwrap();
let grad_u_interpolated = flatten_vertically(&values.grad_u_interpolated).unwrap();
Expand All @@ -325,15 +336,17 @@ fn spatially_indexed_interpolation_tet4() {
{
// For boundary points, we only check values since gradients are discontinuous
// at element interfaces
let values = compute_expected_interpolation_test_values::<_, U1, _>(&space, &interface_points, &u_weights_scalar);
let values =
compute_expected_interpolation_test_values::<_, U1, _>(&space, &interface_points, &u_weights_scalar);
let u_interpolated = flatten_vertically(&values.u_interpolated).unwrap();
let u_expected = flatten_vertically(&values.u_expected).unwrap();
assert_matrix_eq!(u_interpolated, u_expected, comp = abs, tol = 1e-12);
}

{
// Repeat interior quadrature points for vector function
let values = compute_expected_interpolation_test_values::<_, U3, _>(&space, &interior_points, &u_weights_vector);
let values =
compute_expected_interpolation_test_values::<_, U3, _>(&space, &interior_points, &u_weights_vector);
let u_interpolated = flatten_vertically(&values.u_interpolated).unwrap();
let u_expected = flatten_vertically(&values.u_expected).unwrap();
let grad_u_interpolated = flatten_vertically(&values.grad_u_interpolated).unwrap();
Expand All @@ -344,7 +357,8 @@ fn spatially_indexed_interpolation_tet4() {

{
// Repeat interface quadrature points for vector function
let values = compute_expected_interpolation_test_values::<_, U3, _>(&space, &interface_points, &u_weights_vector);
let values =
compute_expected_interpolation_test_values::<_, U3, _>(&space, &interface_points, &u_weights_vector);
let u_interpolated = flatten_vertically(&values.u_interpolated).unwrap();
let u_expected = flatten_vertically(&values.u_expected).unwrap();
assert_matrix_eq!(u_interpolated, u_expected, comp = abs, tol = 1e-12);
Expand All @@ -361,7 +375,7 @@ fn spatially_indexed_tet4_find_closest() {

let indexed_mesh = SpatiallyIndexed::from_space(mesh);
let quadrature = quadrature::total_order::tetrahedron(0).unwrap();
for element_idx in 0 .. indexed_mesh.num_elements() {
for element_idx in 0..indexed_mesh.num_elements() {
for xi_q in quadrature.points() {
let x_q = indexed_mesh.map_element_reference_coords(element_idx, xi_q);
let (closest_element_idx, xi_closest) = indexed_mesh
Expand All @@ -372,4 +386,4 @@ fn spatially_indexed_tet4_find_closest() {
assert_matrix_eq!(xi_closest.coords, xi_q.coords, comp = abs, tol = 1e-12);
}
}
}
}
32 changes: 22 additions & 10 deletions tests/unit_tests/element.rs
View file
Original file line number Diff line number Diff line change
@@ -1,17 +1,20 @@
use itertools::{Itertools, izip};
use fenris::element::{FixedNodesReferenceFiniteElement, Hex20Element, Hex27Element, Quad4d2Element, Quad9d2Element, Tet4Element, Tri3d2Element, Tri6d2Element};
use crate::export_mesh_vtk;
use fenris::connectivity::Tet4Connectivity;
use fenris::element::{
FixedNodesReferenceFiniteElement, Hex20Element, Hex27Element, Quad4d2Element, Quad9d2Element, Tet4Element,
Tri3d2Element, Tri6d2Element,
};
use fenris::mesh::Tet4Mesh;
use fenris_traits::Real;
use itertools::{izip, Itertools};
use nalgebra::{Point2, Point3, Vector3};
use num::clamp;
use numeric_literals::replace_float_literals;
use proptest::array::uniform4;
use proptest::prelude::*;
use proptest::strategy::ValueTree;
use proptest::test_runner::TestRunner;
use fenris::connectivity::Tet4Connectivity;
use fenris::mesh::Tet4Mesh;
use util::assert_approx_matrix_eq;
use crate::export_mesh_vtk;

mod hexahedron;
mod quadrilateral;
Expand Down Expand Up @@ -40,7 +43,7 @@ fn point_in_hex_ref_domain() -> impl Strategy<Value = Point3<f64>> {
}

fn point_in_tet_ref_domain() -> impl Strategy<Value = Point3<f64>> {
uniform4(0.0 ..= 1.0f64)
uniform4(0.0..=1.0f64)
.prop_map(|mut barycentric_coords| {
let mut sum = 0.0;
for lambda_i in &mut barycentric_coords {
Expand All @@ -61,7 +64,9 @@ fn point_in_tet_ref_domain() -> impl Strategy<Value = Point3<f64>> {
if cfg!(debug_assertions) {
let sum = barycentric_coords.iter().sum();
debug_assert!(1.0 - 1e-12 <= sum && sum <= 1.0 + 1e-12);
debug_assert!(barycentric_coords.iter().all(|&lambda_i| 0.0 <= lambda_i && lambda_i <= 1.0));
debug_assert!(barycentric_coords
.iter()
.all(|&lambda_i| 0.0 <= lambda_i && lambda_i <= 1.0));
}
barycentric_coords
})
Expand Down Expand Up @@ -100,8 +105,14 @@ proptest! {
fn sample_points_in_tet_ref_domain() {
let mut runner = TestRunner::deterministic();
let tet_point_strategy = point_in_tet_ref_domain();
let points: Vec<_> = (0 .. 1000)
.map(|_| tet_point_strategy.new_tree(&mut runner).unwrap().current().clone())
let points: Vec<_> = (0..1000)
.map(|_| {
tet_point_strategy
.new_tree(&mut runner)
.unwrap()
.current()
.clone()
})
.collect();

// TODO: Actually export a point cloud and not a "fake" tet mesh
Expand All @@ -111,7 +122,8 @@ fn sample_points_in_tet_ref_domain() {
let reference_tet = Tet4Element::reference();
let reference_tet_mesh = Tet4Mesh::from_vertices_and_connectivity(
reference_tet.vertices().to_vec(),
vec![Tet4Connectivity([0, 1, 2, 3])]);
vec![Tet4Connectivity([0, 1, 2, 3])],
);
export_mesh_vtk("sample_points_in_tet_ref_domain", "reference_tet", &reference_tet_mesh);
}

Expand Down
13 changes: 8 additions & 5 deletions tests/unit_tests/element/tetrahedron.rs
View file
Original file line number Diff line number Diff line change
@@ -1,4 +1,6 @@
use crate::unit_tests::element::{is_definitely_in_tet_ref_interior, is_likely_in_tet_ref_interior, point_in_tet_ref_domain, point_in_tri_ref_domain};
use crate::unit_tests::element::{
is_definitely_in_tet_ref_interior, is_likely_in_tet_ref_interior, point_in_tet_ref_domain, point_in_tri_ref_domain,
};
use fenris::connectivity::{Connectivity, Tet4Connectivity};
use fenris::element::{
ClosestPoint, ClosestPointInElement, ElementConnectivity, FiniteElement, FixedNodesReferenceFiniteElement,
Expand All @@ -12,7 +14,10 @@ use fenris_geometry::Triangle;
use fenris_optimize::calculus::{approximate_jacobian, VectorFunctionBuilder};
use itertools::izip;
use matrixcompare::{assert_scalar_eq, prop_assert_matrix_eq};
use nalgebra::{distance, DVectorView, DimName, Dyn, MatrixView, OMatrix, OPoint, Point3, Vector1, Vector3, U1, U10, U20, U3, U4, point};
use nalgebra::{
distance, point, DVectorView, DimName, Dyn, MatrixView, OMatrix, OPoint, Point3, Vector1, Vector3, U1, U10, U20,
U3, U4,
};
use proptest::array::uniform3;
use proptest::prelude::*;
use util::assert_approx_matrix_eq;
Expand Down Expand Up @@ -70,16 +75,14 @@ fn tet20_lagrange_property() {

#[test]
fn tet4_closest_point_failure_case() {
let vertices = [[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [1.0, 1.0, 0.0], [0.5, 0.5, 0.5]]
.map(Point3::from);
let vertices = [[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [1.0, 1.0, 0.0], [0.5, 0.5, 0.5]].map(Point3::from);
let element = Tet4Element::from_vertices(vertices);
let p = point![0.875, 0.375, 0.375];
let closest = element.closest_point(&p);

let xi_closest = closest.point();
let x_closest = element.map_reference_coords(&xi_closest);
assert_ne!(x_closest, p);

}

proptest! {
Expand Down

0 comments on commit 9060e95

Please sign in to comment.