commit a9871fee9cafad8dc0b1c07163af56dc5b7d8625
parent 5d798a24629b997ca0c1087227f13c112b28c25a
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Wed, 31 Aug 2022 14:17:32 +0200
Small refactoring of rnatm_write_vtk_octree function
Diffstat:
4 files changed, 282 insertions(+), 280 deletions(-)
diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt
@@ -68,7 +68,7 @@ set(RNATM_FILES_SRC
rnatm_octree.c
rnatm_properties.c
rnatm_voxel_partition.c
- rnatm_write_vtk_octree.c)
+ rnatm_write_vtk.c)
set(RNATM_FILES_INC
rnatm_c.h
rnatm_log.h
diff --git a/src/rnatm_c.h b/src/rnatm_c.h
@@ -144,7 +144,7 @@ aerosol_copy_and_release
******************************************************************************/
struct accel_struct {
struct svx_tree* octree;
- fpos_t fpos; /* Position of the serialized data in the storage file */
+ fpos_t fpos; /* Position of the serialized octree data in the storage file */
size_t iband; /* Index of the spectral band */
size_t iquad_pt; /* Index of the quadrature point */
};
diff --git a/src/rnatm_write_vtk.c b/src/rnatm_write_vtk.c
@@ -0,0 +1,280 @@
+/* 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
+octree_data_clear(struct octree_data* data)
+{
+ ASSERT(data);
+ htable_vertex_clear(&data->vertex2id);
+ darray_double_clear(&data->vertices);
+ darray_double_clear(&data->data);
+ darray_size_t_clear(&data->cells);
+}
+
+static void
+register_leaf
+ (const struct svx_voxel* leaf,
+ const size_t ileaf,
+ void* context)
+{
+ 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));
+ }
+
+ /* Register leaf data */
+ kext_min = rnatm_get_k_svx_voxel(ctx->atm, leaf, RNATM_RADCOEF_Kext, RNATM_SVX_OP_MIN);
+ kext_max = rnatm_get_k_svx_voxel(ctx->atm, leaf, RNATM_RADCOEF_Kext, RNATM_SVX_OP_MAX);
+ ASSERT(kext_min <= kext_max);
+ CHK(RES_OK == darray_double_push_back(&ctx->data, &kext_min));
+ CHK(RES_OK == darray_double_push_back(&ctx->data, &kext_max));
+}
+
+static res_T
+write_vtk_octree
+ (struct rnatm* atm,
+ const size_t iaccel_struct,
+ struct octree_data* data,
+ FILE* stream)
+{
+ const struct accel_struct* accel_struct = NULL;
+ const double* leaf_data = NULL;
+ size_t nvertices;
+ size_t ncells;
+ size_t i;
+ res_T res = RES_OK;
+
+ ASSERT(atm && stream);
+ ASSERT(iaccel_struct < darray_accel_struct_size_get(&atm->accel_structs));
+
+ res = make_sure_octree_is_loaded(atm, iaccel_struct);
+ if(res != RES_OK) goto error;
+
+ accel_struct = darray_accel_struct_cdata_get(&atm->accel_structs)+iaccel_struct;
+
+ octree_data_clear(data);
+
+ /* Register leaf data */
+ SVX(tree_for_each_leaf(accel_struct->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(accel_struct->octree, &octree_desc));
+ ASSERT(ncells == octree_desc.nleaves);
+ }
+#endif
+
+ /* Write headers */
+ fprintf(stream, "# vtk DataFile Version 2.0\n");
+ fprintf(stream,
+ "Radiative coefficients for band %lu and quadrature point %lu\n",
+ (unsigned long)accel_struct->iband,
+ (unsigned long)accel_struct->iquad_pt);
+ 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 float 2\n");
+ fprintf(stream, "LOOKUP_TABLE default\n");
+ FOR_EACH(i, 0, ncells) {
+ fprintf(stream, "%g %g\n",
+ (float)leaf_data[i*2+0],
+ (float)leaf_data[i*2+1]);
+ }
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
+/*******************************************************************************
+ * Exported function
+ ******************************************************************************/
+res_T
+rnatm_write_vtk_octrees
+ (struct rnatm* atm,
+ const size_t items[2], /* Range of octrees to write */
+ FILE* stream)
+{
+ char buf[128];
+ struct time t0, t1;
+ struct octree_data data;
+ size_t i;
+ res_T res = RES_OK;
+
+ if(!atm
+ || !stream
+ || items[0] >= darray_accel_struct_size_get(&atm->accel_structs)
+ || items[1] >= darray_accel_struct_size_get(&atm->accel_structs)
+ || items[0] > items[1]) {
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ /* Start time recording */
+ log_info(atm, "write VTK octrees\n");
+ time_current(&t0);
+
+ octree_data_init(atm, &data);
+ FOR_EACH(i, items[0], items[1]+1) {
+ res = write_vtk_octree(atm, items[0]+i, &data, stream);
+ if(res != RES_OK) goto error;
+ }
+ 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, "octrees written in %s\n", buf);
+
+exit:
+ return res;
+error:
+ goto exit;
+}
diff --git a/src/rnatm_write_vtk_octree.c b/src/rnatm_write_vtk_octree.c
@@ -1,278 +0,0 @@
-/* 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
-octree_data_clear(struct octree_data* data)
-{
- ASSERT(data);
- htable_vertex_clear(&data->vertex2id);
- darray_double_clear(&data->vertices);
- darray_double_clear(&data->data);
- darray_size_t_clear(&data->cells);
-}
-
-static void
-register_leaf
- (const struct svx_voxel* leaf,
- const size_t ileaf,
- void* context)
-{
- 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));
- }
-
- /* Register leaf data */
- kext_min = rnatm_get_k_svx_voxel(ctx->atm, leaf, RNATM_RADCOEF_Kext, RNATM_SVX_OP_MIN);
- kext_max = rnatm_get_k_svx_voxel(ctx->atm, leaf, RNATM_RADCOEF_Kext, RNATM_SVX_OP_MAX);
- ASSERT(kext_min <= kext_max);
- CHK(RES_OK == darray_double_push_back(&ctx->data, &kext_min));
- CHK(RES_OK == darray_double_push_back(&ctx->data, &kext_max));
-}
-
-static res_T
-write_vtk_octree
- (struct rnatm* atm,
- const size_t iaccel_struct,
- struct octree_data* data,
- FILE* stream)
-{
- const struct accel_struct* accel_struct = NULL;
- const double* leaf_data = NULL;
- size_t nvertices;
- size_t ncells;
- size_t i;
- res_T res = RES_OK;
-
- ASSERT(atm && stream);
- ASSERT(iaccel_struct < darray_accel_struct_size_get(&atm->accel_structs));
-
- res = make_sure_octree_is_loaded(atm, iaccel_struct);
- if(res != RES_OK) goto error;
-
- accel_struct = darray_accel_struct_cdata_get(&atm->accel_structs)+iaccel_struct;
-
- octree_data_clear(data);
-
- /* Register leaf data */
- SVX(tree_for_each_leaf(accel_struct->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(accel_struct->octree, &octree_desc));
- ASSERT(ncells == octree_desc.nleaves);
- }
-#endif
-
- /* Write headers */
- fprintf(stream, "# vtk DataFile Version 2.0\n");
- fprintf(stream,
- "Radiative coefficients for band %lu and quadrature point %lu\n",
- (unsigned long)accel_struct->iband,
- (unsigned long)accel_struct->iquad_pt);
- 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]);
- }
-
-exit:
- return res;
-error:
- goto exit;
-}
-
-/*******************************************************************************
- * Exported function
- ******************************************************************************/
-res_T
-rnatm_write_vtk_octrees
- (struct rnatm* atm,
- const size_t structs[2], /* Range of acceleration structures to write */
- FILE* stream)
-{
- char buf[128];
- struct time t0, t1;
- struct octree_data data;
- size_t i;
- res_T res = RES_OK;
-
- if(!atm
- || !stream
- || structs[0] >= darray_accel_struct_size_get(&atm->accel_structs)
- || structs[1] >= darray_accel_struct_size_get(&atm->accel_structs)
- || structs[0] > structs[1]) {
- res = RES_BAD_ARG;
- goto error;
- }
-
- /* Start time recording */
- log_info(atm, "write VTK octrees\n");
- time_current(&t0);
-
- octree_data_init(atm, &data);
- FOR_EACH(i, structs[0], structs[1]+1) {
- res = write_vtk_octree(atm, structs[0]+i, &data, stream);
- if(res != RES_OK) goto error;
- }
- 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, "octrees written in %s\n", buf);
-
-exit:
- return res;
-error:
- goto exit;
-}
