Skip to content

Commit

Permalink
Merge pull request #653 from DLR-AMR/addition-cmesh-example-squared-disk
Browse files Browse the repository at this point in the history
Cmesh example with curvilinear geometry: Squared Disk [1/n]
  • Loading branch information
Davknapp authored Sep 4, 2023
2 parents 6fb8ef6 + 20f259a commit ce61347
Show file tree
Hide file tree
Showing 9 changed files with 438 additions and 2 deletions.
2 changes: 2 additions & 0 deletions example/cmesh/Makefile.am
Original file line number Diff line number Diff line change
Expand Up @@ -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
86 changes: 86 additions & 0 deletions example/cmesh/t8_cmesh_geometry_examples.cxx
Original file line number Diff line number Diff line change
@@ -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 <t8.h> /* General t8code header, always include this. */
#include <t8_cmesh.h> /* cmesh definition and basic interface. */
#include <t8_forest/t8_forest_general.h> /* forest definition and basic interface. */
#include <t8_schemes/t8_default/t8_default_cxx.hxx> /* default refinement scheme. */
#include <t8_cmesh_vtk_writer.h> /* write file in vtu file */
#include <t8_forest/t8_forest_io.h>
#include <t8_cmesh/t8_cmesh_examples.h>

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;
}
3 changes: 3 additions & 0 deletions src/Makefile.am
Original file line number Diff line number Diff line change
Expand Up @@ -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 \
Expand Down Expand Up @@ -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 \
Expand Down
83 changes: 83 additions & 0 deletions src/t8_cmesh/t8_cmesh_examples.c
Original file line number Diff line number Diff line change
Expand Up @@ -22,9 +22,12 @@

#include <t8_cmesh.h>
#include <t8_cmesh/t8_cmesh_examples.h>
#include <t8_cmesh/t8_cmesh_helpers.h>
#include <t8_cmesh/t8_cmesh_geometry.h>
#include <t8_geometry/t8_geometry_implementations/t8_geometry_linear.h>
#include <t8_geometry/t8_geometry_implementations/t8_geometry_examples.h>
#include <t8_vec.h>
#include <t8_mat.h>
#include <t8_eclass.h>

/* TODO: In p4est a tree edge is joined with itself to denote a domain boundary.
Expand Down Expand Up @@ -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;
}
8 changes: 8 additions & 0 deletions src/t8_cmesh/t8_cmesh_examples.h
Original file line number Diff line number Diff line change
Expand Up @@ -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 */
122 changes: 122 additions & 0 deletions src/t8_geometry/t8_geometry_implementations/t8_geometry_examples.cxx
Original file line number Diff line number Diff line change
@@ -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 <t8_geometry/t8_geometry_implementations/t8_geometry_examples.hxx>
#include <t8_geometry/t8_geometry_helpers.h>
#include <t8_vec.h>

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 ();
Loading

0 comments on commit ce61347

Please sign in to comment.