Improved quadrangulated disk mesh cmesh example [6/n]

example/cmesh/t8_cmesh_geometry_examples.cxx

@@ -104,13 +104,13 @@ main (int argc, char **argv)
    */
 
   {
-    const char *prefix_cmesh = "t8_squared_disk_cmesh";
-    const char *prefix_forest = "t8_squared_disk_forest";
+    const char *prefix_cmesh = "t8_quadrangulated_disk_cmesh";
+    const char *prefix_forest = "t8_quadrangulated_disk_forest";
 
     const int uniform_level = 5;
     const double radius = 1.0;
 
-    t8_cmesh_t cmesh = t8_cmesh_new_squared_disk (radius, comm);
+    t8_cmesh_t cmesh = t8_cmesh_new_quadrangulated_disk (radius, comm);
 
     t8_forest_t forest = t8_forest_new_uniform (cmesh, t8_scheme_new_default_cxx (), uniform_level, 0, comm);
 
@@ -256,6 +256,26 @@ main (int argc, char **argv)
     t8_forest_unref (&forest);
   }
 
+  {
+    const char *prefix_cmesh = "t8_cubed_sphere_cmesh";
+    const char *prefix_forest = "t8_cubed_sphere_forest";
+
+    const int uniform_level = 2;
+    const double radius = 1.0;
+
+    t8_cmesh_t cmesh = t8_cmesh_new_cubed_sphere (radius, comm);
+
+    t8_forest_t forest = t8_forest_new_uniform (cmesh, t8_scheme_new_default_cxx (), uniform_level, 0, comm);
+
+    t8_cmesh_vtk_write_file (cmesh, prefix_cmesh);
+    t8_global_productionf ("Wrote %s.\n", prefix_cmesh);
+
+    t8_write_forest_to_vtu (forest, prefix_forest);
+    t8_global_productionf ("Wrote %s.\n\n", prefix_forest);
+
+    t8_forest_unref (&forest);
+  }
+
   t8_global_productionf ("Done!\n");
 
   /* Finalize the sc library */

src/t8_cmesh/t8_cmesh_examples.c

@@ -2720,29 +2720,36 @@ t8_cmesh_new_row_of_cubes (t8_locidx_t num_trees, const int set_attributes, cons
 }
 
 t8_cmesh_t
-t8_cmesh_new_squared_disk (const double radius, sc_MPI_Comm comm)
+t8_cmesh_new_quadrangulated_disk (const double radius, sc_MPI_Comm comm)
 {
   /* Initialization of the mesh */
   t8_cmesh_t cmesh;
   t8_cmesh_init (&cmesh);
 
-  const double ri = 0.5 * radius;
+  const double ri = 0.6 * radius;
   const double ro = radius;
 
   const double xi = ri / M_SQRT2;
-  const double yi = ri / M_SQRT2;
+  const double yi = xi;
 
   const double xo = ro / M_SQRT2;
-  const double yo = ro / M_SQRT2;
+  const double yo = xo;
 
-  const int ntrees = 5; /* Number of cmesh elements resp. trees. */
-  const int nverts = 4; /* Number of vertices per cmesh element. */
+  const int nquads = 3;  /* Number of quads in the upper-right quadrant. */
+  const int nyturns = 4; /* Number of turns around y-axis. */
+
+  const int ntrees = nyturns * nquads; /* Number of cmesh elements resp. trees. */
+  const int nverts = t8_eclass_num_vertices[T8_ECLASS_QUAD];
+
+  /* Fine tuning parameter to expand the center squares a bit for more equal
+   * element sizes. */
+  const double s = 1.2;
 
   /* Arrays for the face connectivity computations via vertices. */
   double all_verts[ntrees * T8_ECLASS_MAX_CORNERS * T8_ECLASS_MAX_DIM];
   t8_eclass_t all_eclasses[ntrees];
 
-  t8_geometry_c *geometry = t8_geometry_squared_disk_new ();
+  t8_geometry_c *geometry = t8_geometry_quadrangulated_disk_new ();
   t8_cmesh_register_geometry (cmesh, geometry);
 
   /* Defitition of the tree class. */
@@ -2751,44 +2758,42 @@ t8_cmesh_new_squared_disk (const double radius, sc_MPI_Comm comm)
     all_eclasses[itree] = T8_ECLASS_QUAD;
   }
 
