Merge pull request #977 from DLR-AMR/fix-batch_inside_elem_axis_aligned

update number of the corner of the second vertex

example/cmesh/t8_cmesh_hypercube_pad.cxx

@@ -24,17 +24,18 @@
 #include <t8_cmesh.h>
 #include <t8_cmesh_vtk_writer.h>
 #include <t8_cmesh/t8_cmesh_examples.h>
-
+#include <t8_schemes/t8_default/t8_default_cxx.hxx>
+#include <t8_forest/t8_forest_general.h>
+#include <t8_forest/t8_forest_vtk.h>
 int
 main (int argc, char **argv)
 {
   /* The prefix for our output files. */
-  const char prefix[BUFSIZ] = "t8_cmesh_hypercube_pad";
+  const char prefix[BUFSIZ] = "t8_forest_hypercube_pad";
   t8_locidx_t local_num_trees;
   t8_gloidx_t global_num_trees;
 
-  const double boundary_coords[24] = { 1, 0, 0, 4, 0, 0, 0, 6, 0, 5, 5, 0, -1, -2, 8, 9, 0, 10, 0, 8, 9, 10, 10, 10 };
-
+  const double boundary_coords[24] = { 0, 0, 0, 1, 0, 0, 0, 1, 0, 1, 1, 0, 0, 0, 1, 1, 0, 1, 0, 1, 1, 1, 1, 1 };
   /* Initialize MPI. This has to happen before we initialize sc or t8code. */
   int mpiret = sc_MPI_Init (&argc, &argv);
   /* Error check the MPI return value. */
@@ -44,20 +45,19 @@ main (int argc, char **argv)
   sc_init (sc_MPI_COMM_WORLD, 1, 1, NULL, SC_LP_ESSENTIAL);
   /* Initialize t8code with log level SC_LP_PRODUCTION. See sc.h for more info on the log levels. */
   t8_init (SC_LP_PRODUCTION);
-
   /* Add hypercube with given element class. */
-  t8_cmesh_t cmesh = t8_cmesh_new_hypercube_pad (T8_ECLASS_HEX, sc_MPI_COMM_WORLD, boundary_coords, 3, 3, 3, 0);
+  t8_cmesh_t cmesh = t8_cmesh_new_hypercube_pad (T8_ECLASS_HEX, sc_MPI_COMM_WORLD, boundary_coords, 3, 3, 3, true);
 
   /* Compute local and global number of trees. */
   local_num_trees = t8_cmesh_get_num_local_trees (cmesh);
   global_num_trees = t8_cmesh_get_num_trees (cmesh);
   t8_global_productionf (" [step1] Created coarse mesh.\n");
   t8_global_productionf (" [step1] Local number of trees:\t%i\n", local_num_trees);
   t8_global_productionf (" [step1] Global number of trees:\t%li\n", global_num_trees);
-
-  t8_cmesh_vtk_write_file (cmesh, prefix);
-
-  t8_cmesh_destroy (&cmesh);
+  t8_scheme_cxx *scheme = t8_scheme_new_default_cxx ();
+  t8_forest_t forest = t8_forest_new_uniform (cmesh, scheme, 0, 0, sc_MPI_COMM_WORLD);
+  t8_forest_vtk_write_file (forest, prefix, 1, 1, 1, 1, 0, 0, NULL);
+  t8_forest_unref (&forest);
 
   sc_finalize ();
 

src/t8_geometry/t8_geometry_helpers.c

@@ -211,7 +211,7 @@ t8_geom_compute_linear_axis_aligned_geometry (t8_eclass_t tree_class, const doub
     const size_t offset_domain_dim = i_coord * T8_ECLASS_MAX_DIM;
     for (int i_dim = 0; i_dim < T8_ECLASS_MAX_DIM; ++i_dim) {
       out_coords[offset_domain_dim + i_dim] = tree_vertices[i_dim];
-      out_coords[offset_domain_dim + i_dim] += ref_coords[offset_tree_dim] * vector[i_dim];
+      out_coords[offset_domain_dim + i_dim] += ref_coords[offset_tree_dim + i_dim] * vector[i_dim];
     }
   }
 }

src/t8_geometry/t8_geometry_implementations/t8_geometry_linear_axis_aligned.cxx

