diff --git a/notebooks/convergence/Poisson_MMS.ipynb b/notebooks/convergence/Poisson_MMS.ipynb index 0e1a2df8..cb2757a5 100644 --- a/notebooks/convergence/Poisson_MMS.ipynb +++ b/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, diff --git a/src/mesh/procedural.rs b/src/mesh/procedural.rs index 04c81a8d..99c01e07 100644 --- a/src/mesh/procedural.rs +++ b/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(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( unit_length: T, units_x: usize, @@ -290,8 +293,113 @@ pub fn create_rectangular_uniform_tet_mesh( 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, radius: f64, num_sweeps: usize) -> PolyMesh3d { diff --git a/tests/convergence_tests/poisson_3d_mms.rs b/tests/convergence_tests/poisson_3d_mms.rs index 244dbe66..81dfeb5b 100644 --- a/tests/convergence_tests/poisson_3d_mms.rs +++ b/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(); diff --git a/tests/convergence_tests/poisson_mms_common.rs b/tests/convergence_tests/poisson_mms_common.rs index 4dc862e6..bc236b8f 100644 --- a/tests/convergence_tests/poisson_mms_common.rs +++ b/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( &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); diff --git a/tests/convergence_tests/reference_values/poisson2d_mms_quad4_summary.json b/tests/convergence_tests/reference_values/poisson2d_mms_quad4_summary.json index c8c10546..9e9ce5ad 100644 --- a/tests/convergence_tests/reference_values/poisson2d_mms_quad4_summary.json +++ b/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 ] } \ No newline at end of file diff --git a/tests/convergence_tests/reference_values/poisson2d_mms_quad9_summary.json b/tests/convergence_tests/reference_values/poisson2d_mms_quad9_summary.json index a5577d3d..986b538a 100644 --- a/tests/convergence_tests/reference_values/poisson2d_mms_quad9_summary.json +++ b/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 ] } \ No newline at end of file diff --git a/tests/convergence_tests/reference_values/poisson2d_mms_tri3_summary.json b/tests/convergence_tests/reference_values/poisson2d_mms_tri3_summary.json index c0413993..6da73cf1 100644 --- a/tests/convergence_tests/reference_values/poisson2d_mms_tri3_summary.json +++ b/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 ] } \ No newline at end of file diff --git a/tests/convergence_tests/reference_values/poisson2d_mms_tri6_summary.json b/tests/convergence_tests/reference_values/poisson2d_mms_tri6_summary.json index c8ee0ba7..72a47c60 100644 --- a/tests/convergence_tests/reference_values/poisson2d_mms_tri6_summary.json +++ b/tests/convergence_tests/reference_values/poisson2d_mms_tri6_summary.json @@ -1,27 +1,27 @@ { "element_name": "Tri6", "L2_errors": [ - 0.2950216401236216, + 0.2950216401236215, 0.03632508513027139, - 0.004496286928297505, - 0.0005542209405580082, - 0.00006894024090577562, - 8.606868763085124e-6 + 0.00449628692829757, + 0.0005542209405583764, + 0.00006894024090654412, + 8.606868765783322e-6 ], "H1_seminorm_errors": [ - 1.5170463630635265, + 1.5170463630635262, 0.47436999257492424, 0.1301154241424629, - 0.03343726912784644, - 0.008422416891254678, - 0.002109732968527345 + 0.033437269127846427, + 0.00842241689125471, + 0.0021097329685273373 ], "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 ] } \ No newline at end of file diff --git a/tests/convergence_tests/reference_values/poisson3d_mms_hex20_summary.json b/tests/convergence_tests/reference_values/poisson3d_mms_hex20_summary.json index 71128249..b00447b1 100644 --- a/tests/convergence_tests/reference_values/poisson3d_mms_hex20_summary.json +++ b/tests/convergence_tests/reference_values/poisson3d_mms_hex20_summary.json @@ -2,23 +2,23 @@ "element_name": "Hex20", "L2_errors": [ 0.3535530739635877, - 0.01474560653577477, - 0.001702410217160731, - 0.00021309082427604725, - 0.000026651810179155683 + 0.014745606535774733, + 0.0017024102171607136, + 0.000213090824276002, + 0.000026651810179147626 ], "H1_seminorm_errors": [ 1.9238241709395592, 0.24048791360243685, - 0.04721140317056074, - 0.011222124993446552, - 0.002774143155147405 + 0.04721140317056073, + 0.011222124993446556, + 0.002774143155147397 ], "resolutions": [ - 1.0, - 0.5, - 0.25, - 0.125, - 0.0625 + 1.7320508075688772, + 0.8660254037844386, + 0.4330127018922193, + 0.21650635094610965, + 0.10825317547305482 ] } \ No newline at end of file diff --git a/tests/convergence_tests/reference_values/poisson3d_mms_hex27_summary.json b/tests/convergence_tests/reference_values/poisson3d_mms_hex27_summary.json index a0c2352f..97c8b30a 100644 --- a/tests/convergence_tests/reference_values/poisson3d_mms_hex27_summary.json +++ b/tests/convergence_tests/reference_values/poisson3d_mms_hex27_summary.json @@ -2,23 +2,23 @@ "element_name": "Hex27", "L2_errors": [ 0.02375348764151462, - 0.012106012612151357, - 0.0016658951046437955, - 0.00021209247893351423, - 0.000026621538002910954 + 0.012106012612151367, + 0.0016658951046437979, + 0.00021209247893352607, + 0.000026621538002902395 ], "H1_seminorm_errors": [ 0.2531252878279967, 0.17890810725822845, - 0.04445266581233786, + 0.04445266581233785, 0.011072258366999117, - 0.0027651405468598657 + 0.0027651405468598618 ], "resolutions": [ - 1.0, - 0.5, - 0.25, - 0.125, - 0.0625 + 1.7320508075688772, + 0.8660254037844386, + 0.4330127018922193, + 0.21650635094610965, + 0.10825317547305482 ] } \ No newline at end of file diff --git a/tests/convergence_tests/reference_values/poisson3d_mms_hex8_summary.json b/tests/convergence_tests/reference_values/poisson3d_mms_hex8_summary.json index dc54b5c6..dad83756 100644 --- a/tests/convergence_tests/reference_values/poisson3d_mms_hex8_summary.json +++ b/tests/convergence_tests/reference_values/poisson3d_mms_hex8_summary.json @@ -2,26 +2,26 @@ "element_name": "Hex8", "L2_errors": [ 0.3535530739635877, - 0.09287460149899221, - 0.022983032908891955, - 0.0057456018402315245, - 0.0014366736862290893, - 0.0003591900692691314 + 0.09287460149899229, + 0.022983032908891917, + 0.005745601840231595, + 0.0014366736862290864, + 0.000359190069269035 ], "H1_seminorm_errors": [ 1.9238241709395592, - 0.8876300336370881, - 0.4366606790086766, + 0.8876300336370883, + 0.43666067900867656, 0.2181044665052087, 0.10904521972957422, - 0.054522391133906875 + 0.05452239113390689 ], "resolutions": [ - 1.0, - 0.5, - 0.25, - 0.125, - 0.0625, - 0.03125 + 1.7320508075688772, + 0.8660254037844386, + 0.4330127018922193, + 0.21650635094610965, + 0.10825317547305482, + 0.05412658773652741 ] } \ No newline at end of file diff --git a/tests/convergence_tests/reference_values/poisson3d_mms_tet10_summary.json b/tests/convergence_tests/reference_values/poisson3d_mms_tet10_summary.json index 6a338ef0..376d1c85 100644 --- a/tests/convergence_tests/reference_values/poisson3d_mms_tet10_summary.json +++ b/tests/convergence_tests/reference_values/poisson3d_mms_tet10_summary.json @@ -1,24 +1,21 @@ { "element_name": "Tet10", "L2_errors": [ - 0.22232421782401193, - 0.048384177975326764, - 0.00598532880517876, - 0.0007178741067086166, - 0.00008824613889661203 + 0.0903599855246288, + 0.011524996483904263, + 0.0014117645547182703, + 0.0001721924705989381 ], "H1_seminorm_errors": [ - 1.3058959984844163, - 0.5832316748458222, - 0.17008228628799338, - 0.04506625378562449, - 0.01148027703790675 + 1.0480179096637585, + 0.24329418978653872, + 0.06065557352014549, + 0.014586338228415914 ], "resolutions": [ - 1.0, - 0.5, - 0.25, - 0.125, - 0.0625 + 1.4142135623730951, + 0.7071067811865476, + 0.3535533905932738, + 0.1767766952966369 ] } \ No newline at end of file diff --git a/tests/convergence_tests/reference_values/poisson3d_mms_tet4_summary.json b/tests/convergence_tests/reference_values/poisson3d_mms_tet4_summary.json index 845e4ab4..29f6fbf6 100644 --- a/tests/convergence_tests/reference_values/poisson3d_mms_tet4_summary.json +++ b/tests/convergence_tests/reference_values/poisson3d_mms_tet4_summary.json @@ -1,24 +1,24 @@ { "element_name": "Tet4", "L2_errors": [ - 0.3522698228234682, - 0.25843728396409893, - 0.09756747623016845, - 0.027579551652551715, - 0.0071292211857160665 + 0.17154711789929086, + 0.11007731274055563, + 0.029738674333404543, + 0.007675981971733075, + 0.0019466622689481867 ], "H1_seminorm_errors": [ - 1.9181575922071978, - 1.5488285882830495, - 0.9146645014087618, - 0.47957926282573976, - 0.2428023769718095 + 1.115207809759971, + 0.9602553651703537, + 0.485998412034791, + 0.243648213400367, + 0.12189465613669546 ], "resolutions": [ - 1.0, - 0.5, - 0.25, - 0.125, - 0.0625 + 1.4142135623730951, + 0.7071067811865476, + 0.3535533905932738, + 0.1767766952966369, + 0.08838834764831845 ] } \ No newline at end of file diff --git a/tests/unit_tests/mesh.rs b/tests/unit_tests/mesh.rs index dc9712bb..2ac1b37c 100644 --- a/tests/unit_tests/mesh.rs +++ b/tests/unit_tests/mesh.rs @@ -14,6 +14,7 @@ use proptest::collection::vec; use proptest::prelude::*; use std::cmp::max; +mod procedural; mod refinement; #[test] diff --git a/tests/unit_tests/mesh/procedural.rs b/tests/unit_tests/mesh/procedural.rs new file mode 100644 index 00000000..f1e1e07c --- /dev/null +++ b/tests/unit_tests/mesh/procedural.rs @@ -0,0 +1,191 @@ +use fenris::assembly::global::assemble_scalar; +use fenris::connectivity::Connectivity; +use fenris::element::{ElementConnectivity, FiniteElement, SurfaceFiniteElement}; +use fenris::integrate::{dependency::NoDeps, FnFunction, UFunction}; +use fenris::integrate::{integrate_over_element, volume_form, ElementIntegralAssemblerBuilder}; +use fenris::io::vtk::FiniteElementMeshDataSetBuilder; +use fenris::mesh::procedural::{create_rectangular_uniform_hex_mesh, create_rectangular_uniform_tet_mesh}; +use fenris::quadrature::CanonicalMassQuadrature; +use fenris::quadrature::Quadrature; +use fenris::util::global_vector_from_point_fn; +use fenris_geometry::AxisAlignedBoundingBox3d; +use matrixcompare::prop_assert_scalar_eq; +use nalgebra::coordinates::XYZ; +use nalgebra::{dvector, vector, Point3, Vector1, Vector3, Vector4, U1}; +use proptest::prelude::*; +use std::path::PathBuf; + +#[test] +fn rectangular_uniform_tet_mesh_basics() { + let base_path = PathBuf::from("data/unit_tests/mesh/procedural"); + for res in [1, 2] { + let mesh = create_rectangular_uniform_tet_mesh(1.0, 1, 1, 1, res); + FiniteElementMeshDataSetBuilder::from_mesh(&mesh) + .try_export(base_path.join(format!("basic_uniform_{res}.vtk"))) + .unwrap(); + // The idea here is not to inspect the actual debug output but instead to + // look at the VTK files and check that they look reasonable + insta::assert_debug_snapshot!(format!("mesh_{res}"), mesh); + } +} + +fn empty_tet_mesh_params() -> impl Strategy { + let strategy = prop_oneof![Just(0), 0usize..3]; + [strategy.clone(), strategy.clone(), strategy.clone(), strategy] + .prop_filter("At least one of the parameters must be zero", |params: &[usize; 4]| { + params.iter().product::() == 0 + }) +} + +#[derive(Debug)] +struct RectangularUniformTetMeshParams { + unit_length: f64, + units: [usize; 3], + resolution: usize, +} + +impl Arbitrary for RectangularUniformTetMeshParams { + type Parameters = (); + + fn arbitrary_with(_args: ()) -> Self::Strategy { + let units = [1usize..2, 1..2, 1..2]; + (0.1f64..10.0, units, 1..3usize) + .prop_map(|(unit_length, units, resolution)| Self { + unit_length, + units, + resolution, + }) + .boxed() + } + + type Strategy = BoxedStrategy; +} + +proptest! { + #[test] + fn rectangular_uniform_tet_mesh_empty([units_x, units_y, units_z, cells_per_unit] in empty_tet_mesh_params()) { + let mesh = create_rectangular_uniform_tet_mesh(1.0, units_x, units_y, units_z, cells_per_unit); + prop_assert!(mesh.vertices().is_empty()); + prop_assert!(mesh.connectivity().is_empty()); + } + + #[test] + fn rectangular_uniform_tet_mesh_polynomial_integral(params: RectangularUniformTetMeshParams) { + // Integrate a function f(x, u), where u = u(x) is a linear function + // (so that both linear tets and trilinear hex elements are able to represent it exactly) + // and f overall is linear. + let u = |p: &Point3<_>| vector![4.0 * p.x + 2.0 * p.y - 2.0 * p.z + 3.0]; + let f = FnFunction::new(|p: &Point3<_>, u: &Vector1<_>| vector![-3.0 * p.x + 5.0 * p.y - 2.0 * p.z + 5.0 - 3.0 * u[0]]); + + let RectangularUniformTetMeshParams { unit_length, units, resolution } = params; + let [nx, ny, nz] = units; + let tet_mesh = create_rectangular_uniform_tet_mesh(unit_length, nx, ny, nz, resolution); + let tet_quadrature = tet_mesh.canonical_mass_quadrature(); + let u_tet = global_vector_from_point_fn(tet_mesh.vertices(), u); + + // Use minimal element hex mesh to compute reference value on + let hex_mesh = create_rectangular_uniform_hex_mesh(unit_length, nx, ny, nz, 1); + let hex_quadrature = hex_mesh.canonical_mass_quadrature(); + let u_hex = global_vector_from_point_fn(hex_mesh.vertices(), u); + + let hex_assembler = ElementIntegralAssemblerBuilder::new() + .with_quadrature_table(&hex_quadrature) + .with_space(&hex_mesh) + .with_integrand(f.clone()) + .with_interpolation_weights(&u_hex) + .build_integrator(); + let hex_result = assemble_scalar(&hex_assembler).unwrap(); + + let tet_assembler = ElementIntegralAssemblerBuilder::new() + .with_quadrature_table(&tet_quadrature) + .with_space(&tet_mesh) + .with_integrand(f) + .with_interpolation_weights(&u_tet) + .build_integrator(); + let tet_result = assemble_scalar(&tet_assembler).unwrap(); + + prop_assert_scalar_eq!(tet_result, hex_result, comp = abs, tol = hex_result.abs() * 1e-12); + } + + #[test] + fn rectangular_uniform_tet_mesh_geometric_invariants(params: RectangularUniformTetMeshParams) { + let RectangularUniformTetMeshParams { unit_length, units, resolution } = params; + let [nx, ny, nz] = units; + let tet_mesh = create_rectangular_uniform_tet_mesh(unit_length, nx, ny, nz, resolution); + + let extents = [nx, ny, nz].map(|n| (n as f64) * unit_length); + let aabb = AxisAlignedBoundingBox3d::from_points(tet_mesh.vertices()).unwrap(); + prop_assert_eq!(aabb.min(), &Point3::origin()); + prop_assert_eq!(aabb.max(), &Point3::from(extents)); + + for connectivity in tet_mesh.connectivity() { + let volume_element = connectivity.element(tet_mesh.vertices()).unwrap(); + // Jacobian is constant so can pick any reference point + let j_det = volume_element.reference_jacobian(&Point3::origin()).determinant(); + prop_assert!(j_det > 0.0, "element is inverted"); + } + } + + #[test] + fn rectangular_uniform_tet_mesh_element_wise_divergence(params: RectangularUniformTetMeshParams) { + // We check that the divergence theorem applies to each element. This implicitly + // is intended to verify that all elements have correct orientation. This is however + // somewhat redundant given that checking for inversion already accomplishes the same thing. + // However it acts as an additional sanity check that the implementations for the + // normal computation and face connectivity are sane + + let RectangularUniformTetMeshParams { unit_length, units, resolution } = params; + let [nx, ny, nz] = units; + let tet_mesh = create_rectangular_uniform_tet_mesh(unit_length, nx, ny, nz, resolution); + let tet_quadrature = fenris::quadrature::total_order::tetrahedron(1).unwrap(); + let tri_quadrature = fenris::quadrature::total_order::triangle(2).unwrap(); + + let div_f = |p: &Point3| { + let XYZ { x, y, z } = *p.coords; + vector![ + (6.0 * x + 2.0) + + (4.0 * y + 3.0 - 3.0 * x) + + (2.0 * z - 3.0 - 3.0 * x) + ] + }; + let div_f = FnFunction::new(div_f) + .with_dependencies::>(); + + let f = |p: &Point3| { + let XYZ { x, y, z } = *p.coords; + vector![ + 3.0 * x.powi(2) + 2.0 * x + 4.0 * y * z - 3.0, + 2.0 * y.powi(2) + 3.0 * y - 3.0 * x * y + 2.0, + z.powi(2) - 3.0 * z - 3.0 * x * z + x * y - 4.0 + ] + }; + + for connectivity in tet_mesh.connectivity() { + let volume_element = connectivity.element(tet_mesh.vertices()).unwrap(); + // Interpolation weights are irrelevant here since our function does not depend on u, + // just use zero weights. + let u = Vector4::repeat(0.0); + let div_integral = integrate_over_element( + &div_f, + &volume_element, + &tet_quadrature, + &u, + &mut Default::default())[0]; + + let mut flux_integral = 0.0; + for face_idx in 0 .. connectivity.num_faces() { + let face_conn = connectivity.get_face_connectivity(face_idx).unwrap(); + let face_element = face_conn.element(tet_mesh.vertices()).unwrap(); + + for (w, xi, _) in tri_quadrature.iter() { + let n = face_element.normal(xi); + let x = face_element.map_reference_coords(xi); + let j = face_element.reference_jacobian(xi); + flux_integral += w * volume_form(&j) * n.dot(&f(&x)); + } + } + let tol = 1e-12 * div_integral.abs(); + prop_assert_scalar_eq!(flux_integral, div_integral, comp = abs, tol = tol); + } + } +} diff --git a/tests/unit_tests/mesh/snapshots/unit__unit_tests__mesh__procedural__basic_uniform_tet_mesh.snap b/tests/unit_tests/mesh/snapshots/unit__unit_tests__mesh__procedural__basic_uniform_tet_mesh.snap new file mode 100644 index 00000000..7f35c416 --- /dev/null +++ b/tests/unit_tests/mesh/snapshots/unit__unit_tests__mesh__procedural__basic_uniform_tet_mesh.snap @@ -0,0 +1,151 @@ +--- +source: tests/unit_tests/mesh/procedural.rs +expression: mesh +--- +Mesh { + vertices: [ + [ + 0.0, + 0.0, + 0.0, + ], + [ + 1.0, + 0.0, + 0.0, + ], + [ + 0.0, + 1.0, + 0.0, + ], + [ + 1.0, + 1.0, + 0.0, + ], + [ + 0.0, + 0.0, + 1.0, + ], + [ + 1.0, + 0.0, + 1.0, + ], + [ + 0.0, + 1.0, + 1.0, + ], + [ + 1.0, + 1.0, + 1.0, + ], + [ + 0.5, + 0.5, + 0.5, + ], + ], + connectivity: [ + Tet4Connectivity( + [ + 4, + 6, + 2, + 8, + ], + ), + Tet4Connectivity( + [ + 4, + 2, + 0, + 8, + ], + ), + Tet4Connectivity( + [ + 1, + 3, + 7, + 8, + ], + ), + Tet4Connectivity( + [ + 1, + 7, + 5, + 8, + ], + ), + Tet4Connectivity( + [ + 4, + 5, + 1, + 8, + ], + ), + Tet4Connectivity( + [ + 4, + 1, + 0, + 8, + ], + ), + Tet4Connectivity( + [ + 2, + 3, + 7, + 8, + ], + ), + Tet4Connectivity( + [ + 2, + 7, + 6, + 8, + ], + ), + Tet4Connectivity( + [ + 2, + 3, + 1, + 8, + ], + ), + Tet4Connectivity( + [ + 2, + 1, + 0, + 8, + ], + ), + Tet4Connectivity( + [ + 4, + 5, + 7, + 8, + ], + ), + Tet4Connectivity( + [ + 4, + 7, + 6, + 8, + ], + ), + ], +} diff --git a/tests/unit_tests/mesh/snapshots/unit__unit_tests__mesh__procedural__mesh_1.snap b/tests/unit_tests/mesh/snapshots/unit__unit_tests__mesh__procedural__mesh_1.snap new file mode 100644 index 00000000..00c63d8f --- /dev/null +++ b/tests/unit_tests/mesh/snapshots/unit__unit_tests__mesh__procedural__mesh_1.snap @@ -0,0 +1,151 @@ +--- +source: tests/unit_tests/mesh/procedural.rs +expression: mesh +--- +Mesh { + vertices: [ + [ + 0.0, + 0.0, + 0.0, + ], + [ + 1.0, + 0.0, + 0.0, + ], + [ + 0.0, + 1.0, + 0.0, + ], + [ + 1.0, + 1.0, + 0.0, + ], + [ + 0.0, + 0.0, + 1.0, + ], + [ + 1.0, + 0.0, + 1.0, + ], + [ + 0.0, + 1.0, + 1.0, + ], + [ + 1.0, + 1.0, + 1.0, + ], + [ + 0.5, + 0.5, + 0.5, + ], + ], + connectivity: [ + Tet4Connectivity( + [ + 0, + 2, + 6, + 8, + ], + ), + Tet4Connectivity( + [ + 0, + 6, + 4, + 8, + ], + ), + Tet4Connectivity( + [ + 5, + 7, + 3, + 8, + ], + ), + Tet4Connectivity( + [ + 5, + 3, + 1, + 8, + ], + ), + Tet4Connectivity( + [ + 4, + 5, + 1, + 8, + ], + ), + Tet4Connectivity( + [ + 4, + 1, + 0, + 8, + ], + ), + Tet4Connectivity( + [ + 2, + 3, + 7, + 8, + ], + ), + Tet4Connectivity( + [ + 2, + 7, + 6, + 8, + ], + ), + Tet4Connectivity( + [ + 0, + 1, + 3, + 8, + ], + ), + Tet4Connectivity( + [ + 0, + 3, + 2, + 8, + ], + ), + Tet4Connectivity( + [ + 6, + 7, + 5, + 8, + ], + ), + Tet4Connectivity( + [ + 6, + 5, + 4, + 8, + ], + ), + ], +} diff --git a/tests/unit_tests/mesh/snapshots/unit__unit_tests__mesh__procedural__mesh_2.snap b/tests/unit_tests/mesh/snapshots/unit__unit_tests__mesh__procedural__mesh_2.snap new file mode 100644 index 00000000..0fed4c38 --- /dev/null +++ b/tests/unit_tests/mesh/snapshots/unit__unit_tests__mesh__procedural__mesh_2.snap @@ -0,0 +1,953 @@ +--- +source: tests/unit_tests/mesh/procedural.rs +expression: mesh +--- +Mesh { + vertices: [ + [ + 0.0, + 0.0, + 0.0, + ], + [ + 0.5, + 0.0, + 0.0, + ], + [ + 1.0, + 0.0, + 0.0, + ], + [ + 0.0, + 0.5, + 0.0, + ], + [ + 0.5, + 0.5, + 0.0, + ], + [ + 1.0, + 0.5, + 0.0, + ], + [ + 0.0, + 1.0, + 0.0, + ], + [ + 0.5, + 1.0, + 0.0, + ], + [ + 1.0, + 1.0, + 0.0, + ], + [ + 0.0, + 0.0, + 0.5, + ], + [ + 0.5, + 0.0, + 0.5, + ], + [ + 1.0, + 0.0, + 0.5, + ], + [ + 0.0, + 0.5, + 0.5, + ], + [ + 0.5, + 0.5, + 0.5, + ], + [ + 1.0, + 0.5, + 0.5, + ], + [ + 0.0, + 1.0, + 0.5, + ], + [ + 0.5, + 1.0, + 0.5, + ], + [ + 1.0, + 1.0, + 0.5, + ], + [ + 0.0, + 0.0, + 1.0, + ], + [ + 0.5, + 0.0, + 1.0, + ], + [ + 1.0, + 0.0, + 1.0, + ], + [ + 0.0, + 0.5, + 1.0, + ], + [ + 0.5, + 0.5, + 1.0, + ], + [ + 1.0, + 0.5, + 1.0, + ], + [ + 0.0, + 1.0, + 1.0, + ], + [ + 0.5, + 1.0, + 1.0, + ], + [ + 1.0, + 1.0, + 1.0, + ], + [ + 0.25, + 0.25, + 0.25, + ], + [ + 0.75, + 0.25, + 0.25, + ], + [ + 0.25, + 0.75, + 0.25, + ], + [ + 0.75, + 0.75, + 0.25, + ], + [ + 0.25, + 0.25, + 0.75, + ], + [ + 0.75, + 0.25, + 0.75, + ], + [ + 0.25, + 0.75, + 0.75, + ], + [ + 0.75, + 0.75, + 0.75, + ], + ], + connectivity: [ + Tet4Connectivity( + [ + 27, + 28, + 13, + 10, + ], + ), + Tet4Connectivity( + [ + 27, + 28, + 4, + 13, + ], + ), + Tet4Connectivity( + [ + 27, + 28, + 1, + 4, + ], + ), + Tet4Connectivity( + [ + 27, + 28, + 10, + 1, + ], + ), + Tet4Connectivity( + [ + 0, + 3, + 12, + 27, + ], + ), + Tet4Connectivity( + [ + 0, + 12, + 9, + 27, + ], + ), + Tet4Connectivity( + [ + 27, + 29, + 4, + 3, + ], + ), + Tet4Connectivity( + [ + 27, + 29, + 13, + 4, + ], + ), + Tet4Connectivity( + [ + 27, + 29, + 12, + 13, + ], + ), + Tet4Connectivity( + [ + 27, + 29, + 3, + 12, + ], + ), + Tet4Connectivity( + [ + 9, + 10, + 1, + 27, + ], + ), + Tet4Connectivity( + [ + 9, + 1, + 0, + 27, + ], + ), + Tet4Connectivity( + [ + 27, + 31, + 13, + 12, + ], + ), + Tet4Connectivity( + [ + 27, + 31, + 10, + 13, + ], + ), + Tet4Connectivity( + [ + 27, + 31, + 9, + 10, + ], + ), + Tet4Connectivity( + [ + 27, + 31, + 12, + 9, + ], + ), + Tet4Connectivity( + [ + 0, + 1, + 4, + 27, + ], + ), + Tet4Connectivity( + [ + 0, + 4, + 3, + 27, + ], + ), + Tet4Connectivity( + [ + 11, + 14, + 2, + 28, + ], + ), + Tet4Connectivity( + [ + 14, + 5, + 2, + 28, + ], + ), + Tet4Connectivity( + [ + 28, + 30, + 5, + 4, + ], + ), + Tet4Connectivity( + [ + 28, + 30, + 14, + 5, + ], + ), + Tet4Connectivity( + [ + 28, + 30, + 13, + 14, + ], + ), + Tet4Connectivity( + [ + 28, + 30, + 4, + 13, + ], + ), + Tet4Connectivity( + [ + 10, + 11, + 1, + 28, + ], + ), + Tet4Connectivity( + [ + 11, + 2, + 1, + 28, + ], + ), + Tet4Connectivity( + [ + 28, + 32, + 14, + 13, + ], + ), + Tet4Connectivity( + [ + 28, + 32, + 11, + 14, + ], + ), + Tet4Connectivity( + [ + 28, + 32, + 10, + 11, + ], + ), + Tet4Connectivity( + [ + 28, + 32, + 13, + 10, + ], + ), + Tet4Connectivity( + [ + 1, + 2, + 4, + 28, + ], + ), + Tet4Connectivity( + [ + 2, + 5, + 4, + 28, + ], + ), + Tet4Connectivity( + [ + 29, + 30, + 16, + 13, + ], + ), + Tet4Connectivity( + [ + 29, + 30, + 7, + 16, + ], + ), + Tet4Connectivity( + [ + 29, + 30, + 4, + 7, + ], + ), + Tet4Connectivity( + [ + 29, + 30, + 13, + 4, + ], + ), + Tet4Connectivity( + [ + 3, + 6, + 12, + 29, + ], + ), + Tet4Connectivity( + [ + 6, + 15, + 12, + 29, + ], + ), + Tet4Connectivity( + [ + 6, + 7, + 15, + 29, + ], + ), + Tet4Connectivity( + [ + 7, + 16, + 15, + 29, + ], + ), + Tet4Connectivity( + [ + 29, + 33, + 16, + 15, + ], + ), + Tet4Connectivity( + [ + 29, + 33, + 13, + 16, + ], + ), + Tet4Connectivity( + [ + 29, + 33, + 12, + 13, + ], + ), + Tet4Connectivity( + [ + 29, + 33, + 15, + 12, + ], + ), + Tet4Connectivity( + [ + 3, + 4, + 6, + 29, + ], + ), + Tet4Connectivity( + [ + 4, + 7, + 6, + 29, + ], + ), + Tet4Connectivity( + [ + 14, + 17, + 8, + 30, + ], + ), + Tet4Connectivity( + [ + 14, + 8, + 5, + 30, + ], + ), + Tet4Connectivity( + [ + 7, + 8, + 17, + 30, + ], + ), + Tet4Connectivity( + [ + 7, + 17, + 16, + 30, + ], + ), + Tet4Connectivity( + [ + 30, + 34, + 17, + 16, + ], + ), + Tet4Connectivity( + [ + 30, + 34, + 14, + 17, + ], + ), + Tet4Connectivity( + [ + 30, + 34, + 13, + 14, + ], + ), + Tet4Connectivity( + [ + 30, + 34, + 16, + 13, + ], + ), + Tet4Connectivity( + [ + 4, + 5, + 8, + 30, + ], + ), + Tet4Connectivity( + [ + 4, + 8, + 7, + 30, + ], + ), + Tet4Connectivity( + [ + 31, + 32, + 22, + 19, + ], + ), + Tet4Connectivity( + [ + 31, + 32, + 13, + 22, + ], + ), + Tet4Connectivity( + [ + 31, + 32, + 10, + 13, + ], + ), + Tet4Connectivity( + [ + 31, + 32, + 19, + 10, + ], + ), + Tet4Connectivity( + [ + 9, + 12, + 18, + 31, + ], + ), + Tet4Connectivity( + [ + 12, + 21, + 18, + 31, + ], + ), + Tet4Connectivity( + [ + 31, + 33, + 13, + 12, + ], + ), + Tet4Connectivity( + [ + 31, + 33, + 22, + 13, + ], + ), + Tet4Connectivity( + [ + 31, + 33, + 21, + 22, + ], + ), + Tet4Connectivity( + [ + 31, + 33, + 12, + 21, + ], + ), + Tet4Connectivity( + [ + 18, + 19, + 9, + 31, + ], + ), + Tet4Connectivity( + [ + 19, + 10, + 9, + 31, + ], + ), + Tet4Connectivity( + [ + 21, + 22, + 18, + 31, + ], + ), + Tet4Connectivity( + [ + 22, + 19, + 18, + 31, + ], + ), + Tet4Connectivity( + [ + 20, + 23, + 14, + 32, + ], + ), + Tet4Connectivity( + [ + 20, + 14, + 11, + 32, + ], + ), + Tet4Connectivity( + [ + 32, + 34, + 14, + 13, + ], + ), + Tet4Connectivity( + [ + 32, + 34, + 23, + 14, + ], + ), + Tet4Connectivity( + [ + 32, + 34, + 22, + 23, + ], + ), + Tet4Connectivity( + [ + 32, + 34, + 13, + 22, + ], + ), + Tet4Connectivity( + [ + 19, + 20, + 11, + 32, + ], + ), + Tet4Connectivity( + [ + 19, + 11, + 10, + 32, + ], + ), + Tet4Connectivity( + [ + 22, + 23, + 20, + 32, + ], + ), + Tet4Connectivity( + [ + 22, + 20, + 19, + 32, + ], + ), + Tet4Connectivity( + [ + 33, + 34, + 25, + 22, + ], + ), + Tet4Connectivity( + [ + 33, + 34, + 16, + 25, + ], + ), + Tet4Connectivity( + [ + 33, + 34, + 13, + 16, + ], + ), + Tet4Connectivity( + [ + 33, + 34, + 22, + 13, + ], + ), + Tet4Connectivity( + [ + 12, + 15, + 24, + 33, + ], + ), + Tet4Connectivity( + [ + 12, + 24, + 21, + 33, + ], + ), + Tet4Connectivity( + [ + 15, + 16, + 25, + 33, + ], + ), + Tet4Connectivity( + [ + 15, + 25, + 24, + 33, + ], + ), + Tet4Connectivity( + [ + 24, + 25, + 22, + 33, + ], + ), + Tet4Connectivity( + [ + 24, + 22, + 21, + 33, + ], + ), + Tet4Connectivity( + [ + 23, + 26, + 14, + 34, + ], + ), + Tet4Connectivity( + [ + 26, + 17, + 14, + 34, + ], + ), + Tet4Connectivity( + [ + 16, + 17, + 25, + 34, + ], + ), + Tet4Connectivity( + [ + 17, + 26, + 25, + 34, + ], + ), + Tet4Connectivity( + [ + 25, + 26, + 22, + 34, + ], + ), + Tet4Connectivity( + [ + 26, + 23, + 22, + 34, + ], + ), + ], +}