cburstedde / p4est

The "p4est" forest-of-octrees library
www.p4est.org/
GNU General Public License v2.0
251 stars 114 forks source link

Are corners necessary if corner connectivity will not be used? #112

Open efaulhaber opened 3 years ago

efaulhaber commented 3 years ago

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 cornerto* 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 commented 3 years ago

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 commented 3 years ago

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:

grafik

Balancing twice yields this mesh:

grafik

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 ```C #ifndef P4_TO_P8 #include #include #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 commented 3 years ago

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 commented 2 years ago

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.

cburstedde commented 1 year ago

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 commented 1 year ago

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

What we need to do before closing the issue, is

Any further comments are very welcome.

cburstedde commented 11 months ago

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