diff --git a/example/cmesh/t8_cmesh_geometry_examples.cxx b/example/cmesh/t8_cmesh_geometry_examples.cxx index c11b092b81..462220c0c4 100644 --- a/example/cmesh/t8_cmesh_geometry_examples.cxx +++ b/example/cmesh/t8_cmesh_geometry_examples.cxx @@ -104,13 +104,13 @@ main (int argc, char **argv) */ { - const char *prefix_cmesh = "t8_squared_disk_cmesh"; - const char *prefix_forest = "t8_squared_disk_forest"; + const char *prefix_cmesh = "t8_quadrangulated_disk_cmesh"; + const char *prefix_forest = "t8_quadrangulated_disk_forest"; const int uniform_level = 5; const double radius = 1.0; - t8_cmesh_t cmesh = t8_cmesh_new_squared_disk (radius, comm); + t8_cmesh_t cmesh = t8_cmesh_new_quadrangulated_disk (radius, comm); t8_forest_t forest = t8_forest_new_uniform (cmesh, t8_scheme_new_default_cxx (), uniform_level, 0, comm); @@ -256,6 +256,26 @@ main (int argc, char **argv) t8_forest_unref (&forest); } + { + const char *prefix_cmesh = "t8_cubed_sphere_cmesh"; + const char *prefix_forest = "t8_cubed_sphere_forest"; + + const int uniform_level = 2; + const double radius = 1.0; + + t8_cmesh_t cmesh = t8_cmesh_new_cubed_sphere (radius, comm); + + t8_forest_t forest = t8_forest_new_uniform (cmesh, t8_scheme_new_default_cxx (), uniform_level, 0, comm); + + t8_cmesh_vtk_write_file (cmesh, prefix_cmesh); + t8_global_productionf ("Wrote %s.\n", prefix_cmesh); + + t8_write_forest_to_vtu (forest, prefix_forest); + t8_global_productionf ("Wrote %s.\n\n", prefix_forest); + + t8_forest_unref (&forest); + } + t8_global_productionf ("Done!\n"); /* Finalize the sc library */ diff --git a/src/t8_cmesh/t8_cmesh_examples.c b/src/t8_cmesh/t8_cmesh_examples.c index 62e0cfdf51..3a16725c25 100644 --- a/src/t8_cmesh/t8_cmesh_examples.c +++ b/src/t8_cmesh/t8_cmesh_examples.c @@ -2720,29 +2720,36 @@ t8_cmesh_new_row_of_cubes (t8_locidx_t num_trees, const int set_attributes, cons } t8_cmesh_t -t8_cmesh_new_squared_disk (const double radius, sc_MPI_Comm comm) +t8_cmesh_new_quadrangulated_disk (const double radius, sc_MPI_Comm comm) { /* Initialization of the mesh */ t8_cmesh_t cmesh; t8_cmesh_init (&cmesh); - const double ri = 0.5 * radius; + const double ri = 0.6 * radius; const double ro = radius; const double xi = ri / M_SQRT2; - const double yi = ri / M_SQRT2; + const double yi = xi; const double xo = ro / M_SQRT2; - const double yo = ro / M_SQRT2; + const double yo = xo; - const int ntrees = 5; /* Number of cmesh elements resp. trees. */ - const int nverts = 4; /* Number of vertices per cmesh element. */ + const int nquads = 3; /* Number of quads in the upper-right quadrant. */ + const int nyturns = 4; /* Number of turns around y-axis. */ + + const int ntrees = nyturns * nquads; /* Number of cmesh elements resp. trees. */ + const int nverts = t8_eclass_num_vertices[T8_ECLASS_QUAD]; + + /* Fine tuning parameter to expand the center squares a bit for more equal + * element sizes. */ + const double s = 1.2; /* Arrays for the face connectivity computations via vertices. */ double all_verts[ntrees * T8_ECLASS_MAX_CORNERS * T8_ECLASS_MAX_DIM]; t8_eclass_t all_eclasses[ntrees]; - t8_geometry_c *geometry = t8_geometry_squared_disk_new (); + t8_geometry_c *geometry = t8_geometry_quadrangulated_disk_new (); t8_cmesh_register_geometry (cmesh, geometry); /* Defitition of the tree class. */ @@ -2751,44 +2758,42 @@ t8_cmesh_new_squared_disk (const double radius, sc_MPI_Comm comm) all_eclasses[itree] = T8_ECLASS_QUAD; } - /* Central quad. */ - { - const double vertices[4][3] = { { -xi, -yi, 0.0 }, { xi, -yi, 0.0 }, { -xi, yi, 0.0 }, { xi, yi, 0.0 } }; - - t8_cmesh_set_tree_vertices (cmesh, 0, (double *) vertices, 4); - - /* itree = 0; */ - for (int ivert = 0; ivert < nverts; ivert++) { - for (int icoord = 0; icoord < T8_ECLASS_MAX_DIM; icoord++) { - all_verts[T8_3D_TO_1D (ntrees, T8_ECLASS_MAX_CORNERS, T8_ECLASS_MAX_DIM, 0, ivert, icoord)] - = vertices[ivert][icoord]; - } - } - } - - /* Four quads framing the central quad. */ - { - const double vertices[4][3] = { { -xi, yi, 0.0 }, { xi, yi, 0.0 }, { -xo, yo, 0.0 }, { xo, yo, 0.0 } }; + /* Vertices of upper right quarter of the disk. */ + const double vertices_mid[4][3] = { { 0.0, 0.0, 0.0 }, { s * xi, 0.0, 0.0 }, { 0.0, s * yi, 0.0 }, { xi, yi, 0.0 } }; + const double vertices_top[4][3] = { { 0.0, s * yi, 0.0 }, { xi, yi, 0.0 }, { 0.0, yo, 0.0 }, { xo, yo, 0.0 } }; + const double vertices_bot[4][3] = { { s * xi, 0.0, 0.0 }, { xi, yi, 0.0 }, { xo, 0.0, 0.0 }, { xo, yo, 0.0 } }; - for (int itree = 1; itree < ntrees; itree++) { - double rot_mat[3][3]; - double rot_vertices[4][3]; + int itree = 0; + for (int iturn = 0; iturn < 4; iturn++) { + double rot_mat[3][3]; + double rot_vertices_mid[4][3]; + double rot_vertices_top[4][3]; + double rot_vertices_bot[4][3]; - t8_mat_init_zrot (rot_mat, (itree - 1) * 0.5 * M_PI); + t8_mat_init_zrot (rot_mat, -iturn * 0.5 * M_PI); - for (int i = 0; i < 4; i++) { - t8_mat_mult_vec (rot_mat, &(vertices[i][0]), &(rot_vertices[i][0])); - } + for (int i = 0; i < 4; i++) { + t8_mat_mult_vec (rot_mat, &(vertices_mid[i][0]), &(rot_vertices_mid[i][0])); + t8_mat_mult_vec (rot_mat, &(vertices_top[i][0]), &(rot_vertices_top[i][0])); + t8_mat_mult_vec (rot_mat, &(vertices_bot[i][0]), &(rot_vertices_bot[i][0])); + } - t8_cmesh_set_tree_vertices (cmesh, itree, (double *) rot_vertices, 4); + t8_cmesh_set_tree_vertices (cmesh, itree + 0, (double *) rot_vertices_mid, 4); + t8_cmesh_set_tree_vertices (cmesh, itree + 1, (double *) rot_vertices_top, 4); + t8_cmesh_set_tree_vertices (cmesh, itree + 2, (double *) rot_vertices_bot, 4); - for (int ivert = 0; ivert < nverts; ivert++) { - for (int icoord = 0; icoord < T8_ECLASS_MAX_DIM; icoord++) { - all_verts[T8_3D_TO_1D (ntrees, T8_ECLASS_MAX_CORNERS, T8_ECLASS_MAX_DIM, itree, ivert, icoord)] - = rot_vertices[ivert][icoord]; - } + for (int ivert = 0; ivert < nverts; ivert++) { + for (int icoord = 0; icoord < T8_ECLASS_MAX_DIM; icoord++) { + all_verts[T8_3D_TO_1D (ntrees, T8_ECLASS_MAX_CORNERS, T8_ECLASS_MAX_DIM, itree + 0, ivert, icoord)] + = rot_vertices_mid[ivert][icoord]; + all_verts[T8_3D_TO_1D (ntrees, T8_ECLASS_MAX_CORNERS, T8_ECLASS_MAX_DIM, itree + 1, ivert, icoord)] + = rot_vertices_top[ivert][icoord]; + all_verts[T8_3D_TO_1D (ntrees, T8_ECLASS_MAX_CORNERS, T8_ECLASS_MAX_DIM, itree + 2, ivert, icoord)] + = rot_vertices_bot[ivert][icoord]; } } + + itree = itree + nquads; } /* Face connectivity. */ @@ -3219,3 +3224,120 @@ t8_cmesh_new_cubed_spherical_shell (const double inner_radius, const double shel t8_cmesh_new_quadrangulated_spherical_surface, inner_radius, shell_thickness, num_levels, num_layers, comm); } + +t8_cmesh_t +t8_cmesh_new_cubed_sphere (const double radius, sc_MPI_Comm comm) +{ + /* Initialization of the mesh. */ + t8_cmesh_t cmesh; + t8_cmesh_init (&cmesh); + + const double ri = 0.6 * radius; + const double ro = radius; + + const double SQRT3 = 1.7320508075688772; + const double xi = ri / SQRT3; + const double yi = xi; + const double zi = xi; + + const double xo = ro / SQRT3; + const double yo = xo; + const double zo = xo; + + const int nhexs = 4; /* Number of hexs in the front-upper-right octant. */ + const int nzturns = 4; /* Number of turns around z-axis. */ + const int nyturns = 2; /* Number of turns around y-axis. */ + + const int ntrees = nyturns * nzturns * nhexs; /* Number of cmesh elements resp. trees. */ + const int nverts = t8_eclass_num_vertices[T8_ECLASS_HEX]; + + /* Fine tuning parameter to expand the center hex a bit for more equal + * element sizes. */ + const double s = 1.2; + + /* Arrays for the face connectivity computations via vertices. */ + double all_verts[ntrees * T8_ECLASS_MAX_CORNERS * T8_ECLASS_MAX_DIM]; + t8_eclass_t all_eclasses[ntrees]; + + t8_geometry_c *geometry = t8_geometry_cubed_sphere_new (); + t8_cmesh_register_geometry (cmesh, geometry); + + /* Defitition of the tree class. */ + for (int itree = 0; itree < ntrees; itree++) { + t8_cmesh_set_tree_class (cmesh, itree, T8_ECLASS_HEX); + all_eclasses[itree] = T8_ECLASS_HEX; + } + + const double vertices_mid[8][3] + = { { 0.0, 0.0, 0.0 }, { s * xi, 0.0, 0.0 }, { 0.0, s * yi, 0.0 }, { xi, yi, 0.0 }, + { 0.0, 0.0, s * zi }, { xi, 0.0, zi }, { 0.0, yi, zi }, { xi, yi, zi } }; + const double vertices_top[8][3] = { { 0.0, s * yi, 0.0 }, { xi, yi, 0.0 }, { 0.0, yi, zi }, { xi, yi, zi }, + { 0.0, yo, 0.0 }, { xo, yo, 0.0 }, { 0.0, yo, zo }, { xo, yo, zo } }; + const double vertices_bot[8][3] = { { s * xi, 0.0, 0.0 }, { xi, yi, 0.0 }, { xi, 0.0, zi }, { xi, yi, zi }, + { xo, 0.0, 0.0 }, { xo, yo, 0.0 }, { xo, 0.0, zo }, { xo, yo, zo } }; + const double vertices_zen[8][3] = { { 0.0, 0.0, s * zi }, { xi, 0.0, zi }, { 0.0, yi, zi }, { xi, yi, zi }, + { 0.0, 0.0, zo }, { xo, 0.0, zo }, { 0.0, yo, zo }, { xo, yo, zo } }; + + int itree = 0; + for (int yturn = 0; yturn < nyturns; yturn++) { + double yrot_mat[3][3]; + t8_mat_init_yrot (yrot_mat, yturn * M_PI); + + for (int zturn = 0; zturn < nzturns; zturn++) { + double zrot_mat[3][3]; + + double tmp_vertices_mid[8][3]; + double tmp_vertices_top[8][3]; + double tmp_vertices_bot[8][3]; + double tmp_vertices_zen[8][3]; + + double rot_vertices_mid[8][3]; + double rot_vertices_top[8][3]; + double rot_vertices_bot[8][3]; + double rot_vertices_zen[8][3]; + + t8_mat_init_zrot (zrot_mat, -zturn * 0.5 * M_PI); + + for (int i = 0; i < 8; i++) { + t8_mat_mult_vec (zrot_mat, &(vertices_mid[i][0]), &(tmp_vertices_mid[i][0])); + t8_mat_mult_vec (zrot_mat, &(vertices_top[i][0]), &(tmp_vertices_top[i][0])); + t8_mat_mult_vec (zrot_mat, &(vertices_bot[i][0]), &(tmp_vertices_bot[i][0])); + t8_mat_mult_vec (zrot_mat, &(vertices_zen[i][0]), &(tmp_vertices_zen[i][0])); + } + + for (int i = 0; i < 8; i++) { + t8_mat_mult_vec (yrot_mat, &(tmp_vertices_mid[i][0]), &(rot_vertices_mid[i][0])); + t8_mat_mult_vec (yrot_mat, &(tmp_vertices_top[i][0]), &(rot_vertices_top[i][0])); + t8_mat_mult_vec (yrot_mat, &(tmp_vertices_bot[i][0]), &(rot_vertices_bot[i][0])); + t8_mat_mult_vec (yrot_mat, &(tmp_vertices_zen[i][0]), &(rot_vertices_zen[i][0])); + } + + t8_cmesh_set_tree_vertices (cmesh, itree + 0, (double *) rot_vertices_mid, 8); + t8_cmesh_set_tree_vertices (cmesh, itree + 1, (double *) rot_vertices_top, 8); + t8_cmesh_set_tree_vertices (cmesh, itree + 2, (double *) rot_vertices_bot, 8); + t8_cmesh_set_tree_vertices (cmesh, itree + 3, (double *) rot_vertices_zen, 8); + + for (int ivert = 0; ivert < nverts; ivert++) { + for (int icoord = 0; icoord < T8_ECLASS_MAX_DIM; icoord++) { + all_verts[T8_3D_TO_1D (ntrees, T8_ECLASS_MAX_CORNERS, T8_ECLASS_MAX_DIM, itree + 0, ivert, icoord)] + = rot_vertices_mid[ivert][icoord]; + all_verts[T8_3D_TO_1D (ntrees, T8_ECLASS_MAX_CORNERS, T8_ECLASS_MAX_DIM, itree + 1, ivert, icoord)] + = rot_vertices_top[ivert][icoord]; + all_verts[T8_3D_TO_1D (ntrees, T8_ECLASS_MAX_CORNERS, T8_ECLASS_MAX_DIM, itree + 2, ivert, icoord)] + = rot_vertices_bot[ivert][icoord]; + all_verts[T8_3D_TO_1D (ntrees, T8_ECLASS_MAX_CORNERS, T8_ECLASS_MAX_DIM, itree + 2, ivert, icoord)] + = rot_vertices_zen[ivert][icoord]; + } + } + + itree = itree + nhexs; + } + } + + /* Face connectivity. */ + t8_cmesh_set_join_by_vertices (cmesh, ntrees, all_eclasses, all_verts, NULL, 0); + + /* Commit the mesh */ + t8_cmesh_commit (cmesh, comm); + return cmesh; +} diff --git a/src/t8_cmesh/t8_cmesh_examples.h b/src/t8_cmesh/t8_cmesh_examples.h index 13d5f11077..1750d9a240 100644 --- a/src/t8_cmesh/t8_cmesh_examples.h +++ b/src/t8_cmesh/t8_cmesh_examples.h @@ -313,13 +313,13 @@ t8_cmesh_new_long_brick_pyramid (sc_MPI_Comm comm, int num_cubes); t8_cmesh_t t8_cmesh_new_row_of_cubes (t8_locidx_t num_trees, const int set_attributes, const int do_partition, sc_MPI_Comm comm); -/** Construct a squared disk of given radius. +/** Construct a quadrangulated disk of given radius. * \param [in] radius Radius of the sphere. * \param [in] comm The MPI communicator used to commit the cmesh * \return A cmesh representing the spherical surface. */ t8_cmesh_t -t8_cmesh_new_squared_disk (const double radius, sc_MPI_Comm comm); +t8_cmesh_new_quadrangulated_disk (const double radius, sc_MPI_Comm comm); /** Construct a triangulated spherical surface of given radius: octahedron version. * \param [in] radius Radius of the sphere. @@ -382,6 +382,14 @@ t8_cmesh_t t8_cmesh_new_cubed_spherical_shell (const double inner_radius, const double shell_thickness, const int num_levels, const int num_layers, sc_MPI_Comm comm); +/** Construct a cubed sphere of given radius. + * \param [in] inner_radius Radius of the inner side of the shell. + * \param [in] comm The MPI communicator used to commit the cmesh + * \return A cmesh representing the spherical surface. + */ +t8_cmesh_t +t8_cmesh_new_cubed_sphere (const double radius, sc_MPI_Comm comm); + T8_EXTERN_C_END (); #endif /* !T8_CMESH_EXAMPLES */ diff --git a/src/t8_geometry/t8_geometry_implementations/t8_geometry_examples.cxx b/src/t8_geometry/t8_geometry_implementations/t8_geometry_examples.cxx index 885c81877a..1ac8dc7e58 100644 --- a/src/t8_geometry/t8_geometry_implementations/t8_geometry_examples.cxx +++ b/src/t8_geometry/t8_geometry_implementations/t8_geometry_examples.cxx @@ -25,20 +25,16 @@ #include void -t8_geometry_squared_disk::t8_geom_evaluate (t8_cmesh_t cmesh, t8_gloidx_t gtreeid, const double *ref_coords, - const size_t num_coords, double *out_coords) const +t8_geometry_quadrangulated_disk::t8_geom_evaluate (t8_cmesh_t cmesh, t8_gloidx_t gtreeid, const double *ref_coords, + const size_t num_coords, double *out_coords) const { - if (num_coords != 1) - SC_ABORT ("Error: Batch computation of geometry not yet supported."); - double n[3]; /* Normal vector. */ double r[3]; /* Radial vector. */ double s[3]; /* Radial vector for the corrected coordinates. */ double p[3]; /* Vector on the plane resp. quad. */ - /* Center square. */ - if (gtreeid == 0) { - + /* Center quads. */ + if (gtreeid % 3 == 0) { for (size_t i_coord = 0; i_coord < num_coords; i_coord++) { size_t offset = 3 * i_coord; @@ -52,60 +48,44 @@ t8_geometry_squared_disk::t8_geom_evaluate (t8_cmesh_t cmesh, t8_gloidx_t gtreei return; } - /* Four squares framing the central one. */ - { - const double center_ref[3] = { 0.5, 0.5, 0.0 }; - t8_geom_linear_interpolation (center_ref, active_tree_vertices, 3, 2, n); - - /* Normalize vector `n`. */ - const double norm = sqrt (n[0] * n[0] + n[1] * n[1]); - n[0] = n[0] / norm; - n[1] = n[1] / norm; - } + /* Normal vector along one of the straight edges of the quad. */ + n[0] = active_tree_vertices[0 + 0]; + n[1] = active_tree_vertices[0 + 1]; + n[2] = active_tree_vertices[0 + 2]; + t8_vec_normalize (n); - { - /* Radial vector parallel to one of the tilted edges of the quad. */ - r[0] = active_tree_vertices[0]; - r[1] = active_tree_vertices[1]; - - /* Normalize vector `r`. */ - const double norm = sqrt (r[0] * r[0] + r[1] * r[1]); - r[0] = r[0] / norm; - r[1] = r[1] / norm; - } + /* Radial vector parallel to one of the tilted edges of the quad. */ + r[0] = active_tree_vertices[9 + 0]; + r[1] = active_tree_vertices[9 + 1]; + r[2] = active_tree_vertices[9 + 2]; + t8_vec_normalize (r); for (size_t i_coord = 0; i_coord < num_coords; i_coord++) { size_t offset = 3 * i_coord; - const double x = ref_coords[offset + 0]; - const double y = ref_coords[offset + 1]; + const double x_ref = ref_coords[offset + 0]; + const double y_ref = ref_coords[offset + 1]; { double corr_ref_coords[3]; /* Correction in order to rectify elements near the corners. */ - corr_ref_coords[0] = tan (0.5 * M_PI * (x - 0.5)) * 0.5 + 0.5; - corr_ref_coords[1] = y; + corr_ref_coords[0] = tan (0.25 * M_PI * x_ref); + corr_ref_coords[1] = y_ref; corr_ref_coords[2] = 0.0; /* Compute and normalize vector `s`. */ t8_geom_linear_interpolation (corr_ref_coords, active_tree_vertices, 3, 2, s); - - const double norm = sqrt (s[0] * s[0] + s[1] * s[1]); - s[0] = s[0] / norm; - s[1] = s[1] / norm; + t8_vec_normalize (s); } t8_geom_linear_interpolation (ref_coords + offset, active_tree_vertices, 3, 2, p); /* Compute intersection of line with a plane. */ - const double out_radius = (p[0] * n[0] + p[1] * n[1]) / (r[0] * n[0] + r[1] * n[1]); - - const double blend = y * out_radius; /* y \in [0,1] */ - const double dnelb = 1.0 - y; + const double R = (p[0] * n[0] + p[1] * n[1]) / (r[0] * n[0] + r[1] * n[1]); - out_coords[offset + 0] = dnelb * p[0] + blend * s[0]; - out_coords[offset + 1] = dnelb * p[1] + blend * s[1]; + out_coords[offset + 0] = (1.0 - y_ref) * p[0] + y_ref * R * s[0]; + out_coords[offset + 1] = (1.0 - y_ref) * p[1] + y_ref * R * s[1]; out_coords[offset + 2] = 0.0; } } @@ -349,6 +329,66 @@ t8_geometry_cubed_spherical_shell::t8_geom_evaluate (t8_cmesh_t cmesh, t8_gloidx t8_geom_evaluate_sphere_quad_hex (active_tree_vertices, 3, ref_coords, num_coords, out_coords); } +void +t8_geometry_cubed_sphere::t8_geom_evaluate (t8_cmesh_t cmesh, t8_gloidx_t gtreeid, const double *ref_coords, + const size_t num_coords, double *out_coords) const +{ + double n[3]; /* Normal vector. */ + double r[3]; /* Radial vector. */ + double s[3]; /* Radial vector for the corrected coordinates. */ + double p[3]; /* Vector on the plane resp. quad. */ + + /* Center hex. */ + if (gtreeid % 4 == 0) { + for (size_t i_coord = 0; i_coord < num_coords; i_coord++) { + const size_t offset = 3 * i_coord; + t8_geom_linear_interpolation (ref_coords + offset, active_tree_vertices, 3, 3, p); + t8_vec_copy (p, out_coords + offset); + } + return; + } + + t8_vec_copy (active_tree_vertices, n); + t8_vec_normalize (n); + + t8_vec_copy (active_tree_vertices + 7 * 3, r); + t8_vec_normalize (r); + + const double inv_denominator = 1.0 / t8_vec_dot (r, n); + + for (size_t i_coord = 0; i_coord < num_coords; i_coord++) { + const size_t offset = 3 * i_coord; + + const double x_ref = ref_coords[offset + 0]; + const double y_ref = ref_coords[offset + 1]; + const double z_ref = ref_coords[offset + 2]; + + { + double corr_ref_coords[3]; + + /* Correction in order to rectify elements near the corners. Note, this + * is probably not the most accurate correction but it does a decent enough job. + */ + corr_ref_coords[0] = tan (0.25 * M_PI * x_ref); + corr_ref_coords[1] = tan (0.25 * M_PI * y_ref); + corr_ref_coords[2] = z_ref; + + /* Compute and normalize vector `s`. */ + t8_geom_linear_interpolation (corr_ref_coords, active_tree_vertices, 3, 3, s); + t8_vec_normalize (s); + } + + t8_geom_linear_interpolation (ref_coords + offset, active_tree_vertices, 3, 3, p); + + /* Compute intersection of line with a plane. */ + const double out_radius = t8_vec_dot (p, n) * inv_denominator; + + /* Linear blend from flat to curved: `out_coords = (1.0 - z_ref)*p + z_ref_ * out_radius * s`. */ + t8_vec_axy (p, out_coords + offset, 1.0 - z_ref); + t8_vec_axpy (s, out_coords + offset, z_ref * out_radius); + } +} + T8_EXTERN_C_BEGIN (); void @@ -362,9 +402,9 @@ t8_geometry_destroy (t8_geometry_c **geom) /* Satisfy the C interface from t8_geometry_linear.h. */ t8_geometry_c * -t8_geometry_squared_disk_new () +t8_geometry_quadrangulated_disk_new () { - t8_geometry_squared_disk *geom = new t8_geometry_squared_disk (); + t8_geometry_quadrangulated_disk *geom = new t8_geometry_quadrangulated_disk (); return (t8_geometry_c *) geom; } @@ -396,4 +436,11 @@ t8_geometry_cubed_spherical_shell_new () return (t8_geometry_c *) geom; } +t8_geometry_c * +t8_geometry_cubed_sphere_new () +{ + t8_geometry_cubed_sphere *geom = new t8_geometry_cubed_sphere (); + return (t8_geometry_c *) geom; +} + T8_EXTERN_C_END (); diff --git a/src/t8_geometry/t8_geometry_implementations/t8_geometry_examples.h b/src/t8_geometry/t8_geometry_implementations/t8_geometry_examples.h index d7933da542..8af8d89a90 100644 --- a/src/t8_geometry/t8_geometry_implementations/t8_geometry_examples.h +++ b/src/t8_geometry/t8_geometry_implementations/t8_geometry_examples.h @@ -38,11 +38,11 @@ T8_EXTERN_C_BEGIN (); void t8_geometry_destroy (t8_geometry_c **geom); -/** Create a new squared_disk geometry. +/** Create a new quadrangulated_disk geometry. * \return A pointer to an allocated geometry struct. */ t8_geometry_c * -t8_geometry_squared_disk_new (); +t8_geometry_quadrangulated_disk_new (); /** Create a new triangulated_spherical_surface geometry. * \return A pointer to an allocated geometry struct. @@ -68,6 +68,12 @@ t8_geometry_cubed_spherical_shell_new (); t8_geometry_c * t8_geometry_prismed_spherical_shell_new (); +/** Create a new cubed sphere geometry. + * \return A pointer to an allocated geometry struct. + */ +t8_geometry_c * +t8_geometry_cubed_sphere_new (); + T8_EXTERN_C_END (); #endif /* T8_GEOMETRY_EXAMPLE_H */ diff --git a/src/t8_geometry/t8_geometry_implementations/t8_geometry_examples.hxx b/src/t8_geometry/t8_geometry_implementations/t8_geometry_examples.hxx index 0f979e3e09..32f78200fa 100644 --- a/src/t8_geometry/t8_geometry_implementations/t8_geometry_examples.hxx +++ b/src/t8_geometry/t8_geometry_implementations/t8_geometry_examples.hxx @@ -32,11 +32,11 @@ /** This geometry maps five quads to a disk. */ -struct t8_geometry_squared_disk: public t8_geometry_with_vertices +struct t8_geometry_quadrangulated_disk: public t8_geometry_with_vertices { public: /* Basic constructor that sets the dimension and the name. */ - t8_geometry_squared_disk (): t8_geometry_with_vertices (2, "t8_squared_disk") + t8_geometry_quadrangulated_disk (): t8_geometry_with_vertices (2, "t8_quadrangulated_disk") { } @@ -48,19 +48,18 @@ struct t8_geometry_squared_disk: public t8_geometry_with_vertices * \param [in] num_coords The number of points to map. * \param [out] out_coords The mapped coordinates in physical space of \a ref_coords. The length is \a num_coords * 3. * - * This routine expects an input mesh of five squares looking like this: + * This routine expects an input mesh of 4 * 3 squares looking like this. + * Note, only the upper right quadrant is shown for the sake of clarity. * - * +----------+ - * |\___1____/| - * | | | | - * |4| 0 |2| - * | |______| | - * |/ 3 \| - * +----------+ + * ------+ + * |_1__/| + * | 0 |2| + * |___|_| + * * * The central quad (id = 0) is mapped as is, while the outer edges of * other four quads are stretched onto a circle with a radius determined by - * the four outer corners (denoted by "+") in the schematic above. + * the outer corner (denoted by "+") in the schematic above. * */ void @@ -265,4 +264,39 @@ class t8_geometry_prismed_spherical_shell: public t8_geometry_with_vertices { /* Load tree data is inherited from t8_geometry_with_vertices. */ }; +/** This geometry maps specifically arranged hexahedrons to a sphere. + */ +class t8_geometry_cubed_sphere: public t8_geometry_with_vertices { + public: + /* Basic constructor that sets the dimension and the name. */ + t8_geometry_cubed_sphere (): t8_geometry_with_vertices (3, "t8_geometry_cubed_sphere") + { + } + + /** + * Maps specifically arranged hexahedrons to a sphere. + * \param [in] cmesh The cmesh in which the point lies. + * \param [in] gtreeid The global tree (of the cmesh) in which the reference point is. + * \param [in] ref_coords Array of \a dimension x \a num_coords many entries, specifying a point in /f$ [0,1]^\mathrm{dim} /f$. + * \param [in] num_coords The number of points to map. + * \param [out] out_coords The mapped coordinates in physical space of \a ref_coords. The length is \a num_coords * 3. + * + * This routine expects an input mesh of prism arranged as octahedron or similar. + * + */ + void + t8_geom_evaluate (t8_cmesh_t cmesh, t8_gloidx_t gtreeid, const double *ref_coords, const size_t num_coords, + double out_coords[3]) const; + + /* Jacobian, not implemented. */ + void + t8_geom_evaluate_jacobian (t8_cmesh_t cmesh, t8_gloidx_t gtreeid, const double *ref_coords, const size_t num_coords, + double *jacobian) const + { + SC_ABORT_NOT_REACHED (); + } + + /* Load tree data is inherited from t8_geometry_with_vertices. */ +}; + #endif /* T8_GEOMETRY_EXAMPLES_HXX */ diff --git a/src/t8_vec.h b/src/t8_vec.h index 94010cfeb7..5465ae3df3 100644 --- a/src/t8_vec.h +++ b/src/t8_vec.h @@ -57,6 +57,18 @@ t8_vec_normalize (double vec[3]) } } +/** Make a copy of a vector. + * \param [in] vec_in + * \param [out] vec_out + */ +static inline void +t8_vec_copy (const double vec_in[3], double vec_out[3]) +{ + for (int i = 0; i < 3; i++) { + vec_out[i] = vec_in[i]; + } +} + /** Euclidean distance of X and Y. * \param [in] vec_x A 3D vector. * \param [in] vec_y A 3D vector.