-  /* Central quad. */
-  {
-    const double vertices[4][3] = { { -xi, -yi, 0.0 }, { xi, -yi, 0.0 }, { -xi, yi, 0.0 }, { xi, yi, 0.0 } };
-
-    t8_cmesh_set_tree_vertices (cmesh, 0, (double *) vertices, 4);
-
-    /* itree = 0; */
-    for (int ivert = 0; ivert < nverts; ivert++) {
-      for (int icoord = 0; icoord < T8_ECLASS_MAX_DIM; icoord++) {
-        all_verts[T8_3D_TO_1D (ntrees, T8_ECLASS_MAX_CORNERS, T8_ECLASS_MAX_DIM, 0, ivert, icoord)]
-          = vertices[ivert][icoord];
-      }
-    }
-  }
-
-  /* Four quads framing the central quad. */
-  {
-    const double vertices[4][3] = { { -xi, yi, 0.0 }, { xi, yi, 0.0 }, { -xo, yo, 0.0 }, { xo, yo, 0.0 } };
+  /* Vertices of upper right quarter of the disk. */
+  const double vertices_mid[4][3] = { { 0.0, 0.0, 0.0 }, { s * xi, 0.0, 0.0 }, { 0.0, s * yi, 0.0 }, { xi, yi, 0.0 } };
+  const double vertices_top[4][3] = { { 0.0, s * yi, 0.0 }, { xi, yi, 0.0 }, { 0.0, yo, 0.0 }, { xo, yo, 0.0 } };
+  const double vertices_bot[4][3] = { { s * xi, 0.0, 0.0 }, { xi, yi, 0.0 }, { xo, 0.0, 0.0 }, { xo, yo, 0.0 } };
 
-    for (int itree = 1; itree < ntrees; itree++) {
-      double rot_mat[3][3];
-      double rot_vertices[4][3];
+  int itree = 0;
+  for (int iturn = 0; iturn < 4; iturn++) {
+    double rot_mat[3][3];
+    double rot_vertices_mid[4][3];
+    double rot_vertices_top[4][3];
+    double rot_vertices_bot[4][3];
 
-      t8_mat_init_zrot (rot_mat, (itree - 1) * 0.5 * M_PI);
+    t8_mat_init_zrot (rot_mat, -iturn * 0.5 * M_PI);
 
-      for (int i = 0; i < 4; i++) {
-        t8_mat_mult_vec (rot_mat, &(vertices[i][0]), &(rot_vertices[i][0]));
-      }
+    for (int i = 0; i < 4; i++) {
+      t8_mat_mult_vec (rot_mat, &(vertices_mid[i][0]), &(rot_vertices_mid[i][0]));
+      t8_mat_mult_vec (rot_mat, &(vertices_top[i][0]), &(rot_vertices_top[i][0]));
+      t8_mat_mult_vec (rot_mat, &(vertices_bot[i][0]), &(rot_vertices_bot[i][0]));
+    }
 
-      t8_cmesh_set_tree_vertices (cmesh, itree, (double *) rot_vertices, 4);
+    t8_cmesh_set_tree_vertices (cmesh, itree + 0, (double *) rot_vertices_mid, 4);
+    t8_cmesh_set_tree_vertices (cmesh, itree + 1, (double *) rot_vertices_top, 4);
+    t8_cmesh_set_tree_vertices (cmesh, itree + 2, (double *) rot_vertices_bot, 4);
 
-      for (int ivert = 0; ivert < nverts; ivert++) {
-        for (int icoord = 0; icoord < T8_ECLASS_MAX_DIM; icoord++) {
-          all_verts[T8_3D_TO_1D (ntrees, T8_ECLASS_MAX_CORNERS, T8_ECLASS_MAX_DIM, itree, ivert, icoord)]
-            = rot_vertices[ivert][icoord];
-        }
+    for (int ivert = 0; ivert < nverts; ivert++) {
+      for (int icoord = 0; icoord < T8_ECLASS_MAX_DIM; icoord++) {
+        all_verts[T8_3D_TO_1D (ntrees, T8_ECLASS_MAX_CORNERS, T8_ECLASS_MAX_DIM, itree + 0, ivert, icoord)]
+          = rot_vertices_mid[ivert][icoord];
+        all_verts[T8_3D_TO_1D (ntrees, T8_ECLASS_MAX_CORNERS, T8_ECLASS_MAX_DIM, itree + 1, ivert, icoord)]
+          = rot_vertices_top[ivert][icoord];
+        all_verts[T8_3D_TO_1D (ntrees, T8_ECLASS_MAX_CORNERS, T8_ECLASS_MAX_DIM, itree + 2, ivert, icoord)]
+          = rot_vertices_bot[ivert][icoord];
       }
     }
+
+    itree = itree + nquads;
   }
 
   /* Face connectivity. */
@@ -3219,3 +3224,120 @@ t8_cmesh_new_cubed_spherical_shell (const double inner_radius, const double shel
                                        t8_cmesh_new_quadrangulated_spherical_surface, inner_radius, shell_thickness,
                                        num_levels, num_layers, comm);
 }
+
+t8_cmesh_t
+t8_cmesh_new_cubed_sphere (const double radius, sc_MPI_Comm comm)
+{
+  /* Initialization of the mesh. */
+  t8_cmesh_t cmesh;
+  t8_cmesh_init (&cmesh);
+
+  const double ri = 0.6 * radius;
+  const double ro = radius;
+
+  const double SQRT3 = 1.7320508075688772;
+  const double xi = ri / SQRT3;
+  const double yi = xi;
+  const double zi = xi;
+
+  const double xo = ro / SQRT3;
+  const double yo = xo;
+  const double zo = xo;
+
+  const int nhexs = 4;   /* Number of hexs in the front-upper-right octant. */
+  const int nzturns = 4; /* Number of turns around z-axis. */
+  const int nyturns = 2; /* Number of turns around y-axis. */
+
+  const int ntrees = nyturns * nzturns * nhexs; /* Number of cmesh elements resp. trees. */
+  const int nverts = t8_eclass_num_vertices[T8_ECLASS_HEX];
+
+  /* Fine tuning parameter to expand the center hex a bit for more equal
+   * element sizes. */
+  const double s = 1.2;
+
+  /* Arrays for the face connectivity computations via vertices. */
+  double all_verts[ntrees * T8_ECLASS_MAX_CORNERS * T8_ECLASS_MAX_DIM];
+  t8_eclass_t all_eclasses[ntrees];
+
+  t8_geometry_c *geometry = t8_geometry_cubed_sphere_new ();
+  t8_cmesh_register_geometry (cmesh, geometry);
+
+  /* Defitition of the tree class. */
+  for (int itree = 0; itree < ntrees; itree++) {
+    t8_cmesh_set_tree_class (cmesh, itree, T8_ECLASS_HEX);
+    all_eclasses[itree] = T8_ECLASS_HEX;
+  }
+
+  const double vertices_mid[8][3]
+    = { { 0.0, 0.0, 0.0 },    { s * xi, 0.0, 0.0 }, { 0.0, s * yi, 0.0 }, { xi, yi, 0.0 },
+        { 0.0, 0.0, s * zi }, { xi, 0.0, zi },      { 0.0, yi, zi },      { xi, yi, zi } };
+  const double vertices_top[8][3] = { { 0.0, s * yi, 0.0 }, { xi, yi, 0.0 }, { 0.0, yi, zi }, { xi, yi, zi },
+                                      { 0.0, yo, 0.0 },     { xo, yo, 0.0 }, { 0.0, yo, zo }, { xo, yo, zo } };
+  const double vertices_bot[8][3] = { { s * xi, 0.0, 0.0 }, { xi, yi, 0.0 }, { xi, 0.0, zi }, { xi, yi, zi },
+                                      { xo, 0.0, 0.0 },     { xo, yo, 0.0 }, { xo, 0.0, zo }, { xo, yo, zo } };
+  const double vertices_zen[8][3] = { { 0.0, 0.0, s * zi }, { xi, 0.0, zi }, { 0.0, yi, zi }, { xi, yi, zi },
+                                      { 0.0, 0.0, zo },     { xo, 0.0, zo }, { 0.0, yo, zo }, { xo, yo, zo } };
+
+  int itree = 0;
+  for (int yturn = 0; yturn < nyturns; yturn++) {
+    double yrot_mat[3][3];
+    t8_mat_init_yrot (yrot_mat, yturn * M_PI);
+
+    for (int zturn = 0; zturn < nzturns; zturn++) {
+      double zrot_mat[3][3];
+
+      double tmp_vertices_mid[8][3];
+      double tmp_vertices_top[8][3];
+      double tmp_vertices_bot[8][3];
+      double tmp_vertices_zen[8][3];
+
+      double rot_vertices_mid[8][3];
+      double rot_vertices_top[8][3];
+      double rot_vertices_bot[8][3];
+      double rot_vertices_zen[8][3];
+
+      t8_mat_init_zrot (zrot_mat, -zturn * 0.5 * M_PI);
+
+      for (int i = 0; i < 8; i++) {
+        t8_mat_mult_vec (zrot_mat, &(vertices_mid[i][0]), &(tmp_vertices_mid[i][0]));
+        t8_mat_mult_vec (zrot_mat, &(vertices_top[i][0]), &(tmp_vertices_top[i][0]));
+        t8_mat_mult_vec (zrot_mat, &(vertices_bot[i][0]), &(tmp_vertices_bot[i][0]));
+        t8_mat_mult_vec (zrot_mat, &(vertices_zen[i][0]), &(tmp_vertices_zen[i][0]));
+      }
+
+      for (int i = 0; i < 8; i++) {
+        t8_mat_mult_vec (yrot_mat, &(tmp_vertices_mid[i][0]), &(rot_vertices_mid[i][0]));
+        t8_mat_mult_vec (yrot_mat, &(tmp_vertices_top[i][0]), &(rot_vertices_top[i][0]));
+        t8_mat_mult_vec (yrot_mat, &(tmp_vertices_bot[i][0]), &(rot_vertices_bot[i][0]));
+        t8_mat_mult_vec (yrot_mat, &(tmp_vertices_zen[i][0]), &(rot_vertices_zen[i][0]));
+      }
+
+      t8_cmesh_set_tree_vertices (cmesh, itree + 0, (double *) rot_vertices_mid, 8);
+      t8_cmesh_set_tree_vertices (cmesh, itree + 1, (double *) rot_vertices_top, 8);
+      t8_cmesh_set_tree_vertices (cmesh, itree + 2, (double *) rot_vertices_bot, 8);
+      t8_cmesh_set_tree_vertices (cmesh, itree + 3, (double *) rot_vertices_zen, 8);
+
+      for (int ivert = 0; ivert < nverts; ivert++) {
+        for (int icoord = 0; icoord < T8_ECLASS_MAX_DIM; icoord++) {
+          all_verts[T8_3D_TO_1D (ntrees, T8_ECLASS_MAX_CORNERS, T8_ECLASS_MAX_DIM, itree + 0, ivert, icoord)]
+            = rot_vertices_mid[ivert][icoord];
+          all_verts[T8_3D_TO_1D (ntrees, T8_ECLASS_MAX_CORNERS, T8_ECLASS_MAX_DIM, itree + 1, ivert, icoord)]
+            = rot_vertices_top[ivert][icoord];
+          all_verts[T8_3D_TO_1D (ntrees, T8_ECLASS_MAX_CORNERS, T8_ECLASS_MAX_DIM, itree + 2, ivert, icoord)]
+            = rot_vertices_bot[ivert][icoord];
+          all_verts[T8_3D_TO_1D (ntrees, T8_ECLASS_MAX_CORNERS, T8_ECLASS_MAX_DIM, itree + 2, ivert, icoord)]
+            = rot_vertices_zen[ivert][icoord];
+        }
+      }
+
+      itree = itree + nhexs;
+    }
+  }
+
+  /* Face connectivity. */
+  t8_cmesh_set_join_by_vertices (cmesh, ntrees, all_eclasses, all_verts, NULL, 0);
+
+  /* Commit the mesh */
+  t8_cmesh_commit (cmesh, comm);
+  return cmesh;
+}

src/t8_cmesh/t8_cmesh_examples.h

@@ -313,13 +313,13 @@ t8_cmesh_new_long_brick_pyramid (sc_MPI_Comm comm, int num_cubes);
 t8_cmesh_t
 t8_cmesh_new_row_of_cubes (t8_locidx_t num_trees, const int set_attributes, const int do_partition, sc_MPI_Comm comm);
 
-/** Construct a squared disk of given radius.
+/** Construct a quadrangulated disk of given radius.
  * \param [in] radius        Radius of the sphere.
  * \param [in] comm          The MPI communicator used to commit the cmesh
  * \return                   A cmesh representing the spherical surface.
  */
 t8_cmesh_t
-t8_cmesh_new_squared_disk (const double radius, sc_MPI_Comm comm);
+t8_cmesh_new_quadrangulated_disk (const double radius, sc_MPI_Comm comm);
 
 /** Construct a triangulated spherical surface of given radius: octahedron version.
  * \param [in] radius        Radius of the sphere.
@@ -382,6 +382,14 @@ t8_cmesh_t
 t8_cmesh_new_cubed_spherical_shell (const double inner_radius, const double shell_thickness, const int num_levels,
                                     const int num_layers, sc_MPI_Comm comm);
 
+/** Construct a cubed sphere of given radius.
+ * \param [in] inner_radius       Radius of the inner side of the shell.
+ * \param [in] comm               The MPI communicator used to commit the cmesh
+ * \return                        A cmesh representing the spherical surface.
+ */
+t8_cmesh_t
+t8_cmesh_new_cubed_sphere (const double radius, sc_MPI_Comm comm);
+
 T8_EXTERN_C_END ();
 
 #endif /* !T8_CMESH_EXAMPLES */