Merge pull request #67 from Andlon/tet_closest_point

ClosestPoint impl for Tet4Element, enables interpolation together with SpatiallyIndexed

fenris-solid/Cargo.toml

@@ -14,7 +14,7 @@ rustdoc-args = [ "--html-in-header", "assets/doc-header.html",
                  "--html-before-content", "assets/doc-header.html" ]
 
 [dependencies]
-fenris = { workspace = true, path = ".." }
+fenris = { workspace = true }
 serde = "1.0.126"
 numeric_literals = "0.2.0"
 

src/element/tetrahedron.rs

@@ -1,12 +1,18 @@
 use numeric_literals::replace_float_literals;
+use std::cmp::Ordering;
 
-use crate::connectivity::{Tet10Connectivity, Tet20Connectivity, Tet4Connectivity};
-use crate::element::{ElementConnectivity, FiniteElement, FixedNodesReferenceFiniteElement};
+use crate::connectivity::{Connectivity, Tet10Connectivity, Tet20Connectivity, Tet4Connectivity};
+use crate::element::{
+    BoundsForElement, ClosestPoint, ClosestPointInElement, ElementConnectivity, FiniteElement,
+    FixedNodesReferenceFiniteElement,
+};
 use crate::nalgebra::{
     distance, Matrix1x4, Matrix3, Matrix3x4, OMatrix, OPoint, Point3, Scalar, Vector3, U1, U10, U20, U3, U4,
 };
 use crate::Real;
+use fenris_geometry::AxisAlignedBoundingBox;
 use itertools::Itertools;
+use nalgebra::distance_squared;
 
 impl<T> ElementConnectivity<T> for Tet4Connectivity
 where
@@ -76,6 +82,10 @@ where
     }
 }
 
+/// A Tet10 element, see documentation of [`Tet10Connectivity`](crate::connectivity::Tet10Connectivity).
+///
+/// We currently assume that the reference-to-physical transformation is affine, meaning
+/// that the geometry of the element is assumed to be affine.
 #[derive(Copy, Clone, Debug, PartialEq, Eq)]
 pub struct Tet10Element<T>
 where
@@ -558,6 +568,12 @@ where
     }
 }
 
+#[replace_float_literals(T::from_f64(literal).unwrap())]
+fn is_likely_in_tet_ref_interior<T: Real>(xi: &Point3<T>) -> bool {
+    let eps = 4.0 * T::default_epsilon();
+    xi.x >= -1.0 - eps && xi.y >= -1.0 - eps && xi.z >= -1.0 - eps && xi.x + xi.y + xi.z <= -1.0 + eps
+}
+
 impl<T> FiniteElement<T> for Tet4Element<T>
 where
     T: Real,
@@ -590,3 +606,67 @@ where
             .fold(T::zero(), |a, b| a.max(b.clone()))
     }
 }
+
+impl<T: Real> BoundsForElement<T> for Tet4Element<T> {
+    fn element_bounds(&self) -> AxisAlignedBoundingBox<T, Self::GeometryDim> {
+        AxisAlignedBoundingBox::from_points(self.vertices()).unwrap()
+    }
+}
+
+impl<T: Real> ClosestPointInElement<T> for Tet4Element<T> {
+    #[allow(non_snake_case)]
+    fn closest_point(&self, p: &OPoint<T, Self::GeometryDim>) -> ClosestPoint<T, Self::ReferenceDim> {
+        let xi_interior = {
+            // Transformation is affine, so Jacobian is constant:
+            //  p = A xi + p0
+            // for some p0 which we can determine by evaluating at xi = 0
+            let A = self.reference_jacobian(&Point3::origin());
+            A.try_inverse()
+                .map(|a_inv| {
+                    let p0 = self.map_reference_coords(&Point3::origin());
+                    Point3::from(a_inv * (p - p0))
+                })
+                // If the inverse transformation doesn't lead to a point clearly inside
+                // the reference domain, we assume that the closest point is on the boundary
+                .filter(is_likely_in_tet_ref_interior)
+        };
+
+        let conn = Tet4Connectivity([0, 1, 2, 3]);
+        let face_elements_iter = (0..4)
+            .map(|face_idx| conn.get_face_connectivity(face_idx).unwrap())
+            .map(|face_conn| face_conn.element(self.vertices()).unwrap());
+
+        let (face_idx, _, xi_face, dist2_face) = face_elements_iter
+            .enumerate()
+            .map(|(face_idx, tri_element)| {
+                let xi_closest = tri_element.closest_point(p).point().clone();
+                let x_closest = tri_element.map_reference_coords(&xi_closest);
+                let dist2 = distance_squared(&x_closest, p);
+                (face_idx, tri_element, xi_closest, dist2)
+            })
+            .min_by(|(_, _, _, d1), (_, _, _, d2)| d1.partial_cmp(d2).unwrap_or(Ordering::Less))
+            .expect("Always have 4 > 0 faces");
+
+        if let Some(xi_interior) = xi_interior {
+            let x_interior = self.map_reference_coords(&xi_interior);
+            let dist2_interior = distance_squared(p, &x_interior);
+            if dist2_interior < dist2_face {
+                return ClosestPoint::InElement(xi_interior);
+            }
+        }
+
+        // Next, we need to obtain the coordinates for the point on the face in the
+        // tetrahedron reference element.
+        // We can do this by considering the corresponding tetrahedron reference element
+        // and use the same face index to map into "physical space" which will in fact
+        // be the reference coordinates that we need
+        let reference_element = Self::reference();
+        let xi = conn
+            .get_face_connectivity(face_idx)
+            .unwrap()
+            .element(reference_element.vertices())
+            .unwrap()
+            .map_reference_coords(&xi_face);
+        ClosestPoint::ClosestPoint(xi)
+    }
+}

