Are corners necessary if corner connectivity will not be used?

[efaulhaber]

I am using p4est in a DG code, so I only need face connectivity and no corner connectivity.
In the docs, it says

The corners are only stored when they connect trees.

and

If num_corners == 0, tree_to_corner and corner_to_* arrays are set to NULL.

I figured I could just use num_corners == 0 since I don't need any corner connectivity. This worked perfectly fine until I refined one quadrant to level 4. In this case, I needed to call p4est_balance twice to obtain a balanced forest.

Is it fine to use num_corners == 0 if corner connectivity is not needed, and is this a bug?
Or did I use this incorrectly, and should I have specified the corner connectivity even if I don't use it myself?

[cburstedde]

I figured I could just use num_corners == 0 since I don't need any corner connectivity. This worked perfectly fine until I refined one quadrant to level 4. In this case, I needed to call p4est_balance twice to obtain a balanced forest.

You should be able to omit connecting corners and edges as well in your case.

Make sure you call balance only for faces, there is a parameter to the function.
It should never be necessary to call balance twice. If this still occurs, please report.

Is it fine to use num_corners == 0 if corner connectivity is not needed, and is this a bug?
Or did I use this incorrectly, and should I have specified the corner connectivity even if I don't use it myself?

You would not need corners or edges. This should all work -- please follow up with more information.

[efaulhaber]

I created a periodic rectangle connectivity, similar to p4est_connectivity_new_periodic but with num_corners = 0. Then, I refined the bottom left quadrant to level 4.
Balancing once yields this mesh:

Balancing twice yields this mesh:

This does not happen when I use p4est_connectivity_new_periodic (with corner connectivity).

Here's a quick and dirty MWE to demonstrate this behavior. It yields the following output showing that the first balancing produces 25 quadrants while the second balancing creates 3 more.

[p4est] Into p4est_new with min quadrants 0 level 0 uniform 1
[p4est]  New p4est with 1 trees on 1 processors
[p4est] Done p4est_new with 1 total quadrants
[p4est] Into p4est_refine with 1 total quadrants, allowed level 29
[p4est] Done p4est_refine with 13 total quadrants
[p4est] Into p4est_balance FACE with 13 total quadrants
[p4est] Done p4est_balance with 25 total quadrants
[p4est] Into p4est_balance FACE with 25 total quadrants
[p4est] Done p4est_balance with 28 total quadrants

MWE code

#ifndef P4_TO_P8

#include <p4est_vtk.h>
#include <p4est_extended.h>
#include "hw32.h"


static int
refine_fn (p4est_t * p4est, p4est_topidx_t which_tree,
           p4est_quadrant_t * quadrant)
{
  // Refine bottom left quadrant to level 4
  if (quadrant->x == 0 && quadrant->y == 0 && quadrant->level < 4) {
    return 1;
  }

  return 0;
}


p4est_connectivity_t* p4est_connectivity_new_test(void)
{
  // This is exactly copied from p4est_connectivity_new_periodic
  // but with num_corners = 0
  const p4est_topidx_t num_vertices = P4EST_CHILDREN;
  const p4est_topidx_t num_trees = 1;
  const p4est_topidx_t num_corners = 0;
  const double        vertices[P4EST_CHILDREN * 3] = {
    0, 0, 0,
    1, 0, 0,
    0, 1, 0,
    1, 1, 0,
  };
  const p4est_topidx_t tree_to_vertex[1 * P4EST_CHILDREN] = {
    0, 1, 2, 3,
  };

  const p4est_topidx_t tree_to_tree[1 * P4EST_FACES] = {
    0, 0, 0, 0,
  };
  const int8_t        tree_to_face[1 * P4EST_FACES] = {
    1, 0, 3, 2,
  };
  // const p4est_topidx_t tree_to_corner[1 * P4EST_CHILDREN] = {
  //   0, 0, 0, 0,
  // };
  const p4est_topidx_t* tree_to_corner = NULL;

  const p4est_topidx_t ctt_offset[1] = {
    0, //P4EST_CHILDREN,
  };
  // const p4est_topidx_t corner_to_tree[P4EST_CHILDREN] = {
  //   0, 0, 0, 0,
  // };
  const p4est_topidx_t* corner_to_tree = NULL;

  // const int8_t        corner_to_corner[P4EST_CHILDREN] = {
  //   0, 1, 2, 3,
  // };
  const int8_t* corner_to_corner = NULL;

  return p4est_connectivity_new_copy (num_vertices, num_trees,
                                      num_corners, vertices, tree_to_vertex,
                                      tree_to_tree, tree_to_face,
                                      tree_to_corner, ctt_offset,
                                      corner_to_tree, corner_to_corner);
}


