diff --git a/example/cmesh/Makefile.am b/example/cmesh/Makefile.am index 605c62037b..d5ad9154bd 100644 --- a/example/cmesh/Makefile.am +++ b/example/cmesh/Makefile.am @@ -7,10 +7,12 @@ bin_PROGRAMS += \ example/cmesh/t8_cmesh_create_partitioned \ example/cmesh/t8_cmesh_refine \ example/cmesh/t8_cmesh_set_join_by_vertices \ + example/cmesh/t8_cmesh_geometry_examples \ example/cmesh/t8_cmesh_hypercube_pad example_cmesh_t8_cmesh_partition_SOURCES = example/cmesh/t8_cmesh_partition.cxx example_cmesh_t8_cmesh_refine_SOURCES = example/cmesh/t8_cmesh_refine.cxx example_cmesh_t8_cmesh_set_join_by_vertices_SOURCES = example/cmesh/t8_cmesh_set_join_by_vertices.cxx +example_cmesh_t8_cmesh_geometry_examples_SOURCES = example/cmesh/t8_cmesh_geometry_examples.cxx example_cmesh_t8_cmesh_create_partitioned_SOURCES = example/cmesh/t8_cmesh_create_partitioned.cxx example_cmesh_t8_cmesh_hypercube_pad_SOURCES = example/cmesh/t8_cmesh_hypercube_pad.cxx diff --git a/example/cmesh/t8_cmesh_geometry_examples.cxx b/example/cmesh/t8_cmesh_geometry_examples.cxx new file mode 100644 index 0000000000..8f9f09f6d1 --- /dev/null +++ b/example/cmesh/t8_cmesh_geometry_examples.cxx @@ -0,0 +1,86 @@ +/* + 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 types in parallel. + + Copyright (C) 2023 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. +*/ + +/* Show-case several cmesh examples with curvilinear geometries. */ + +#include /* General t8code header, always include this. */ +#include /* cmesh definition and basic interface. */ +#include /* forest definition and basic interface. */ +#include /* default refinement scheme. */ +#include /* write file in vtu file */ +#include +#include + +int +main (int argc, char **argv) +{ + /* + * Initialization. + */ + + /* Initialize MPI. This has to happen before we initialize sc or t8code. */ + int mpiret = sc_MPI_Init (&argc, &argv); + /* Error check the MPI return value. */ + SC_CHECK_MPI (mpiret); + + /* Initialize the sc library, has to happen before we initialize t8code. */ + sc_init (sc_MPI_COMM_WORLD, 1, 1, NULL, SC_LP_PRODUCTION); + /* Initialize t8code with log level SC_LP_PRODUCTION. See sc.h for more info on the log levels. */ + t8_init (SC_LP_PRODUCTION); + + /* We will use MPI_COMM_WORLD as a communicator. */ + sc_MPI_Comm comm = sc_MPI_COMM_WORLD; + + /* + * Creation of several meshes and storing them to disk. + */ + + { + const char *prefix_cmesh = "t8_squared_disk_cmesh"; + const char *prefix_forest = "t8_squared_disk_forest"; + + const int uniform_level = 5; + const double radius = 1.0; + + t8_cmesh_t cmesh = t8_cmesh_new_squared_disk (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, 1.0); + t8_global_productionf ("Wrote %s.\n", prefix_cmesh); + + t8_forest_write_vtk_ext (forest, prefix_forest, 1, 1, 1, 1, 0, 1, 0, 0, NULL); + t8_global_productionf ("Wrote %s.\n\n", prefix_forest); + + t8_forest_unref (&forest); + } + + /* More examples will be added soon. */ + + /* Finalize the sc library */ + sc_finalize (); + + mpiret = sc_MPI_Finalize (); + SC_CHECK_MPI (mpiret); + + return 0; +} diff --git a/src/Makefile.am b/src/Makefile.am index db738ed554..28f686305e 100644 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -65,10 +65,12 @@ libt8_installed_headers_geometry_impl = \ src/t8_geometry/t8_geometry_implementations/t8_geometry_linear.h \ src/t8_geometry/t8_geometry_implementations/t8_geometry_linear_axis_aligned.h \ src/t8_geometry/t8_geometry_implementations/t8_geometry_analytic.hxx \ + src/t8_geometry/t8_geometry_implementations/t8_geometry_examples.h \ src/t8_geometry/t8_geometry_implementations/t8_geometry_occ.h \ src/t8_geometry/t8_geometry_implementations/t8_geometry_occ.hxx \ src/t8_geometry/t8_geometry_implementations/t8_geometry_linear.hxx \ src/t8_geometry/t8_geometry_implementations/t8_geometry_linear_axis_aligned.hxx \ + src/t8_geometry/t8_geometry_implementations/t8_geometry_examples.hxx \ src/t8_geometry/t8_geometry_implementations/t8_geometry_zero.hxx libt8_installed_headers_vtk = \ src/t8_vtk/t8_vtk_reader.hxx \ @@ -124,6 +126,7 @@ libt8_compiled_sources = \ src/t8_geometry/t8_geometry_implementations/t8_geometry_linear.cxx \ src/t8_geometry/t8_geometry_implementations/t8_geometry_linear_axis_aligned.cxx \ src/t8_geometry/t8_geometry_implementations/t8_geometry_zero.cxx \ + src/t8_geometry/t8_geometry_implementations/t8_geometry_examples.cxx \ src/t8_forest/t8_forest_partition.cxx src/t8_forest/t8_forest_cxx.cxx \ src/t8_forest/t8_forest_private.c src/t8_forest/t8_forest_vtk.cxx \ src/t8_forest/t8_forest_ghost.cxx src/t8_forest/t8_forest_iterate.cxx \ diff --git a/src/t8_cmesh/t8_cmesh_examples.c b/src/t8_cmesh/t8_cmesh_examples.c index 17f4f0c926..32a5fef373 100644 --- a/src/t8_cmesh/t8_cmesh_examples.c +++ b/src/t8_cmesh/t8_cmesh_examples.c @@ -22,9 +22,12 @@ #include #include +#include #include #include +#include #include +#include #include /* TODO: In p4est a tree edge is joined with itself to denote a domain boundary. @@ -2655,3 +2658,83 @@ t8_cmesh_new_row_of_cubes (t8_locidx_t num_trees, const int set_attributes, cons t8_cmesh_commit (cmesh, comm); return cmesh; } + +t8_cmesh_t +t8_cmesh_new_squared_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 ro = radius; + + const double xi = ri / M_SQRT2; + const double yi = ri / M_SQRT2; + + const double xo = ro / M_SQRT2; + const double yo = ro / M_SQRT2; + + const int ntrees = 5; /* Number of cmesh elements resp. trees. */ + const int nverts = 4; /* Number of vertices per cmesh element. */ + + /* 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_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_QUAD); + 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 } }; + + for (int itree = 1; itree < ntrees; itree++) { + double rot_mat[3][3]; + double rot_vertices[4][3]; + + t8_mat_init_zrot (rot_mat, (itree - 1) * 0.5 * M_PI); + + for (int i = 0; i < 4; i++) { + t8_mat_mult_vec (rot_mat, &(vertices[i][0]), &(rot_vertices[i][0])); + } + + t8_cmesh_set_tree_vertices (cmesh, itree, (double *) rot_vertices, 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]; + } + } + } + } + + /* 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 9282498ec9..986966e01b 100644 --- a/src/t8_cmesh/t8_cmesh_examples.h +++ b/src/t8_cmesh/t8_cmesh_examples.h @@ -314,6 +314,14 @@ 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. + * \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_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 new file mode 100644 index 0000000000..a211b8b51a --- /dev/null +++ b/src/t8_geometry/t8_geometry_implementations/t8_geometry_examples.cxx @@ -0,0 +1,122 @@ +/* + 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) 2023 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 + +void +t8_geometry_squared_disk::t8_geom_evaluate (t8_cmesh_t cmesh, t8_gloidx_t gtreeid, const double *ref_coords, + double out_coords[3]) 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. */ + + const double x = ref_coords[0]; + const double y = ref_coords[1]; + + t8_locidx_t ltreeid = t8_cmesh_get_local_id (cmesh, gtreeid); + double *tree_vertices = t8_cmesh_get_tree_vertices (cmesh, ltreeid); + + t8_geom_linear_interpolation (ref_coords, tree_vertices, 3, 2, p); + + /* Center square. */ + if (gtreeid == 0) { + out_coords[0] = p[0]; + out_coords[1] = p[1]; + out_coords[2] = 0.0; + + return; + } + + /* Four squares framing the central one. */ + { + const double center_ref[3] = { 0.5, 0.5, 0.0 }; + t8_geom_linear_interpolation (center_ref, 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; + } + + { + /* Radial vector parallel to one of the tilted edges of the quad. */ + r[0] = tree_vertices[0]; + r[1] = 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; + } + + { + 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[2] = 0.0; + + /* Compute and normalize vector `s`. */ + t8_geom_linear_interpolation (corr_ref_coords, 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; + } + + { + 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; + + out_coords[0] = dnelb * p[0] + blend * s[0]; + out_coords[1] = dnelb * p[1] + blend * s[1]; + out_coords[2] = 0.0; + } +} + +T8_EXTERN_C_BEGIN (); + +void +t8_geometry_destroy (t8_geometry_c **geom) +{ + T8_ASSERT (geom != NULL); + + delete *geom; + *geom = NULL; +} + +/* Satisfy the C interface from t8_geometry_linear.h. */ +t8_geometry_c * +t8_geometry_squared_disk_new () +{ + t8_geometry_squared_disk *geom = new t8_geometry_squared_disk (); + 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 new file mode 100644 index 0000000000..eba1db0458 --- /dev/null +++ b/src/t8_geometry/t8_geometry_implementations/t8_geometry_examples.h @@ -0,0 +1,49 @@ +/* + 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) 2023 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. +*/ + +/** \file t8_geometry_examples.h + * This header provides the C interface to create a cubed_sphere geometry. + */ + +#ifndef T8_GEOMETRY_EXAMPLES_H +#define T8_GEOMETRY_EXAMPLES_H + +#include +#include + +T8_EXTERN_C_BEGIN (); + +/** Destroy a geometry object. + * \param [in,out] geom A pointer to a geometry object. Set to NULL on output. + */ +void +t8_geometry_destroy (t8_geometry_c **geom); + +/** Create a new squared_disk geometry. + * \return A pointer to an allocated geometry struct. + */ +t8_geometry_c * +t8_geometry_squared_disk_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 new file mode 100644 index 0000000000..b069c8969d --- /dev/null +++ b/src/t8_geometry/t8_geometry_implementations/t8_geometry_examples.hxx @@ -0,0 +1,83 @@ +/* + 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) 2023 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. +*/ + +/** \file t8_geometry_example.hxx + * Various mappings for several cmesh examples. + */ + +#ifndef T8_GEOMETRY_EXAMPLES_HXX +#define T8_GEOMETRY_EXAMPLES_HXX + +#include +#include + +/** This geometry maps five quads to a disk. + */ +class t8_geometry_squared_disk: public t8_geometry { + public: + /* Basic constructor that sets the dimension and the name. */ + t8_geometry_squared_disk (): t8_geometry (2, "t8_squared_disk") + { + } + + /** + * Map five quads to a disk. + * \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 many entries, specifying a point in [0,1]^dimension. + * \param [out] out_coords The mapped coordinates in physical space of \a ref_coords. + * + * This routine expects an input mesh of five squares looking like this: + * + * +----------+ + * |\___1____/| + * | | | | + * |4| 0 |2| + * | |______| | + * |/ 3 \| + * +----------+ + * + * 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. + * + */ + void + t8_geom_evaluate (t8_cmesh_t cmesh, t8_gloidx_t gtreeid, const double *ref_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, double *jacobian) const + { + SC_ABORT_NOT_REACHED (); + } + + /* Load tree data is empty since we have no tree data. + * We need to provide an implementation anyways. */ + void + t8_geom_load_tree_data (t8_cmesh_t cmesh, t8_gloidx_t gtreeid) + { + SC_ABORT_NOT_REACHED (); + } +}; + +#endif /* T8_GEOMETRY_EXAMPLES_HXX */ diff --git a/src/t8_mat.h b/src/t8_mat.h index 373d6cfdc3..dfbc2ddf37 100644 --- a/src/t8_mat.h +++ b/src/t8_mat.h @@ -104,7 +104,7 @@ t8_mat_init_zrot (double mat[3][3], const double angle) * \param [in,out] b 3-vector. */ static inline void -t8_mat_mult_vec (const double mat[3][3], const double a[3], double b[3]) +t8_mat_mult_vec (double mat[3][3], const double a[3], double b[3]) { b[0] = mat[0][0] * a[0] + mat[0][1] * a[1] + mat[0][2] * a[2]; b[1] = mat[1][0] * a[0] + mat[1][1] * a[1] + mat[1][2] * a[2]; @@ -117,7 +117,7 @@ t8_mat_mult_vec (const double mat[3][3], const double a[3], double b[3]) * \param [int,out] C 3x3-matrix. */ static inline void -t8_mat_mult_mat (const double A[3][3], const double B[3][3], double C[3][3]) +t8_mat_mult_mat (double A[3][3], double B[3][3], double C[3][3]) { for (int i = 0; i < 3; i++) { for (int j = 0; j < 3; j++) {