Feature edge-hanging corners

src/fclaw2d_convenience.c

@@ -324,16 +324,24 @@ fclaw2d_domain_new (p4est_wrap_t * wrap, sc_keyvalue_t * attributes)
 }
 
 fclaw2d_domain_t *
-fclaw2d_domain_new_p4est (p4est_t *p4est)
+fclaw2d_domain_new_p4est (p4est_t * p4est)
 {
     FCLAW_ASSERT (p4est != NULL);
     FCLAW_ASSERT (p4est->user_pointer == NULL);
 
     p4est_wrap_t *wrap;
-
-    /* create p4est_wrap from the given p4est */
-    wrap = p4est_wrap_new_p4est (p4est, 0, P4EST_CONNECT_FULL, NULL, NULL);
-
+    p4est_wrap_params_t params;
+
+    /* wrap the p4est with edgehanging corners neighbors enabled in the mesh */
+    p4est_wrap_params_init (&params);
+    params.hollow = 0;
+    params.mesh_params.btype = P4EST_CONNECT_FULL;
+    params.mesh_params.compute_level_lists = 1;
+    params.mesh_params.compute_tree_index = 1;
+#ifdef P4_TO_P8
+    params.mesh_params.edgehanging_corners = 1;
+#endif
+    wrap = p4est_wrap_new_p4est_params (p4est, &params);
     FCLAW_ASSERT (wrap->p4est->data_size == 0);
 
     /* attributes of the created domain is initialized by sc_keyvalue_new */
@@ -352,52 +360,43 @@ fclaw2d_check_initial_level (sc_MPI_Comm mpicomm, int initial_level)
                      "Initial level %d too fine for p4est", initial_level);
 }
 
+#ifndef P4_TO_P8
 fclaw2d_domain_t *
 fclaw2d_domain_new_unitsquare (sc_MPI_Comm mpicomm, int initial_level)
 {
-    fclaw2d_check_initial_level (mpicomm, initial_level);
-    return fclaw2d_domain_new (p4est_wrap_new_unitsquare (mpicomm,
-                                                          initial_level),
-                               NULL);
+    return fclaw2d_domain_new_conn (mpicomm, initial_level,
+                                    p4est_connectivity_new_unitsquare ());
 }
 
-#ifndef P4_TO_P8
-
 fclaw2d_domain_t *
 fclaw2d_domain_new_torus (sc_MPI_Comm mpicomm, int initial_level)
 {
-    fclaw2d_check_initial_level (mpicomm, initial_level);
-    return
-        fclaw2d_domain_new (p4est_wrap_new_periodic (mpicomm, initial_level),
-                            NULL);
+    return fclaw2d_domain_new_conn (mpicomm, initial_level,
+                                    p4est_connectivity_new_periodic ());
 }
 
 fclaw2d_domain_t *
 fclaw2d_domain_new_twosphere (sc_MPI_Comm mpicomm, int initial_level)
 {
-    fclaw2d_check_initial_level (mpicomm, initial_level);
-    return
-        fclaw2d_domain_new (p4est_wrap_new_pillow (mpicomm, initial_level),
-                            NULL);
+    return fclaw2d_domain_new_conn (mpicomm, initial_level,
+                                    p4est_connectivity_new_pillow ());
 }
 
 fclaw2d_domain_t *
 fclaw2d_domain_new_cubedsphere (sc_MPI_Comm mpicomm, int initial_level)
 {
-    fclaw2d_check_initial_level (mpicomm, initial_level);
-    return fclaw2d_domain_new (p4est_wrap_new_cubed (mpicomm, initial_level),
-                               NULL);
+    return fclaw2d_domain_new_conn (mpicomm, initial_level,
+                                    p4est_connectivity_new_cubed ());
 }
 
 fclaw2d_domain_t *
 fclaw2d_domain_new_disk (sc_MPI_Comm mpicomm,
                          int periodic_in_x, int periodic_in_y,
                          int initial_level)
 {
-    fclaw2d_check_initial_level (mpicomm, initial_level);
-    return fclaw2d_domain_new
-        (p4est_wrap_new_disk (mpicomm, periodic_in_x, periodic_in_y,
-                              initial_level), NULL);
+    return fclaw2d_domain_new_conn (mpicomm, initial_level,
+                                    p4est_connectivity_new_disk (periodic_in_x,
+                                                                 periodic_in_y));
 }
 
 #endif /* P4_TO_P8 */
@@ -414,33 +413,43 @@ fclaw2d_domain_new_brick (sc_MPI_Comm mpicomm,
 #endif
                           int initial_level)
 {
-    p4est_wrap_t *wrap;
-
-    fclaw2d_check_initial_level (mpicomm, initial_level);
-    wrap = p4est_wrap_new_brick (mpicomm, blocks_in_x, blocks_in_y,
+    return fclaw2d_domain_new_conn (mpicomm, initial_level,
+                                    p4est_connectivity_new_brick (blocks_in_x,
+                                                                  blocks_in_y,
 #ifdef P4_TO_P8
-                                 blocks_in_z,
+                                                                  blocks_in_z,
 #endif
-                                 periodic_in_x, periodic_in_y,
+                                                                  periodic_in_x,
+                                                                  periodic_in_y
 #ifdef P4_TO_P8
-                                 periodic_in_z,
+                                                                  ,
+                                                                  periodic_in_z
 #endif
-                                 initial_level);
-    return fclaw2d_domain_new (wrap, NULL);
+                                    ));
 }
 
 fclaw2d_domain_t *
 fclaw2d_domain_new_conn (sc_MPI_Comm mpicomm, int initial_level,
                          p4est_connectivity_t * conn)
 {
     p4est_wrap_t *wrap;
-    fclaw2d_domain_t *domain;
+    p4est_wrap_params_t params;
 
     fclaw2d_check_initial_level (mpicomm, initial_level);
-    wrap = p4est_wrap_new_conn (mpicomm, conn, initial_level);
-    domain = fclaw2d_domain_new (wrap, NULL);
 
-    return domain;
+    /* Wrap the connectivity. This does the same as p4est_wrap_new_conn,
+     * except for setting edgehanging_corners to 1. */
+    p4est_wrap_params_init (&params);
+    params.hollow = 0;
+    params.mesh_params.btype = P4EST_CONNECT_FULL;
+    params.mesh_params.compute_level_lists = 1;
+    params.mesh_params.compute_tree_index = 1;
+#ifdef P4_TO_P8
+    params.mesh_params.edgehanging_corners = 1;
+#endif
+    wrap = p4est_wrap_new_params (mpicomm, conn, initial_level, &params);
+
+    return fclaw2d_domain_new (wrap, NULL);
 }
 
 void

