diff --git a/src/p4est_extended.h b/src/p4est_extended.h index 62046e6ff..3ee30f7d7 100644 --- a/src/p4est_extended.h +++ b/src/p4est_extended.h @@ -382,6 +382,8 @@ p4est_t *p4est_new_ext (sc_MPI_Comm mpicomm, void *user_pointer); /** Create a new mesh. + * This function sets a subset of the mesh creation parameters. For full control + * use \ref p4est_mesh_new_params. * \param [in] p4est A forest that is fully 2:1 balanced. * \param [in] ghost The ghost layer created from the * provided p4est. @@ -389,8 +391,8 @@ p4est_t *p4est_new_ext (sc_MPI_Comm mpicomm, * compute the quad_to_tree list. * \param [in] compute_level_lists Boolean to decide whether to compute the * level lists in quad_level. - * \param [in] btype Currently ignored, only face neighbors - * are stored. + * \param [in] btype Flag indicating the connection types (face, + corner) stored in the mesh. * \return A fully allocated mesh structure. */ p4est_mesh_t *p4est_mesh_new_ext (p4est_t * p4est, diff --git a/src/p4est_mesh.c b/src/p4est_mesh.c index e6fd6d05e..35b6aca9b 100644 --- a/src/p4est_mesh.c +++ b/src/p4est_mesh.c @@ -145,6 +145,15 @@ mesh_edge_process_inter_tree_edges (p8est_iter_edge_info_t * info, int goodones = 0; int edgeid_offset = mesh->local_num_quadrants + mesh->ghost_num_quadrants; + /* variables needed for adding edge-hanging corner information */ + int add_hedges = 0; + int cornerid; + int cid; + int cgoodones = 0; + int8_t *ccorners; + p4est_locidx_t *cquads; + int8_t *pccorner; + p4est_locidx_t *pcquad; /* sanity check for subedge_id */ P4EST_ASSERT (-1 <= subedge_id && subedge_id < 2); @@ -162,6 +171,16 @@ mesh_edge_process_inter_tree_edges (p8est_iter_edge_info_t * info, equads = P4EST_ALLOC (p4est_locidx_t, nAdjacentQuads); eedges = P4EST_ALLOC (int8_t, nAdjacentQuads); + ccorners = NULL; + cquads = NULL; + if (mesh->params.edgehanging_corners && + mesh->params.btype >= P8EST_CONNECT_CORNER) { + add_hedges = 1; + /* we have at most one hanging corner neighbor per edge-neighboring tree */ + cquads = P4EST_ALLOC (p4est_locidx_t, cz - 1); + ccorners = P4EST_ALLOC (int8_t, cz - 1); + } + P4EST_ASSERT (0 <= side1->treeid && side1->treeid < info->p4est->connectivity->num_trees); tree1 = p4est_tree_array_index (info->p4est->trees, side1->treeid); @@ -217,6 +236,18 @@ mesh_edge_process_inter_tree_edges (p8est_iter_edge_info_t * info, equads[goodones] = qid2; eedges[goodones] = P8EST_EDGES * localOri + (int) side2->edge; ++goodones; + if (add_hedges) { + /* store edge-hanging corner neighbor */ + int8_t subEdgeIdx = ((subedge_id ^ 1) + localOri) % 2; + qid2 = + side2->is.hanging.quadid[subEdgeIdx] + + (side2->is.hanging.is_ghost[subEdgeIdx] ? + mesh->local_num_quadrants : tree2->quadrants_offset); + cquads[cgoodones] = qid2; + ccorners[cgoodones] = + p8est_edge_corners[side2->edge][subEdgeIdx ^ 1]; + ++cgoodones; + } } else if (!side1->is_hanging && side2->is_hanging) { /* check if we have to swap hanging quads in order to store @@ -265,6 +296,23 @@ mesh_edge_process_inter_tree_edges (p8est_iter_edge_info_t * info, /* we have excluded between 0 and all quadrants. */ P4EST_ASSERT (0 <= goodones && goodones < nAdjacentQuads); + if (add_hedges) { + if (cgoodones > 0) { + /* Allocate and fill corner information in the mesh structure */ + cornerid = mesh_corner_allocate (mesh, cgoodones, &pcquad, &pccorner); + /* "link" to arrays encoding inter-tree corner-neighborhood */ + cid = p8est_edge_corners[side1->edge][subedge_id ^ 1]; + P4EST_ASSERT (mesh->quad_to_corner[P8EST_CHILDREN * qid1 + cid] == -1); + mesh->quad_to_corner[P8EST_CHILDREN * qid1 + cid] = + edgeid_offset + cornerid; + /* populate allocated memory */ + memcpy (pcquad, cquads, cgoodones * sizeof (p4est_locidx_t)); + memcpy (pccorner, ccorners, cgoodones * sizeof (int8_t)); + } + P4EST_FREE (cquads); + P4EST_FREE (ccorners); + } + if (goodones > 0) { /* Allocate and fill edge information in the mesh structure */ edgeid = mesh_edge_allocate (mesh, goodones, &pequad, &peedge); @@ -554,56 +602,22 @@ mesh_iter_edge (p8est_iter_edge_info_t * info, void *user_data) return; } if (cz == 2) { + P4EST_ASSERT (info->tree_boundary); for (iz = 0; iz < cz; ++iz) { side1 = (p8est_iter_edge_side_t *) sc_array_index (&info->sides, iz); P4EST_ASSERT (0 <= side1->treeid && side1->treeid < info->p4est->connectivity->num_trees); P4EST_ASSERT (0 <= side1->edge && side1->edge < P8EST_EDGES); - - if (info->tree_boundary) { - if (side1->is_hanging) { - for (j = 0; j < 2; ++j) { - if (!side1->is.hanging.is_ghost[j]) { - mesh_edge_process_inter_tree_edges (info, side1, j, mesh, cz, - iz); - } - } - } - else { - if (!side1->is.full.is_ghost) { - mesh_edge_process_inter_tree_edges (info, side1, -1, mesh, cz, - iz); + if (side1->is_hanging) { + for (j = 0; j < 2; ++j) { + if (!side1->is.hanging.is_ghost[j]) { + mesh_edge_process_inter_tree_edges (info, side1, j, mesh, cz, iz); } } } else { - if (!side1->is_hanging) { - if (!side1->is.full.is_ghost) { - tree1 = - p4est_tree_array_index (info->p4est->trees, side1->treeid); - qid1 = side1->is.full.quadid + tree1->quadrants_offset; - - P4EST_ASSERT (0 <= qid1 && qid1 < mesh->local_num_quadrants); - P4EST_ASSERT - (mesh->quad_to_edge[P8EST_EDGES * qid1 + side1->edge] == -1); - - mesh->quad_to_edge[P8EST_EDGES * qid1 + side1->edge] = -3; - } - } - else { - for (j = 0; j < 2; ++j) { - if (!side1->is.hanging.is_ghost[j]) { - tree1 = - p4est_tree_array_index (info->p4est->trees, side1->treeid); - qid1 = side1->is.hanging.quadid[j] + tree1->quadrants_offset; - - P4EST_ASSERT (0 <= qid1 && qid1 < mesh->local_num_quadrants); - - P4EST_ASSERT - (mesh->quad_to_edge[P8EST_EDGES * qid1 + side1->edge] == -1); - mesh->quad_to_edge[P8EST_EDGES * qid1 + side1->edge] = -3; - } - } + if (!side1->is.full.is_ghost) { + mesh_edge_process_inter_tree_edges (info, side1, -1, mesh, cz, iz); } } } @@ -647,7 +661,7 @@ mesh_iter_edge (p8est_iter_edge_info_t * info, void *user_data) int8_t *peedge; p4est_locidx_t *pequad; - P4EST_ASSERT (!info->tree_boundary); + P4EST_ASSERT (!info->tree_boundary && cz == P4EST_HALF); side1 = (p8est_iter_edge_side_t *) sc_array_index (&info->sides, 0); P4EST_ASSERT (0 <= side1->treeid && @@ -657,8 +671,8 @@ mesh_iter_edge (p8est_iter_edge_info_t * info, void *user_data) memset (visited, 0, P4EST_HALF * sizeof (int8_t)); - /* TODO: someone has time double-checking this loop bound? */ - for (iz = 0; iz < ((cz + 1) >> 1); ++iz) { + /* search cz/2 == P4EST_HALF/2 pairs of opposing edges */ + for (iz = 0; iz < (cz >> 1); ++iz) { side1 = side2 = NULL; qid1 = -1; eid1 = -1; @@ -887,6 +901,68 @@ mesh_iter_edge (p8est_iter_edge_info_t * info, void *user_data) mesh->quad_to_edge[in_qtoe] = qid1; } } + if (mesh->params.edgehanging_corners && + mesh->params.btype >= P8EST_CONNECT_CORNER) { + /* add corner neighbor information across edge-hanging corner */ + int notk; + int cid; + p4est_locidx_t in_qtoc; + + for (k = 0; k < 2; ++k) { + notk = k ^ 1; + + /* do not record anything if both sides are ghost */ + if (side1->is.hanging.is_ghost[k] && + side2->is.hanging.is_ghost[notk]) { + continue; + } + + /* get qid1 and qid2 */ + if (!side1->is.hanging.is_ghost[k]) { + qid1 = side1->is.hanging.quadid[k] + qoffset; + P4EST_ASSERT (0 <= qid1 + && qid1 < mesh->local_num_quadrants); + } + else { + P4EST_ASSERT (side1->is.hanging.quad[k] != NULL); + P4EST_ASSERT (0 <= side1->is.hanging.quadid[k] && + side1->is.hanging.quadid[k] < + mesh->ghost_num_quadrants); + qid1 = + mesh->local_num_quadrants + side1->is.hanging.quadid[k]; + } + if (!side2->is.hanging.is_ghost[notk]) { + qid2 = side2->is.hanging.quadid[notk] + qoffset; + P4EST_ASSERT (0 <= qid2 + && qid2 < mesh->local_num_quadrants); + } + else { + P4EST_ASSERT (side2->is.hanging.quad[notk] != NULL); + P4EST_ASSERT (0 <= side2->is.hanging.quadid[notk] && + side2->is.hanging.quadid[notk] < + mesh->ghost_num_quadrants); + qid2 = + mesh->local_num_quadrants + + side2->is.hanging.quadid[notk]; + } + + /* write values */ + /* The corner lies on the edge, so we can use + * p8est_edge_corners as lookup. k indexes a proper corner, + * notk indexes the hanging corner we are interested in */ + cid = p8est_edge_corners[side1->edge][notk]; + if (!side1->is.hanging.is_ghost[k]) { + in_qtoc = P4EST_CHILDREN * qid1 + cid; + P4EST_ASSERT (mesh->quad_to_corner[in_qtoc] == -1); + mesh->quad_to_corner[in_qtoc] = qid2; + } + if (!side2->is.hanging.is_ghost[notk]) { + in_qtoc = P4EST_CHILDREN * qid2 + (cid ^ 7); + P4EST_ASSERT (mesh->quad_to_corner[in_qtoc] == -1); + mesh->quad_to_corner[in_qtoc] = qid1; + } + } + } } visited[j] = 1; side1 = side2 = NULL; @@ -1136,6 +1212,19 @@ p4est_mesh_memory_used (p4est_mesh_t * mesh) return all_memory; } +void +p4est_mesh_params_init (p4est_mesh_params_t * params) +{ + memset (params, 0, sizeof (p4est_mesh_params_t)); + + params->compute_level_lists = 0; + params->compute_tree_index = 0; + params->btype = P4EST_CONNECT_FACE; +#ifdef P4_TO_P8 + params->edgehanging_corners = 0; +#endif +} + p4est_mesh_t * p4est_mesh_new (p4est_t * p4est, p4est_ghost_t * ghost, p4est_connect_type_t btype) @@ -1147,6 +1236,24 @@ p4est_mesh_t * p4est_mesh_new_ext (p4est_t * p4est, p4est_ghost_t * ghost, int compute_tree_index, int compute_level_lists, p4est_connect_type_t btype) +{ + p4est_mesh_params_t params; + + /* initialize parameter struct to pass to mesh_new_params */ + p4est_mesh_params_init (¶ms); + params.btype = btype; + params.compute_level_lists = compute_level_lists; + params.compute_tree_index = compute_tree_index; +#ifdef P4_TO_P8 + params.edgehanging_corners = 0; +#endif + + return p4est_mesh_new_params (p4est, ghost, ¶ms); +} + +p4est_mesh_t * +p4est_mesh_new_params (p4est_t * p4est, p4est_ghost_t * ghost, + p4est_mesh_params_t * params) { int do_corner = 0; #ifdef P4_TO_P8 @@ -1159,27 +1266,36 @@ p4est_mesh_new_ext (p4est_t * p4est, p4est_ghost_t * ghost, p4est_mesh_t *mesh; /* check whether input condition for p4est is met */ - P4EST_ASSERT (p4est_is_balanced (p4est, btype)); + P4EST_ASSERT (p4est_is_balanced (p4est, params->btype)); mesh = P4EST_ALLOC_ZERO (p4est_mesh_t, 1); + /* store mesh creation parameters in mesh */ + if (params != NULL) { + mesh->params = *params; + } + else { + p4est_mesh_params_init (&mesh->params); + } + /* number of local quadrants and number of local ghost cells */ lq = mesh->local_num_quadrants = p4est->local_num_quadrants; ng = mesh->ghost_num_quadrants = (p4est_locidx_t) ghost->ghosts.elem_count; /* decide which callback function have to be activated */ #ifdef P4_TO_P8 - if (btype >= P8EST_CONNECT_EDGE) { + if (mesh->params.btype >= P8EST_CONNECT_EDGE) { do_edge = 1; } #endif - if (btype >= P4EST_CONNECT_FULL) { + if (mesh->params.btype >= P4EST_CONNECT_FULL) { do_corner = 1; } - do_volume = compute_tree_index || compute_level_lists; + do_volume = mesh->params.compute_tree_index + || mesh->params.compute_level_lists; /* Optional map of tree index for each quadrant */ - if (compute_tree_index) { + if (mesh->params.compute_tree_index) { mesh->quad_to_tree = P4EST_ALLOC (p4est_topidx_t, lq); } @@ -1190,7 +1306,7 @@ p4est_mesh_new_ext (p4est_t * p4est, p4est_ghost_t * ghost, mesh->quad_to_half = sc_array_new (P4EST_HALF * sizeof (p4est_locidx_t)); /* Allocate optional per-level lists of quadrants */ - if (compute_level_lists) { + if (mesh->params.compute_level_lists) { mesh->quad_level = P4EST_ALLOC (sc_array_t, P4EST_QMAXLEVEL + 1); for (jl = 0; jl <= P4EST_QMAXLEVEL; ++jl) { @@ -1923,6 +2039,7 @@ get_corner_neighbors (p4est_t * p4est, p4est_ghost_t * ghost, return 0; } +/*** OUTDATED FUNCTION ***/ p4est_locidx_t p4est_mesh_get_neighbors (p4est_t * p4est, p4est_ghost_t * ghost, p4est_mesh_t * mesh, p4est_locidx_t curr_quad_id, diff --git a/src/p4est_mesh.h b/src/p4est_mesh.h index ee4030b60..2df4213d0 100644 --- a/src/p4est_mesh.h +++ b/src/p4est_mesh.h @@ -24,7 +24,19 @@ /** \file p4est_mesh.h * - * forest topology in a conventional mesh format + * Forest topology in a conventional mesh format. + * + * A typical workflow starts with \ref p4est_mesh_params_init to initialize a + * \ref p4est_mesh_params_t, followed by eventual user-dependent changes to the + * parameters. + * + * Next a \ref p4est_mesh_t is created with \ref p4est_mesh_new_params. + * + * Now, the user can create a \ref p4est_mesh_face_neighbor_t with + * \ref p4est_mesh_face_neighbor_init and loop over a quadrants face neighbors + * by repeated calls to \ref p4est_mesh_face_neighbor_next. + * + * Once done, the mesh has to be destroyed with \ref p4est_mesh_destroy. * * \ingroup p4est */ @@ -36,17 +48,34 @@ SC_EXTERN_C_BEGIN; +/** This structure contains the different parameters of mesh creation. + * A default instance can be initialized by calling \ref p4est_mesh_params_init + * and used for mesh creation by calling \ref p4est_mesh_new_params. */ +typedef struct +{ + int compute_tree_index; /**< Boolean to decide whether to + allocate and compute the + quad_to_tree list. */ + int compute_level_lists; /**< Boolean to decide whether to + compute the level lists in + quad_level. */ + p4est_connect_type_t btype; /**< Flag indicating the + connection types (face, edge, + corner) stored in the mesh. */ +} +p4est_mesh_params_t; + /** This structure contains complete mesh information on a 2:1 balanced forest. * It stores the locally relevant neighborhood, that is, all locally owned * quadrants and one layer of adjacent ghost quadrants and their owners. * * For each local quadrant, its tree number is stored in quad_to_tree. * The quad_to_tree array is NULL by default and can be enabled using - * \ref p4est_mesh_new_ext. + * \ref p4est_mesh_new_ext or \ref p4est_mesh_new_params. * For each ghost quadrant, its owner rank is stored in ghost_to_proc. * For each level, an array of local quadrant numbers is stored in quad_level. * The quad_level array is NULL by default and can be enabled using - * \ref p4est_mesh_new_ext. + * \ref p4est_mesh_new_ext or \ref p4est_mesh_new_params. * * The quad_to_quad list stores one value for each local quadrant's face. * This value is in 0..local_num_quadrants-1 for local quadrants, or in @@ -93,18 +122,18 @@ SC_EXTERN_C_BEGIN; * only happens on the domain boundary, which is necessarily a tree boundary. * Corner-neighbors for hanging nodes are assigned the value -1. * - * TODO: In case of an inter-tree corner neighbor relation in a brick-like - * situation (exactly one neighbor, diagonally opposite corner number), - * use the same encoding as for corners within a tree. + * The params struct describes the parameters the mesh was created with. + * For full control over the parameters, use \ref p8est_mesh_new_params for + * mesh creation. */ typedef struct { - p4est_locidx_t local_num_quadrants; - p4est_locidx_t ghost_num_quadrants; + p4est_locidx_t local_num_quadrants; /**< number of process-local quadrants */ + p4est_locidx_t ghost_num_quadrants; /**< number of ghost-layer quadrants */ p4est_topidx_t *quad_to_tree; /**< tree index for each local quad. - Is NULL by default, but may be - enabled by \ref p4est_mesh_new_ext. */ + Is NULL if compute_tree_index in + params is 0. */ int *ghost_to_proc; /**< processor for each ghost quad */ p4est_locidx_t *quad_to_quad; /**< one index for each of the 4 faces */ @@ -115,15 +144,21 @@ typedef struct The array has entries indexed by 0..P4EST_QMAXLEVEL inclusive that are arrays of local quadrant ids. - Is NULL by default, but may be - enabled by \ref p4est_mesh_new_ext. */ - - /* These members are NULL if corners are not requested in \ref p4est_mesh_new. */ - p4est_locidx_t local_num_corners; /* tree-boundary corners */ - p4est_locidx_t *quad_to_corner; /* 4 indices for each local quad */ - sc_array_t *corner_offset; /* local_num_corners + 1 entries */ - sc_array_t *corner_quad; /* corner_offset indexes into this */ - sc_array_t *corner_corner; /* and this one too (type int8_t) */ + Is NULL if compute_level_lists in + params is 0. */ + + /* These members are NULL if btype in params is < P8EST_CONNECT_CORNER and can + * be requested in \ref p4est_mesh_new. */ + p4est_locidx_t local_num_corners; /**< tree-boundary corners */ + p4est_locidx_t *quad_to_corner; /**< 4 indices for each local quad */ + sc_array_t *corner_offset; /**< local_num_corners + 1 entries */ + sc_array_t *corner_quad; /**< corner_offset indexes into this */ + sc_array_t *corner_corner; /**< and this one too (type int8_t) */ + + p4est_mesh_params_t params; /**< parameters the mesh was created + with, e.g. by passing them to + \ref p4est_mesh_new_ext or + \ref p4est_mesh_new_params */ } p4est_mesh_t; @@ -133,21 +168,22 @@ p4est_mesh_t; typedef struct { /* forest information */ - p4est_t *p4est; - p4est_ghost_t *ghost; - p4est_mesh_t *mesh; + p4est_t *p4est; /**< the forest */ + p4est_ghost_t *ghost; /**< the ghost layer of the forest */ + p4est_mesh_t *mesh; /**< a mesh derived from the forest */ /* quadrant information */ - p4est_topidx_t which_tree; - p4est_locidx_t quadrant_id; /* tree-local quadrant index */ - p4est_locidx_t quadrant_code; /* 4 * (quadrant_id + tree_offset) */ + p4est_topidx_t which_tree; /**< the current tree index */ + p4est_locidx_t quadrant_id; /**< tree-local quadrant index */ + p4est_locidx_t quadrant_code; /**< 4 * (quadrant_id + tree_offset) */ /* neighbor information */ - int face; /* Face number in 0..3. */ - int subface; /* Hanging neighbor number in 0..1. */ + int face; /**< Face number in 0..3. */ + int subface; /**< Hanging neighbor number in 0..1. */ /* internal information */ - p4est_locidx_t current_qtq; + p4est_locidx_t current_qtq; /**< track index of current neighboring + quadrant */ } p4est_mesh_face_neighbor_t; @@ -157,9 +193,14 @@ p4est_mesh_face_neighbor_t; */ size_t p4est_mesh_memory_used (p4est_mesh_t * mesh); +/** Initialize a default p4est_mesh_params_t structure. + * The parameters are set to create the most basic mesh structure, without + * tree index and level lists and considering only face connections. */ +void p4est_mesh_params_init (p4est_mesh_params_t * params); + /** Create a p4est_mesh structure. * This function does not populate the quad_to_tree and quad_level fields. - * To populate them, use \ref p4est_mesh_new_ext. + * To populate them, use \ref p4est_mesh_new_params. * \param [in] p4est A forest that is fully 2:1 balanced. * \param [in] ghost The ghost layer created from the provided p4est. * \param [in] btype Determines the highest codimension of neighbors. @@ -169,6 +210,17 @@ p4est_mesh_t *p4est_mesh_new (p4est_t * p4est, p4est_ghost_t * ghost, p4est_connect_type_t btype); +/** Create a new mesh. + * \param [in] p4est A forest that is fully 2:1 balanced. + * \param [in] ghost The ghost layer created from the provided p4est. + * \param [in] params The mesh creation parameters. If NULL, the function + * defaults to the parameters of \ref p4est_mesh_params_init. + * \return A fully allocated mesh structure. + */ +p4est_mesh_t *p4est_mesh_new_params (p4est_t * p4est, + p4est_ghost_t * ghost, + p4est_mesh_params_t * params); + /** Destroy a p4est_mesh structure. * \param [in] mesh Mesh structure previously created by p4est_mesh_new. */ @@ -187,6 +239,7 @@ p4est_quadrant_t *p4est_mesh_get_quadrant (p4est_t * p4est, p4est_mesh_t * mesh, p4est_locidx_t qid); +/*** OUTDATED FUNCTION ***/ /** Lookup neighboring quads of quadrant in a specific direction. * \param [in] p4est Forest to be worked with. * \param [in] ghost Ghost layer. @@ -263,6 +316,9 @@ p4est_quadrant_t *p4est_mesh_quadrant_cumulative (p4est_t * p4est, /** Initialize a mesh neighbor iterator by quadrant index. * \param [out] mfn A p4est_mesh_face_neighbor_t to be initialized. + * \param [in] p4est Forest to be worked with. + * \param [in] ghost Ghost layer of the forest. + * \param [in] mesh A mesh derived from the forest. * \param [in] which_tree Tree of quadrant whose neighbors are looped over. * \param [in] quadrant_id Index relative to which_tree of quadrant. */ @@ -276,6 +332,9 @@ void p4est_mesh_face_neighbor_init2 (p4est_mesh_face_neighbor_t /** Initialize a mesh neighbor iterator by quadrant pointer. * \param [out] mfn A p4est_mesh_face_neighbor_t to be initialized. + * \param [in] p4est Forest to be worked with. + * \param [in] ghost Ghost layer of the forest. + * \param [in] mesh A mesh derived from the forest. * \param [in] which_tree Tree of quadrant whose neighbors are looped over. * \param [in] quadrant Pointer to quadrant contained in which_tree. */ @@ -292,7 +351,8 @@ void p4est_mesh_face_neighbor_init (p4est_mesh_face_neighbor_t * \param [out] ntree If not NULL, the tree number of the neighbor. * \param [out] nquad If not NULL, the quadrant number within tree. * For ghosts instead the number in ghost layer. - * \param [out] nface If not NULL, neighbor's face as in p4est_mesh_t. + * \param [out] nface If not NULL, neighbor's face encoding as in + quad_to_face array of p4est_mesh_t. * \param [out] nrank If not NULL, the owner process of the neighbor. * \return Either a real quadrant or one from the ghost layer. * Returns NULL when the iterator is done. diff --git a/src/p4est_to_p8est.h b/src/p4est_to_p8est.h index 07df7980a..719fcea82 100644 --- a/src/p4est_to_p8est.h +++ b/src/p4est_to_p8est.h @@ -162,6 +162,7 @@ #define p4est_iter_corner_t p8est_iter_corner_t #define p4est_iter_corner_side_t p8est_iter_corner_side_t #define p4est_iter_corner_info_t p8est_iter_corner_info_t +#define p4est_mesh_params_t p8est_mesh_params_t #define p4est_search_query_t p8est_search_query_t #define p4est_search_local_t p8est_search_local_t #define p4est_search_reorder_t p8est_search_reorder_t @@ -283,6 +284,8 @@ #define p4est_quadrant_set_morton_ext128 p8est_quadrant_set_morton_ext128 #define p4est_new_ext p8est_new_ext #define p4est_mesh_new_ext p8est_mesh_new_ext +#define p4est_mesh_new_params p8est_mesh_new_params +#define p4est_mesh_params_init p8est_mesh_params_init #define p4est_copy_ext p8est_copy_ext #define p4est_refine_ext p8est_refine_ext #define p4est_coarsen_ext p8est_coarsen_ext diff --git a/src/p8est_extended.h b/src/p8est_extended.h index 94adbf505..680cab5a7 100644 --- a/src/p8est_extended.h +++ b/src/p8est_extended.h @@ -386,6 +386,8 @@ p8est_t *p8est_new_ext (sc_MPI_Comm mpicomm, void *user_pointer); /** Create a new mesh. + * This function sets a subset of the mesh creation parameters. For full control + * use \ref p8est_mesh_new_params. * \param [in] p8est A forest that is fully 2:1 balanced. * \param [in] ghost The ghost layer created from the * provided p4est. @@ -393,11 +395,11 @@ p8est_t *p8est_new_ext (sc_MPI_Comm mpicomm, * compute the quad_to_tree list. * \param [in] compute_level_lists Boolean to decide whether to compute the * level lists in quad_level. - * \param [in] btype Currently ignored, only face neighbors - * are stored. + * \param [in] btype Flag indicating the connection types (face, + edge, corner) stored in the mesh. * \return A fully allocated mesh structure. */ -p8est_mesh_t *p8est_mesh_new_ext (p8est_t * p4est, +p8est_mesh_t *p8est_mesh_new_ext (p8est_t * p8est, p8est_ghost_t * ghost, int compute_tree_index, int compute_level_lists, diff --git a/src/p8est_mesh.h b/src/p8est_mesh.h index a087be352..09af56b9f 100644 --- a/src/p8est_mesh.h +++ b/src/p8est_mesh.h @@ -24,7 +24,19 @@ /** \file p8est_mesh.h * - * forest topology in a conventional mesh format + * Forest topology in a conventional mesh format. + * + * A typical workflow starts with \ref p8est_mesh_params_init to initialize a + * \ref p8est_mesh_params_t, followed by eventual user-dependent changes to the + * parameters. + * + * Next a \ref p8est_mesh_t is created with \ref p8est_mesh_new_params. + * + * Now, the user can create a \ref p8est_mesh_face_neighbor_t with + * \ref p8est_mesh_face_neighbor_init and loop over a quadrants face neighbors + * by repeated calls to \ref p8est_mesh_face_neighbor_next. + * + * Once done, the mesh has to be destroyed with \ref p8est_mesh_destroy. * * \ingroup p8est */ @@ -36,17 +48,37 @@ SC_EXTERN_C_BEGIN; +/** This structure contains the different parameters of mesh creation. + * A default instance can be initialzed by calling \ref p8est_mesh_params_init + * and used for mesh creation by calling \ref p8est_mesh_new_params. */ +typedef struct +{ + int compute_tree_index; /**< Boolean to decide whether to + allocate and compute the + quad_to_tree list. */ + int compute_level_lists; /**< Boolean to decide whether to + compute the level lists in + quad_level. */ + p8est_connect_type_t btype; /**< Flag indicating the + connection types (face, edge, + corner) stored in the mesh. */ + int edgehanging_corners; /**< Boolean to decide whether to + add corner neighbors across + coarse edges. */ +} +p8est_mesh_params_t; + /** This structure contains complete mesh information on a 2:1 balanced forest. * It stores the locally relevant neighborhood, that is, all locally owned * quadrants and one layer of adjacent ghost quadrants and their owners. * * For each local quadrant, its tree number is stored in quad_to_tree. * The quad_to_tree array is NULL by default and can be enabled using - * \ref p8est_mesh_new_ext. + * \ref p8est_mesh_new_ext or \ref p8est_mesh_new_params. * For each ghost quadrant, its owner rank is stored in ghost_to_proc. * For each level, an array of local quadrant numbers is stored in quad_level. * The quad_level array is NULL by default and can be enabled using - * \ref p8est_mesh_new_ext. + * \ref p8est_mesh_new_ext or \ref p8est_mesh_new_params. * * The quad_to_quad list stores one value for each local quadrant's face. * This value is in 0..local_num_quadrants-1 for local quadrants, or in @@ -125,20 +157,23 @@ SC_EXTERN_C_BEGIN; * * Corners with no diagonal neighbor at all are assigned the value -3. This * only happens on the domain boundary, which is necessarily a tree boundary. - * Corner-neighbors for face- and edge-hanging nodes are assigned the value -1. + * If edgehanging_corners in params is 0, all corner-neighbors for face- and + * edge-hanging nodes are assigned the value -1. + * If it is 1, we check for corner neighbors across coarse edges and assign -1 + * for the remaining face- and edge-hanging nodes. * - * TODO: In case of an inter-tree neighbor relation in a brick-like - * situation (one same-size neighbor, diagonally opposite edge/corner), - * use the same encoding as for edges/corners within a tree. + * The params struct describes the parameters the mesh was created with. + * For full control over the parameters, use \ref p8est_mesh_new_params for + * mesh creation. */ typedef struct { - p4est_locidx_t local_num_quadrants; - p4est_locidx_t ghost_num_quadrants; + p4est_locidx_t local_num_quadrants; /**< number of process-local quadrants */ + p4est_locidx_t ghost_num_quadrants; /**< number of ghost-layer quadrants */ p4est_topidx_t *quad_to_tree; /**< tree index for each local quad. - Is NULL by default, but may be - enabled by \ref p8est_mesh_new_ext. */ + Is NULL if compute_tree_index in + params is 0. */ int *ghost_to_proc; /**< processor for each ghost quad */ p4est_locidx_t *quad_to_quad; /**< one index for each of the 6 faces */ @@ -148,22 +183,29 @@ typedef struct The array has entries indexed by 0..P4EST_QMAXLEVEL inclusive that are arrays of local quadrant ids. - Is NULL by default, but may be - enabled by \ref p8est_mesh_new_ext. */ + Is NULL if compute_level_lists in + params is 0. */ - /* These members are NULL if edges are not requested in \ref p8est_mesh_new. */ + /* These members are NULL if btype in params is < P8EST_CONNECT_EDGE and can + * be requested in \ref p8est_mesh_new. */ p4est_locidx_t local_num_edges; /**< unsame-size and tree-boundary edges */ p4est_locidx_t *quad_to_edge; /**< 12 indices for each local quad */ sc_array_t *edge_offset; /**< local_num_edges + 1 entries */ sc_array_t *edge_quad; /**< edge_offset indexes into this */ sc_array_t *edge_edge; /**< and this one too (type int8_t) */ - /* These members are NULL if corners are not requested in \ref p8est_mesh_new. */ - p4est_locidx_t local_num_corners; /* tree-boundary corners */ - p4est_locidx_t *quad_to_corner; /* 8 indices for each local quad */ - sc_array_t *corner_offset; /* local_num_corners + 1 entries */ - sc_array_t *corner_quad; /* corner_offset indexes into this */ - sc_array_t *corner_corner; /* and this one too (type int8_t) */ + /* These members are NULL if btype in params is < P8EST_CONNECT_CORNER and can + * be requested in \ref p8est_mesh_new. */ + p4est_locidx_t local_num_corners; /**< tree-boundary corners */ + p4est_locidx_t *quad_to_corner; /**< 8 indices for each local quad */ + sc_array_t *corner_offset; /**< local_num_corners + 1 entries */ + sc_array_t *corner_quad; /**< corner_offset indexes into this */ + sc_array_t *corner_corner; /**< and this one too (type int8_t) */ + + p8est_mesh_params_t params; /**< parameters the mesh was created + with, e.g. by passing them to + \ref p8est_mesh_new_ext or + \ref p8est_mesh_new_params */ } p8est_mesh_t; @@ -173,21 +215,22 @@ p8est_mesh_t; typedef struct { /* forest information */ - p8est_t *p4est; - p8est_ghost_t *ghost; - p8est_mesh_t *mesh; + p8est_t *p4est; /**< the forest */ + p8est_ghost_t *ghost; /**< the ghost layer of the forest */ + p8est_mesh_t *mesh; /**< a mesh derived from the forest */ /* quadrant information */ - p4est_topidx_t which_tree; - p4est_locidx_t quadrant_id; /* tree-local quadrant index */ - p4est_locidx_t quadrant_code; /* 6 * (quadrant_id + tree_offset) */ + p4est_topidx_t which_tree; /**< the current tree index */ + p4est_locidx_t quadrant_id; /**< tree-local quadrant index */ + p4est_locidx_t quadrant_code; /**< 6 * (quadrant_id + tree_offset) */ /* neighbor information */ - int face; /* Face number in 0..5. */ - int subface; /* Hanging neighbor number in 0..3. */ + int face; /**< Face number in 0..5. */ + int subface; /**< Hanging neighbor number in 0..3. */ /* internal information */ - p4est_locidx_t current_qtq; + p4est_locidx_t current_qtq; /**< track index of current neighboring + quadrant */ } p8est_mesh_face_neighbor_t; @@ -197,9 +240,15 @@ p8est_mesh_face_neighbor_t; */ size_t p8est_mesh_memory_used (p8est_mesh_t * mesh); +/** Initialize a default p8est_mesh_params_t structure. + * The parameters are set to create the most basic mesh structure, without + * tree index and level lists and considering only face connections. */ +void p8est_mesh_params_init (p8est_mesh_params_t * params); + /** Create a p8est_mesh structure. - * This function does not populate the quad_to_tree and quad_level fields. - * To populate them, use \ref p8est_mesh_new_ext. + * This function does not populate the quad_to_tree and quad_level fields and + * ignores corner neighbors across edge-hanging corners. + * To populate them, use \ref p8est_mesh_new_params. * \param [in] p8est A forest that is fully 2:1 balanced. * \param [in] ghost The ghost layer created from the provided p4est. * \param [in] btype Determines the highest codimension of neighbors. @@ -208,6 +257,17 @@ size_t p8est_mesh_memory_used (p8est_mesh_t * mesh); p8est_mesh_t *p8est_mesh_new (p8est_t * p8est, p8est_ghost_t * ghost, p8est_connect_type_t btype); +/** Create a new mesh. + * \param [in] p8est A forest that is fully 2:1 balanced. + * \param [in] ghost The ghost layer created from the provided p4est. + * \param [in] params The mesh creation parameters. If NULL, the function + * defaults to the parameters of \ref p8est_mesh_params_init. + * \return A fully allocated mesh structure. + */ +p8est_mesh_t *p8est_mesh_new_params (p8est_t * p8est, + p8est_ghost_t * ghost, + p8est_mesh_params_t * params); + /** Destroy a p8est_mesh structure. * \param [in] mesh Mesh structure previously created by p8est_mesh_new. */ @@ -226,6 +286,7 @@ p8est_quadrant_t *p8est_mesh_get_quadrant (p8est_t * p4est, p8est_mesh_t * mesh, p4est_locidx_t qid); +/*** OUTDATED FUNCTION ***/ /** Lookup neighboring quads of quadrant in a specific direction * \param [in] p4est Forest to be worked with. * \param [in] ghost Ghost quadrants. @@ -238,7 +299,7 @@ p8est_quadrant_t *p8est_mesh_get_quadrant (p8est_t * p4est, * 18 .. 25 neighbor(-s) across c_{i-18} * \param [out] neighboring_quads Array containing neighboring quad(-s) * Needs to be empty, contains - * p4est_quadrant_t*. May be NULL, then \ref + * p4est_quadrant_t*. May be NULL, then \b * neighboring_qids must not be NULL. * \param [out] neighboring_encs Array containing encodings for neighboring * quads as described below @@ -304,6 +365,9 @@ p8est_quadrant_t *p8est_mesh_quadrant_cumulative (p8est_t * p8est, /** Initialize a mesh neighbor iterator by quadrant index. * \param [out] mfn A p8est_mesh_face_neighbor_t to be initialized. + * \param [in] p8est Forest to be worked with. + * \param [in] ghost Ghost layer of the forest. + * \param [in] mesh A mesh derived from the forest. * \param [in] which_tree Tree of quadrant whose neighbors are looped over. * \param [in] quadrant_id Index relative to which_tree of quadrant. */ @@ -317,6 +381,9 @@ void p8est_mesh_face_neighbor_init2 (p8est_mesh_face_neighbor_t /** Initialize a mesh neighbor iterator by quadrant pointer. * \param [out] mfn A p8est_mesh_face_neighbor_t to be initialized. + * \param [in] p8est Forest to be worked with. + * \param [in] ghost Ghost layer of the forest. + * \param [in] mesh A mesh derived from the forest. * \param [in] which_tree Tree of quadrant whose neighbors are looped over. * \param [in] quadrant Pointer to quadrant contained in which_tree. */ @@ -333,7 +400,8 @@ void p8est_mesh_face_neighbor_init (p8est_mesh_face_neighbor_t * \param [out] ntree If not NULL, the tree number of the neighbor. * \param [out] nquad If not NULL, the quadrant number within tree. * For ghosts instead the number in ghost layer. - * \param [out] nface If not NULL, neighbor's face as in p8est_mesh_t. + * \param [out] nface If not NULL, neighbor's face encoding as in + quad_to_face array of p8est_mesh_t. * \param [out] nrank If not NULL, the owner process of the neighbor. * \return Either a real quadrant or one from the ghost layer. * Returns NULL when the iterator is done.