commit bb92a7fe3a83bc25f8c741e5f5c37ea5c4137952
parent e78ad3f15077b773595450e12bf0a14fe8290853
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Tue, 26 Jul 2022 18:39:02 +0200
Add rnatm_write_vtk_octree function
Diffstat:
3 files changed, 236 insertions(+), 1 deletion(-)
diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt
@@ -67,7 +67,8 @@ set(RNATM_FILES_SRC
rnatm_mesh.c
rnatm_octree.c
rnatm_properties.c
- rnatm_voxel_partition.c)
+ rnatm_voxel_partition.c
+ rnatm_write_vtk_octree.c)
set(RNATM_FILES_INC
rnatm_c.h
rnatm_log.h
diff --git a/src/rnatm.h b/src/rnatm.h
@@ -144,4 +144,9 @@ RNATM_API res_T
rnatm_validate
(const struct rnatm* atm);
+RNATM_API res_T
+rnatm_write_vtk_octree
+ (const struct rnatm* atm,
+ FILE* stream);
+
#endif /* RNATM_H */
diff --git a/src/rnatm_write_vtk_octree.c b/src/rnatm_write_vtk_octree.c
@@ -0,0 +1,229 @@
+/* Copyright (C) 2022 Centre National de la Recherche Scientifique
+ * Copyright (C) 2022 Institut de Physique du Globe de Paris
+ * Copyright (C) 2022 |Méso|Star> (contact@meso-star.com)
+ * Copyright (C) 2022 Université de Reims Champagne-Ardenne
+ * Copyright (C) 2022 Université de Versaille Saint-Quentin
+ * Copyright (C) 2022 Université Paul Sabatier (contact@laplace.univ-tlse.fr)
+ *
+ * This program 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 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program 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 this program. If not, see <http://www.gnu.org/licenses/>. */
+
+#include "rnatm.h"
+#include "rnatm_c.h"
+#include "rnatm_log.h"
+#include "rnatm_voxel.h"
+
+#include <rsys/clock_time.h>
+#include <rsys/dynamic_array_double.h>
+#include <rsys/dynamic_array_size_t.h>
+#include <rsys/hash_table.h>
+
+#include <star/svx.h>
+
+struct vertex {
+ double x;
+ double y;
+ double z;
+};
+
+static char
+vertex_eq(const struct vertex* v0, const struct vertex* v1)
+{
+ return eq_eps(v0->x, v1->x, 1.e-6)
+ && eq_eps(v0->y, v1->y, 1.e-6)
+ && eq_eps(v0->z, v1->z, 1.e-6);
+}
+
+/* Generate the htable_vertex data structure */
+#define HTABLE_NAME vertex
+#define HTABLE_KEY struct vertex
+#define HTABLE_DATA size_t
+#define HTABLE_KEY_FUNCTOR_EQ vertex_eq
+#include <rsys/hash_table.h>
+
+/* Temporary data structure used to dump the octree data into a VTK file */
+struct octree_data {
+ struct htable_vertex vertex2id; /* Map a coordinate to its vertex id */
+ struct darray_double vertices; /* Array of unique vertices */
+ struct darray_double data; /* List of registered leaf data */
+ struct darray_size_t cells; /* Ids of the cell vertices */
+ const struct rnatm* atm;
+};
+
+/*******************************************************************************
+ * Helper functions
+ ******************************************************************************/
+static void
+octree_data_init
+ (const struct rnatm* atm,
+ struct octree_data* data)
+{
+ ASSERT(data);
+ htable_vertex_init(atm->allocator, &data->vertex2id);
+ darray_double_init(atm->allocator, &data->vertices);
+ darray_double_init(atm->allocator, &data->data);
+ darray_size_t_init(atm->allocator, &data->cells);
+ data->atm = atm;
+}
+
+static void
+octree_data_release(struct octree_data* data)
+{
+ ASSERT(data);
+ htable_vertex_release(&data->vertex2id);
+ darray_double_release(&data->vertices);
+ darray_double_release(&data->data);
+ darray_size_t_release(&data->cells);
+}
+
+static void
+register_leaf
+ (const struct svx_voxel* leaf,
+ const size_t ileaf,
+ void* context)
+{
+ struct octree_data* ctx = context;
+ const float* vx = NULL;
+ struct vertex v[8];
+ double kext_min;
+ double kext_max;
+ int i;
+ ASSERT(leaf && ctx);
+ (void)ileaf;
+
+ /* Compute the leaf vertices */
+ v[0].x = leaf->lower[0]; v[0].y = leaf->lower[1]; v[0].z = leaf->lower[2];
+ v[1].x = leaf->upper[0]; v[1].y = leaf->lower[1]; v[1].z = leaf->lower[2];
+ v[2].x = leaf->lower[0]; v[2].y = leaf->upper[1]; v[2].z = leaf->lower[2];
+ v[3].x = leaf->upper[0]; v[3].y = leaf->upper[1]; v[3].z = leaf->lower[2];
+ v[4].x = leaf->lower[0]; v[4].y = leaf->lower[1]; v[4].z = leaf->upper[2];
+ v[5].x = leaf->upper[0]; v[5].y = leaf->lower[1]; v[5].z = leaf->upper[2];
+ v[6].x = leaf->lower[0]; v[6].y = leaf->upper[1]; v[6].z = leaf->upper[2];
+ v[7].x = leaf->upper[0]; v[7].y = leaf->upper[1]; v[7].z = leaf->upper[2];
+
+ FOR_EACH(i, 0, 8) {
+ size_t *pid = htable_vertex_find(&ctx->vertex2id, v+i);
+ size_t id;
+ if(pid) {
+ id = *pid;
+ } else { /* Register the leaf vertex */
+ id = darray_double_size_get(&ctx->vertices)/3;
+ CHK(RES_OK == htable_vertex_set(&ctx->vertex2id, v+i, &id));
+ CHK(RES_OK == darray_double_push_back(&ctx->vertices, &v[i].x));
+ CHK(RES_OK == darray_double_push_back(&ctx->vertices, &v[i].y));
+ CHK(RES_OK == darray_double_push_back(&ctx->vertices, &v[i].z));
+ }
+ /* Add the vertex id to the leaf cell */
+ CHK(RES_OK == darray_size_t_push_back(&ctx->cells, &id));
+ }
+
+ /* Register leaf data */
+ vx = leaf->data;
+ kext_min = vx[voxel_idata(RNATM_RADCOEF_Kext, RNATM_SVX_OP_MIN)];
+ kext_max = vx[voxel_idata(RNATM_RADCOEF_Kext, RNATM_SVX_OP_MAX)];
+ CHK(RES_OK == darray_double_push_back(&ctx->data, &kext_min));
+ CHK(RES_OK == darray_double_push_back(&ctx->data, &kext_max));
+}
+
+
+/*******************************************************************************
+ * Exported function
+ ******************************************************************************/
+res_T
+rnatm_write_vtk_octree(const struct rnatm* atm, FILE* stream)
+{
+ char buf[128];
+ struct time t0, t1;
+ struct octree_data data;
+ const double* leaf_data;
+ size_t nvertices;
+ size_t ncells;
+ size_t i;
+ res_T res = RES_OK;
+
+ if(!atm || !stream) {
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ /* Start time recording */
+ log_info(atm, "write VTK octree\n");
+ time_current(&t0);
+
+ octree_data_init(atm, &data);
+
+ /* Register leaf data */
+ SVX(tree_for_each_leaf(atm->octree, register_leaf, &data));
+ nvertices = darray_double_size_get(&data.vertices)/3/*#coords per vertex*/;
+ ncells = darray_size_t_size_get(&data.cells)/8/*#ids per cell*/;
+#ifndef NDEBUG
+ {
+ struct svx_tree_desc octree_desc = SVX_TREE_DESC_NULL;
+ SVX(tree_get_desc(atm->octree, &octree_desc));
+ ASSERT(ncells == octree_desc.nleaves);
+ }
+#endif
+
+ /* Write headers */
+ fprintf(stream, "# vtk DataFile Version 2.0\n");
+ fprintf(stream, "Radiative coefficients\n");
+ fprintf(stream, "ASCII\n");
+ fprintf(stream, "DATASET UNSTRUCTURED_GRID\n");
+
+ /* Write vertex coordinates */
+ fprintf(stream, "POINTS %lu float\n", (unsigned long)nvertices);
+ FOR_EACH(i, 0, nvertices) {
+ fprintf(stream, "%g %g %g\n",
+ SPLIT3(darray_double_cdata_get(&data.vertices) + i*3));
+ }
+
+ /* Write the cells */
+ fprintf(stream, "CELLS %lu %lu\n",
+ (unsigned long)ncells,
+ (unsigned long)(ncells*(8/*#verts per cell*/ + 1/*1st field of a cell*/)));
+ FOR_EACH(i, 0, ncells) {
+ fprintf(stream, "8 %lu %lu %lu %lu %lu %lu %lu %lu\n",
+ (unsigned long)darray_size_t_cdata_get(&data.cells)[i*8+0],
+ (unsigned long)darray_size_t_cdata_get(&data.cells)[i*8+1],
+ (unsigned long)darray_size_t_cdata_get(&data.cells)[i*8+2],
+ (unsigned long)darray_size_t_cdata_get(&data.cells)[i*8+3],
+ (unsigned long)darray_size_t_cdata_get(&data.cells)[i*8+4],
+ (unsigned long)darray_size_t_cdata_get(&data.cells)[i*8+5],
+ (unsigned long)darray_size_t_cdata_get(&data.cells)[i*8+6],
+ (unsigned long)darray_size_t_cdata_get(&data.cells)[i*8+7]);
+ }
+
+ /* Write the cell type */
+ fprintf(stream, "CELL_TYPES %lu\n", (unsigned long)ncells);
+ FOR_EACH(i, 0, ncells) fprintf(stream, "11\n");
+
+ /* Write the cell data */
+ leaf_data = darray_double_cdata_get(&data.data);
+ fprintf(stream, "CELL_DATA %lu\n", (unsigned long)ncells);
+ fprintf(stream, "SCALARS Kext double 2\n");
+ fprintf(stream, "LOOKUP_TABLE default\n");
+ FOR_EACH(i, 0, ncells) {
+ fprintf(stream, "%g %g\n", leaf_data[i*2+0], leaf_data[i*2+1]);
+ }
+ octree_data_release(&data);
+
+ /* Log elapsed time */
+ time_sub(&t0, time_current(&t1), &t0);
+ time_dump(&t0, TIME_ALL, NULL, buf, sizeof(buf));
+ log_info(atm, "octree written in %s\n", buf);
+
+exit:
+ return res;
+error:
+ goto exit;
+}
