commit 93a23f206fd9258a1c3304c2e05dadd0d09865f5
parent 1f74f76ba128ff6593160cf05132a4ebff9b17ce
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Mon, 25 Jan 2021 21:35:27 +0100
Add the atrstm_dump_svx_octree function
Diffstat:
3 files changed, 245 insertions(+), 3 deletions(-)
diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt
@@ -55,6 +55,7 @@ set(VERSION ${VERSION_MAJOR}.${VERSION_MINOR}.${VERSION_PATCH})
set(ATRSTM_FILES_SRC
atrstm.c
atrstm_cache.c
+ atrstm_dump_svx_octree.c
atrstm_log.c
atrstm_partition.c
atrstm_radcoefs.c
diff --git a/src/atrstm.h b/src/atrstm.h
@@ -217,6 +217,9 @@ struct atrstm_dump_svx_octree_args {
static const struct atrstm_dump_svx_octree_args
ATRSTM_DUMP_SVX_OCTREE_ARGS_DEFAULT = ATRSTM_DUMP_SVX_OCTREE_ARGS_DEFAULT__;
+typedef double atrstm_radcoefs_T[ATRSTM_RADCOEFS_COUNT__];
+typedef double atrstm_radcoefs_svx_T[ATRSTM_RADCOEFS_COUNT__][ATRSTM_SVX_OPS_COUNT__];
+
/* Forward declaration of extern data types */
struct logger;
struct mem_allocator;
@@ -252,19 +255,19 @@ ATRSTM_API res_T
atrstm_fetch_radcoefs
(const struct atrstm* atrstm,
const struct atrstm_fetch_radcoefs_args* args,
- double radcoefs[ATRSTM_RADCOEFS_COUNT__]);/*In m^-1*/
+ atrstm_radcoefs_T radcoefs); /* In m^-1 */
ATRSTM_API res_T
atrstm_fetch_radcoefs_svx
(const struct atrstm* atrstm,
const struct atrstm_fetch_radcoefs_svx_args* args,
- double radcoefs[ATRSTM_RADCOEFS_COUNT__][ATRSTM_SVX_OPS_COUNT__]);/*In m^-1*/
+ atrstm_radcoefs_svx_T radcoefs); /* In m^-1 */
ATRSTM_API res_T
atrstm_fetch_radcoefs_svx_voxel
(const struct atrstm* atrstm,
const struct atrstm_fetch_radcoefs_svx_voxel_args* args,
- double radcoefs[ATRSTM_RADCOEFS_COUNT__][ATRSTM_SVX_OPS_COUNT__]);/*In m^-1*/
+ atrstm_radcoefs_svx_T radcoefs); /* In m^-1 */
ATRSTM_API res_T
atrstm_dump_svx_octree
diff --git a/src/atrstm_dump_svx_octree.c b/src/atrstm_dump_svx_octree.c
@@ -0,0 +1,238 @@
+/* Copyright (C) 2020 CNRS
+ *
+ * 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 "atrstm.h"
+#include "atrstm_c.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 atrstm* atrstm;
+};
+
+/*******************************************************************************
+ * Helper functions
+ ******************************************************************************/
+static int
+check_dump_octree_args
+ (const struct atrstm* atrstm,
+ const struct atrstm_dump_svx_octree_args* args)
+{
+ ASSERT(atrstm && args);
+ (void)atrstm;
+ return args
+ && args->iband == ATRSTM_DUMP_SVX_OCTREE_ARGS_DEFAULT.iband /* Not use yet */
+ && args->iquad == ATRSTM_DUMP_SVX_OCTREE_ARGS_DEFAULT.iquad; /* Not use yet */
+}
+
+static void
+octree_data_init
+ (const struct atrstm* atrstm,
+ struct octree_data* data)
+{
+ ASSERT(data);
+ htable_vertex_init(atrstm->allocator, &data->vertex2id);
+ darray_double_init(atrstm->allocator, &data->vertices);
+ darray_double_init(atrstm->allocator, &data->data);
+ darray_size_t_init(atrstm->allocator, &data->cells);
+ data->atrstm = atrstm;
+}
+
+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 atrstm_fetch_radcoefs_svx_voxel_args args =
+ ATRSTM_FETCH_RADCOEFS_SVX_VOXEL_ARGS_DEFAULT;
+ atrstm_radcoefs_svx_T radcoefs;
+ struct octree_data* ctx = context;
+ 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));
+ }
+
+ /* Fetch leaf data */
+ args.voxel = *leaf;
+ args.components_mask = ATRSTM_CPNTS_MASK_ALL;
+ args.radcoefs_mask = ATRSTM_RADCOEF_FLAG_Kext;
+ args.operations_mask = ATRSTM_SVX_OPS_MASK_ALL;
+ ATRSTM(fetch_radcoefs_svx_voxel(ctx->atrstm, &args, radcoefs));
+
+ /* Register leaf data */
+ kext_min = radcoefs[ATRSTM_RADCOEF_Kext][ATRSTM_SVX_OP_MIN];
+ kext_max = radcoefs[ATRSTM_RADCOEF_Kext][ATRSTM_SVX_OP_MAX];
+ if(kext_min > kext_max) {/* Empty voxel */
+ kext_min = kext_max = 0;
+ }
+
+ CHK(RES_OK == darray_double_push_back(&ctx->data, &kext_min));
+ CHK(RES_OK == darray_double_push_back(&ctx->data, &kext_max));
+}
+
+/*******************************************************************************
+ * Exported functions
+ ******************************************************************************/
+res_T
+atrstm_dump_svx_octree
+ (const struct atrstm* atrstm,
+ const struct atrstm_dump_svx_octree_args* args,
+ FILE* stream)
+{
+ struct octree_data data;
+ const double* leaf_data;
+ size_t nvertices;
+ size_t ncells;
+ size_t i;
+ res_T res = RES_OK;
+
+ if(!atrstm || !check_dump_octree_args(atrstm, args)) {
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ octree_data_init(atrstm, &data);
+
+ /* Register leaf data */
+ SVX(tree_for_each_leaf(atrstm->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(atrstm->octree, &octree_desc));
+ ASSERT(ncells == octree_desc.nleaves);
+ }
+#endif
+
+ /* Write headers */
+ fprintf(stream, "# vtk DataFile Version 2.0\n");
+ fprintf(stream, "Radiative coefficients in [%g, %g]\n",
+ atrstm->wlen_range[0], atrstm->wlen_range[1]);
+ 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);
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