int main(int argc, char **argv)
{
  int                 mpiret;
  sc_MPI_Comm         mpicomm;
  p4est_t            *p4est;
  p4est_connectivity_t *conn;

  /* Initialize MPI; see sc_mpi.h.
   * If configure --enable-mpi is given these are true MPI calls.
   * Else these are dummy functions that simulate a single-processor run. */
  mpiret = sc_MPI_Init (&argc, &argv);
  SC_CHECK_MPI (mpiret);
  mpicomm = sc_MPI_COMM_WORLD;

  /* These functions are optional.  If called they store the MPI rank as a
   * static variable so subsequent global p4est log messages are only issued
   * from processor zero.  Here we turn off most of the logging; see sc.h. */
  sc_init (mpicomm, 1, 1, NULL, SC_LP_ESSENTIAL);
  p4est_init (NULL, SC_LP_PRODUCTION);

  conn = p4est_connectivity_new_test();
  p4est = p4est_new_ext (mpicomm, conn, 0, 0, 1, 0, NULL, NULL);

  p4est_tree_t *tree = (p4est_tree_t *) (p4est->trees->array);
  sc_array_t *quad_array = &tree->quadrants;

  // Refine bottom left element to level 4
  p4est_refine (p4est, 1, refine_fn, NULL);

  p4est_balance(p4est, P4EST_CONNECT_FACE, NULL);
  p4est_balance(p4est, P4EST_CONNECT_FACE, NULL);

  /* Destroy the p4est and the connectivity structure. */
  p4est_destroy (p4est);
  p4est_connectivity_destroy (conn);

  /* Verify that allocations internal to p4est and sc do not leak memory.
   * This should be called if sc_init () has been called earlier. */
  sc_finalize ();

  /* This is standard MPI programs.  Without --enable-mpi, this is a dummy. */
  mpiret = sc_MPI_Finalize ();
  SC_CHECK_MPI (mpiret);
  return 0;
}

#endif

[cburstedde]

Thanks, this is really helpful! I'll copy @tisaac, since it may be related to one of the later evolutions of the balance algorithm.

(It may also be related to #101, which is so far not part of the official branches but ideally should be in the future.)

[tisaac]

I'm sorry about missing this issue last year!

My interpretation is that omitting corners from this periodic mesh should be considered incorrect. When you omit the corners data structures, you are telling the topology traversal routines, for example, "the only tree-corners neighboring (tree 0, corner 0) are the ones that can be found across faces (and edges in 3D)". In this case, this gets interpreted as "the only tree-corners neighboring (tree 0, corner 0) are (tree 0, corner 1) [the neighbor across the x face] and (tree 0, corner 2) [the neighbor across the y face]." That means that in balance the affect of (tree 0, corner 0) on (tree 0, corner 3) are not considered.

• We should update the documentation to make it clearer that corners can only be omitted if all corners neighbors can be found (non-recursively) in the face and edge neighbors
• At this point I consider corner-omission and edge-omission to be an unsafe optimization with very limited upside in terms of the size of the connectivity and we should discourage users from doing it.

[cburstedde]

Let me come back to this. Thanks again to @efaulhaber on your pictures!

Your top picture, even without corner information, should be face balancing the top right quadrant one level deeper. This is apparently not happening. So I'd suspect that corner information has some effect backwards on face connection, which is somewhat against the usual flow of information in p4est.

This may be related to #238, which is being investigated.

[cburstedde]

After discussing the issue some more, we basically concur with @tisaac's view. Technically, even for face- or edge-only balance, we require (parts of) the full 3x3 neighborhood of a quadrant (3x3x3 of an octant) to be known. Omitting corner (or edge) connections removes certain blocks from the neighborhood at tree boundaries, turning it into a cross-shape instead of a full cube, which leads to incomplete balance. The principle that should be maintained can be stated as follows:

Connectivity completeness: If a 3D connectivity contains natural connections between trees that are edge neighbors without being face neighbors, these edges shall be encoded explicitly in the connectivity structure. If a connectivity implies natural connections between trees that are corner neighbors without being face (or edge) neighbors, these corners shall be encoded explicitly.

A connectivity can be

• valid if it is free of contradictions, even when being incomplete according to the above principle. The function p4est_connectivity_is_valid checks for some, yet not all, of potential contradictions.
• complete if it is valid and fully connected at edges and corners as implied by iterating through the encoded tree neighbor faces. While a connectivity must always be valid, it may be incomplete for specific use cases.

What we need to do before closing the issue, is

• as stated by @tisaac updating the documentation in p4est_connectivity.h,
• finish and document the existing function p4est_connectivity_complete,
• provide the new function p4est_connectivity_is_complete.

Any further comments are very welcome.

[cburstedde]

The documentation of the 2D connectivity has been extended in #236.