Merge pull request #747 from DLR-AMR/enhance_hypercube_pad

Add an option to use axis-aligned geometry for hex/quad/lines in hype…

example/cmesh/t8_cmesh_hypercube_pad.cxx

@@ -46,7 +46,7 @@ main (int argc, char **argv)
   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);
+  t8_cmesh_t cmesh = t8_cmesh_new_hypercube_pad (T8_ECLASS_HEX, sc_MPI_COMM_WORLD, boundary_coords, 3, 3, 3, 0);
 
   /* Compute local and global number of trees. */
   local_num_trees = t8_cmesh_get_num_local_trees (cmesh);

example/geometry/t8_example_geometries.cxx

@@ -76,7 +76,8 @@ typedef enum {
  * 
  * This geometry does not provide a jacobian.
  */
-class t8_geometry_sincos: public t8_geometry {
+struct t8_geometry_sincos: public t8_geometry
+{
  public:
   /* Basic constructor that sets the dimension and the name. */
   t8_geometry_sincos (): t8_geometry (2, "t8_sincos_geometry")
@@ -143,7 +144,8 @@ class t8_geometry_sincos: public t8_geometry {
  * It inherits from the w_vertices geometry since we use the tree's vertex coordinates.
  * This geometry does not provide a jacobian.
  */
-class t8_geometry_moebius: public t8_geometry_with_vertices {
+struct t8_geometry_moebius: public t8_geometry_with_vertices
+{
  public:
   /* Basic constructor that sets the dimension and the name. */
   t8_geometry_moebius (): t8_geometry_with_vertices (2, "t8_moebius_geometry")
@@ -208,7 +210,8 @@ class t8_geometry_moebius: public t8_geometry_with_vertices {
  * 
  * This geometry does not provide a jacobian.
  */
-class t8_geometry_cylinder: public t8_geometry {
+struct t8_geometry_cylinder: public t8_geometry
+{
  public:
   /* Basic constructor that sets the dimension and the name. */
   t8_geometry_cylinder (): t8_geometry (2, "t8_cylinder_geometry")
@@ -269,7 +272,8 @@ class t8_geometry_cylinder: public t8_geometry {
  * 
  * This geometry does not provide a jacobian.
  */
-class t8_geometry_circle: public t8_geometry_with_vertices {
+struct t8_geometry_circle: public t8_geometry_with_vertices
+{
  public:
   /* Basic constructor that sets the dimension and the name. */
   t8_geometry_circle (): t8_geometry_with_vertices (2, "t8_circle_geometry")
@@ -339,7 +343,8 @@ class t8_geometry_circle: public t8_geometry_with_vertices {
  * The geometry can only be used with single tree cmeshes (unit square).
  */
 
-class t8_geometry_moving: public t8_geometry {
+struct t8_geometry_moving: public t8_geometry
+{
  public:
   /* Basic constructor that sets the dimension the name and the time pointer. */
   t8_geometry_moving (const double *time): t8_geometry (2, "t8_moving_geometry"), ptime (time)
@@ -415,7 +420,8 @@ class t8_geometry_moving: public t8_geometry {
  * in z direction.
  * Can be used with 1 tree unit cube cmesh only.
  */
-class t8_geometry_cube_zdistorted: public t8_geometry {
+struct t8_geometry_cube_zdistorted: public t8_geometry
+{
  public:
   /* Basic constructor that sets the dimension and the name. */
   t8_geometry_cube_zdistorted (): t8_geometry (3, "t8_cube_zdistorted_geometry")

example/remove/t8_example_empty_trees.cxx

@@ -25,6 +25,7 @@
 #include <t8_forest/t8_forest_types.h>
 #include <t8_cmesh/t8_cmesh_examples.h>
 #include <t8_schemes/t8_default/t8_default_cxx.hxx>
+#include <t8_geometry/t8_geometry_implementations/t8_geometry_linear_axis_aligned.h>
 #include <sc_options.h>
 #include <string>
 
@@ -49,8 +50,10 @@ t8_strip_of_quads (t8_gloidx_t num_trees, t8_gloidx_t empty_tree, const char **v
 {
 
   const double boundary_coords[12] = { 0, 0, 0, 1, 0, 0, 0, 1, 0, 1, 1, 0 };
+  const t8_geometry_c *geometry = t8_geometry_linear_axis_aligned_new (2);
 
-  t8_cmesh_t cmesh = t8_cmesh_new_hypercube_pad (T8_ECLASS_QUAD, sc_MPI_COMM_WORLD, boundary_coords, num_trees, 1, 0);
+  t8_cmesh_t cmesh
+    = t8_cmesh_new_hypercube_pad (T8_ECLASS_QUAD, sc_MPI_COMM_WORLD, boundary_coords, num_trees, 1, 0, geometry);
 
   t8_forest_t forest = t8_forest_new_uniform (cmesh, t8_scheme_new_default_cxx (), 0, 0, sc_MPI_COMM_WORLD);
 

src/t8_cmesh/t8_cmesh_examples.c

@@ -24,7 +24,9 @@
 #include <t8_cmesh/t8_cmesh_examples.h>
 #include <t8_cmesh/t8_cmesh_helpers.h>
 #include <t8_cmesh/t8_cmesh_geometry.h>
+#include <t8_geometry/t8_geometry_base.h>
 #include <t8_geometry/t8_geometry_implementations/t8_geometry_linear.h>
+#include <t8_geometry/t8_geometry_implementations/t8_geometry_linear_axis_aligned.h>
 #include <t8_geometry/t8_geometry_implementations/t8_geometry_examples.h>
 #include <t8_vec.h>
 #include <t8_mat.h>
@@ -926,12 +928,13 @@ t8_resize_box (const int dim, double *box_corners, const double *box_dir, const
  * \param [in] boundary     The boundary vertices of \a cmesh.
  * \param [in] quads_x      The number of quads along the x-axis.
  * \param [in] quads_y      The number of quads along the y-axis.
+ * \param [in] use_axis_aligned_geom Flag if cmesh uses the axis_aligned_geometry. Only available for T8_ECLASS_QUAD
  * \note each quad of \a quads_x * \a quads_y quads in \a boundary contains one
  * tree of \a eclass T8_ECLASS_QUAD or two of T8_ECLASS_TRIANGLE.
  */
 static void
 t8_cmesh_set_vertices_2D (t8_cmesh_t cmesh, const t8_eclass_t eclass, const double *boundary, const t8_locidx_t quads_x,
-                          const t8_locidx_t quads_y)
+                          const t8_locidx_t quads_y, const int use_axis_aligned_geom)
 {
   T8_ASSERT (!t8_cmesh_is_committed (cmesh));
   T8_ASSERT (eclass == T8_ECLASS_QUAD || eclass == T8_ECLASS_TRIANGLE);
@@ -982,7 +985,7 @@ t8_cmesh_set_vertices_2D (t8_cmesh_t cmesh, const t8_eclass_t eclass, const doub
   for (t8_locidx_t quad_y_id = 0; quad_y_id < quads_y; quad_y_id++) {
     for (t8_locidx_t quad_x_id = 0; quad_x_id < quads_x; quad_x_id++) {
 
-      t8_vec_axy (box, vertices, 1.0);                /* Vertex 0 */
+      memcpy (vertices, box, 3 * sizeof (double));    /* Vertex 0 */
       t8_vec_axpyz (box, box_dir, vertices + 6, 1.0); /* Vertex 2 */
 
       /* Reduce box along x axis */
@@ -991,12 +994,20 @@ t8_cmesh_set_vertices_2D (t8_cmesh_t cmesh, const t8_eclass_t eclass, const doub
 
       t8_vec_axy (box, vertices + 3, 1.0);            /* Vertex 1 */
       t8_vec_axpyz (box, box_dir, vertices + 9, 1.0); /* Vertex 3 */
+      if (use_axis_aligned_geom && eclass == T8_ECLASS_QUAD) {
+        /* Copy vertex 3 into the place of vertex 1. The box-procedure has to be done to compute 
+         * vertex 3 correctly. */
+        memcpy (vertices + 3, vertices + 9, 3 * sizeof (double));
+      }
+      else if (use_axis_aligned_geom && eclass != T8_ECLASS_QUAD) {
+        SC_ABORTF ("Axis aligned geometry is not available for eclass %s!\n", t8_eclass_to_string[eclass]);
+      }
 
       /* Map vertices of current quad on to respective trees inside. */
       if (eclass == T8_ECLASS_QUAD) {
         /* No mapping is required. */
         const t8_gloidx_t tree_id = quad_y_id * quads_x + quad_x_id;
-        t8_cmesh_set_tree_vertices (cmesh, tree_id, vertices, 4);
+        t8_cmesh_set_tree_vertices (cmesh, tree_id, vertices, use_axis_aligned_geom ? 2 : 4);
       }
       else {
         T8_ASSERT (eclass == T8_ECLASS_TRIANGLE);
@@ -1039,12 +1050,13 @@ t8_cmesh_set_vertices_2D (t8_cmesh_t cmesh, const t8_eclass_t eclass, const doub
  * \param [in] hexs_x       The number of hexs along the x-axis.
  * \param [in] hexs_y       The number of hexs along the y-axis.
  * \param [in] hexs_z       The number of hexs along the z-axis.
+ * \param [in] use_axis_aligned_geom Flag if cmesh uses the axis aligned_geometry. Only available for T8_ECLASS_QUAD
  * \note each hex of \a hexs_x * \a hexs_y * \a hexs_z hexs in \a boundary 
  * contains several trees of class \a eclass.
  */
 static void
 t8_cmesh_set_vertices_3D (t8_cmesh_t cmesh, const t8_eclass_t eclass, const double *boundary, const t8_locidx_t hexs_x,
-                          const t8_locidx_t hexs_y, const t8_locidx_t hexs_z)
+                          const t8_locidx_t hexs_y, const t8_locidx_t hexs_z, const int use_axis_aligned_geom)
 {
   T8_ASSERT (!t8_cmesh_is_committed (cmesh));
   /* x axes */
@@ -1120,8 +1132,7 @@ t8_cmesh_set_vertices_3D (t8_cmesh_t cmesh, const t8_eclass_t eclass, const doub
   for (t8_locidx_t hex_z_id = 0; hex_z_id < hexs_z; hex_z_id++) {
     for (t8_locidx_t hex_y_id = 0; hex_y_id < hexs_y; hex_y_id++) {
       for (t8_locidx_t hex_x_id = 0; hex_x_id < hexs_x; hex_x_id++) {
-
-        t8_vec_axy (box, vertices, 1.0);                     /* Vertex 0 */
+        memcpy (vertices, box, 3 * sizeof (double));         /* Vertex 0 */
         t8_vec_axpyz (box, box_dir + 12, vertices + 6, 1.0); /* Vertex 2 */
 
         /* Reduce box along z axis and face 4. */
@@ -1145,11 +1156,20 @@ t8_cmesh_set_vertices_3D (t8_cmesh_t cmesh, const t8_eclass_t eclass, const doub
         t8_vec_axy (box, vertices + 3, 1.0);                 /* Vertex 1 */
         t8_vec_axpyz (box, box_dir + 12, vertices + 9, 1.0); /* Vertex 3 */
 
+        if (use_axis_aligned_geom && eclass == T8_ECLASS_HEX) {
+          /* Copy vertex 7 into the place of vertex 1. The box-procedure has to be done to compute 
+         * vertex 7 correctly. */
+          memcpy (vertices + 3, vertices + 21, 3 * sizeof (double));
+        }
+        else if (use_axis_aligned_geom && eclass != T8_ECLASS_HEX) {
+          SC_ABORTF ("Axis aligned geometry is not available for eclass %s!\n", t8_eclass_to_string[eclass]);
+        }
+
         /* Map vertices of current hex on to respective trees inside. */
         const t8_gloidx_t hex_id = hex_z_id * hexs_y * hexs_x + hex_y_id * hexs_x + hex_x_id;
         if (eclass == T8_ECLASS_HEX) {
           /* No mapping is required. */
-          t8_cmesh_set_tree_vertices (cmesh, hex_id, vertices, 8);
+          t8_cmesh_set_tree_vertices (cmesh, hex_id, vertices, use_axis_aligned_geom ? 2 : 8);
         }
         else if (eclass == T8_ECLASS_TET) {
           const t8_gloidx_t tree_id_0 = hex_id * 6;
@@ -1236,7 +1256,7 @@ t8_cmesh_set_vertices_3D (t8_cmesh_t cmesh, const t8_eclass_t eclass, const doub
 
 t8_cmesh_t
 t8_cmesh_new_hypercube_pad (const t8_eclass_t eclass, sc_MPI_Comm comm, const double *boundary, t8_locidx_t polygons_x,
-                            t8_locidx_t polygons_y, t8_locidx_t polygons_z)
+                            t8_locidx_t polygons_y, t8_locidx_t polygons_z, const t8_geometry_c *geometry)
 {
   SC_CHECK_ABORT (eclass != T8_ECLASS_PYRAMID, "Pyramids are not yet supported.");
   const int dim = t8_eclass_to_dimension[eclass];
@@ -1257,9 +1277,9 @@ t8_cmesh_new_hypercube_pad (const t8_eclass_t eclass, sc_MPI_Comm comm, const do
   t8_cmesh_t cmesh;
   t8_cmesh_init (&cmesh);
 
-  /* We use standard linear geometry */
-  const t8_geometry_c *linear_geom = t8_geometry_linear_new (dim);
-  t8_cmesh_register_geometry (cmesh, linear_geom);
+  const int is_axis_aligned = t8_geom_get_type (geometry) == T8_GEOMETRY_TYPE_LINEAR_AXIS_ALIGNED;
+
+  t8_cmesh_register_geometry (cmesh, geometry);
 
   /* Number of trees inside each polygon of given eclass. */
   const t8_locidx_t num_trees_for_single_hypercube[T8_ECLASS_COUNT] = { 1, 1, 1, 2, 1, 6, 2, -1 };
@@ -1273,11 +1293,11 @@ t8_cmesh_new_hypercube_pad (const t8_eclass_t eclass, sc_MPI_Comm comm, const do
   /* Set the vertices of all trees. */
   if (dim == 3) {
     T8_ASSERT (eclass == T8_ECLASS_HEX || eclass == T8_ECLASS_TET || eclass == T8_ECLASS_PRISM);
-    t8_cmesh_set_vertices_3D (cmesh, eclass, boundary, polygons_x, polygons_y, polygons_z);
+    t8_cmesh_set_vertices_3D (cmesh, eclass, boundary, polygons_x, polygons_y, polygons_z, is_axis_aligned);
   }
   else if (dim == 2) {
     T8_ASSERT (eclass == T8_ECLASS_QUAD || eclass == T8_ECLASS_TRIANGLE);
-    t8_cmesh_set_vertices_2D (cmesh, eclass, boundary, polygons_x, polygons_y);
+    t8_cmesh_set_vertices_2D (cmesh, eclass, boundary, polygons_x, polygons_y, is_axis_aligned);
   }
   else if (dim == 1) {
     T8_ASSERT (eclass == T8_ECLASS_LINE);

src/t8_cmesh/t8_cmesh_examples.h

@@ -30,6 +30,7 @@
 #include <t8_cmesh.h>
 #include <p4est_connectivity.h>
 #include <p8est_connectivity.h>
+#include <t8_cmesh/t8_cmesh_geometry.h>
 
 T8_EXTERN_C_BEGIN ();
 
@@ -103,6 +104,7 @@ t8_cmesh_new_hypercube (t8_eclass_t eclass, sc_MPI_Comm comm, int do_bcast, int
  *                          Only required if \a eclass is 2D or 3D.
  * \param [in] polygons_z   The number of polygons along the z-axis.
  *                          Only required if \a eclass is 3D.
+ * \param [in] geometry     The geometry to use. If the geometry is axis_aligned only two points per tree are stored
  * \return                  A committed t8_cmesh structure with 
  *                          \a polygons_x * \a polygons_z * \a polygons_y many 
  *                          sub-hypercubes of class \a eclass.
@@ -129,7 +131,7 @@ t8_cmesh_new_hypercube (t8_eclass_t eclass, sc_MPI_Comm comm, int do_bcast, int
  */
 t8_cmesh_t
 t8_cmesh_new_hypercube_pad (const t8_eclass_t eclass, sc_MPI_Comm comm, const double *boundary, t8_locidx_t polygons_x,
-                            t8_locidx_t polygons_y, t8_locidx_t polygons_z);
+                            t8_locidx_t polygons_y, t8_locidx_t polygons_z, const t8_geometry_c *geometry);
 
 /** Hybercube with 6 Tets, 6 Prism, 4 Hex. 
  * \param [in]  comm            The mpi communicator to be used.

src/t8_geometry/t8_geometry.cxx

@@ -368,7 +368,5 @@ t8_geometry_get_type (t8_cmesh_t cmesh, t8_gloidx_t gtreeid)
   t8_geom_handler_update_tree (geom_handler, cmesh, gtreeid);
 
   /* Return the type. */
-  /* *INDENT-OFF* */
   return geom_handler->active_geometry->t8_geom_get_type ();
-  /* *INDENT-ON* */
 }

src/t8_geometry/t8_geometry_base.cxx

@@ -50,3 +50,15 @@ t8_geom_get_name (const t8_geometry_c *geom)
 
   return geom->t8_geom_get_name ();
 }
+
+/** Get the type of the geometry.
+ * \param [in]  geom  A geometry.
+ * \return            The type of \a geom.
+ */
+t8_geometry_type_t
+t8_geom_get_type (const t8_geometry_c *geom)
+{
+  T8_ASSERT (geom != NULL);
+
+  return geom->t8_geom_get_type ();
+}

src/t8_geometry/t8_geometry_base.hxx

@@ -21,7 +21,7 @@
 */
 
 /** \file t8_geometry_base.hxx
- * Implements the base pure virtual class t8_geometry which
+ * Implements the base pure virtual struct t8_geometry which
  * provides a general template for all geometries.
  */
 
@@ -30,6 +30,7 @@
 
 #include <t8.h>
 #include <t8_cmesh.h>
+#include <t8_geometry/t8_geometry.h>
 
 T8_EXTERN_C_BEGIN ();
 

src/t8_geometry/t8_geometry_implementations/t8_geometry_examples.hxx

@@ -32,7 +32,8 @@
 
 /** This geometry maps five quads to a disk.
  */
-class t8_geometry_squared_disk: public t8_geometry_with_vertices {
+struct t8_geometry_squared_disk: public t8_geometry_with_vertices
+{
  public:
   /* Basic constructor that sets the dimension and the name. */
   t8_geometry_squared_disk (): t8_geometry_with_vertices (2, "t8_squared_disk")
@@ -89,7 +90,8 @@ class t8_geometry_squared_disk: public t8_geometry_with_vertices {
 
 /** This geometry maps the faces of an oktaeder to a spherical surface.
  */
-class t8_geometry_triangulated_spherical_surface: public t8_geometry_with_vertices {
+struct t8_geometry_triangulated_spherical_surface: public t8_geometry_with_vertices
+{
  public:
   /* Basic constructor that sets the dimension and the name. */
   t8_geometry_triangulated_spherical_surface (): t8_geometry_with_vertices (2, "t8_triangulated_spherical_surface")

src/t8_geometry/t8_geometry_with_vertices.hxx

@@ -21,7 +21,7 @@
 */
 
 /** \file t8_geometry_with_vertices.hxx
- * Implements the inherited class t8_geometry_with_vertices, which can be 
+ * Implements the inherited struct t8_geometry_with_vertices, which can be 
  * used for geometries that use vertex coordinate information of the cmesh.
  */
 
@@ -35,7 +35,8 @@
 
 T8_EXTERN_C_BEGIN ();
 
-class t8_geometry_with_vertices: public t8_geometry {
+struct t8_geometry_with_vertices: public t8_geometry
+{
  public:
   /* Basic constructor that sets the dimension, the name, and the name for the attribute. */
   t8_geometry_with_vertices (int dimension, const char *name, const char *attribute_name = NULL)

test/t8_cmesh/t8_gtest_cmesh_set_join_by_vertices.cxx

@@ -27,6 +27,7 @@
 #include <t8_cmesh/t8_cmesh_examples.h>
 #include <t8_cmesh/t8_cmesh_helpers.h>
 #include <t8_cmesh/t8_cmesh_testcases.h>
+#include <t8_geometry/t8_geometry_implementations/t8_geometry_linear.hxx>
 #include <test/t8_gtest_macros.hxx>
 
 #include <p8est_geometry.h>
@@ -242,7 +243,8 @@ TEST (t8_cmesh_set_join_by_vertices, test_cmesh_set_join_by_vertices)
     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 };
 
     t8_eclass_t eclass = T8_ECLASS_HEX;
-    t8_cmesh_t cmesh = t8_cmesh_new_hypercube_pad (eclass, comm, boundary_coords, 2, 2, 2);
+    t8_geometry_c *geometry = new t8_geometry_linear (3);
+    t8_cmesh_t cmesh = t8_cmesh_new_hypercube_pad (eclass, comm, boundary_coords, 2, 2, 2, geometry);
     test_with_cmesh (cmesh);
     t8_cmesh_destroy (&cmesh);
   }