src/element/triangle.rs

@@ -1,9 +1,3 @@
-use fenris_geometry::AxisAlignedBoundingBox;
-use itertools::Itertools;
-use nalgebra::distance_squared;
-use numeric_literals::replace_float_literals;
-use std::cmp::Ordering;
-
 use crate::connectivity::{Tri3d2Connectivity, Tri3d3Connectivity, Tri6d2Connectivity};
 use crate::element::{
     BoundsForElement, ClosestPoint, ClosestPointInElement, ElementConnectivity, FiniteElement,
@@ -15,6 +9,11 @@ use crate::nalgebra::{
     Vector2, Vector3, U2, U3, U6,
 };
 use crate::Real;
+use fenris_geometry::{AxisAlignedBoundingBox, LineSegment3d};
+use itertools::Itertools;
+use nalgebra::distance_squared;
+use numeric_literals::replace_float_literals;
+use std::cmp::Ordering;
 
 /// A finite element representing linear basis functions on a triangle, in two dimensions.
 ///
@@ -445,7 +444,7 @@ where
 #[replace_float_literals(T::from_f64(literal).unwrap())]
 fn is_likely_in_tri_ref_interior<T: Real>(xi: &Point2<T>) -> bool {
     let eps = 4.0 * T::default_epsilon();
-    xi.x >= -1.0 + eps && xi.y >= -1.0 + eps && xi.x + xi.y <= eps
+    xi.x >= -1.0 - eps && xi.y >= -1.0 - eps && xi.x + xi.y <= eps
 }
 
 impl<T: Real> ClosestPointInElement<T> for Tri3d2Element<T> {
@@ -532,3 +531,73 @@ impl<T: Real> BoundsForElement<T> for Tri3d2Element<T> {
         AxisAlignedBoundingBox::from_points(self.vertices()).expect("Never fails since we always have > 0 vertices")
     }
 }
+
+#[allow(non_snake_case)]
+impl<T: Real> ClosestPointInElement<T> for Tri3d3Element<T> {
+    fn closest_point(&self, p: &OPoint<T, Self::GeometryDim>) -> ClosestPoint<T, Self::ReferenceDim> {
+        // This follows the same idea as for the implementation for Tri3d2Element,
+        // except that we need to project the point into the triangle of the plane rather
+        // than directly inverting the mapping
+        let xi_interior = {
+            // Want to project the point into the plane by solving
+            //  min_xi 0.5 * ||A xi + p0 - p||^2,
+            // with solution
+            //  A^TA xi = A^T (p - p0)
+            let ref A = self.reference_jacobian(&Point2::origin());
+            let p0 = self.map_reference_coords(&Point2::origin());
+            let ATA = A.transpose() * A;
+            ATA.try_inverse()
+                .map(|ATA_inv| ATA_inv * (A.transpose() * (p - p0)))
+                .map(Point2::from)
+                .filter(is_likely_in_tri_ref_interior)
+        };
+
+        // Compute the closest point on each edge and take the point corresponding to the
+        // smallest distance (squared)
+        let [a, b, c] = self.vertices();
+        let edges = [(a, b), (b, c), (c, a)];
+        let (idx, t, dist2_edge) = edges
+            .into_iter()
+            .map(|(x1, x2)| LineSegment3d::from_end_points(x1.clone(), x2.clone()))
+            .enumerate()
+            .map(|(idx, segment)| {
+                // Parameter is [0, 1]
+                let t = segment.closest_point_parametric(p);
+                let point = segment.point_from_parameter(t);
+                let dist2 = distance_squared(p, &point);
+                (idx, t, dist2)
+            })
+            .min_by(|(_, _, dist2_a), (_, _, dist2_b)| {
+                dist2_a
+                    .partial_cmp(&dist2_b)
+                    // TODO: This is an arbitrary choice, see comment in 2D tri element impl
+                    .unwrap_or(Ordering::Less)
+            })
+            .expect("We always have exactly 3 items in the iterator");
+
+        // Use parameter representation to transfer result to 2D reference element
+        let reference_element = Tri3d2Element::reference();
+        let a = reference_element.vertices()[(idx + 0) % 3];
+        let b = reference_element.vertices()[(idx + 1) % 3];
+        let edge_ref_coords = LineSegment2d::from_end_points(a, b).point_from_parameter(t);
+
+        if let Some(xi_interior) = xi_interior {
+            let x_interior = self.map_reference_coords(&xi_interior);
+            let dist2_interior = distance_squared(p, &x_interior);
+            if dist2_interior < dist2_edge {
+                // We will essentially never find a point that's "inside" the element, since
+                // we here have a 2D element embedded in 3D, and with floating point it's hardly
+                // meaningful to talk about being "in the element" in this case
+                return ClosestPoint::ClosestPoint(xi_interior);
+            }
+        }
+
+        ClosestPoint::ClosestPoint(edge_ref_coords)
+    }
+}
+
+impl<T: Real> BoundsForElement<T> for Tri3d3Element<T> {
+    fn element_bounds(&self) -> AxisAlignedBoundingBox<T, Self::GeometryDim> {
+        AxisAlignedBoundingBox::from_points(self.vertices()).expect("AABB is always well defined")
+    }
+}

src/proptest.rs

@@ -1,11 +1,10 @@
-use crate::element::{Tet4Element, Tri3d2Element};
-use crate::geometry::proptest::Triangle3dParams;
+use crate::element::{Tet4Element, Tri3d2Element, Tri3d3Element};
 use crate::geometry::Orientation::Counterclockwise;
 use crate::mesh::procedural::create_rectangular_uniform_quad_mesh_2d;
 use crate::mesh::QuadMesh2d;
 use ::proptest::prelude::*;
 use fenris_geometry::proptest::Triangle2dParams;
-use fenris_geometry::{Triangle2d, Triangle3d};
+use fenris_geometry::Triangle2d;
 use nalgebra::{Point2, Point3, Vector2};
 use std::cmp::max;
 
@@ -36,32 +35,33 @@ impl Arbitrary for Tri3d2Element<f64> {
     }
 }
 
