Merge pull request #66 from Andlon/bcc_uniform_tet_mesh

Reimplement uniform tet mesh generation in terms of BCC lattice

notebooks/convergence/Poisson_MMS.ipynb

@@ -126,14 +126,14 @@
     "\n",
     "# L2 error\n",
     "fig = prepare_convergence_plot(summaries_2d, error_type='L2')\n",
-    "draw_convergence_line(fig.gca(), p=3, x0=0.25, y0=1e-3, label_xoffset = 0.1)\n",
-    "draw_convergence_line(fig.gca(), p=2, x0=0.25, y0=2e-1, label_xoffset = -0.1)\n",
+    "draw_convergence_line(fig.gca(), p=3, x0=0.25, y0=5e-4, label_xoffset = 0.1)\n",
+    "draw_convergence_line(fig.gca(), p=2, x0=0.25, y0=1e-1, label_xoffset = -0.1)\n",
     "plt.show(fig)\n",
     "\n",
     "# H1 seminorm errors\n",
     "fig = prepare_convergence_plot(summaries_2d, error_type='H1_seminorm')\n",
-    "draw_convergence_line(fig.gca(), p=1, x0=0.25, y0=1.25, label_xoffset = -0.1)\n",
-    "draw_convergence_line(fig.gca(), p=2, x0=0.25, y0=3e-2, label_xoffset = 0.1)\n",
+    "draw_convergence_line(fig.gca(), p=1, x0=0.25, y0=1.0, label_xoffset = -0.1)\n",
+    "draw_convergence_line(fig.gca(), p=2, x0=0.25, y0=2e-2, label_xoffset = 0.1)\n",
     "plt.show(fig)\n"
    ]
   },
@@ -169,13 +169,13 @@
     "summaries_3d = sorted(load_summaries(data_folder_3d), key = summary_key_3d)\n",
     "\n",
     "fig = prepare_convergence_plot(summaries_3d, error_type='L2')\n",
-    "draw_convergence_line(fig.gca(), p=3, x0=0.25, y0=1e-3, label_xoffset = 0.1)\n",
-    "draw_convergence_line(fig.gca(), p=2, x0=0.25, y0=4e-2, label_xoffset = -0.06)\n",
+    "draw_convergence_line(fig.gca(), p=3, x0=0.25, y0=2e-4, label_xoffset = 0.1)\n",
+    "draw_convergence_line(fig.gca(), p=2, x0=0.25, y0=2e-2, label_xoffset = -0.06)\n",
     "plt.show(fig)\n",
     "\n",
     "fig = prepare_convergence_plot(summaries_3d, error_type='H1_seminorm')\n",
-    "draw_convergence_line(fig.gca(), p=1, x0=0.25, y0=1.25, label_xoffset = -0.1)\n",
-    "draw_convergence_line(fig.gca(), p=2, x0=0.25, y0=3e-2, label_xoffset = 0.1)\n",
+    "draw_convergence_line(fig.gca(), p=1, x0=0.25, y0=0.5, label_xoffset = -0.1)\n",
+    "draw_convergence_line(fig.gca(), p=2, x0=0.25, y0=1e-2, label_xoffset = 0.1)\n",
     "plt.show(fig)"
    ]
   },
@@ -196,7 +196,7 @@
  ],
  "metadata": {
   "kernelspec": {
-   "display_name": "Python 3",
+   "display_name": "Python 3 (ipykernel)",
    "language": "python",
    "name": "python3"
   },
@@ -210,7 +210,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.8.3"
+   "version": "3.7.4"
   }
  },
  "nbformat": 4,

src/mesh/procedural.rs

@@ -1,11 +1,12 @@
 //! Basic procedural mesh generation routines.
-use crate::connectivity::{Hex8Connectivity, Quad4d2Connectivity};
+use crate::connectivity::{Hex8Connectivity, Quad4d2Connectivity, Tet4Connectivity};
 use crate::geometry::polymesh::PolyMesh3d;
 use crate::geometry::sdf::BoundedSdf;
 use crate::geometry::{AxisAlignedBoundingBox2d, HalfSpace};
 use crate::mesh::{HexMesh, Mesh, QuadMesh2d, Tet4Mesh, TriangleMesh2d};
 use crate::Real;
