diff --git a/src/t8_forest/t8_forest.cxx b/src/t8_forest/t8_forest.cxx index 3bef5f8e89..c23d2920cd 100644 --- a/src/t8_forest/t8_forest.cxx +++ b/src/t8_forest/t8_forest.cxx @@ -418,6 +418,14 @@ t8_forest_element_coordinate (t8_forest_t forest, t8_locidx_t ltree_id, const t8 const t8_gloidx_t gtreeid = t8_forest_global_tree_id (forest, ltree_id); /* Get the cmesh */ const t8_cmesh_t cmesh = t8_forest_get_cmesh (forest); +#if T8_ENABLE_DEBUG + if (tree_class == T8_ECLASS_TET) { + T8_ASSERT (vertex_coords[1] >= 0.0); + T8_ASSERT (vertex_coords[2] >= vertex_coords[1]); + T8_ASSERT (vertex_coords[0] >= vertex_coords[2]); + T8_ASSERT (vertex_coords[0] <= 1.0); + } +#endif /* Evaluate the geometry */ t8_geometry_evaluate (cmesh, gtreeid, vertex_coords, 1, coordinates); } diff --git a/src/t8_schemes/t8_scheme.hxx b/src/t8_schemes/t8_scheme.hxx index a115d58c67..8e21cfb7eb 100644 --- a/src/t8_schemes/t8_scheme.hxx +++ b/src/t8_schemes/t8_scheme.hxx @@ -103,7 +103,11 @@ struct t8_scheme t8_standalone_scheme, t8_standalone_scheme, t8_standalone_scheme, + t8_standalone_scheme, t8_standalone_scheme, + t8_standalone_scheme, + t8_standalone_scheme, + t8_standalone_scheme, /* Subelement schemes */ t8_subelem_scheme_hanging_nodes_quad, t8_subelem_scheme_hanging_nodes_tri diff --git a/src/t8_schemes/t8_standalone/t8_standalone.cxx b/src/t8_schemes/t8_standalone/t8_standalone.cxx index d1b81d76c0..31b6bd594b 100644 --- a/src/t8_schemes/t8_standalone/t8_standalone.cxx +++ b/src/t8_schemes/t8_standalone/t8_standalone.cxx @@ -38,11 +38,11 @@ t8_scheme_new_standalone (void) builder.add_eclass_scheme> (); builder.add_eclass_scheme> (); builder.add_eclass_scheme> (); - builder.add_eclass_scheme (); + builder.add_eclass_scheme> (); builder.add_eclass_scheme> (); - builder.add_eclass_scheme (); - builder.add_eclass_scheme (); - builder.add_eclass_scheme (); + builder.add_eclass_scheme> (); + builder.add_eclass_scheme> (); + builder.add_eclass_scheme> (); // builder.add_eclass_scheme> (); //ADD if PYRAMID in "element_get_shape" return builder.build_scheme (); } diff --git a/src/t8_schemes/t8_standalone/t8_standalone_elements.hxx b/src/t8_schemes/t8_standalone/t8_standalone_elements.hxx index d2f46e4e97..c0fc4371e0 100644 --- a/src/t8_schemes/t8_standalone/t8_standalone_elements.hxx +++ b/src/t8_schemes/t8_standalone/t8_standalone_elements.hxx @@ -60,10 +60,148 @@ constexpr uint8_t T8_ELEMENT_NUM_FACES[T8_ECLASS_COUNT] = { 0, 2, 4, 3, 6, 4, 5, /** Number of face children of the standalone element types */ constexpr uint8_t T8_ELEMENT_MAX_NUM_FACECHILDREN[T8_ECLASS_COUNT] = { 0, 1, 2, 2, 4, 4, 4, 4 }; +/** Number of faces per corner of the standalone element types */ +constexpr uint8_t T8_ELEMENT_NUM_CORNER_FACES[T8_ECLASS_COUNT] = { 0, 1, 2, 2, 3, 3, 3, 4 }; + +/** Number of corners per face of the standalone element types */ +constexpr uint8_t T8_ELEMENT_NUM_FACE_CORNERS[T8_ECLASS_COUNT] = { 0, 2, 2, 2, 4, 3, 4, 4 }; + /** Number of equations of the standalone element types */ constexpr uint8_t T8_ELEMENT_NUM_EQUATIONS[T8_ECLASS_COUNT] = { 0, 0, 0, 1, 0, 3, 1, 2 }; -/** Type for one entry of the element coordinates. */ +/* The following lookup tables (LUTs) are only specialized for element classes with type equations + * (triangle, tet, prism, pyramid), see t8_standalone_lut/. Unless stated otherwise, the first index + * is the element type interpreted as integer (type.to_ulong ()). */ + +/**PARENT CHILD BIJECTION*/ +/** Type of a child. Indexed by [parent type][child id (Iloc)]. */ +template +constexpr int8_t t8_element_type_Iloc_to_childtype[1 << T8_ELEMENT_NUM_EQUATIONS[TEclass]] + [T8_ELEMENT_NUM_CHILDREN[TEclass]]; +/** Cube id of a child, i.e. its position inside the parent's cube. Indexed by [parent type][child id (Iloc)]. */ +template +constexpr int8_t t8_element_type_Iloc_to_childcubeid[1 << T8_ELEMENT_NUM_EQUATIONS[TEclass]] + [T8_ELEMENT_NUM_CHILDREN[TEclass]]; +/** Type of the parent. Indexed by [child type][child cube id]. */ +template +constexpr int8_t t8_element_type_cubeid_to_parenttype[1 << T8_ELEMENT_NUM_EQUATIONS[TEclass]] + [1 << T8_ELEMENT_DIM[TEclass]]; +/** Child id (Iloc) of an element relative to its parent. Indexed by [child type][child cube id]. */ +template +constexpr int8_t t8_element_type_cubeid_to_Iloc[1 << T8_ELEMENT_NUM_EQUATIONS[TEclass]][1 << T8_ELEMENT_DIM[TEclass]]; + +/**TYPE EQUATIONS*/ +/** The two coordinate dimensions compared by each type equation. Indexed by [equation][0|1]. + * Bit i of the element type encodes the ordering of the coordinates in these two dimensions. */ +template +constexpr int8_t t8_type_edge_equations[T8_ELEMENT_NUM_EQUATIONS[TEclass]][2]; + +/**VERTEX*/ +/** Offset (0 or 1, in units of the element length) of a vertex relative to the anchor node. + * Indexed by [type][vertex][dim]. */ +template +constexpr int8_t t8_type_vertex_dim_to_binary[1 << T8_ELEMENT_NUM_EQUATIONS[TEclass]][T8_ELEMENT_NUM_CORNERS[TEclass]] + [T8_ELEMENT_DIM[TEclass]]; + +/**FACE*/ +/** 1 if the face lies inside the element's cube (neighbor differs only in its type), 0 if it lies on the cube boundary. + * Indexed by [type][face]. */ +template +constexpr int8_t t8_standalone_lut_face_internal[1 << T8_ELEMENT_NUM_EQUATIONS[TEclass]][T8_ELEMENT_NUM_FACES[TEclass]]; +/** Type bit that is flipped to obtain the neighbor across an internal face, -1 for cube boundary faces. + * Indexed by [type][face]. */ +template +constexpr int8_t t8_standalone_lut_type_face_to_typebit[1 << T8_ELEMENT_NUM_EQUATIONS[TEclass]] + [T8_ELEMENT_NUM_FACES[TEclass]]; +/** Orientation (+1/-1) of the outward face normal along \ref t8_standalone_lut_type_face_to_facenormal_dim, + * 0 for internal faces. Indexed by [type][face]. */ +template +constexpr int8_t t8_standalone_lut_type_face_to_sign[1 << T8_ELEMENT_NUM_EQUATIONS[TEclass]] + [T8_ELEMENT_NUM_FACES[TEclass]]; +/** Coordinate dimension of the face normal for cube boundary faces, -1 for internal faces. Indexed by [type][face]. */ +template +constexpr int8_t t8_standalone_lut_type_face_to_facenormal_dim[1 << T8_ELEMENT_NUM_EQUATIONS[TEclass]] + [T8_ELEMENT_NUM_FACES[TEclass]]; +/** Face number of the face neighbor's face that coincides with the given face. Indexed by [type][face]. */ +template +constexpr int8_t t8_standalone_lut_type_face_to_neighface[1 << T8_ELEMENT_NUM_EQUATIONS[TEclass]] + [T8_ELEMENT_NUM_FACES[TEclass]]; +/** Face of the parent on which the given face of the child lies, -1 if it lies inside the parent. + * Indexed by [child type][child cube id][child face]. */ +template +constexpr int8_t t8_standalone_lut_type_cubeid_face_to_parentface[1 << T8_ELEMENT_NUM_EQUATIONS[TEclass]] + [1 << T8_ELEMENT_DIM[TEclass]] + [T8_ELEMENT_NUM_FACES[TEclass]]; +/** 1 if the face lies on the upper boundary (coordinate = anchor + length) of the face normal dimension, + * 0 if it lies on the lower boundary. Indexed by [type][face]. */ +template +constexpr int8_t t8_standalone_lut_type_face_to_is_1_boundary[1 << T8_ELEMENT_NUM_EQUATIONS[TEclass]] + [T8_ELEMENT_NUM_FACES[TEclass]]; +/** Cube id of the last descendant touching the face, used to compute the last face descendant. + * Indexed by [type][face]. */ +template +constexpr int8_t t8_standalone_lut_type_face_to_last_facechilds_cubeid[1 << T8_ELEMENT_NUM_EQUATIONS[TEclass]] + [T8_ELEMENT_NUM_FACES[TEclass]]; +/** Child id of the i-th child touching the face. Indexed by [parent type][parent face][face child id]. */ +template +constexpr int8_t t8_standalone_lut_type_face_facechildid_to_childid[1 << T8_ELEMENT_NUM_EQUATIONS[TEclass]] + [T8_ELEMENT_NUM_FACES[TEclass]] + [T8_ELEMENT_MAX_NUM_FACECHILDREN[TEclass]]; +/** Face of the child that lies on the given parent face, -1 if none. Indexed by [parent type][child id][parent face]. */ +template +constexpr int8_t t8_standalone_lut_type_childid_face_to_childface[1 << T8_ELEMENT_NUM_EQUATIONS[TEclass]] + [T8_ELEMENT_NUM_CHILDREN[TEclass]] + [T8_ELEMENT_MAX_NUM_FACECHILDREN[TEclass]]; +/** Root tree face on which the element face can lie, -1 if an element of this type + * cannot touch the tree boundary with this face. Indexed by [type][face]. */ +template +constexpr int8_t t8_standalone_lut_type_face_to_tree_face[1 << T8_ELEMENT_NUM_EQUATIONS[TEclass]] + [T8_ELEMENT_NUM_FACES[TEclass]]; +/** Coordinate dimension of the boundary face element corresponding to a dimension of the volume element, + * -1 if the dimension is dropped. Indexed by [root face][volume dim]. */ +template +constexpr int8_t t8_standalone_lut_rootface_dim_to_facedim[T8_ELEMENT_NUM_FACES[TEclass]][T8_ELEMENT_DIM[TEclass]]; +/** Type equation of the boundary face element corresponding to a type equation of the volume element, + * -1 if the equation is dropped. Indexed by [root face][volume equation]. */ +template +constexpr int8_t t8_standalone_lut_rootface_eq_to_faceeq[T8_ELEMENT_NUM_FACES[TEclass]] + [T8_ELEMENT_NUM_EQUATIONS[TEclass]]; +/** Element face that lies on the given root tree face, -1 if an element of this type cannot touch it. + * Indexed by [type][root face]. */ +template +constexpr int8_t t8_standalone_lut_type_rootface_to_face[1 << T8_ELEMENT_NUM_EQUATIONS[TEclass]] + [T8_ELEMENT_NUM_FACES[TEclass]]; +/** Element class of each root tree face. Indexed by [root face]. */ +template +constexpr t8_eclass t8_standalone_lut_rootface_to_eclass[T8_ELEMENT_NUM_FACES[TEclass]]; + +/** Faces adjacent to a corner. Indexed by [type][corner][i]. */ +template +constexpr int8_t t8_standalone_lut_cornerface[1 << T8_ELEMENT_NUM_EQUATIONS[TEclass]][T8_ELEMENT_NUM_CORNERS[TEclass]] + [T8_ELEMENT_NUM_CORNER_FACES[TEclass]]; + +/** Corners of a face. Indexed by [type][face][i]. */ +template +constexpr int8_t t8_standalone_lut_facecorner[1 << T8_ELEMENT_NUM_EQUATIONS[TEclass]][T8_ELEMENT_NUM_FACES[TEclass]] + [T8_ELEMENT_NUM_FACE_CORNERS[TEclass]]; +/** Matrix mapping reference coordinates of the reference element to the element of the given type. + * Indexed by [type][out dim][in dim]. */ +template +constexpr int8_t t8_standalone_lut_transform_coords[1 << T8_ELEMENT_NUM_EQUATIONS[TEclass]][T8_ELEMENT_DIM[TEclass]] + [T8_ELEMENT_DIM[TEclass]]; + +/** Matrix mapping the internal coordinate frame back to the reference coordinates of the tree. + * Indexed by [type][out dim][in dim]; currently only type 0 is used. */ +template +constexpr int8_t t8_standalone_lut_backtransform_coords[1 << T8_ELEMENT_NUM_EQUATIONS[TEclass]][T8_ELEMENT_DIM[TEclass]] + [T8_ELEMENT_DIM[TEclass]]; + +#include "t8_standalone_lut/t8_standalone_lut_triangle.hxx" +#include "t8_standalone_lut/t8_standalone_lut_prism.hxx" +#include "t8_standalone_lut/t8_standalone_lut_pyra.hxx" +#include "t8_standalone_lut/t8_standalone_lut_tet.hxx" + +/** Type for the element coordinates. */ typedef int32_t t8_element_coord; /** Type for the element level. */ typedef uint8_t t8_element_level; @@ -92,6 +230,8 @@ struct t8_standalone_element t8_element_coords coords; /** The refinement level of the element relative to the root at level 0. */ t8_element_level level; + /** Bit array: which inequality is fulfilled at which level. */ + t8_element_type type; }; #endif /* T8_STANDALONE_ELEMENTS_HXX */ diff --git a/src/t8_schemes/t8_standalone/t8_standalone_implementation.hxx b/src/t8_schemes/t8_standalone/t8_standalone_implementation.hxx index 5255538575..ec6aa6ba45 100644 --- a/src/t8_schemes/t8_standalone/t8_standalone_implementation.hxx +++ b/src/t8_schemes/t8_standalone/t8_standalone_implementation.hxx @@ -148,7 +148,12 @@ struct t8_standalone_scheme: public t8_scheme_helpers *el = (const t8_standalone_element *) elem; + if constexpr (TEclass == T8_ECLASS_PYRAMID) { + if (el->type == T8_DPYRAMID_FIRST_PYRA_TYPE || el->type == T8_DPYRAMID_SECOND_PYRA_TYPE) { + return T8_ECLASS_PYRAMID; + } + else { + return T8_ECLASS_TET; + } + } return TEclass; } @@ -213,9 +232,15 @@ struct t8_standalone_scheme: public t8_scheme_helpers *el = (t8_standalone_element *) element; + return t8_standalone_lut_facecorner[el->type.to_ulong ()][face][corner]; + } } /** Return the face numbers of the faces sharing an element's corner. @@ -236,7 +261,15 @@ struct t8_standalone_scheme: public t8_scheme_helpers> face & 1) + 2 * face; + T8_ASSERTF (t8_standalone_lut_facecorner[0][1][1] == 2, "LUT is wrong!"); + + if constexpr (!T8_ELEMENT_NUM_EQUATIONS[TEclass]) { + return (corner >> face & 1) + 2 * face; + } + else { + t8_standalone_element *el = (t8_standalone_element *) element; + return t8_standalone_lut_cornerface[el->type.to_ulong ()][corner][face]; + } } /** Compute the shape of the face of an element. @@ -252,6 +285,7 @@ struct t8_standalone_scheme: public t8_scheme_helpers *el = (const t8_standalone_element *) elem; switch (TEclass) { case T8_ECLASS_VERTEX: SC_ABORT ("Vertices do not have faces.\n"); @@ -259,8 +293,29 @@ struct t8_standalone_scheme: public t8_scheme_helperstype[0] != el->type[1]) { + return T8_ECLASS_TRIANGLE; + } + else { + if (face == 4) // TODO: Check + return T8_ECLASS_QUAD; + else + return T8_ECLASS_TRIANGLE; + } default: SC_ABORT ("This function is not implemented yet.\n"); } @@ -305,9 +360,7 @@ struct t8_standalone_scheme: public t8_scheme_helperscoords[idim] != el2->coords[idim]) return 0; } - /* return el1->type == el2->type; - ToDo-Type */ - return 1; + return el1->type == el2->type; } // ################################################____ACCESSOR____################################################ @@ -336,9 +389,7 @@ struct t8_standalone_scheme: public t8_scheme_helperscoords[idim] = 0; } - /* el->type = 0; - ToDo-Type */ - return; + el->type = 0; } /** Compute the parent of a given element \b elem and store it in \b parent. @@ -364,7 +415,8 @@ struct t8_standalone_scheme: public t8_scheme_helperslevel > 0); if constexpr (T8_ELEMENT_NUM_EQUATIONS[TEclass]) { - SC_ABORT ("Only implemented for hypercubes.\n"); + const t8_cube_id cube_id = compute_cubeid (el, el->level); + parent_elem->type = t8_element_type_cubeid_to_parenttype[el->type.to_ulong ()][cube_id]; } const t8_element_coord length = element_get_len ((el->level)); @@ -394,7 +446,9 @@ struct t8_standalone_scheme: public t8_scheme_helperslevel && el->level <= T8_ELEMENT_MAXLEVEL[TEclass]); /* To get the number siblings, we first get the parent and then get the number of children of that parent*/ if constexpr (refines_irregular ()) { - SC_ABORT ("This function is not implemented yet.\n"); + t8_standalone_element parent; + element_get_parent ((t8_element_t *) el, (t8_element_t *) &parent); + return element_get_num_children ((t8_element_t *) &parent); } else { return T8_ELEMENT_NUM_CHILDREN[TEclass]; @@ -447,7 +501,8 @@ struct t8_standalone_scheme: public t8_scheme_helpers[el->type.to_ulong ()][childid]; + c->type = t8_element_type_Iloc_to_childtype[el->type.to_ulong ()][childid]; } else { cube_id = childid; @@ -470,7 +525,12 @@ struct t8_standalone_scheme: public t8_scheme_helperslevel); t8_child_id child_id; if constexpr (T8_ELEMENT_NUM_EQUATIONS[TEclass]) { - SC_ABORT ("Only implemented for hypercubes.\n"); + child_id = t8_element_type_cubeid_to_Iloc[el->type.to_ulong ()][cube_id]; } else { child_id = cube_id; @@ -579,29 +639,14 @@ struct t8_standalone_scheme: public t8_scheme_helpers parent, compare; - /* Take the parent of the first element as baseline to compare against */ - element_get_parent ((const t8_element_t *) fam[0], (t8_element_t *) &parent); - const int num_children = element_get_num_children ((const t8_element_t *) &parent); - for (int childid = 0; childid < num_children; childid++) { - /* check whether each element has the same parent */ - element_get_parent ((const t8_element_t *) fam[childid], (t8_element_t *) &compare); - if (element_compare ((const t8_element_t *) &parent, (const t8_element_t *) &compare)) { - return 0; - } - - /* check whether each element is the correct child of the collective parent */ - /* Could be replaced by type comparison as level is already checked in parent comparison */ - element_get_child ((const t8_element_t *) &parent, childid, (t8_element_t *) &compare); - if (element_compare ((const t8_element_t *) fam[childid], (const t8_element_t *) &compare)) { + const int num_siblings = element_get_num_siblings (fam[0]); + const int level = element_get_level (fam[0]); + for (int ielement = 1; ielement < num_siblings; ielement++) { + if (element_get_level (fam[ielement]) != level) { return 0; } } @@ -675,14 +720,18 @@ struct t8_standalone_scheme: public t8_scheme_helpers *el1 = (const t8_standalone_element *) elem1; const t8_standalone_element *el2 = (const t8_standalone_element *) elem2; /* get the first possible level of the nca*/ - int cube_ancestor_level = element_get_cube_nca_level (el1, el2); - int real_level; + const int cube_ancestor_level = element_get_cube_nca_level (el1, el2); + int real_level = cube_ancestor_level; if constexpr (T8_ELEMENT_NUM_EQUATIONS[TEclass]) { - SC_ABORT ("Only implemented for hypercubes.\n"); - } - else { - real_level = cube_ancestor_level; + t8_element_type ancestor_type, ancestor2_type; + do { + ancestor_type = element_compute_type_at_level (el1, real_level); + ancestor2_type = element_compute_type_at_level (el2, real_level); + real_level--; + } while (ancestor_type != ancestor2_type); + real_level++; /* we subtracted once too much*/ } + /* get the ancestor at the calculated level*/ element_get_ancestor (el1, real_level, (t8_standalone_element *) nca); T8_ASSERT (element_is_valid (nca)); @@ -739,7 +788,7 @@ struct t8_standalone_scheme: public t8_scheme_helperscoords[idim] |= coord_offset; } - T8_ASSERT (element_is_valid (desc)); + T8_ASSERT (element_is_valid ((t8_element_t *) d)); } // ################################################____FACE REFINEMENT____################################################ @@ -785,10 +834,18 @@ struct t8_standalone_scheme: public t8_scheme_helpers[el->type.to_ulong ()][face][ifacechild]; + } } for (int ifacechild = num_children - 1; ifacechild >= 0; ifacechild--) { element_get_child ((const t8_element_t *) el, child_indices[ifacechild], @@ -820,7 +877,15 @@ struct t8_standalone_scheme: public t8_scheme_helpers *el = (const t8_standalone_element *) elem; + const int child_id + = t8_standalone_lut_type_face_facechildid_to_childid[el->type.to_ulong ()][face][face_child]; + return t8_standalone_lut_type_childid_face_to_childface[el->type.to_ulong ()][child_id][face]; + } } /** Given a face of an element return the face number @@ -851,8 +916,8 @@ struct t8_standalone_scheme: public t8_scheme_helperslevel); + return t8_standalone_lut_type_cubeid_face_to_parentface[el->type.to_ulong ()][cube_id][face]; } } @@ -876,16 +941,29 @@ struct t8_standalone_scheme: public t8_scheme_helperslevel = level; if constexpr (T8_ELEMENT_NUM_EQUATIONS[TEclass]) { - SC_ABORT ("Only implemented for hypercubes.\n"); + first_descendant->type = el->type; } std::copy (el->coords.begin (), el->coords.end (), first_descendant->coords.begin ()); - const bool face_is_1_boundary = face % 2; + int face_is_1_boundary; - if (face_is_1_boundary) { //the face is a xi=1 boundary - const int facenormal_dim = face / 2; + if constexpr (!T8_ELEMENT_NUM_EQUATIONS[TEclass]) { + face_is_1_boundary = face % 2; + } + else { + face_is_1_boundary = t8_standalone_lut_type_face_to_is_1_boundary[el->type.to_ulong ()][face]; + } - const t8_element_coord coord_offset = element_get_len (el->level) - element_get_len (level); + if (face_is_1_boundary) { //the face is a xi=1 boundary + int facenormal_dim; + if constexpr (!T8_ELEMENT_NUM_EQUATIONS[TEclass]) { + facenormal_dim = face / 2; + } + else { + facenormal_dim = t8_standalone_lut_type_face_to_facenormal_dim[el->type.to_ulong ()][face]; + T8_ASSERT (facenormal_dim != -1); + } + t8_element_coord coord_offset = element_get_len (el->level) - element_get_len (level); first_descendant->coords[facenormal_dim] += coord_offset; } @@ -913,11 +991,18 @@ struct t8_standalone_scheme: public t8_scheme_helperslevel) - element_get_len (level); if constexpr (T8_ELEMENT_NUM_EQUATIONS[TEclass]) { - SC_ABORT ("Only implemented for hypercubes.\n"); + last_descendant->type = el->type; /**TODO: Check if this is always true! */ } for (int idim = 0; idim < T8_ELEMENT_DIM[TEclass]; idim++) { - const int multiplier = (idim == face / 2) ? face % 2 : 1; + int multiplier = 1; + if constexpr (!T8_ELEMENT_NUM_EQUATIONS[TEclass]) { + multiplier = (idim == face / 2) ? face % 2 : 1; + } + else { + t8_cube_id cube_id = t8_standalone_lut_type_face_to_last_facechilds_cubeid[el->type.to_ulong ()][face]; + multiplier = (cube_id & (1 << idim)) >> idim; // = cubeid[idim]; + } last_descendant->coords[idim] = el->coords[idim] + multiplier * coord_offset; } } @@ -947,7 +1032,12 @@ struct t8_standalone_scheme: public t8_scheme_helpers= x_j-a_j or x_j-a_j <= x_d-a_d if constexpr (T8_ELEMENT_NUM_EQUATIONS[TEclass]) { - SC_ABORT ("Only implemented for hypercubes.\n"); + for (int ieq = 0; ieq < T8_ELEMENT_NUM_EQUATIONS[TEclass]; ++ieq) { + if ((t8_type_edge_equations[ieq][0] == dim && el->type[ieq]) + || (t8_type_edge_equations[ieq][1] == dim && !el->type[ieq])) { + return 0; + } + } } } else { @@ -958,7 +1048,12 @@ struct t8_standalone_scheme: public t8_scheme_helpers= x_d-a_d if constexpr (T8_ELEMENT_NUM_EQUATIONS[TEclass]) { - SC_ABORT ("Only implemented for hypercubes.\n"); + for (int ieq = 0; ieq < T8_ELEMENT_NUM_EQUATIONS[TEclass]; ieq++) { + if ((t8_type_edge_equations[ieq][0] == dim && !el->type[ieq]) + || (t8_type_edge_equations[ieq][1] == dim && el->type[ieq])) { + return 0; + } + } } } } @@ -967,8 +1062,13 @@ struct t8_standalone_scheme: public t8_scheme_helpers[el->type.to_ulong ()][face]; + if (el->coords[t8_type_edge_equations[ieq][0]] + != el->coords[t8_type_edge_equations[ieq][1]]) { + return 0; + } } + else { SC_ABORT ("Cubes should not have internal faces!\n"); } @@ -993,8 +1093,13 @@ struct t8_standalone_scheme: public t8_scheme_helpers *el = (const t8_standalone_element *) elem; + return t8_standalone_lut_type_face_to_tree_face[el->type.to_ulong ()][face]; + } } /** Construct the face neighbor of a given element if this face neighbor @@ -1025,17 +1130,47 @@ struct t8_standalone_scheme: public t8_scheme_helpers *neighbor = (t8_standalone_element *) neigh; if constexpr (T8_ELEMENT_NUM_EQUATIONS[TEclass]) { - SC_ABORT ("Only implemented for hypercubes.\n"); + const int internal_face = t8_standalone_lut_face_internal[el->type.to_ulong ()][face]; + if (internal_face) { + /**Determine typebit*/ + int typebit = t8_standalone_lut_type_face_to_typebit[el->type.to_ulong ()][face]; + /**Change typebit*/ + neighbor->type.flip (typebit); + *neigh_face = t8_standalone_lut_type_face_to_neighface[el->type.to_ulong ()][face]; + return element_is_inside_root (neighbor); + } + } + + int facenormal_dim, sign; + + if constexpr (T8_ELEMENT_NUM_EQUATIONS[TEclass]) { + sign = t8_standalone_lut_type_face_to_sign[el->type.to_ulong ()][face]; + facenormal_dim = t8_standalone_lut_type_face_to_facenormal_dim[el->type.to_ulong ()][face]; + T8_ASSERT (facenormal_dim != -1); + } + else { + sign = face % 2 ? 1 : -1; + facenormal_dim = face / 2; } - const int facenormal_dim = face / 2; - const int sign = face % 2 ? 1 : -1; /**Adapt coordinates*/ const t8_element_coord length = element_get_len (el->level); - neighbor->coords[facenormal_dim] += length * sign; - *neigh_face = face ^ 1; + /**Adapt typebits*/ + if constexpr (T8_ELEMENT_NUM_EQUATIONS[TEclass]) { + for (int ieq = 0; ieq < T8_ELEMENT_NUM_EQUATIONS[TEclass]; ieq++) { + /**For all neighboring typebits, change typebit*/ + if (t8_type_edge_equations[ieq][0] == facenormal_dim + || t8_type_edge_equations[ieq][1] == facenormal_dim) { + neighbor->type.flip (ieq); /*ASSERT that flip is in correct direction */ + } + } + *neigh_face = t8_standalone_lut_type_face_to_neighface[el->type.to_ulong ()][face]; + } + else { + *neigh_face = face ^ 1; + } T8_ASSERT (element_is_valid ((t8_element_t *) neighbor)); /**check inside root*/ @@ -1091,7 +1226,7 @@ struct t8_standalone_scheme: public t8_scheme_helpers::element_is_valid (face)); return extrude_face ((t8_standalone_element *) face, elem, root_face); break; + case T8_ECLASS_TRIANGLE: + T8_ASSERT (t8_standalone_scheme::element_is_valid (face)); + return extrude_face ((t8_standalone_element *) face, elem, root_face); + break; default: - if constexpr (T8_ELEMENT_NUM_EQUATIONS[TEclass]) { - SC_ABORT ("Only implemented for hypercubes.\n"); - } + SC_ABORT ("Only useful for two-dimensional eclasses.\n"); return 0; break; } @@ -1133,7 +1270,7 @@ struct t8_standalone_scheme: public t8_scheme_helpers (elem, root_face, (t8_standalone_element *) boundary); return; + case T8_ECLASS_TRIANGLE: + compute_boundary_face (elem, root_face, + (t8_standalone_element *) boundary); + return; default: - if constexpr (T8_ELEMENT_NUM_EQUATIONS[TEclass]) { - SC_ABORT ("Only implemented for hypercubes.\n"); - } + SC_ABORT ("Only useful for two-dimensional eclasses.\n"); return; } } @@ -1323,7 +1462,7 @@ struct t8_standalone_scheme: public t8_scheme_helpers= 0); if constexpr (TEclass == T8_ECLASS_PYRAMID) { - SC_ABORT ("Not implemented yet.\n"); + const t8_linearidx_t two_to_l = 1LL << level; + const t8_linearidx_t eight_to_l = 1LL << (3 * level); + return ((eight_to_l << 2) - two_to_l) / 3; } return 1LL << (level * T8_ELEMENT_DIM[TEclass]); } @@ -1407,20 +1548,56 @@ struct t8_standalone_scheme: public t8_scheme_helpers *el = (const t8_standalone_element *) elem; double *current_ref_coords = (double *) ref_coords; double *current_out_coords = out_coords; t8_element_coord length = element_get_len (element_get_level (elem)); - for (size_t coord = 0; coord < num_coords; ++coord) { - for (int dim = 0; dim < T8_ELEMENT_DIM[TEclass]; ++dim) { - current_out_coords[dim] - = ((t8_standalone_element *) elem)->coords[dim] + current_ref_coords[dim] * length; + if constexpr (T8_ELEMENT_NUM_EQUATIONS[TEclass]) { + double tmp_coords[T8_ELEMENT_DIM[TEclass]] = { 0.0 }; + for (int idim = 0; idim < T8_ELEMENT_DIM[TEclass]; ++idim) { + for (int jdim = 0; jdim < T8_ELEMENT_DIM[TEclass]; ++jdim) { + tmp_coords[idim] + += t8_standalone_lut_transform_coords[el->type.to_ulong ()][idim][jdim] * current_ref_coords[jdim]; + } + if (TEclass == T8_ECLASS_PYRAMID && el->type == T8_DPYRAMID_SECOND_PYRA_TYPE) { + tmp_coords[idim] += 1.0; + } + } - current_out_coords[dim] /= (double) get_root_len (); + for (size_t coord = 0; coord < num_coords; ++coord) { + double tmp_out_coords[T8_ELEMENT_DIM[TEclass]] = {}; + for (int dim = 0; dim < T8_ELEMENT_DIM[TEclass]; ++dim) { + current_out_coords[dim] = el->coords[dim] + tmp_coords[dim] * length; + + current_out_coords[dim] /= (double) get_root_len (); + } + for (int dim = 0; dim < T8_ELEMENT_DIM[TEclass]; ++dim) { + for (int jdim = 0; jdim < T8_ELEMENT_DIM[TEclass]; ++jdim) { + tmp_out_coords[dim] + += t8_standalone_lut_backtransform_coords[0][dim][jdim] * current_out_coords[jdim]; + } + } + for (int dim = 0; dim < T8_ELEMENT_DIM[TEclass]; ++dim) { + current_out_coords[dim] = tmp_out_coords[dim]; + } + // current_out_coords = tmp_out_coords; + current_ref_coords += T8_ECLASS_MAX_DIM; + current_out_coords += T8_ELEMENT_DIM[TEclass]; } + } - current_ref_coords += T8_ECLASS_MAX_DIM; - current_out_coords += T8_ELEMENT_DIM[TEclass]; + else { + for (size_t coord = 0; coord < num_coords; ++coord) { + for (int dim = 0; dim < T8_ELEMENT_DIM[TEclass]; ++dim) { + current_out_coords[dim] = el->coords[dim] + current_ref_coords[dim] * length; + + current_out_coords[dim] /= (double) get_root_len (); + } + + current_ref_coords += T8_ECLASS_MAX_DIM; + current_out_coords += T8_ELEMENT_DIM[TEclass]; + } } } @@ -1574,13 +1751,12 @@ struct t8_standalone_scheme: public t8_scheme_helpers *el = (const t8_standalone_element *) elem; t8_debugf ("level: %i\n", el->level); - for (int i = 0; i < T8_ELEMENT_DIM[TEclass]; i++) { - t8_debugf ("x_%i: %i \n", i, el->coords[i]); + for (int idim = 0; idim < T8_ELEMENT_DIM[TEclass]; idim++) { + t8_debugf ("x_%i: %i \n", idim, el->coords[idim]); + } + for (int e_num = 0; e_num < T8_ELEMENT_NUM_EQUATIONS[TEclass]; e_num++) { + t8_debugf ("t_%i: %i \n", e_num, (int) el->type.test (e_num)); } - /** for (int e = 0; e < T8_ELEMENT_NUM_EQUATIONS[TEclass]; e++) { - * t8_debugf ("t_%i: %i \n", e, el->type[e]); - *} - * ToDo-Type */ } #endif @@ -1626,6 +1802,10 @@ struct t8_standalone_scheme: public t8_scheme_helperslevel, 1, sc_MPI_INT8_T, send_buffer, buffer_size, position, comm); SC_CHECK_MPI (mpiret); + + int8_t type = els[ielem]->type.to_ulong (); + mpiret = sc_MPI_Pack (&type, 1, sc_MPI_INT8_T, send_buffer, buffer_size, position, comm); + SC_CHECK_MPI (mpiret); } } @@ -1646,10 +1826,10 @@ struct t8_standalone_scheme: public t8_scheme_helperslevel), 1, sc_MPI_INT8_T, comm); SC_CHECK_MPI (mpiret); + + int8_t type; + mpiret = sc_MPI_Unpack (recvbuf, buffer_size, position, &type, 1, sc_MPI_INT8_T, comm); + SC_CHECK_MPI (mpiret); + els[ielem]->type = type; } } @@ -1734,7 +1919,7 @@ struct t8_standalone_scheme: public t8_scheme_helperstype = element_compute_type_at_level (ancestor, level); } /* The coordinates and the type of the ancestor are defined by the level. */ @@ -1855,7 +2040,14 @@ struct t8_standalone_scheme: public t8_scheme_helpers *el = (t8_standalone_element *) elem; + const int8_t type = el->type.to_ulong (); + int tmp_out_coords[T8_ELEMENT_DIM[TEclass]] = {}; + for (int idim = 0; idim < T8_ELEMENT_DIM[TEclass]; idim++) { + coords[idim] + = el->coords[idim] + t8_type_vertex_dim_to_binary[type][vertex][idim] * element_get_len (el->level); + } + for (int idim = 0; idim < T8_ELEMENT_DIM[TEclass]; ++idim) { + for (int jdim = 0; jdim < T8_ELEMENT_DIM[TEclass]; ++jdim) { + tmp_out_coords[idim] += t8_standalone_lut_backtransform_coords[0][idim][jdim] * coords[jdim]; + } + } + for (int idim = 0; idim < T8_ELEMENT_DIM[TEclass]; ++idim) { + coords[idim] = tmp_out_coords[idim]; + } + //coords = tmp_out_coords; } else { //Hypercubes @@ -1896,7 +2103,8 @@ struct t8_standalone_scheme: public t8_scheme_helpers[elem->type.to_ulong ()][face]; } } @@ -1936,10 +2144,9 @@ struct t8_standalone_scheme: public t8_scheme_helpers[elem->type.to_ulong ()][face]; } } - /** Given a face of an element that is also a hypercube face, determine if it is the boundary x_i == 1. * \param [in] elem The input element. * \param [in] face The input face. Needs to be a face of the hypercube the element is embedded in. @@ -1954,7 +2161,7 @@ struct t8_standalone_scheme: public t8_scheme_helpers[elem->type.to_ulong ()][face]; } } @@ -1970,16 +2177,16 @@ struct t8_standalone_scheme: public t8_scheme_helpers[root_type][root_face]; } } /** Get the eclass of the face for the element eclass * \return The eclass of the face. - * Note: Only implemented for hypercubes */ static constexpr t8_eclass_t - get_face_eclass () noexcept + get_root_face_eclass (const int root_face) noexcept { switch (TEclass) { case T8_ECLASS_VERTEX: @@ -1993,7 +2200,7 @@ struct t8_standalone_scheme: public t8_scheme_helpers[root_face]; } return T8_ECLASS_INVALID; } @@ -2043,11 +2250,16 @@ struct t8_standalone_scheme: public t8_scheme_helperstype = 0; + for (int ieq = 0; ieq < T8_ELEMENT_NUM_EQUATIONS[TEclass]; ieq++) { + const int ifaceeq = t8_standalone_lut_rootface_eq_to_faceeq[root_face][ieq]; + if (ifaceeq != -1) { + boundary->type[ifaceeq] = el->type[ieq]; + } + } + T8_ASSERT (t8_standalone_scheme::element_is_valid ((t8_element_t *) boundary)); } - T8_ASSERT (t8_standalone_scheme::element_is_valid ((t8_element_t *) boundary)); } - /** Given a boundary face inside a root tree's face construct * the element inside the root tree that has the given face as a * face. @@ -2093,6 +2305,7 @@ struct t8_standalone_scheme: public t8_scheme_helpers root_type = 0; + el->type = root_type; + for (int ieq = 0; ieq < T8_ELEMENT_NUM_EQUATIONS[TEclass]; ieq++) { + const int ifaceeq = t8_standalone_lut_rootface_eq_to_faceeq[root_face][ieq]; + if (ifaceeq != -1) { + el->type[ieq] = face->type[ifaceeq]; + } + } + /** Set those typebits, that are connected to the face_normaldim of root_face*/ + for (int ieq = 0; ieq < T8_ELEMENT_NUM_EQUATIONS[TEclass]; ieq++) { + const int facenormal_dim + = t8_standalone_lut_type_face_to_facenormal_dim[root_type.to_ulong ()][root_face]; + if (t8_type_edge_equations[ieq][0] == facenormal_dim) { + el->type[ieq] = !t8_standalone_lut_type_face_to_is_1_boundary[root_type.to_ulong ()][root_face]; + } + else if (t8_type_edge_equations[ieq][1] == facenormal_dim) { + el->type[ieq] = t8_standalone_lut_type_face_to_is_1_boundary[root_type.to_ulong ()][root_face]; + } + } + return t8_standalone_lut_type_rootface_to_face[el->type.to_ulong ()][root_face]; } } @@ -2152,9 +2384,38 @@ struct t8_standalone_scheme: public t8_scheme_helpers[root_face][idim]; } - return 0; + } + + /* For each typebit, consider the coordinate information between level and p->level |10...11|xxxx|0...0| + * of both inequality defining dimensions */ + + static inline t8_element_type + element_compute_type_at_level (const t8_standalone_element *el, const int level) + { + t8_element_type type = 0; + T8_ASSERT (0 <= el->level && el->level <= T8_ELEMENT_MAXLEVEL[TEclass]); + + for (int e = 0; e < T8_ELEMENT_NUM_EQUATIONS[TEclass]; e++) { + t8_element_coord coord_v0 = el->coords[t8_type_edge_equations[e][0]]; + t8_element_coord coord_v1 = el->coords[t8_type_edge_equations[e][1]]; + + coord_v0 = (coord_v0 << level) & ((1 << T8_ELEMENT_MAXLEVEL[TEclass]) - 1); + coord_v1 = (coord_v1 << level) & ((1 << T8_ELEMENT_MAXLEVEL[TEclass]) - 1); + + if (coord_v0 == coord_v1) { + type[e] = el->type[e]; + } + else if (coord_v0 < coord_v1) { + type |= (1 << e); + } + else { + T8_ASSERT (coord_v0 > coord_v1); + T8_ASSERT (!(type & (t8_element_type) (1 << e)).all ()); + } + } + return type; } }; @@ -2347,7 +2608,76 @@ t8_standalone_scheme::element_transform_face ([[maybe_unused [[maybe_unused]] const int sign, [[maybe_unused]] const int is_smaller_face) noexcept { - SC_ABORT ("Only implemented for hypercubes.\n"); + T8_ASSERT (0 <= orientation && orientation <= 2); + const t8_standalone_element *el1 = (const t8_standalone_element *) elem1; + t8_standalone_element *el2 = (t8_standalone_element *) elem2; + + const t8_element_coord length = element_get_len (el1->level); + + // copy elem1 in tmp so that results can directly be set in elem2 + t8_standalone_element tmp; + element_copy (elem1, (t8_element_t *) &tmp); + + /* + * The corners of the triangle are enumerated like this + * type 0 type 1 + * also root tree + * v_2 v_1 v_2 + * x x--x + * /| | / + * / | |/ + * x--x x + * v_0 v_1 v_0 + * + */ + if (sign) { + /* The tree faces have the same topological orientation, and + * thus we have to perform a coordinate switch. */ + if (el1->type == 0) { + tmp.coords[1] = el1->coords[0] - el1->coords[1]; + } + else { + tmp.coords[1] = el1->coords[0] - el1->coords[1] - length; + } + } + + if (!is_smaller_face && orientation != 0 && !sign) { + /* Translate orientation if triangle1 is not on the smaller face. + * sign = 0 sign = 1 + * 0 -> 0 0 -> 0 + * 1 -> 2 1 -> 1 + * 2 -> 1 2 -> 2 + */ + orientation = 3 - orientation; + } + + el2->level = tmp.level; + el2->type = tmp.type; + switch (orientation) { + case 0: + element_copy ((t8_element_t *) &tmp, (t8_element_t *) el2); + break; + case 1: + el2->coords[0] = get_root_len () - length - tmp.coords[1]; + if (tmp.type == 0) { + el2->coords[1] = tmp.coords[0] - tmp.coords[1]; + } + else { + el2->coords[1] = tmp.coords[0] - tmp.coords[1] - length; + } + break; + case 2: + if (tmp.type == 0) { + el2->coords[0] = get_root_len () - length + tmp.coords[1] - tmp.coords[0]; + } + else { + el2->coords[0] = get_root_len () + tmp.coords[1] - tmp.coords[0]; + } + el2->coords[1] = get_root_len () - length - tmp.coords[0]; + break; + default: + SC_ABORT_NOT_REACHED (); + } }; #endif /* T8_STANDALONE_IMPLEMENTATION_HXX */ diff --git a/src/t8_schemes/t8_standalone/t8_standalone_lut/t8_standalone_lut_prism.hxx b/src/t8_schemes/t8_standalone/t8_standalone_lut/t8_standalone_lut_prism.hxx new file mode 100644 index 0000000000..37ce9e5dc4 --- /dev/null +++ b/src/t8_schemes/t8_standalone/t8_standalone_lut/t8_standalone_lut_prism.hxx @@ -0,0 +1,293 @@ +/* + This file is part of t8code. + t8code is a C library to manage a collection (a forest) of multiple + connected adaptive space-trees of general element classes in parallel. + + Copyright (C) 2025 the developers + + t8code is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + t8code is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with t8code; if not, write to the Free Software Foundation, Inc., + 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. +*/ + +#ifndef T8_SELE_LUT_PRISM_HXX +#define T8_SELE_LUT_PRISM_HXX +// PRISMS +// clang-format off +template<> +inline constexpr int8_t t8_element_type_Iloc_to_childtype[1< +inline constexpr int8_t t8_element_type_Iloc_to_childcubeid[1< +inline constexpr int8_t t8_element_type_cubeid_to_parenttype[1< +inline constexpr int8_t t8_element_type_cubeid_to_Iloc[1< +inline constexpr int8_t t8_type_edge_equations[T8_ELEMENT_NUM_EQUATIONS[T8_ECLASS_PRISM]][2] = { + {0, 1} +}; + +template<> +inline constexpr int8_t t8_type_vertex_dim_to_binary[1 << T8_ELEMENT_NUM_EQUATIONS[T8_ECLASS_PRISM]] + [T8_ELEMENT_NUM_CORNERS[T8_ECLASS_PRISM]][T8_ELEMENT_DIM[T8_ECLASS_PRISM]] = { + { + {0, 0, 0}, + {1, 0, 0}, + {1, 1, 0}, + {0, 0, 1}, + {1, 0, 1}, + {1, 1, 1} + }, + { + /* PAY ATTENTION TO ORDERING OF ELEMENTS! What is the ordering? **/ + {0, 0, 0}, + {0, 1, 0}, + {1, 1, 0}, + {0, 0, 1}, + {0, 1, 1}, + {1, 1, 1} + } +}; + + +template<> +inline constexpr int8_t t8_standalone_lut_face_internal[1< +inline constexpr int8_t t8_standalone_lut_type_face_to_typebit[1< +inline constexpr int8_t t8_standalone_lut_type_face_to_sign[1< +inline constexpr int8_t t8_standalone_lut_type_face_to_facenormal_dim[1< +inline constexpr int8_t t8_standalone_lut_type_face_to_tree_face[1< +inline constexpr int8_t t8_standalone_lut_type_face_to_neighface[1< +inline constexpr int8_t t8_standalone_lut_type_cubeid_face_to_parentface[1< +inline constexpr int8_t t8_standalone_lut_type_childid_face_to_childface[1< +inline constexpr int8_t t8_standalone_lut_type_face_to_is_1_boundary[1< +inline constexpr int8_t t8_standalone_lut_type_face_to_last_facechilds_cubeid[1< +inline constexpr int8_t t8_standalone_lut_type_face_facechildid_to_childid[1< +inline constexpr int8_t t8_standalone_lut_rootface_dim_to_facedim[T8_ELEMENT_NUM_FACES[T8_ECLASS_PRISM]][T8_ELEMENT_DIM[T8_ECLASS_PRISM]]={ + {-1,0,1}, + {0,0,1}, + {0,-1,1}, + {0,1,-1}, + {0,1,-1} +}; +template<> +inline constexpr int8_t t8_standalone_lut_rootface_eq_to_faceeq[T8_ELEMENT_NUM_FACES[T8_ECLASS_PRISM]][T8_ELEMENT_NUM_EQUATIONS[T8_ECLASS_PRISM]]={ + {-1}, + {-1}, + {-1}, + {0}, + {0} +}; + +template<> +inline constexpr int8_t t8_standalone_lut_type_rootface_to_face[1< +inline constexpr int8_t t8_standalone_lut_cornerface[1< +inline constexpr int8_t t8_standalone_lut_facecorner[1< +inline constexpr int8_t t8_standalone_lut_transform_coords[1 << T8_ELEMENT_NUM_EQUATIONS[T8_ECLASS_PRISM]][T8_ELEMENT_DIM[T8_ECLASS_PRISM]][T8_ELEMENT_DIM[T8_ECLASS_PRISM]]={ + { + {1, 0, 0}, + {0, 1, 0}, + {0, 0, 1} + }, + { + {0, 1, 0}, + {1, 0, 0}, + {0, 0, 1} + } +}; + +template<> +inline constexpr int8_t t8_standalone_lut_backtransform_coords[1 << T8_ELEMENT_NUM_EQUATIONS[T8_ECLASS_PRISM]][T8_ELEMENT_DIM[T8_ECLASS_PRISM]][T8_ELEMENT_DIM[T8_ECLASS_PRISM]]={ + { + {1, 0, 0}, + {0, 1, 0}, + {0, 0, 1} + }, + { + {0, 1, 0}, + {1, 0, 0}, + {0, 0, 1} + } +}; + +template<> +inline constexpr t8_eclass_t t8_standalone_lut_rootface_to_eclass[T8_ELEMENT_NUM_FACES[T8_ECLASS_PRISM]]={ + T8_ECLASS_QUAD, + T8_ECLASS_QUAD, + T8_ECLASS_QUAD, + T8_ECLASS_TRIANGLE, + T8_ECLASS_TRIANGLE +}; +// clang-format on +#endif /* T8_SELE_LUT_PRISM_HXX */ diff --git a/src/t8_schemes/t8_standalone/t8_standalone_lut/t8_standalone_lut_pyra.hxx b/src/t8_schemes/t8_standalone/t8_standalone_lut/t8_standalone_lut_pyra.hxx new file mode 100644 index 0000000000..12ab0a5e1f --- /dev/null +++ b/src/t8_schemes/t8_standalone/t8_standalone_lut/t8_standalone_lut_pyra.hxx @@ -0,0 +1,450 @@ +/* + This file is part of t8code. + t8code is a C library to manage a collection (a forest) of multiple + connected adaptive space-trees of general element classes in parallel. + + Copyright (C) 2025 the developers + + t8code is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + t8code is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with t8code; if not, write to the Free Software Foundation, Inc., + 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. +*/ + +#ifndef T8_SELE_LUT_PYRA_HXX +#define T8_SELE_LUT_PYRA_HXX +// PYRAS +// clang-format off + +/** The first type of pyramids in the shape of a pyramid*/ +#define T8_DPYRAMID_FIRST_PYRA_TYPE 0 + +/** The second type of pyramids in the shape of a pyramid*/ +#define T8_DPYRAMID_SECOND_PYRA_TYPE 3 + + +template<> +inline constexpr int8_t t8_element_type_Iloc_to_childtype[1< +inline constexpr int8_t t8_element_type_Iloc_to_childcubeid[1< +inline constexpr int8_t t8_element_type_cubeid_to_parenttype[1< +inline constexpr int8_t t8_element_type_cubeid_to_Iloc[1< +inline constexpr int8_t t8_type_edge_equations[T8_ELEMENT_NUM_EQUATIONS[T8_ECLASS_PYRAMID]][2] = { + {1, 2}, + {0, 2} +}; + +template<> +inline constexpr int8_t t8_type_vertex_dim_to_binary[1 << T8_ELEMENT_NUM_EQUATIONS[T8_ECLASS_PYRAMID]] + [T8_ELEMENT_NUM_CORNERS[T8_ECLASS_PYRAMID]][T8_ELEMENT_DIM[T8_ECLASS_PYRAMID]] = { + { + {0, 0, 0}, + {1, 0, 0}, + {0, 1, 0}, + {1, 1, 0}, + {1, 1, 1} + }, + { + {0, 0, 0}, + {1, 0, 0}, + {1, 0, 1}, + {1, 1, 1}, + {-1, -1, -1} + }, + { + {0, 0, 0}, + {0, 1, 0}, + {0, 1, 1}, + {1, 1, 1}, + {-1, -1, -1} + }, + { + /* PAY ATTENTION TO ORDERING OF ELEMENTS! **/ + // {0, 0, 1}, + // {1, 0, 1}, + // {0, 1, 1}, + // {1, 1, 1}, + // {0, 0, 0} + {1, 1, 1}, + {0, 1, 1}, + {1, 0, 1}, + {0, 0, 1}, + {0, 0, 0} + } +}; + +template<> +inline constexpr int8_t t8_standalone_lut_face_internal[1< +inline constexpr int8_t t8_standalone_lut_type_face_to_typebit[1< +inline constexpr int8_t t8_standalone_lut_type_face_to_sign[1< +inline constexpr int8_t t8_standalone_lut_type_face_to_facenormal_dim[1< +inline constexpr int8_t t8_standalone_lut_type_face_to_neighface[1< +inline constexpr int8_t t8_standalone_lut_type_face_to_tree_face[1< +inline constexpr int8_t t8_standalone_lut_type_cubeid_face_to_parentface[1< +inline constexpr int8_t t8_standalone_lut_type_childid_face_to_childface[1< +inline constexpr int8_t t8_standalone_lut_type_face_to_is_1_boundary[1< +inline constexpr int8_t t8_standalone_lut_type_face_to_last_facechilds_cubeid[1< +inline constexpr int8_t t8_standalone_lut_type_face_facechildid_to_childid[1< +inline constexpr int8_t t8_standalone_lut_rootface_dim_to_facedim[T8_ELEMENT_NUM_FACES[T8_ECLASS_PYRAMID]][T8_ELEMENT_DIM[T8_ECLASS_PYRAMID]]={ + { 1,0,1}, + {-1,0,1}, + {0,1,1}, + {0,-1,1}, + {0,1,-1} +}; +template<> +inline constexpr int8_t t8_standalone_lut_rootface_eq_to_faceeq[T8_ELEMENT_NUM_FACES[T8_ECLASS_PYRAMID]][T8_ELEMENT_NUM_EQUATIONS[T8_ECLASS_PYRAMID]]={ + {0,-1}, + {0,-1}, + {-1,0}, + {-1,0}, + {-1,-1} +}; + +template<> +inline constexpr int8_t t8_standalone_lut_type_rootface_to_face[1< +inline constexpr int8_t t8_standalone_lut_cornerface[1< +inline constexpr int8_t t8_standalone_lut_facecorner[1< +inline constexpr int8_t t8_standalone_lut_transform_coords[1 << T8_ELEMENT_NUM_EQUATIONS[T8_ECLASS_PYRAMID]][T8_ELEMENT_DIM[T8_ECLASS_PYRAMID]][T8_ELEMENT_DIM[T8_ECLASS_PYRAMID]]={ + { + {1, 0, 0}, + {0, 1, 0}, + {0, 0, 1} + }, + { + {1, 0, 0}, + {0, 1, 0}, + {0, 0, 1} + }, + { + {0, 1, 0}, + {1, 0, 0}, + {0, 0, 1} + }, + { + {-1, 0, 0}, + {0, -1, 0}, + {0, 0, -1} + } +}; + +template<> +inline constexpr int8_t t8_standalone_lut_backtransform_coords[1 << T8_ELEMENT_NUM_EQUATIONS[T8_ECLASS_PYRAMID]][T8_ELEMENT_DIM[T8_ECLASS_PYRAMID]][T8_ELEMENT_DIM[T8_ECLASS_PYRAMID]]={ + { + {1, 0, 0}, + {0, 1, 0}, + {0, 0, 1} + }, + { + {1, 0, 0}, + {0, 1, 0}, + {0, 0, 1} + }, + { + {0, 1, 0}, + {1, 0, 0}, + {0, 0, 1} + }, + { + {-1, 0, 0}, + {0, -1, 0}, + {0, 0, -1} + } +}; + +template<> +inline constexpr t8_eclass_t t8_standalone_lut_rootface_to_eclass[T8_ELEMENT_NUM_FACES[T8_ECLASS_PYRAMID]]={ + T8_ECLASS_TRIANGLE, + T8_ECLASS_TRIANGLE, + T8_ECLASS_TRIANGLE, + T8_ECLASS_TRIANGLE, + T8_ECLASS_QUAD +}; +// clang-format on +#endif /* T8_SELE_LUT_PYRA_HXX */ diff --git a/src/t8_schemes/t8_standalone/t8_standalone_lut/t8_standalone_lut_tet.hxx b/src/t8_schemes/t8_standalone/t8_standalone_lut/t8_standalone_lut_tet.hxx new file mode 100644 index 0000000000..360c248cd8 --- /dev/null +++ b/src/t8_schemes/t8_standalone/t8_standalone_lut/t8_standalone_lut_tet.hxx @@ -0,0 +1,688 @@ +/* + This file is part of t8code. + t8code is a C library to manage a collection (a forest) of multiple + connected adaptive space-trees of general element classes in parallel. + + Copyright (C) 2025 the developers + + t8code is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + t8code is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with t8code; if not, write to the Free Software Foundation, Inc., + 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. +*/ + +#ifndef T8_SELE_LUT_TET_HXX +#define T8_SELE_LUT_TET_HXX +// TETS + +// clang-format off +template<> +inline constexpr int8_t t8_element_type_Iloc_to_childtype[1< +inline constexpr int8_t t8_element_type_Iloc_to_childcubeid[1< +inline constexpr int8_t t8_element_type_cubeid_to_parenttype[1< +inline constexpr int8_t t8_element_type_cubeid_to_Iloc[1<= y >= z*/ +template<> +inline constexpr int8_t t8_type_edge_equations[T8_ELEMENT_NUM_EQUATIONS[T8_ECLASS_TET]][2] = { + {0, 1}, + {1, 2}, + {0, 2} +}; + +template<> +inline constexpr int8_t t8_type_vertex_dim_to_binary[1 << T8_ELEMENT_NUM_EQUATIONS[T8_ECLASS_TET]] + [T8_ELEMENT_NUM_CORNERS[T8_ECLASS_TET]][T8_ELEMENT_DIM[T8_ECLASS_TET]] = { + { + {0, 0, 0}, + {1, 0, 0}, + {1, 1, 0}, + {1, 1, 1} + }, + { + {0, 0, 0}, + {0, 1, 0}, + {1, 1, 0}, + {1, 1, 1} + }, + { + {0, 0, 0}, + {1, 0, 0}, + {1, 0, 1}, + {1, 1, 1} + }, + { + {-1, -1, -1}, + {-1, -1, -1}, + {-1, -1, -1}, + {-1, -1, -1} + }, + { + {-1, -1, -1}, + {-1, -1, -1}, + {-1, -1, -1}, + {-1, -1, -1} + }, + { + {0, 0, 0}, + {0, 1, 0}, + {0, 1, 1}, + {1, 1, 1} + }, + { + {0, 0, 0}, + {0, 0, 1}, + {1, 0, 1}, + {1, 1, 1} + }, + { + {0, 0, 0}, + {0, 0, 1}, + {0, 1, 1}, + {1, 1, 1} + } +}; +template<> +inline constexpr int8_t t8_standalone_lut_face_internal[1< +inline constexpr int8_t t8_standalone_lut_type_face_to_typebit[1< +inline constexpr int8_t t8_standalone_lut_type_face_to_sign[1< +inline constexpr int8_t t8_standalone_lut_type_face_to_facenormal_dim[1< +inline constexpr int8_t t8_standalone_lut_type_face_to_is_1_boundary[1< +inline constexpr int8_t t8_standalone_lut_type_face_to_neighface[1< +inline constexpr int8_t t8_standalone_lut_type_face_to_tree_face[1< +inline constexpr int8_t t8_standalone_lut_type_cubeid_face_to_parentface[1< +inline constexpr int8_t t8_standalone_lut_type_childid_face_to_childface[1< +inline constexpr int8_t t8_standalone_lut_type_face_to_last_facechilds_cubeid[1< +inline constexpr int8_t t8_standalone_lut_type_face_facechildid_to_childid[1< +inline constexpr int8_t t8_standalone_lut_rootface_dim_to_facedim[T8_ELEMENT_NUM_FACES[T8_ECLASS_TET]][T8_ELEMENT_DIM[T8_ECLASS_TET]]={ + {-1,0,1}, + {0,0,1}, + {0,1,1}, + {0,1,-1} +}; +template<> +inline constexpr int8_t t8_standalone_lut_rootface_eq_to_faceeq[T8_ELEMENT_NUM_FACES[T8_ECLASS_TET]][T8_ELEMENT_NUM_EQUATIONS[T8_ECLASS_TET]]={ + {-1, 0,-1}, + {-1, 0, 0}, + { 0,-1, 0}, + { 0,-1,-1} +}; + +template<> +inline constexpr int8_t t8_standalone_lut_type_rootface_to_face[1< +inline constexpr int8_t t8_standalone_lut_cornerface[1< +inline constexpr int8_t t8_standalone_lut_facecorner[1< +inline constexpr int8_t t8_standalone_lut_transform_coords[1 << T8_ELEMENT_NUM_EQUATIONS[T8_ECLASS_TET]][T8_ELEMENT_DIM[T8_ECLASS_TET]][T8_ELEMENT_DIM[T8_ECLASS_TET]]={ + { + {1, 0, 0}, + {0, 0, 1}, + {0, 1, 0} + }, + { + {0, 0, 1}, + {1, 0, 0}, + {0, 1, 0} + }, + { + {1, 0, 0}, + {0, 1, 0}, + {0, 0, 1} + }, + { + {0, 0, 0}, + {0, 0, 0}, + {0, 0, 0} + }, + { + {0, 0, 0}, + {0, 0, 0}, + {0, 0, 0} + }, + { + {0, 1, 0}, + {1, 0, 0}, + {0, 0, 1} + }, + { + {0, 0, 1}, + {0, 1, 0}, + {1, 0, 0} + }, + { + {0, 1, 0}, + {0, 0, 1}, + {1, 0, 0} + } +}; + +template<> +inline constexpr int8_t t8_standalone_lut_backtransform_coords[1 << T8_ELEMENT_NUM_EQUATIONS[T8_ECLASS_TET]][T8_ELEMENT_DIM[T8_ECLASS_TET]][T8_ELEMENT_DIM[T8_ECLASS_TET]] = { + { + {1, 0, 0}, + {0, 0, 1}, + {0, 1, 0} + }, + { + {0, 1, 0}, + {0, 0, 1}, + {1, 0, 0} + }, + { + {1, 0, 0}, + {0, 1, 0}, + {0, 0, 1} + }, + { + {0, 0, 0}, + {0, 0, 0}, + {0, 0, 0} + }, + { + {0, 0, 0}, + {0, 0, 0}, + {0, 0, 0} + }, + { + {0, 1, 0}, + {1, 0, 0}, + {0, 0, 1} + }, + { + {0, 0, 1}, + {0, 1, 0}, + {1, 0, 0} + }, + { + {0, 0, 1}, + {1, 0, 0}, + {0, 1, 0} + } +}; + +template<> +inline constexpr t8_eclass_t t8_standalone_lut_rootface_to_eclass[T8_ELEMENT_NUM_FACES[T8_ECLASS_TET]]={ + T8_ECLASS_TRIANGLE, + T8_ECLASS_TRIANGLE, + T8_ECLASS_TRIANGLE, + T8_ECLASS_TRIANGLE +}; +// clang-format on +#endif /* T8_SELE_LUT_TET_HXX */ diff --git a/src/t8_schemes/t8_standalone/t8_standalone_lut/t8_standalone_lut_triangle.hxx b/src/t8_schemes/t8_standalone/t8_standalone_lut/t8_standalone_lut_triangle.hxx new file mode 100644 index 0000000000..078010cc24 --- /dev/null +++ b/src/t8_schemes/t8_standalone/t8_standalone_lut/t8_standalone_lut_triangle.hxx @@ -0,0 +1,251 @@ +/* + This file is part of t8code. + t8code is a C library to manage a collection (a forest) of multiple + connected adaptive space-trees of general element classes in parallel. + + Copyright (C) 2025 the developers + + t8code is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + t8code is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with t8code; if not, write to the Free Software Foundation, Inc., + 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. +*/ + +#ifndef T8_SELE_LUT_TRIANGLE_HXX +#define T8_SELE_LUT_TRIANGLE_HXX +// TRIANGLES +// clang-format off +template<> +inline constexpr int8_t t8_element_type_Iloc_to_childtype[1< +inline constexpr int8_t t8_element_type_Iloc_to_childcubeid[1< +inline constexpr int8_t t8_element_type_cubeid_to_parenttype[1< +inline constexpr int8_t t8_element_type_cubeid_to_Iloc[1< +inline constexpr int8_t t8_type_edge_equations[T8_ELEMENT_NUM_EQUATIONS[T8_ECLASS_TRIANGLE]][2] = { + {0, 1} +}; + +template<> +inline constexpr int8_t t8_type_vertex_dim_to_binary[1 << T8_ELEMENT_NUM_EQUATIONS[T8_ECLASS_TRIANGLE]] + [T8_ELEMENT_NUM_CORNERS[T8_ECLASS_TRIANGLE]][T8_ELEMENT_DIM[T8_ECLASS_TRIANGLE]] = { + { + {0, 0}, + {1, 0}, + {1, 1} + }, + { + /* PAY ATTENTION TO ORDERING OF ELEMENTS! What is the ordering? **/ + {0, 0}, + {0, 1}, + {1, 1} + } +}; + + +template<> +inline constexpr int8_t t8_standalone_lut_face_internal[1< +inline constexpr int8_t t8_standalone_lut_type_face_to_typebit[1< +inline constexpr int8_t t8_standalone_lut_type_face_to_sign[1< +inline constexpr int8_t t8_standalone_lut_type_face_to_is_1_boundary[1< +inline constexpr int8_t t8_standalone_lut_type_face_to_facenormal_dim[1< +inline constexpr int8_t t8_standalone_lut_type_face_to_neighface[1< +inline constexpr int8_t t8_standalone_lut_type_face_to_tree_face[1< +inline constexpr int8_t t8_standalone_lut_type_cubeid_face_to_parentface[1< +inline constexpr int8_t t8_standalone_lut_type_childid_face_to_childface[1< +inline constexpr int8_t t8_standalone_lut_type_face_to_last_facechilds_cubeid[1< +inline constexpr int8_t t8_standalone_lut_type_face_facechildid_to_childid[1< +inline constexpr int8_t t8_standalone_lut_rootface_dim_to_facedim[T8_ELEMENT_NUM_FACES[T8_ECLASS_TRIANGLE]][T8_ELEMENT_DIM[T8_ECLASS_TRIANGLE]]={ + {-1,0}, + {0,0}, + {0,-1} +}; +template<> +inline constexpr int8_t t8_standalone_lut_rootface_eq_to_faceeq[T8_ELEMENT_NUM_FACES[T8_ECLASS_TRIANGLE]][T8_ELEMENT_NUM_EQUATIONS[T8_ECLASS_TRIANGLE]]={ + {-1}, + {-1}, + {-1} +}; + +template<> +inline constexpr int8_t t8_standalone_lut_type_rootface_to_face[1< +inline constexpr int8_t t8_standalone_lut_cornerface[1< +inline constexpr int8_t t8_standalone_lut_facecorner[1< +inline constexpr int8_t t8_standalone_lut_transform_coords[1 << T8_ELEMENT_NUM_EQUATIONS[T8_ECLASS_TRIANGLE]][T8_ELEMENT_DIM[T8_ECLASS_TRIANGLE]][T8_ELEMENT_DIM[T8_ECLASS_TRIANGLE]]={ + { + {1, 0}, + {0, 1} + }, + { + {0, 1}, + {1, 0} + } +}; + +template<> +inline constexpr int8_t t8_standalone_lut_backtransform_coords[1 << T8_ELEMENT_NUM_EQUATIONS[T8_ECLASS_TRIANGLE]][T8_ELEMENT_DIM[T8_ECLASS_TRIANGLE]][T8_ELEMENT_DIM[T8_ECLASS_TRIANGLE]]={ + { + {1, 0}, + {0, 1} + }, + { + {0, 1}, + {1, 0} + } +}; + +template<> +inline constexpr t8_eclass_t t8_standalone_lut_rootface_to_eclass[T8_ELEMENT_NUM_FACES[T8_ECLASS_TRIANGLE]]={ + T8_ECLASS_LINE, + T8_ECLASS_LINE, + T8_ECLASS_LINE +}; + +// clang-format on +#endif /* T8_SELE_LUT_TRIANGLE_HXX */ diff --git a/test/t8_forest/t8_gtest_element_volume.cxx b/test/t8_forest/t8_gtest_element_volume.cxx index dbce03fc9f..54b0d19a7a 100644 --- a/test/t8_forest/t8_gtest_element_volume.cxx +++ b/test/t8_forest/t8_gtest_element_volume.cxx @@ -3,7 +3,7 @@ t8code is a C library to manage a collection (a forest) of multiple connected adaptive space-trees of general element classes in parallel. - Copyright (C) 2015 the developers + Copyright (C) 2025 the developers t8code is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by @@ -67,29 +67,37 @@ struct t8_forest_volume: public testing::TestWithParampyramid.level) * 3); - /* Ancestors switch the shape. A tetrahedron has a 1/16th of its parents volume. - * For all levels we already divided the control-volume by 8, hence we - * divide it by 2 once. */ - if (pyra->switch_shape_at_level > 0) { - control_volume /= 2; + int level = scheme->element_get_level (eclass, pyra); + for (int ilvl = level; ilvl > 0; ilvl--) { + t8_element_shape_t shape = scheme->element_get_shape (eclass, pyra); + scheme->element_get_parent (eclass, pyra, pyra); + t8_element_shape_t parent_shape = scheme->element_get_shape (eclass, pyra); + control_volume = control_volume * shape_factor[shape][parent_shape]; } return control_volume; @@ -112,7 +120,7 @@ TEST_P (t8_forest_volume, volume_check) const t8_element_t *element = t8_forest_get_leaf_element_in_tree (forest, itree, ielement); const double volume = t8_forest_element_volume (forest, itree, element); if (eclass == T8_ECLASS_PYRAMID) { - const double shape_volume = pyramid_control_volume ((t8_dpyramid_t *) element); + const double shape_volume = pyramid_control_volume ((t8_element_t *) element, eclass, scheme); EXPECT_NEAR (volume, shape_volume, T8_PRECISION_SQRT_EPS); } else { diff --git a/test/t8_schemes/t8_gtest_default_pyramid_count_leaves.cxx b/test/t8_schemes/t8_gtest_default_pyramid_count_leaves.cxx new file mode 100644 index 0000000000..22820a96aa --- /dev/null +++ b/test/t8_schemes/t8_gtest_default_pyramid_count_leaves.cxx @@ -0,0 +1,50 @@ +/* + This file is part of t8code. + t8code is a C library to manage a collection (a forest) of multiple + connected adaptive space-trees of general element classes in parallel. + + Copyright (C) 2015 the developers + + t8code is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + t8code is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with t8code; if not, write to the Free Software Foundation, Inc., + 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. +*/ + +#include +#include +#include + +TEST (class_element_leaves, test_default_pyramid_count_leaves_root) +{ + + const t8_scheme *scheme; + scheme = t8_scheme_new_default (); + t8_eclass_t eclass = T8_ECLASS_PYRAMID; + + const int maxlevel = scheme->get_maxlevel (eclass); + t8_gloidx_t compare_value = 1; + t8_gloidx_t sum1 = 1; + t8_gloidx_t sum2 = 1; + + for (int level = 0; level <= maxlevel; ++level) { + const t8_gloidx_t leaf_count = scheme->count_leaves_from_root (eclass, level); + ASSERT_EQ (leaf_count, compare_value) + << "Incorrect leaf count " << leaf_count << " at eclass " << t8_eclass_to_string[eclass] << " and level " << level + << " (expecting " << compare_value << ")"; + //Calculate the copare value for the default pyramid + sum1 *= 8; + sum2 *= 6; + compare_value = 2 * sum1 - sum2; + } + scheme->unref (); +} diff --git a/test/t8_schemes/t8_gtest_element_count_leaves.cxx b/test/t8_schemes/t8_gtest_element_count_leaves.cxx index 20be78302d..e676043468 100644 --- a/test/t8_schemes/t8_gtest_element_count_leaves.cxx +++ b/test/t8_schemes/t8_gtest_element_count_leaves.cxx @@ -23,6 +23,9 @@ #include #include #include +#include +#include +#include /* * In this file we test whether the t8_element_count_leaves{_from_root} @@ -44,38 +47,60 @@ struct class_element_leaves: public testing::TestWithParam (GetParam ()); eclass = scheme->get_eclass_scheme_eclass (eclass); + + scheme->element_new (eclass, 1, &element); } void TearDown () override { + scheme->element_destroy (eclass, 1, &element); scheme->unref (); } + t8_element_t *element; const t8_scheme *scheme; t8_eclass_t eclass; }; +static int +adapt_all ([[maybe_unused]] t8_forest_t forest, [[maybe_unused]] t8_forest_t forest_from, + [[maybe_unused]] t8_locidx_t which_tree, [[maybe_unused]] const t8_eclass_t tree_class, + [[maybe_unused]] t8_locidx_t lelement_id, [[maybe_unused]] const t8_scheme_c *scheme, + [[maybe_unused]] const int is_family, [[maybe_unused]] const int num_elements, + [[maybe_unused]] t8_element_t *elements[]) +{ + return 1; +} + TEST_P (class_element_leaves, test_element_count_leaves_root) { - const int maxlevel = scheme->get_maxlevel (eclass); +#if T8_TEST_LEVEL_INT >= 1 + const int maxlevel = 3; +#else + const int maxlevel = 5; +#endif t8_gloidx_t compare_value = 1; - t8_gloidx_t sum1 = 1; - t8_gloidx_t sum2 = 1; - + t8_gloidx_t test_value = 1; + t8_cmesh_t cmesh; + t8_cmesh_init (&cmesh); + t8_cmesh_new_from_class (cmesh, eclass, sc_MPI_COMM_WORLD); + t8_cmesh_ref (cmesh); + t8_forest_t forest = t8_forest_new_uniform (cmesh, scheme, 0, 0, sc_MPI_COMM_WORLD); + scheme->ref (); for (int level = 0; level <= maxlevel; ++level) { const t8_gloidx_t leaf_count = scheme->count_leaves_from_root (eclass, level); ASSERT_EQ (leaf_count, compare_value) << "Incorrect leaf count " << leaf_count << " at eclass " << t8_eclass_to_string[eclass] << " and level " << level << " (expecting " << compare_value << ")"; - /* Multiply the compare_value with 2^dim (= number of children per element) */ - if (eclass == T8_ECLASS_PYRAMID) { - sum1 *= 8; - sum2 *= 6; - compare_value = 2 * sum1 - sum2; - } - else { - compare_value *= 1 << t8_eclass_to_dimension[eclass]; + forest = t8_forest_new_adapt (forest, adapt_all, 0, 0, NULL); + compare_value = t8_forest_get_local_num_leaf_elements (forest); + + if (!scheme->refines_irregular (eclass)) { + test_value *= 1 << t8_eclass_to_dimension[eclass]; + EXPECT_EQ (test_value, compare_value); } } + t8_cmesh_unref (&cmesh); + t8_forest_unref (&forest); } /* Tests whether the leaf count for the same level is equal to 1