atrstm_dump_svx_octree.c (7966B)
1 /* Copyright (C) 2022, 2023 |Méso|Star> (contact@meso-star.com) 2 * Copyright (C) 2020, 2021 Centre National de la Recherche Scientifique 3 * 4 * This program is free software: you can redistribute it and/or modify 5 * it under the terms of the GNU General Public License as published by 6 * the Free Software Foundation, either version 3 of the License, or 7 * (at your option) any later version. 8 * 9 * This program is distributed in the hope that it will be useful, 10 * but WITHOUT ANY WARRANTY; without even the implied warranty of 11 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 12 * GNU General Public License for more details. 13 * 14 * You should have received a copy of the GNU General Public License 15 * along with this program. If not, see <http://www.gnu.org/licenses/>. */ 16 17 #include "atrstm.h" 18 #include "atrstm_c.h" 19 20 #include <rsys/dynamic_array_double.h> 21 #include <rsys/dynamic_array_size_t.h> 22 #include <rsys/hash_table.h> 23 24 #include <star/svx.h> 25 26 struct vertex { 27 double x; 28 double y; 29 double z; 30 }; 31 32 static char 33 vertex_eq(const struct vertex* v0, const struct vertex* v1) 34 { 35 return eq_eps(v0->x, v1->x, 1.e-6) 36 && eq_eps(v0->y, v1->y, 1.e-6) 37 && eq_eps(v0->z, v1->z, 1.e-6); 38 } 39 40 /* Generate the htable_vertex data structure */ 41 #define HTABLE_NAME vertex 42 #define HTABLE_KEY struct vertex 43 #define HTABLE_DATA size_t 44 #define HTABLE_KEY_FUNCTOR_EQ vertex_eq 45 #include <rsys/hash_table.h> 46 47 /* Temporary data structure used to dump the octree data into a VTK file */ 48 struct octree_data { 49 struct htable_vertex vertex2id; /* Map a coordinate to its vertex id */ 50 struct darray_double vertices; /* Array of unique vertices */ 51 struct darray_double data; /* List of registered leaf data */ 52 struct darray_size_t cells; /* Ids of the cell vertices */ 53 const struct atrstm* atrstm; 54 }; 55 56 /******************************************************************************* 57 * Helper functions 58 ******************************************************************************/ 59 static int 60 check_dump_octree_args 61 (const struct atrstm* atrstm, 62 const struct atrstm_dump_svx_octree_args* args) 63 { 64 ASSERT(atrstm); 65 (void)atrstm; 66 return args 67 && args->iband == ATRSTM_DUMP_SVX_OCTREE_ARGS_DEFAULT.iband /* Not use yet */ 68 && args->iquad == ATRSTM_DUMP_SVX_OCTREE_ARGS_DEFAULT.iquad; /* Not use yet */ 69 } 70 71 static void 72 octree_data_init 73 (const struct atrstm* atrstm, 74 struct octree_data* data) 75 { 76 ASSERT(data); 77 htable_vertex_init(atrstm->allocator, &data->vertex2id); 78 darray_double_init(atrstm->allocator, &data->vertices); 79 darray_double_init(atrstm->allocator, &data->data); 80 darray_size_t_init(atrstm->allocator, &data->cells); 81 data->atrstm = atrstm; 82 } 83 84 static void 85 octree_data_release(struct octree_data* data) 86 { 87 ASSERT(data); 88 htable_vertex_release(&data->vertex2id); 89 darray_double_release(&data->vertices); 90 darray_double_release(&data->data); 91 darray_size_t_release(&data->cells); 92 } 93 94 static void 95 register_leaf 96 (const struct svx_voxel* leaf, 97 const size_t ileaf, 98 void* context) 99 { 100 struct atrstm_fetch_radcoefs_svx_voxel_args args = 101 ATRSTM_FETCH_RADCOEFS_SVX_VOXEL_ARGS_DEFAULT; 102 atrstm_radcoefs_svx_T radcoefs; 103 struct octree_data* ctx = context; 104 struct vertex v[8]; 105 double kext_min; 106 double kext_max; 107 int i; 108 ASSERT(leaf && ctx); 109 (void)ileaf; 110 111 /* Compute the leaf vertices */ 112 v[0].x = leaf->lower[0]; v[0].y = leaf->lower[1]; v[0].z = leaf->lower[2]; 113 v[1].x = leaf->upper[0]; v[1].y = leaf->lower[1]; v[1].z = leaf->lower[2]; 114 v[2].x = leaf->lower[0]; v[2].y = leaf->upper[1]; v[2].z = leaf->lower[2]; 115 v[3].x = leaf->upper[0]; v[3].y = leaf->upper[1]; v[3].z = leaf->lower[2]; 116 v[4].x = leaf->lower[0]; v[4].y = leaf->lower[1]; v[4].z = leaf->upper[2]; 117 v[5].x = leaf->upper[0]; v[5].y = leaf->lower[1]; v[5].z = leaf->upper[2]; 118 v[6].x = leaf->lower[0]; v[6].y = leaf->upper[1]; v[6].z = leaf->upper[2]; 119 v[7].x = leaf->upper[0]; v[7].y = leaf->upper[1]; v[7].z = leaf->upper[2]; 120 121 FOR_EACH(i, 0, 8) { 122 size_t *pid = htable_vertex_find(&ctx->vertex2id, v+i); 123 size_t id; 124 if(pid) { 125 id = *pid; 126 } else { /* Register the leaf vertex */ 127 id = darray_double_size_get(&ctx->vertices)/3; 128 CHK(RES_OK == htable_vertex_set(&ctx->vertex2id, v+i, &id)); 129 CHK(RES_OK == darray_double_push_back(&ctx->vertices, &v[i].x)); 130 CHK(RES_OK == darray_double_push_back(&ctx->vertices, &v[i].y)); 131 CHK(RES_OK == darray_double_push_back(&ctx->vertices, &v[i].z)); 132 } 133 /* Add the vertex id to the leaf cell */ 134 CHK(RES_OK == darray_size_t_push_back(&ctx->cells, &id)); 135 } 136 137 /* Fetch leaf data */ 138 args.voxel = *leaf; 139 args.components_mask = ATRSTM_CPNTS_MASK_ALL; 140 args.radcoefs_mask = ATRSTM_RADCOEF_FLAG_Kext; 141 args.operations_mask = ATRSTM_SVX_OPS_MASK_ALL; 142 ATRSTM(fetch_radcoefs_svx_voxel(ctx->atrstm, &args, radcoefs)); 143 144 /* Register leaf data */ 145 kext_min = radcoefs[ATRSTM_RADCOEF_Kext][ATRSTM_SVX_OP_MIN]; 146 kext_max = radcoefs[ATRSTM_RADCOEF_Kext][ATRSTM_SVX_OP_MAX]; 147 CHK(RES_OK == darray_double_push_back(&ctx->data, &kext_min)); 148 CHK(RES_OK == darray_double_push_back(&ctx->data, &kext_max)); 149 } 150 151 /******************************************************************************* 152 * Exported functions 153 ******************************************************************************/ 154 res_T 155 atrstm_dump_svx_octree 156 (const struct atrstm* atrstm, 157 const struct atrstm_dump_svx_octree_args* args, 158 FILE* stream) 159 { 160 struct octree_data data; 161 const double* leaf_data; 162 size_t nvertices; 163 size_t ncells; 164 size_t i; 165 res_T res = RES_OK; 166 167 if(!atrstm || !check_dump_octree_args(atrstm, args)) { 168 res = RES_BAD_ARG; 169 goto error; 170 } 171 172 octree_data_init(atrstm, &data); 173 174 /* Register leaf data */ 175 SVX(tree_for_each_leaf(atrstm->octree, register_leaf, &data)); 176 nvertices = darray_double_size_get(&data.vertices)/3/*#coords per vertex*/; 177 ncells = darray_size_t_size_get(&data.cells)/8/*#ids per cell*/; 178 #ifndef NDEBUG 179 { 180 struct svx_tree_desc octree_desc = SVX_TREE_DESC_NULL; 181 SVX(tree_get_desc(atrstm->octree, &octree_desc)); 182 ASSERT(ncells == octree_desc.nleaves); 183 } 184 #endif 185 186 /* Write headers */ 187 fprintf(stream, "# vtk DataFile Version 2.0\n"); 188 fprintf(stream, "Radiative coefficients in [%g, %g]\n", 189 atrstm->wlen_range[0], atrstm->wlen_range[1]); 190 fprintf(stream, "ASCII\n"); 191 fprintf(stream, "DATASET UNSTRUCTURED_GRID\n"); 192 193 /* Write vertex coordinates */ 194 fprintf(stream, "POINTS %lu float\n", (unsigned long)nvertices); 195 FOR_EACH(i, 0, nvertices) { 196 fprintf(stream, "%g %g %g\n", 197 SPLIT3(darray_double_cdata_get(&data.vertices) + i*3)); 198 } 199 200 /* Write the cells */ 201 fprintf(stream, "CELLS %lu %lu\n", 202 (unsigned long)ncells, 203 (unsigned long)(ncells*(8/*#verts per cell*/ + 1/*1st field of a cell*/))); 204 FOR_EACH(i, 0, ncells) { 205 fprintf(stream, "8 %lu %lu %lu %lu %lu %lu %lu %lu\n", 206 (unsigned long)darray_size_t_cdata_get(&data.cells)[i*8+0], 207 (unsigned long)darray_size_t_cdata_get(&data.cells)[i*8+1], 208 (unsigned long)darray_size_t_cdata_get(&data.cells)[i*8+2], 209 (unsigned long)darray_size_t_cdata_get(&data.cells)[i*8+3], 210 (unsigned long)darray_size_t_cdata_get(&data.cells)[i*8+4], 211 (unsigned long)darray_size_t_cdata_get(&data.cells)[i*8+5], 212 (unsigned long)darray_size_t_cdata_get(&data.cells)[i*8+6], 213 (unsigned long)darray_size_t_cdata_get(&data.cells)[i*8+7]); 214 } 215 216 /* Write the cell type */ 217 fprintf(stream, "CELL_TYPES %lu\n", (unsigned long)ncells); 218 FOR_EACH(i, 0, ncells) fprintf(stream, "11\n"); 219 220 /* Write the cell data */ 221 leaf_data = darray_double_cdata_get(&data.data); 222 fprintf(stream, "CELL_DATA %lu\n", (unsigned long)ncells); 223 fprintf(stream, "SCALARS Kext double 2\n"); 224 fprintf(stream, "LOOKUP_TABLE default\n"); 225 FOR_EACH(i, 0, ncells) { 226 fprintf(stream, "%g %g\n", leaf_data[i*2+0], leaf_data[i*2+1]); 227 } 228 octree_data_release(&data); 229 230 exit: 231 return res; 232 error: 233 goto exit; 234 } 235