src/fclaw2d_to_3d.h

@@ -124,7 +124,6 @@ OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 #define fclaw2d_file_read_array         fclaw3d_file_read_array
 #define fclaw2d_file_error_string       fclaw3d_file_error_string
 #define fclaw2d_file_close              fclaw3d_file_close
-#define fclaw2d_domain_new_unitsquare   fclaw3d_domain_new_unitcube
 #define fclaw2d_domain_new_brick        fclaw3d_domain_new_brick
 #define fclaw2d_domain_new_conn         fclaw3d_domain_new_conn
 #define fclaw2d_domain_new_p4est        fclaw3d_domain_new_p8est

src/fclaw3d_convenience.c

@@ -25,3 +25,10 @@ OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 
 #include <fclaw2d_to_3d.h>
 #include "fclaw2d_convenience.c"
+
+fclaw2d_domain_t *
+fclaw3d_domain_new_unitcube (sc_MPI_Comm mpicomm, int initial_level)
+{
+    return fclaw2d_domain_new_conn (mpicomm, initial_level,
+                                    p8est_connectivity_new_unitcube ());
+}

src/forestclaw2d.c

@@ -1052,76 +1052,6 @@ fclaw2d_patch_corner_neighbors (fclaw2d_domain_t * domain,
             *rcorner = cornerno ^ (P4EST_CHILDREN - 1);
         }
     }
-#ifdef P4_TO_P8
-    else if(qid == -1)
-    {
-        /* workaround for hanging edge corners */
-        fclaw3d_block_t *block = &domain->blocks[blockno];
-        fclaw3d_patch_t *patch = &block->patches[patchno];
-        int childid = fclaw3d_patch_childid (patch);
-        int edge;
-        int face;
-        const int* corner_edges = p8est_corner_edges[cornerno];
-        const int* corner_faces = p8est_corner_faces[cornerno];
-
-        switch (childid ^ cornerno)
-        {
-        case 1:
-            edge = corner_edges[0]; /* hanging on edge parallel to x */
-            face = corner_faces[0];
-            break;
-        case 2:
-            edge = corner_edges[1]; /* hanging on edge parallel to y */
-            face = corner_faces[1];
-            break;
-        case 4:
-            edge = corner_edges[2]; /* hanging on edge parallel to z */
-            face = corner_faces[2];
-            break;
-        default:
-            /* 3, 5, 6 is hanging on face; 0 or 7 is not hanging at all */
-            edge = -1;
-        }
-
-        if(edge != -1)
-        {
-            /* traverse across same-size face neighbor */
-            p4est_locidx_t f_qid = mesh->quad_to_quad[P4EST_FACES*local_num+face];
-            int v = mesh->quad_to_face[P4EST_FACES*local_num+face];
-
-            /* In the hanging edge case, the face neighboring quadrant we are trying
-               to traverse will need to be local. If there is ever a partitioning
-               where this isn't the case, then the workaround can't be performed. */
-            if(f_qid >= 0 && f_qid < mesh->local_num_quadrants && v >= 0 && v <= 23)
-            {
-                p4est_locidx_t qte = mesh->quad_to_edge[P8EST_EDGES * f_qid + edge];
-                if (qte >= 0 && qte < mesh->local_num_quadrants + mesh->ghost_num_quadrants)
-                {
-                    /* same size, same tree */
-                    qid = qte;
-                    *rcorner = cornerno ^ (P8EST_CHILDREN - 1);
-                }
-                else if(qte >= 0)
-                {
-                    /* possibly same-size different block */
-                    qte -= mesh->local_num_quadrants + mesh->ghost_num_quadrants;
-                    p4est_locidx_t cstart = fclaw2d_array_index_locidx (mesh->edge_offset, qte);
-                    p4est_locidx_t cend = fclaw2d_array_index_locidx (mesh->edge_offset, qte + 1);
-                    int8_t e = *(int8_t *) sc_array_index_int (mesh->edge_edge, (int) cstart);
-                    if(cstart + 1 == cend && e >= 0 && e <= 24)
-                    {
-                        /* TODO this can be rotated 
-                           check that it is not for now */
-                        FCLAW_ASSERT(((e%12)^3) == edge);
-                        /* same-size different block */
-                        qid = fclaw2d_array_index_locidx (mesh->edge_quad, (int) cstart);
-                        *rcorner = cornerno ^ (P8EST_CHILDREN - 1);
-                    }
-                }
-            }
-        }
-    }
-#endif
 
     if (qid < 0)
     {