-use nalgebra::{convert, try_convert, Point2, Point3, Unit, Vector2, Vector3};
+use itertools::{iproduct, Itertools};
+use nalgebra::{convert, point, try_convert, vector, Point2, Point3, Unit, Vector2, Vector3};
 use numeric_literals::replace_float_literals;
 use ordered_float::NotNan;
 use std::cmp::min;
@@ -37,8 +38,7 @@ pub fn create_unit_box_uniform_tet_mesh_3d<T>(cells_per_dim: usize) -> Tet4Mesh<
 where
     T: Real,
 {
-    let hex_mesh = create_unit_box_uniform_hex_mesh_3d(cells_per_dim);
-    Tet4Mesh::from(&hex_mesh)
+    create_rectangular_uniform_tet_mesh(T::one(), 1, 1, 1, cells_per_dim)
 }
 
 /// Generates an axis-aligned rectangular uniform mesh given a unit length,
@@ -278,8 +278,11 @@ where
 
 /// Creates a rectangular uniform tetrahedral mesh.
 ///
-/// Currently the same as [`create_rectangular_uniform_hex_mesh`], except that the hexahedral
-/// mesh is tetrahedralized.
+/// The implementation uses a BCC lattice, where each pair of adjacent cell centers
+/// along each coordinate direction is connected by an octahedron, which is further subdivided
+/// into four tetrahedra along the edge between the two cell centers. The boundaries of the
+/// cuboidal domain are finally filled with pyramids that are subdivided into two tetrahedra.
+#[replace_float_literals(T::from_f64(literal).unwrap())]
 pub fn create_rectangular_uniform_tet_mesh<T>(
     unit_length: T,
     units_x: usize,
@@ -290,8 +293,113 @@ pub fn create_rectangular_uniform_tet_mesh<T>(
 where
     T: Real,
 {
-    let hex_mesh = create_rectangular_uniform_hex_mesh(unit_length, units_x, units_y, units_z, cells_per_unit);
-    Tet4Mesh::from(&hex_mesh)
+    if units_x == 0 || units_y == 0 || units_z == 0 || cells_per_unit == 0 {
+        return Mesh::from_vertices_and_connectivity(vec![], vec![]);
+    }
+
+    let cell_size = unit_length / T::from_usize(cells_per_unit).unwrap();
+    // Cell, vertex counts along each dimension
+    let [cx, cy, cz] = [units_x, units_y, units_z].map(|units| units * cells_per_unit);
+    let [vx, vy, vz] = [cx, cy, cz].map(|num_cells| num_cells + 1);
+
+    // Construct all vertices first: first all vertices of the (implicit) uniform hex mesh,
+    // then all vertices that correspond to cell centers.
+    let mut vertices = Vec::new();
+    for (k, j, i) in iproduct!(0..vz, 0..vy, 0..vx) {
+        vertices.push(point![
+            cell_size * T::from_usize(i).unwrap(),
+            cell_size * T::from_usize(j).unwrap(),
+            cell_size * T::from_usize(k).unwrap()
+        ]);
+    }
+    let cell_center_offset = vertices.len();
+    for (k, j, i) in iproduct!(0..cz, 0..cy, 0..cx) {
+        vertices.push(point![
+            cell_size * (0.5 + T::from_usize(i).unwrap()),
+            cell_size * (0.5 + T::from_usize(j).unwrap()),
+            cell_size * (0.5 + T::from_usize(k).unwrap())
+        ])
+    }
+
+    let vertex_to_global_idx = |[i, j, k]: [usize; 3]| (vx * vy) * k + vx * j + i;
+    let cell_to_global_midpoint_idx = |[i, j, k]: [usize; 3]| (cx * cy) * k + cx * j + i + cell_center_offset;
+
+    let mut connectivity = Vec::new();
+
+    // Offsets to [i, j, k] coordinates to obtain the vertices of the face connecting [i, j, k] and
+    // its neighbor in the positive direction along each axis 0, 1, 2
+    let positive_face_deltas_for_each_axis = [
+        [[1, 0, 1], [1, 1, 1], [1, 1, 0], [1, 0, 0]].map(Vector3::from),
+        [[0, 1, 0], [1, 1, 0], [1, 1, 1], [0, 1, 1]].map(Vector3::from),
+        [[0, 1, 1], [1, 1, 1], [1, 0, 1], [0, 0, 1]].map(Vector3::from),
+    ];
+
+    let connect_centers_with_tets = |connectivity: &mut Vec<_>, [i, j, k]: [usize; 3], axis: usize| {
+        // Make four tets connecting (i, j, k) and (i + di, j + dj, k + dk).
+        // The octahedron formed by the two cell centers and the common face vertices
+        // is split into four tetrahedra along the edge between the cell centers.
+        let cell = Vector3::from([i, j, k]);
+        let cell_delta = Vector3::from_fn(|idx, _| (idx == axis) as usize);
+
+        let shared_face_vertices = positive_face_deltas_for_each_axis[axis]
+            .map(|delta| cell + delta)
+            .map(|v| vertex_to_global_idx(v.into()));
+        let c1 = cell_to_global_midpoint_idx([i, j, k]);
+        let c2 = cell_to_global_midpoint_idx((cell + cell_delta).into());
+        for (v1, v2) in shared_face_vertices
+            .into_iter()
+            .cycle()
+            .take(5)
+            .tuple_windows()
+        {
+            connectivity.push(Tet4Connectivity([c1, c2, v2, v1]));
+        }
+    };
+
+    let make_pyramid = |connectivity: &mut Vec<_>, [i, j, k]: [usize; 3], axis: usize, positive_dir: bool| {
+        let positive_face_deltas = positive_face_deltas_for_each_axis[axis];
+        let mut face_vertices = positive_face_deltas.map(|delta_coord| delta_coord + vector![i, j, k]);
+        if !positive_dir {
+            // Face vertices are oriented such that they are only correct for the positive
+            // direction, need to flip otherwise.
+            face_vertices.reverse();
+            // Pick the faces one coordinate unit lower
+            for coord in &mut face_vertices {
+                coord[axis] -= 1;
+            }
+        }
+
+        let [a, b, c, d] = face_vertices.map(|v| vertex_to_global_idx(v.into()));
+        let center = cell_to_global_midpoint_idx([i, j, k]);
+
+        // Ensure that the diagonal choice alternates along the boundary,
+        // to prevent excessive diagonal bias along the surface
+        if (i + j + k) % 2 == 0 {
+            connectivity.push(Tet4Connectivity([a, b, c, center]));
+            connectivity.push(Tet4Connectivity([a, c, d, center]));
+        } else {
+            connectivity.push(Tet4Connectivity([a, b, d, center]));
+            connectivity.push(Tet4Connectivity([b, c, d, center]));
+        }
+    };
+
+    for (k, j, i) in iproduct!(0..cz, 0..cy, 0..cx) {
+        let cell = [i, j, k];
+        let num_cells = [cx, cy, cz];
+        for axis in [0, 1, 2] {
+            if cell[axis] + 1 < num_cells[axis] {
+                connect_centers_with_tets(&mut connectivity, cell, axis);
+            }
+            if cell[axis] == 0 {
+                make_pyramid(&mut connectivity, cell, axis, false);
+            }
+            if cell[axis] + 1 == num_cells[axis] {
+                make_pyramid(&mut connectivity, cell, axis, true);
+            }
+        }
+    }
+
+    Mesh::from_vertices_and_connectivity(vertices, connectivity)
 }
 
 pub fn create_simple_stupid_sphere(center: &Point3<f64>, radius: f64, num_sweeps: usize) -> PolyMesh3d<f64> {

tests/convergence_tests/poisson_3d_mms.rs

@@ -119,7 +119,7 @@ fn poisson_3d_tet4() {
 
 #[test]
 fn poisson_3d_tet10() {
-    let resolutions = [1, 2, 4, 8, 16];
+    let resolutions = [1, 2, 4, 8];
     let mesh_producer = |res| Tet10Mesh::from(&create_unit_box_uniform_tet_mesh_3d(res));
     // TODO: Use "correct" quadrature
     let quadrature = quadrature::total_order::tetrahedron(2).unwrap();

tests/convergence_tests/poisson_mms_common.rs

@@ -8,7 +8,7 @@ use fenris::assembly::local::{
     ElementEllipticAssemblerBuilder, ElementSourceAssemblerBuilder, SourceFunction, UniformQuadratureTable,
 };
 use fenris::assembly::operators::LaplaceOperator;
-use fenris::element::ElementConnectivity;
+use fenris::element::{ElementConnectivity, FiniteElement};
 use fenris::error::{estimate_H1_seminorm_error, estimate_L2_error};
 use fenris::io::vtk::{FiniteElementMeshDataSetBuilder, VtkCellConnectivity};
 use fenris::mesh::Mesh;
@@ -241,9 +241,14 @@ pub fn solve_and_produce_output<C, D, Source>(
             &u_exact_grad,
         );
 
-        // Resolution measures number of cells per unit-length, and the unit square is one unit
-        // long.
-        let h = 1.0 / resolution as f64;
+        // We use the maximum diameter as a measure of resolution
+        let h = mesh
+            .connectivity()
+            .iter()
+            .map(|conn| conn.element(mesh.vertices()).unwrap())
+            .map(|element| element.diameter())
+            .max_by(f64::total_cmp)
+            .unwrap();
         summary.resolutions.push(h);
         summary.L2_errors.push(result.L2_error);
         summary.H1_seminorm_errors.push(result.H1_seminorm_error);

tests/convergence_tests/reference_values/poisson2d_mms_quad4_summary.json

@@ -3,25 +3,25 @@
   "L2_errors": [
     0.49999970147862477,
     0.11900335666838031,
-    0.030180169603514433,
-    0.007587213803594591,
-    0.0018997045680722057,
-    0.0004751116684116985
+    0.030180169603514505,
+    0.007587213803594678,
+    0.001899704568072491,
+    0.0004751116684120406
   ],
   "H1_seminorm_errors": [
     2.221441469078788,
     0.9965116176145944,
-    0.5013692054047112,
+    0.5013692054047113,
     0.2515137803184585,
-    0.12587387281638515,
-    0.06295197000212865
+    0.12587387281638512,
+    0.06295197000212867
   ],
   "resolutions": [
-    1.0,
-    0.5,
-    0.25,
-    0.125,
-    0.0625,
-    0.03125
+    1.4142135623730951,
+    0.7071067811865476,
+    0.3535533905932738,
+    0.1767766952966369,
+    0.08838834764831845,
+    0.04419417382415922
   ]
 }

tests/convergence_tests/reference_values/poisson2d_mms_quad9_summary.json

@@ -2,26 +2,26 @@
   "element_name": "Quad9",
   "L2_errors": [
     0.1274022764904238,
-    0.012885935334998632,
-    0.0018756562014574107,
-    0.00024331327849208738,
-    0.00003068953043125697,
-    3.844775193267211e-6
+    0.012885935334998653,
+    0.0018756562014574667,
+    0.0002433132784921437,
+    0.00003068953043123288,
+    3.844775193310826e-6
   ],
   "H1_seminorm_errors": [
     0.5975135375836107,
     0.20744142717034694,
-    0.05125983229930842,
-    0.012778779118497884,
-    0.0031924802336273845,
-    0.0007979824507208034
+    0.05125983229930836,
+    0.012778779118497887,
+    0.003192480233627371,
+    0.0007979824507207673
   ],
   "resolutions": [
-    1.0,
-    0.5,
-    0.25,
-    0.125,
-    0.0625,
-    0.03125
+    1.4142135623730951,
+    0.7071067811865476,
+    0.3535533905932738,
+    0.1767766952966369,
+    0.08838834764831845,
+    0.04419417382415922
   ]
 }

tests/convergence_tests/reference_values/poisson2d_mms_tri3_summary.json

@@ -2,26 +2,26 @@
   "element_name": "Tri3",
   "L2_errors": [
     0.49666268112934125,
-    0.2912097235354121,
-    0.09341017643520341,
-    0.02501754539089316,
-    0.006369056084509816,
-    0.0015996524117294227
+    0.291209723535412,
+    0.0934101764352033,
+    0.025017545390893087,
+    0.006369056084510224,
+    0.0015996524117291757
   ],
   "H1_seminorm_errors": [
     2.2349292649333714,
     1.5332420495006416,
-    0.8421072409316398,
-    0.4322326169347917,
-    0.217590327161479,
-    0.10898216346345377
+    0.8421072409316397,
+    0.4322326169347916,
+    0.21759032716147897,
+    0.10898216346345381
   ],
   "resolutions": [
-    1.0,
-    0.5,
-    0.25,
-    0.125,
-    0.0625,
-    0.03125
+    1.4142135623730951,
+    0.7071067811865476,
+    0.3535533905932738,
+    0.1767766952966369,
+    0.08838834764831845,
+    0.04419417382415922
   ]
 }