@@ -77,9 +77,14 @@ t8_geometry_linear_axis_aligned::t8_geom_point_batch_inside_element (t8_forest_t
   double v_min[3];
   double v_max[3];
 
+  const t8_eclass_t elem_class = t8_forest_get_tree_class (forest, ltreeid);
+  T8_ASSERT (elem_class == T8_ECLASS_LINE || elem_class == T8_ECLASS_HEX || elem_class == T8_ECLASS_QUAD);
+
+  const int max_corner_id = (elem_class == T8_ECLASS_LINE) ? 1 : (elem_class == T8_ECLASS_QUAD) ? 3 : 7;
+
   /*Geometry is fully described by v_min and v_max*/
   t8_forest_element_coordinate (forest, ltreeid, element, 0, v_min);
-  t8_forest_element_coordinate (forest, ltreeid, element, 1, v_max);
+  t8_forest_element_coordinate (forest, ltreeid, element, max_corner_id, v_max);
 #if T8_ENABLE_DEBUG
   const double coords[6] = { v_min[0], v_min[1], v_min[2], v_max[0], v_max[1], v_max[2] };
   T8_ASSERT (correct_point_order (coords));

test/t8_geometry/t8_geometry_implementations/t8_gtest_geometry_linear.cxx

@@ -25,88 +25,143 @@
  */
 
 #include <test/t8_gtest_macros.hxx>
-#include <test/t8_geometry/t8_gtest_geometry_macros.hxx>
+#include <test/t8_gtest_custom_assertion.hxx>
 #include <gtest/gtest.h>
 #include <t8_eclass.h>
 #include <t8_cmesh.hxx>
 #include <t8_cmesh/t8_cmesh_examples.h>
 #include <t8_geometry/t8_geometry.h>
 #include <t8_geometry/t8_geometry_implementations/t8_geometry_linear.hxx>
+#include <t8_geometry/t8_geometry_implementations/t8_geometry_linear_axis_aligned.hxx>
 #include <t8_element.h>
 
-t8_cmesh_t
-t8_geometry_testing_tree_from_class (const t8_eclass_t eclass)
-{
-  t8_cmesh_t cmesh;
-  t8_cmesh_init (&cmesh);
-  const int num_vertices = t8_eclass_num_vertices[eclass];
-  double *vertices = T8_ALLOC (double, T8_ECLASS_MAX_DIM *num_vertices);
-  for (int i_vertex = 0; i_vertex < num_vertices; ++i_vertex) {
-    for (int dim = 0; dim < T8_ECLASS_MAX_DIM; ++dim) {
-      vertices[i_vertex * T8_ECLASS_MAX_DIM + dim] = t8_element_corner_ref_coords[eclass][i_vertex][dim];
+class geometry_test: public testing::TestWithParam<std::tuple<int, t8_eclass>> {
+ public:
+  static void
+  SetUpTestSuite ()
+  {
+    seed = time (NULL); /* RNG seed */
+  }
+  static int seed;
+
+ protected:
+  void
+  SetUp () override
+  {
+    const int geom_int = std::get<0> (GetParam ());
+    eclass = std::get<1> (GetParam ());
+    t8_cmesh_init (&cmesh);
+    if (geom_int == T8_GEOMETRY_TYPE_LINEAR_AXIS_ALIGNED
+        && !(eclass == T8_ECLASS_LINE || eclass == T8_ECLASS_QUAD || eclass == T8_ECLASS_HEX)) {
+      GTEST_SKIP ();
+    }
+
+    const int num_vertices = t8_eclass_num_vertices[eclass];
+    t8_cmesh_set_tree_class (cmesh, 0, eclass);
+    double *vertices = T8_ALLOC_ZERO (double, num_vertices *T8_ECLASS_MAX_DIM);
+    for (int i_vertex = 0; i_vertex < num_vertices; ++i_vertex) {
+      for (int dim = 0; dim < T8_ECLASS_MAX_DIM; ++dim) {
+        vertices[i_vertex * T8_ECLASS_MAX_DIM + dim] = t8_element_corner_ref_coords[eclass][i_vertex][dim];
+      }
     }
+    switch (geom_int) {
+    case T8_GEOMETRY_TYPE_LINEAR:
+      geom = new t8_geometry_linear (t8_eclass_to_dimension[eclass]);
+      t8_cmesh_register_geometry<t8_geometry_linear> (cmesh, t8_eclass_to_dimension[eclass]);
+      break;
+    case T8_GEOMETRY_TYPE_LINEAR_AXIS_ALIGNED:
+      geom = new t8_geometry_linear_axis_aligned (t8_eclass_to_dimension[eclass]);
+      /* Copy last vertex to the second position*/
+      vertices[3] = vertices[3 * (num_vertices - 1)];
+      vertices[4] = vertices[3 * (num_vertices - 1) + 1];
+      vertices[5] = vertices[3 * (num_vertices - 1) + 2];
+
+      t8_cmesh_register_geometry<t8_geometry_linear_axis_aligned> (cmesh, t8_eclass_to_dimension[eclass]);
+      break;
+    default:
+      break;
+    }
+    t8_cmesh_set_tree_vertices (cmesh, 0, vertices,
+                                geom_int == T8_GEOMETRY_TYPE_LINEAR ? t8_eclass_num_vertices[eclass] : 2);
+    t8_cmesh_commit (cmesh, sc_MPI_COMM_WORLD);
+    T8_FREE (vertices);
   }
-  t8_cmesh_set_tree_class (cmesh, 0, eclass);
-  t8_cmesh_register_geometry<t8_geometry_linear> (cmesh, t8_eclass_to_dimension[eclass]);
-  t8_cmesh_set_tree_vertices (cmesh, 0, vertices, t8_eclass_num_vertices[eclass]);
-  t8_cmesh_commit (cmesh, sc_MPI_COMM_WORLD);
-  T8_FREE (vertices);
-  return cmesh;
-}
+  void
+  TearDown () override
+  {
+    t8_cmesh_unref (&cmesh);
+  }
+  t8_eclass_t eclass;
+  t8_geometry_with_vertices *geom;
+  t8_cmesh_t cmesh;
+#ifdef T8_ENABLE_LESS_TESTS
+  const int num_points = 1000;
+#else
+  const int num_points = 10000;
+#endif
+};
 
-/* Check whether the linear geometry map is correct.
- * We create a cmesh of one tree and unit element
+int geometry_test::seed;
+
+/* Check whether the linear axis aligned geometry map is correct.
+ * We create a cmesh of one tree and unit line, quad or hex
  * geometry. We then create random points in a reference tree and
  * check whether the evaluation is correct. */
-TEST_P (geometry_test, cmesh_geometry_linear)
+TEST_P (geometry_test, cmesh_geometry)
 {
   /* TODO: Add a test for the jacobian, as soon as its implemented. */
 
-  t8_debugf ("Testing linear geometry evaluation.\n");
-
   /* Create random points in [0,1]^d and check if they are mapped correctly. */
-  t8_geometry_linear linear_geom (t8_eclass_to_dimension[eclass]);
-  t8_cmesh_t cmesh;
 
-  double point[3];
   double point_mapped[3];
-  const int seed = 0; /* RNG seed */
   const t8_geometry_c *cmesh_geom;
 
-  cmesh = t8_geometry_testing_tree_from_class (eclass);
-
-  /* Double check that the geometry is the linear geometry. */
+  /* Double check that the geometry is the linear axis aligned geometry. */
   cmesh_geom = t8_cmesh_get_tree_geometry (cmesh, 0);
   ASSERT_TRUE (cmesh_geom != NULL) << "Could not get cmesh's geometry.";
-  int has_same_name = cmesh_geom->t8_geom_get_name () == linear_geom.t8_geom_get_name ();
-  ASSERT_EQ (has_same_name, 1) << "cmesh's geometry is not the linear geometry.";
+  ASSERT_EQ (cmesh_geom->t8_geom_get_hash (), geom->t8_geom_get_hash ())
+    << "cmesh's geometry is not the expected geometry.";
 
   srand (seed);
   for (int ipoint = 0; ipoint < num_points; ++ipoint) {
+    double point[3] = { 0 };
     /* Compute random coordinates in [0,1].
      * These are seen as reference coordinates in the single
      * cmesh tree. Our geometry will map them into the physical
      * space. Since this space is also [0,1] and the cmesh only
      * has one tree, the mapped coordinates must be the same as the 
      * reference coordinates. */
-    point[0] = (double) rand () / RAND_MAX;
-    point[1] = (double) rand () / RAND_MAX;
-    point[2] = (double) rand () / RAND_MAX;
+    for (int idim = 0; idim < t8_eclass_to_dimension[eclass]; ++idim) {
+      point[idim] = (double) rand () / RAND_MAX;
+    }
 
     /* Evaluate the geometry */
     t8_geometry_evaluate (cmesh, 0, point, 1, point_mapped);
     /* Check that the first dim coordinates are the same */
-    int idim;
-    for (idim = 0; idim < t8_eclass_to_dimension[eclass]; ++idim) {
-      ASSERT_NEAR (point[idim], point_mapped[idim], T8_PRECISION_SQRT_EPS) << "Linear geometry computed wrong value.";
-    }
-    /* Check that the remaining entries are 0. */
-    for (; idim < 3; ++idim) {
-      ASSERT_EQ (point_mapped[idim], 0) << "Linear geometry computed wrong value.";
-    }
+    EXPECT_VEC3_EQ (point, point_mapped, T8_PRECISION_SQRT_EPS);
   }
-  /* Destroy the cmesh */
-  t8_cmesh_destroy (&cmesh);
 }
 
-INSTANTIATE_TEST_SUITE_P (t8_gtest_geometry, geometry_test, AllEclasses);
+auto print_test = [] (const testing::TestParamInfo<std::tuple<int, t8_eclass>> &info) {
+  std::string name;
+  const int geom_int = std::get<0> (info.param);
+  const t8_eclass_t eclass = std::get<1> (info.param);
+  switch (geom_int) {
+  case T8_GEOMETRY_TYPE_LINEAR:
+    name = std::string ("linear_geometry_");
+    break;
+  case T8_GEOMETRY_TYPE_LINEAR_AXIS_ALIGNED:
+    name = std::string ("linear_axis_aligned_geometry_");
+    break;
+  default:
+    name = std::string ("geometry_not_allowed_");
+    break;
+  }
+  name += std::string (t8_eclass_to_string[eclass]);
+  return name;
+};
+
+INSTANTIATE_TEST_SUITE_P (
+  t8_gtest_geometry, geometry_test,
+  ::testing::Combine (::testing::Values (T8_GEOMETRY_TYPE_LINEAR, T8_GEOMETRY_TYPE_LINEAR_AXIS_ALIGNED), AllEclasses),
+  print_test);

test/t8_geometry/t8_gtest_geometry_macros.hxx

@@ -1,47 +0,0 @@
-/*
-  This file is part of t8code.
-  t8code is a C library to manage a collection (a forest) of multiple
-  connected adaptive space-trees of general element classes in parallel.
-
-  Copyright (C) 2024 the developers
-
-  t8code is free software; you can redistribute it and/or modify
-  it under the terms of the GNU General Public License as published by
-  the Free Software Foundation; either version 2 of the License, or
-  (at your option) any later version.
-
-  t8code is distributed in the hope that it will be useful,
-  but WITHOUT ANY WARRANTY; without even the implied warranty of
-  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-  GNU General Public License for more details.
-
-  You should have received a copy of the GNU General Public License
-  along with t8code; if not, write to the Free Software Foundation, Inc.,
-  51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
-*/
-
-/** \file t8_gtest_geometry_macros.hxx
- * This file contains macros and base classes for geometry tests.
- */
-
-#include <gtest/gtest.h>
-#include <t8_eclass.h>
-
-class geometry_test: public testing::TestWithParam<t8_eclass_t> {
- protected:
-  void
-  SetUp () override
-  {
-    eclass = GetParam ();
-  }
-  void
-  TearDown () override
-  {
-  }
-  t8_eclass_t eclass;
-#ifdef T8_ENABLE_LESS_TESTS
-  const int num_points = 1000;
-#else
-  const int num_points = 10000;
-#endif
-};

test/t8_geometry/t8_gtest_point_inside.cxx

@@ -251,10 +251,12 @@ TEST_P (geometry_point_inside, test_point_inside)
         num_steps = 2;
       }
       /* Corrected number of points due to possible rounding errors in pow */
+      const int min_points_outside = 6;
       const int num_points = sc_intpow (num_steps, num_corners - 1);
-      double *test_point = T8_ALLOC_ZERO (double, num_points * 3);
-      int *point_is_inside = T8_ALLOC (int, num_points);
-      int *point_is_recognized_as_inside = T8_ALLOC (int, num_points);
+      const int total_points = num_points + min_points_outside;
+      double *test_point = T8_ALLOC_ZERO (double, total_points * 3);
+      int *point_is_inside = T8_ALLOC (int, total_points);
+      int *point_is_recognized_as_inside = T8_ALLOC (int, total_points);
       double step = (barycentric_range_upper_bound - barycentric_range_lower_bound) / (num_steps - 1);
       //t8_debugf ("step size %g, steps %i, points %i (corners %i)\n", step,
       //           num_steps, num_points, num_corners);
@@ -301,6 +303,15 @@ TEST_P (geometry_point_inside, test_point_inside)
 
         num_in += point_is_inside ? 1 : 0;
       }
+      /* Create a set of points that are outside of the cube. 
+       * The coordinates are given by the corner points except for one coordinate. For each side of the cube
+       * we place one point outside of it. */
+      for (int side = 0; side < min_points_outside; ++side) {
+        test_point[3 * (num_points + side)] = (side == 0) ? -2.0 : (side == 1) ? 2.0 : -0.1;
+        test_point[3 * (num_points + side) + 1] = (side == 2) ? -2.0 : (side == 3) ? 2.0 : 0.3;
+        test_point[3 * (num_points + side) + 2] = (side == 4) ? -2.0 : (side == 5) ? 2.0 : 0.15;
+        point_is_inside[num_points + side] = false;
+      }
       /* We now check whether the point inside function correctly sees whether
          * the point is inside the element or not. */
       t8_forest_element_points_inside (forest, 0, element, test_point, num_points, point_is_recognized_as_inside,