From 8dfa831b3f725b7fdf12de1907fcdd589bb6b07b Mon Sep 17 00:00:00 2001 From: Carsten Burstedde Date: Sun, 21 Jul 2024 14:29:46 +0200 Subject: [PATCH] p4est_geometry: compute lnodes coordinates --- src/p4est_geometry.c | 446 +++++++++++++++++++++++++++++++++++++++++++ src/p4est_geometry.h | 52 ++++- src/p4est_to_p8est.h | 1 + src/p8est_geometry.h | 52 ++++- 4 files changed, 549 insertions(+), 2 deletions(-) diff --git a/src/p4est_geometry.c b/src/p4est_geometry.c index 979cab60f..2a0c18e9a 100644 --- a/src/p4est_geometry.c +++ b/src/p4est_geometry.c @@ -29,8 +29,10 @@ */ #ifndef P4_TO_P8 +#include #include #else +#include #include #endif @@ -619,3 +621,447 @@ p4est_geometry_new_sphere2d (p4est_connectivity_t * conn, double R) } /* p4est_geometry_new_sphere2d */ #endif /* !P4_TO_P8 */ + +static void +p4est_geometry_Q2_node (const p4est_quadrant_t *quad, +#ifdef P4_TO_P8 + int k, +#endif + int j, int i, p4est_quadrant_t *cnode) +{ + const p4est_qcoord_t h2 = P4EST_QUADRANT_LEN (quad->level + 1); + + P4EST_ASSERT (p4est_quadrant_is_valid (quad)); + P4EST_ASSERT (0 <= i && i <= 2); + P4EST_ASSERT (0 <= j && j <= 2); +#ifdef P4_TO_P8 + P4EST_ASSERT (0 <= k && k <= 2); +#endif + + cnode->x = quad->x + i * h2; + cnode->y = quad->y + j * h2; +#ifdef P4_TO_P8 + cnode->z = quad->z + k * h2; +#endif + cnode->level = P4EST_MAXLEVEL; +} + +/** A geometry coordinate tuple with tree and node information. */ +typedef struct p4est_geometry_node_coordinate +{ + p4est_locidx_t coord_index; /**< Index of node coordinate. */ + p4est_locidx_t local_node; /**< Local node index of point. */ + + /** Tree boundary index in [0, \ref P4EST_INSUL) of a tree node. + * To ensure correct node coordinates even in the case of periodic meshes, + * we hash the tree number, the local node, and this tree boundary index. + */ + int8_t tree_bound; +} +p4est_geometry_node_coordinate_t; + +static unsigned int +p4est_geometry_node_hash (const void *v, const void *u) +{ + uint32_t utt, uln, c; + const p4est_geometry_node_coordinate_t *nt = + (const p4est_geometry_node_coordinate_t *) v; + + P4EST_ASSERT (nt != NULL); + P4EST_ASSERT (0 <= nt->tree_bound && nt->tree_bound < P4EST_INSUL); + + /* execute primitive hash function */ + utt = (uint32_t) 0x7FFFFFFF; + uln = (uint32_t) nt->local_node; + c = (uint32_t) nt->tree_bound; + sc_hash_mix (utt, uln, c); + sc_hash_final (utt, uln, c); + + return (unsigned int) c; +} + +static int +p4est_geometry_node_equal (const void *v1, const void *v2, const void *u) +{ + const p4est_geometry_node_coordinate_t *nt1 = + (const p4est_geometry_node_coordinate_t *) v1; + const p4est_geometry_node_coordinate_t *nt2 = + (const p4est_geometry_node_coordinate_t *) v2; + + P4EST_ASSERT (nt1 != NULL); + P4EST_ASSERT (nt2 != NULL); + P4EST_ASSERT (0 <= nt1->tree_bound && nt1->tree_bound < P4EST_INSUL); + P4EST_ASSERT (0 <= nt2->tree_bound && nt2->tree_bound < P4EST_INSUL); + + return nt1->local_node == nt2->local_node && + nt1->tree_bound == nt2->tree_bound; +} + +static double * +p4est_geometry_coordinates_volume (p4est_locidx_t *volcoord, + sc_array_t *coordinates, + p4est_locidx_t *result_index, + p4est_geometry_node_coordinate_t *tn) +{ + p4est_locidx_t vval; + + /* basic checks */ + P4EST_ASSERT (volcoord != NULL); + P4EST_ASSERT (coordinates != NULL && + coordinates->elem_size == 3 * sizeof (double)); + P4EST_ASSERT (result_index != NULL); + P4EST_ASSERT (tn != NULL && tn->tree_bound == p4est_volume_point); + + /* determine whether this local node has already been processed */ + if ((vval = volcoord[tn->local_node]) >= 0) { + + /* this coordinate node is already computed */ + P4EST_ASSERT (0 <= vval && (size_t) vval < coordinates->elem_count); + *result_index = vval; + + /* no more to be done for this element node */ + return NULL; + } + P4EST_ASSERT (vval == -1); + + /* install a new element node coordinate in output array */ + *result_index = volcoord[tn->local_node] = + (p4est_locidx_t) coordinates->elem_count; + + /* we will write to this element node coordinate */ + return (double *) sc_array_push (coordinates); +} + +static double * +p4est_geometry_coordinates_insert (sc_hash_t *hash, sc_mempool_t *pool, + sc_array_t *coordinates, + p4est_locidx_t *result_index, + p4est_geometry_node_coordinate_t *tn) +{ + void **found; + p4est_geometry_node_coordinate_t *inserted; + + /* basic checks */ + P4EST_ASSERT (hash != NULL); + P4EST_ASSERT (pool != NULL && pool->elem_size == sizeof (*tn)); + P4EST_ASSERT (coordinates != NULL && + coordinates->elem_size == 3 * sizeof (double)); + P4EST_ASSERT (result_index != NULL); + P4EST_ASSERT (tn != NULL && tn->tree_bound != p4est_volume_point); + P4EST_ASSERT (0 <= tn->tree_bound && tn->tree_bound < P4EST_INSUL); + + /* determine whether this local node has already been processed */ + if (!sc_hash_insert_unique (hash, tn, &found)) { + + /* this coordinate node is already computed */ + inserted = *(p4est_geometry_node_coordinate_t **) found; + P4EST_ASSERT (p4est_geometry_node_equal (inserted, tn, hash->user_data)); + P4EST_ASSERT (0 <= inserted->coord_index && (size_t) + inserted->coord_index < coordinates->elem_count); + *result_index = inserted->coord_index; + + /* no more to be done for this element node */ + return NULL; + } + + /* install a new element node coordinate in output array */ + inserted = *(p4est_geometry_node_coordinate_t **) found = + (p4est_geometry_node_coordinate_t *) sc_mempool_alloc (pool); + *result_index = inserted->coord_index = + (p4est_locidx_t) coordinates->elem_count; + + /* remember node coordinate key */ + inserted->local_node = tn->local_node; + inserted->tree_bound = tn->tree_bound; + + /* we will write to this element node coordinate */ + return (double *) sc_array_push (coordinates); +} + +void +p4est_geometry_coordinates_lnodes (p4est_t *p4est, + p4est_lnodes_t *lnodes, + const double *refloc, + p4est_geometry_t *geom, + sc_array_t *coordinates, + sc_array_t *element_coordinates) +{ + static const double irlen = 1. / P4EST_ROOT_LEN; + int vno, vd, deg; + int i, kjixo, j, kjxo, kji, kxo; + int cid, fcd; + int n; +#ifdef P4_TO_P8 + int k; + int l; + int e; +#endif + int dtb[P4EST_DIM], dth[P4EST_DIM], dts; + size_t numenodes; + size_t volquery, treequery; + size_t collected, duplicates; + double abc[3], *xyz, *ret; + p4est_topidx_t tt; + p4est_locidx_t el, ne; + p4est_locidx_t *enodes, *ecoords, *volcoord; + p4est_lnodes_code_t *fcodes, fc; + p4est_quadrant_t sparent, *parent = &sparent, *quad; + p4est_quadrant_t scnode, *cnode = &scnode, *thequad; + p4est_tree_t *tree; + sc_mempool_t *pool; + sc_hash_t *hash; + p4est_geometry_node_coordinate_t stn, *tn = &stn; + + /* basic checks */ + P4EST_ASSERT (p4est != NULL); + P4EST_ASSERT (p4est->connectivity != NULL); + P4EST_ASSERT (geom == NULL || geom->X != NULL); + P4EST_ASSERT (lnodes != NULL); + P4EST_ASSERT (lnodes->degree == 1 || lnodes->degree == 2); + P4EST_ASSERT ((lnodes->degree == 1 && lnodes->vnodes == P4EST_CHILDREN) || + (lnodes->degree == 2 && lnodes->vnodes == P4EST_INSUL)); + P4EST_ASSERT (p4est->local_num_quadrants == lnodes->num_local_elements); + P4EST_ASSERT (coordinates != NULL && + coordinates->elem_size == 3 * sizeof (double)); + P4EST_ASSERT (element_coordinates == NULL || + element_coordinates->elem_size == sizeof (p4est_locidx_t)); + + /* prepare lookup information */ + pool = sc_mempool_new (sizeof (p4est_geometry_node_coordinate_t)); + volcoord = P4EST_ALLOC (p4est_locidx_t, lnodes->num_local_nodes); + memset (volcoord, -1, sizeof (p4est_locidx_t) * lnodes->num_local_nodes); + + /* loop through p4est qaudrants in natural order */ + memset (tn, 0, sizeof (stn)); + vd = (deg = lnodes->degree) + 1; + vno = lnodes->vnodes; + P4EST_ASSERT (vno == P4EST_DIM_POW (vd)); + numenodes = (size_t) lnodes->num_local_elements * (size_t) vno; + fcodes = lnodes->face_code; + enodes = lnodes->element_nodes; + if (element_coordinates == NULL) { + sc_array_resize (coordinates, lnodes->num_local_nodes); + sc_array_memset (coordinates, -1); + ecoords = NULL; + } + else { + sc_array_resize (coordinates, 0); + sc_array_resize (element_coordinates, numenodes); + + /* don't use sc_array_index since the array may have no elements */ + ecoords = (p4est_locidx_t *) element_coordinates->array; + } + volquery = treequery = 0; + collected = duplicates = 0; + for (tt = p4est->first_local_tree; tt <= p4est->last_local_tree; ++tt) { + + /* hashing is done anew for each tree */ + hash = sc_hash_new (p4est_geometry_node_hash, + p4est_geometry_node_equal, NULL, NULL); + + /* access tree and its local elements */ + tree = p4est_tree_array_index (p4est->trees, tt); + ne = (p4est_locidx_t) tree->quadrants.elem_count; + for (el = 0; el < ne; ++el, enodes += vno) { + + /* access quadrant in forest */ + quad = p4est_quadrant_array_index (&tree->quadrants, (size_t) el); + fc = *(fcodes++); + if (fc) { + /* this element is hanging */ + cid = fc & (P4EST_CHILDREN - 1); + P4EST_ASSERT (cid == p4est_quadrant_child_id (quad)); + fc >>= P4EST_DIM; + P4EST_ASSERT (fc); + p4est_quadrant_parent (quad, parent); + } + else { + /* non-hanging element */ + cid = 0; + P4EST_QUADRANT_INIT (parent); + } + + /* iterate through the local nodes referenced by the element */ + kji = 0; +#ifndef P4_TO_P8 + kxo = 0; +#else + for (k = 0; k < vd; ++k) { + kxo = (((cid & 4) ? deg - k : k) != 0) << 2; +#if 0 + } +#endif +#endif + for (j = 0; j < vd; ++j) { + kjxo = kxo + ((((cid & 2) ? deg - j : j) != 0) << 1); + for (i = 0; i < vd; ++i, ++kji) { + kjixo = kjxo + (((cid & 1) ? deg - i : i) != 0); + P4EST_ASSERT (0 <= kjixo && kjixo < P4EST_CHILDREN); + P4EST_ASSERT (deg != 1 || (kjixo ^ cid) == kji); + tn->local_node = enodes[kji]; + + /* if we identify coordinates with local nodes, shortcut */ + if (ecoords == NULL) { + ++volquery; + if (volcoord[tn->local_node] >= 0) { + /* this local node has been computed already */ + ++duplicates; + continue; + } + P4EST_ASSERT (volcoord[tn->local_node] == -1); + volcoord[tn->local_node] = 0; + } + + /* compute relevant quadrant to determine node coordinates */ + thequad = quad; + if (fc) { + fcd = p4est_lnodes_corner_hanging[kjixo]; + if (fcd >= 0 && (fc & (1 << fcd))) { + thequad = parent; + } + } + + /* compute coordinates of quadrant node */ + if (deg == 1) { + p4est_quadrant_corner_node (thequad, kji, cnode); + } + else { + P4EST_ASSERT (deg == 2); + p4est_geometry_Q2_node (thequad, +#ifdef P4_TO_P8 + k, +#endif + j, i, cnode); + } + + /* another shortcut for identified coordinates */ + if (ecoords == NULL) { + xyz = (double *) sc_array_index (coordinates, tn->local_node); + goto tnodes_shortcut; + } + + /* compute tree boundary status of quadrant node */ + dts = +#ifdef P4_TO_P8 + (dtb[2] = + (dth[2] = (cnode->z == P4EST_ROOT_LEN)) || cnode->z == 0) + +#endif + (dtb[1] = + (dth[1] = (cnode->y == P4EST_ROOT_LEN)) || cnode->y == 0) + + (dtb[0] = + (dth[0] = (cnode->x == P4EST_ROOT_LEN)) || cnode->x == 0); + switch (dts) { + case 0: + /* the most frequent case comes first */ + tn->tree_bound = p4est_volume_point; + break; + case 1: + /* the node sits inside a tree face */ + for (n = 0; n < P4EST_DIM; ++n) { + if (dtb[n]) { + tn->tree_bound = p4est_face_points[(n << 1) + dth[n]]; + break; + } + } + P4EST_ASSERT (n < P4EST_DIM); + break; +#ifdef P4_TO_P8 + case 2: + /* the node sits inside a tree edge */ + e = 0; + l = 0; + for (n = 0; n < P4EST_DIM; ++n) { + if (!dtb[n]) { + e += n << 2; + } + else { + e += dth[n] ? (1 << l) : 0; + ++l; + } + } + P4EST_ASSERT (l == 2); + tn->tree_bound = p8est_edge_points[e]; + break; +#endif + case P4EST_DIM: + /* the node sits on a tree corner equal a quadrant corner */ + tn->tree_bound = (deg == 1) ? p4est_corner_points[kji] : kji; + break; + default: + SC_ABORT_NOT_REACHED (); + } + P4EST_ASSERT (0 <= tn->tree_bound && tn->tree_bound < P4EST_INSUL); + + /* determine whether this local node has already been processed */ + P4EST_ASSERT ((dts == 0) == (tn->tree_bound == p4est_volume_point)); + if (dts == 0) { + /* volume points, never on a tree boundary, need no hashing */ + xyz = p4est_geometry_coordinates_volume (volcoord, coordinates, + ecoords + kji, tn); + ++volquery; + } + else { + /* codimension points may be on a tree boundary */ + xyz = p4est_geometry_coordinates_insert (hash, pool, coordinates, + ecoords + kji, tn); + ++treequery; + } + + /* evaluate reference coordinates if this point is a new one */ + if (xyz == NULL) { + ++duplicates; + } + else { +tnodes_shortcut: + ++collected; + + /* apply geometry transformation */ + ret = (geom == NULL) ? xyz : abc; + ret[0] = cnode->x * irlen; + ret[1] = cnode->y * irlen; +#ifndef P4_TO_P8 + ret[2] = 0.; +#else + ret[2] = cnode->z * irlen; +#endif + if (geom != NULL) { + /* transform tree reference into geometry image */ + geom->X (geom, tt, abc, xyz); + } + } + } /* i loop */ + } /* j loop */ +#ifdef P4_TO_P8 +#if 0 + { +#endif + } /* k loop */ +#endif + if (ecoords != NULL) { + P4EST_ASSERT (element_coordinates != NULL); + ecoords += vno; + } + } /* element loop */ + + /* free per-tree working memory */ + sc_hash_destroy (hash); + sc_mempool_truncate (pool); + } /* tree loop */ + P4EST_ASSERT (fcodes - lnodes->face_code == + (ptrdiff_t) lnodes->num_local_elements); + P4EST_ASSERT (enodes - lnodes->element_nodes == (ptrdiff_t) numenodes); + P4EST_ASSERT (volquery + treequery == numenodes); + P4EST_ASSERT (collected + duplicates == numenodes); + P4EST_ASSERT (collected == coordinates->elem_count); + P4EST_ASSERT (collected >= (size_t) lnodes->num_local_nodes); + + /* log some statistics */ + P4EST_INFOF (P4EST_STRING "_geometry_coordinates_lnodes " + "%lld nodes %lld coordinates\n", + (long long) lnodes->num_local_nodes, + (long long) collected); + + /* free temporary memory */ + P4EST_FREE (volcoord); + sc_mempool_destroy (pool); +} diff --git a/src/p4est_geometry.h b/src/p4est_geometry.h index a609335ab..20927eb14 100644 --- a/src/p4est_geometry.h +++ b/src/p4est_geometry.h @@ -41,7 +41,7 @@ #ifndef P4EST_GEOMETRY_H #define P4EST_GEOMETRY_H -#include +#include SC_EXTERN_C_BEGIN; @@ -173,6 +173,56 @@ p4est_geometry_t *p4est_geometry_new_disk2d (p4est_connectivity_t * conn, p4est_geometry_t *p4est_geometry_new_sphere2d (p4est_connectivity_t * conn, double R); +/** Compute node coordinates for a \ref p4est_lnodes structure. + * Presently we allow for an lnodes degree of 1 or 2. Cubic + * or higher degrees may be transparently enabled in the future. + * + * The simple mode assigns one tree reference coordinate to each lnode. + * This may not be suitable for visualizing periodic connectivities. + * + * In a more advanced mode indicated by NULL \c element_coordinates input, + * the coordinates are made unique by reference location: If a tree is + * periodic, for example, its corners reference the same lnode but will + * generate separate coordinate entries for proper visualization. + * There will be more coordinates generated than there are lnodes. + * + * The coordinate numbers generated by the present version of the function + * are partition-dependent. This may be seen as a flaw. Looking into it. + * + * \param [in] p4est A valid forest structure. + * \param [in] lnodes A valid \ref p4est_lnodes structure of degree + * 1 or 2. Higher degrees not presently allowed. + * Must be derived from the \c p4est. + * \param [in] refloc Eventually used for cubic and upwards degrees. + * We will expect degree + 1 many values for the + * one-dimensional reference node spacing. Out of + * these, the indices from 1 to (degree - 1) / 2 + * inclusive will be accessed by this function. + * The others default by symmetry considerations. + * \param [in] geom May be NULL for generating the tree-reference + * coordinates, or a valid geometry object for + * transforming the reference into mapped space. + * \param [in,out] coordinates On input, an array with entries of + * 3 double variables each. Resized in this + * function and populated with coordinate tuples. + * With a NULL geometry, these are in [0, 1]**2. + * Otherwise, they are mapped by the geometry. + * \param [in,out] element_coordinates This may be NULL, in which case + * we generate one coordinate tuple for each lnode. + * Otherwise, this must be an array with entries of + * type \ref p4est_locidx_t. Is resized to the same + * number of entries as \c lnodes->element_nodes. + * Its entries point into the coordinates array. + * The tree index of any given entry is implicit + * in that this array is derived from a p4est, + * where sets of (degree + 1)**2 entries each + * correspond to the forest elements in order. + */ +void p4est_geometry_coordinates_lnodes + (p4est_t *p4est, p4est_lnodes_t *lnodes, + const double *refloc, p4est_geometry_t *geom, + sc_array_t *coordinates, sc_array_t *element_coordinates); + SC_EXTERN_C_END; #endif /* !P4EST_GEOMETRY_H */ diff --git a/src/p4est_to_p8est.h b/src/p4est_to_p8est.h index c34efcd09..92391976f 100644 --- a/src/p4est_to_p8est.h +++ b/src/p4est_to_p8est.h @@ -521,6 +521,7 @@ #define p4est_geometry_destroy p8est_geometry_destroy #define p4est_geometry_new_connectivity p8est_geometry_new_connectivity #define p4est_geometry_connectivity_X p8est_geometry_connectivity_X +#define p4est_geometry_coordinates_lnodes p8est_geometry_coordinates_lnodes /* functions in p4est_vtk */ #define p4est_vtk_context_new p8est_vtk_context_new diff --git a/src/p8est_geometry.h b/src/p8est_geometry.h index 34158c350..158d74f72 100644 --- a/src/p8est_geometry.h +++ b/src/p8est_geometry.h @@ -38,7 +38,7 @@ #ifndef P8EST_GEOMETRY_H #define P8EST_GEOMETRY_H -#include +#include SC_EXTERN_C_BEGIN; @@ -151,6 +151,56 @@ p8est_geometry_t *p8est_geometry_new_torus (p8est_connectivity_t * conn, double R0, double R1, double R2); +/** Compute node coordinates for a \ref p8est_lnodes structure. + * Presently we allow for an lnodes degree of 1 or 2. Cubic + * or higher degrees may be transparently enabled in the future. + * + * The simple mode assigns one tree reference coordinate to each lnode. + * This may not be suitable for visualizing periodic connectivities. + * + * In a more advanced mode indicated by NULL \c element_coordinates input, + * the coordinates are made unique by reference location: If a tree is + * periodic, for example, its corners reference the same lnode but will + * generate separate coordinate entries for proper visualization. + * There will be more coordinates generated than there are lnodes. + * + * The coordinate numbers generated by the present version of the function + * are partition-dependent. This may be seen as a flaw. Looking into it. + * + * \param [in] p8est A valid forest structure. + * \param [in] lnodes A valid \ref p4est_lnodes structure of degree + * 1 or 2. Higher degrees not presently allowed. + * Must be derived from the \c p8est. + * \param [in] refloc Eventually used for cubic and upwards degrees. + * We will expect degree + 1 many values for the + * one-dimensional reference node spacing. Out of + * these, the indices from 1 to (degree - 1) / 2 + * inclusive will be accessed by this function. + * The others default by symmetry considerations. + * \param [in] geom May be NULL for generating the tree-reference + * coordinates, or a valid geometry object for + * transforming the reference into mapped space. + * \param [in,out] coordinates On input, an array with entries of + * 3 double variables each. Resized in this + * function and populated with coordinate tuples. + * With a NULL geometry, these are in [0, 1]**3. + * Otherwise, they are mapped by the geometry. + * \param [in,out] element_coordinates This may be NULL, in which case + * we generate one coordinate tuple for each lnode. + * Otherwise, this must be an array with entries of + * type \ref p4est_locidx_t. Is resized to the same + * number of entries as \c lnodes->element_nodes. + * Its entries point into the coordinates array. + * The tree index of any given entry is implicit + * in that this array is derived from a p4est, + * where sets of (degree + 1)**3 entries each + * correspond to the forest elements in order. + */ +void p8est_geometry_coordinates_lnodes + (p8est_t *p8est, p8est_lnodes_t *lnodes, + const double *refloc, p8est_geometry_t *geom, + sc_array_t *coordinates, sc_array_t *element_coordinates); + SC_EXTERN_C_END; #endif /* !P8EST_GEOMETRY_H */