+impl Arbitrary for Tri3d3Element<f64> {
+    type Parameters = ();
+    type Strategy = BoxedStrategy<Self>;
+
+    fn arbitrary_with(_args: Self::Parameters) -> Self::Strategy {
+        let vertices: [_; 3] = std::array::from_fn(|_| point3());
+        vertices
+            .prop_map(|vertices| Tri3d3Element::from_vertices(vertices))
+            .boxed()
+    }
+}
+
 impl Arbitrary for Tet4Element<f64> {
     // TODO: Reasonable parameters?
     type Parameters = ();
     type Strategy = BoxedStrategy<Self>;
 
     fn arbitrary_with(_args: Self::Parameters) -> Self::Strategy {
-        any_with::<Triangle3d<f64>>(Triangle3dParams::default().with_orientation(Counterclockwise))
-            .prop_flat_map(|triangle| {
-                // To create an arbitrary tetrahedron element, we take a counter-clockwise oriented
-                // triangle, and pick a point somewhere on the "positive" side of the
-                // triangle plane. We do this by associating a parameter with each
-                // tangent vector defined by the sides of the triangle,
-                // plus a non-negative parameter that scales along the normal direction
-                let range = -10.0..10.0;
-                let tangent_params = [range.clone(), range.clone(), range.clone()];
-                let normal_param = 0.0..=10.0;
-                (Just(triangle), tangent_params, normal_param)
-            })
-            .prop_map(|(triangle, tangent_params, normal_param)| {
-                let mut tangent_pos = triangle.centroid();
-
-                for (side, param) in triangle.sides().iter().zip(&tangent_params) {
-                    tangent_pos.coords += *param * side;
-                }
-                let coord = tangent_pos + normal_param * triangle.normal_dir().normalize();
-                Tet4Element::from_vertices([triangle.0[0], triangle.0[1], triangle.0[2], coord])
+        let l = 5.0;
+        (Tri3d3Element::arbitrary(), -l..l, -l..l, 0.0..l)
+            .prop_map(|(tri_element, t1_scale, t2_scale, n_scale)| {
+                let [a, b, c] = tri_element.vertices().clone();
+                let t1 = b - a;
+                let t2 = c - a;
+                let n = t1.cross(&t2);
+                let d = a + t1_scale * t1 + t2_scale * t2 + n_scale * n;
+                Tet4Element::from_vertices([a, b, c, d])
             })
             .boxed()
     }