Skip to content

Commit

Permalink
Merge pull request ForestClaw#309 from hannesbrandt/feature-hanging-c…
Browse files Browse the repository at this point in the history 
…orners

Feature edge-hanging corners
  • Loading branch information
@scottaiton
scottaiton authored Feb 6, 2024
2 parents 0d2ec2d + 90e45fe commit 388c6a6
Show file tree
Hide file tree
Showing 4 changed files with 55 additions and 110 deletions.
87 changes: 48 additions & 39 deletions src/fclaw2d_convenience.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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 */
Expand All @@ -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 */
Expand All @@ -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
Expand Down
1 change: 0 additions & 1 deletion src/fclaw2d_to_3d.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
7 changes: 7 additions & 0 deletions src/fclaw3d_convenience.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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 ());
}
70 changes: 0 additions & 70 deletions src/forestclaw2d.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
{
Expand Down

0 comments on commit 388c6a6

Please sign in to comment.