rnatm_write_vtk.c (9242B)
1 /* Copyright (C) 2022, 2023, 2025 Centre National de la Recherche Scientifique 2 * Copyright (C) 2022, 2023, 2025 Institut Pierre-Simon Laplace 3 * Copyright (C) 2022, 2023, 2025 Institut de Physique du Globe de Paris 4 * Copyright (C) 2022, 2023, 2025 |Méso|Star> (contact@meso-star.com) 5 * Copyright (C) 2022, 2023, 2025 Observatoire de Paris 6 * Copyright (C) 2022, 2023, 2025 Université de Reims Champagne-Ardenne 7 * Copyright (C) 2022, 2023, 2025 Université de Versaille Saint-Quentin 8 * Copyright (C) 2022, 2023, 2025 Université Paul Sabatier 9 * 10 * This program is free software: you can redistribute it and/or modify 11 * it under the terms of the GNU General Public License as published by 12 * the Free Software Foundation, either version 3 of the License, or 13 * (at your option) any later version. 14 * 15 * This program is distributed in the hope that it will be useful, 16 * but WITHOUT ANY WARRANTY; without even the implied warranty of 17 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 18 * GNU General Public License for more details. 19 * 20 * You should have received a copy of the GNU General Public License 21 * along with this program. If not, see <http://www.gnu.org/licenses/>. */ 22 23 #include "rnatm.h" 24 #include "rnatm_c.h" 25 #include "rnatm_log.h" 26 #include "rnatm_voxel.h" 27 28 #include <rsys/clock_time.h> 29 #include <rsys/dynamic_array_double.h> 30 #include <rsys/dynamic_array_size_t.h> 31 #include <rsys/hash_table.h> 32 33 #include <star/svx.h> 34 35 struct vertex { 36 double x; 37 double y; 38 double z; 39 }; 40 41 static char 42 vertex_eq(const struct vertex* v0, const struct vertex* v1) 43 { 44 return eq_eps(v0->x, v1->x, 1.e-6) 45 && eq_eps(v0->y, v1->y, 1.e-6) 46 && eq_eps(v0->z, v1->z, 1.e-6); 47 } 48 49 /* Generate the htable_vertex data structure */ 50 #define HTABLE_NAME vertex 51 #define HTABLE_KEY struct vertex 52 #define HTABLE_DATA size_t 53 #define HTABLE_KEY_FUNCTOR_EQ vertex_eq 54 #include <rsys/hash_table.h> 55 56 /* Temporary data structure used to dump the octree data into a VTK file */ 57 struct octree_data { 58 struct htable_vertex vertex2id; /* Map a coordinate to its vertex id */ 59 struct darray_double vertices; /* Array of unique vertices */ 60 struct darray_double data; /* List of registered leaf data */ 61 struct darray_size_t cells; /* Ids of the cell vertices */ 62 const struct rnatm* atm; 63 }; 64 65 /******************************************************************************* 66 * Helper functions 67 ******************************************************************************/ 68 static void 69 octree_data_init 70 (const struct rnatm* atm, 71 struct octree_data* data) 72 { 73 ASSERT(data); 74 htable_vertex_init(atm->allocator, &data->vertex2id); 75 darray_double_init(atm->allocator, &data->vertices); 76 darray_double_init(atm->allocator, &data->data); 77 darray_size_t_init(atm->allocator, &data->cells); 78 data->atm = atm; 79 } 80 81 static void 82 octree_data_release(struct octree_data* data) 83 { 84 ASSERT(data); 85 htable_vertex_release(&data->vertex2id); 86 darray_double_release(&data->vertices); 87 darray_double_release(&data->data); 88 darray_size_t_release(&data->cells); 89 } 90 91 static void 92 octree_data_clear(struct octree_data* data) 93 { 94 ASSERT(data); 95 htable_vertex_clear(&data->vertex2id); 96 darray_double_clear(&data->vertices); 97 darray_double_clear(&data->data); 98 darray_size_t_clear(&data->cells); 99 } 100 101 static void 102 register_leaf 103 (const struct svx_voxel* leaf, 104 const size_t ileaf, 105 void* context) 106 { 107 struct octree_data* ctx = context; 108 struct vertex v[8]; 109 double kext_min; 110 double kext_max; 111 int i; 112 ASSERT(leaf && ctx); 113 (void)ileaf; 114 115 /* Compute the leaf vertices */ 116 v[0].x = leaf->lower[0]; v[0].y = leaf->lower[1]; v[0].z = leaf->lower[2]; 117 v[1].x = leaf->upper[0]; v[1].y = leaf->lower[1]; v[1].z = leaf->lower[2]; 118 v[2].x = leaf->lower[0]; v[2].y = leaf->upper[1]; v[2].z = leaf->lower[2]; 119 v[3].x = leaf->upper[0]; v[3].y = leaf->upper[1]; v[3].z = leaf->lower[2]; 120 v[4].x = leaf->lower[0]; v[4].y = leaf->lower[1]; v[4].z = leaf->upper[2]; 121 v[5].x = leaf->upper[0]; v[5].y = leaf->lower[1]; v[5].z = leaf->upper[2]; 122 v[6].x = leaf->lower[0]; v[6].y = leaf->upper[1]; v[6].z = leaf->upper[2]; 123 v[7].x = leaf->upper[0]; v[7].y = leaf->upper[1]; v[7].z = leaf->upper[2]; 124 125 FOR_EACH(i, 0, 8) { 126 size_t *pid = htable_vertex_find(&ctx->vertex2id, v+i); 127 size_t id; 128 if(pid) { 129 id = *pid; 130 } else { /* Register the leaf vertex */ 131 id = darray_double_size_get(&ctx->vertices)/3; 132 CHK(RES_OK == htable_vertex_set(&ctx->vertex2id, v+i, &id)); 133 CHK(RES_OK == darray_double_push_back(&ctx->vertices, &v[i].x)); 134 CHK(RES_OK == darray_double_push_back(&ctx->vertices, &v[i].y)); 135 CHK(RES_OK == darray_double_push_back(&ctx->vertices, &v[i].z)); 136 } 137 /* Add the vertex id to the leaf cell */ 138 CHK(RES_OK == darray_size_t_push_back(&ctx->cells, &id)); 139 } 140 141 /* Register leaf data */ 142 kext_min = rnatm_get_k_svx_voxel(ctx->atm, leaf, RNATM_RADCOEF_Kext, RNATM_SVX_OP_MIN); 143 kext_max = rnatm_get_k_svx_voxel(ctx->atm, leaf, RNATM_RADCOEF_Kext, RNATM_SVX_OP_MAX); 144 ASSERT(kext_min <= kext_max); 145 CHK(RES_OK == darray_double_push_back(&ctx->data, &kext_min)); 146 CHK(RES_OK == darray_double_push_back(&ctx->data, &kext_max)); 147 } 148 149 static res_T 150 write_vtk_octree 151 (struct rnatm* atm, 152 const size_t iaccel_struct, 153 struct octree_data* data, 154 FILE* stream) 155 { 156 const struct accel_struct* accel_struct = NULL; 157 const double* leaf_data = NULL; 158 size_t nvertices; 159 size_t ncells; 160 size_t i; 161 res_T res = RES_OK; 162 163 ASSERT(atm && stream); 164 ASSERT(iaccel_struct < darray_accel_struct_size_get(&atm->accel_structs)); 165 166 res = make_sure_octree_is_loaded(atm, iaccel_struct); 167 if(res != RES_OK) goto error; 168 169 accel_struct = darray_accel_struct_cdata_get(&atm->accel_structs)+iaccel_struct; 170 171 octree_data_clear(data); 172 173 /* Register leaf data */ 174 SVX(tree_for_each_leaf(accel_struct->octree, register_leaf, data)); 175 nvertices = darray_double_size_get(&data->vertices)/3/*#coords per vertex*/; 176 ncells = darray_size_t_size_get(&data->cells)/8/*#ids per cell*/; 177 #ifndef NDEBUG 178 { 179 struct svx_tree_desc octree_desc = SVX_TREE_DESC_NULL; 180 SVX(tree_get_desc(accel_struct->octree, &octree_desc)); 181 ASSERT(ncells == octree_desc.nleaves); 182 } 183 #endif 184 185 /* Write headers */ 186 fprintf(stream, "# vtk DataFile Version 2.0\n"); 187 fprintf(stream, 188 "Radiative coefficients for band %lu and quadrature point %lu\n", 189 (unsigned long)accel_struct->iband, 190 (unsigned long)accel_struct->iquad_pt); 191 fprintf(stream, "ASCII\n"); 192 fprintf(stream, "DATASET UNSTRUCTURED_GRID\n"); 193 194 /* Write vertex coordinates */ 195 fprintf(stream, "POINTS %lu float\n", (unsigned long)nvertices); 196 FOR_EACH(i, 0, nvertices) { 197 fprintf(stream, "%g %g %g\n", 198 SPLIT3(darray_double_cdata_get(&data->vertices) + i*3)); 199 } 200 201 /* Write the cells */ 202 fprintf(stream, "CELLS %lu %lu\n", 203 (unsigned long)ncells, 204 (unsigned long)(ncells*(8/*#verts per cell*/ + 1/*1st field of a cell*/))); 205 FOR_EACH(i, 0, ncells) { 206 fprintf(stream, "8 %lu %lu %lu %lu %lu %lu %lu %lu\n", 207 (unsigned long)darray_size_t_cdata_get(&data->cells)[i*8+0], 208 (unsigned long)darray_size_t_cdata_get(&data->cells)[i*8+1], 209 (unsigned long)darray_size_t_cdata_get(&data->cells)[i*8+2], 210 (unsigned long)darray_size_t_cdata_get(&data->cells)[i*8+3], 211 (unsigned long)darray_size_t_cdata_get(&data->cells)[i*8+4], 212 (unsigned long)darray_size_t_cdata_get(&data->cells)[i*8+5], 213 (unsigned long)darray_size_t_cdata_get(&data->cells)[i*8+6], 214 (unsigned long)darray_size_t_cdata_get(&data->cells)[i*8+7]); 215 } 216 217 /* Write the cell type */ 218 fprintf(stream, "CELL_TYPES %lu\n", (unsigned long)ncells); 219 FOR_EACH(i, 0, ncells) fprintf(stream, "11\n"); 220 221 /* Write the cell data */ 222 leaf_data = darray_double_cdata_get(&data->data); 223 fprintf(stream, "CELL_DATA %lu\n", (unsigned long)ncells); 224 fprintf(stream, "SCALARS Kext float 2\n"); 225 fprintf(stream, "LOOKUP_TABLE default\n"); 226 FOR_EACH(i, 0, ncells) { 227 fprintf(stream, "%g %g\n", 228 (float)leaf_data[i*2+0], 229 (float)leaf_data[i*2+1]); 230 } 231 232 exit: 233 return res; 234 error: 235 goto exit; 236 } 237 238 /******************************************************************************* 239 * Exported function 240 ******************************************************************************/ 241 res_T 242 rnatm_write_vtk_octrees 243 (struct rnatm* atm, 244 const size_t items[2], /* Range of octrees to write */ 245 FILE* stream) 246 { 247 char buf[128]; 248 struct time t0, t1; 249 struct octree_data data; 250 size_t i; 251 res_T res = RES_OK; 252 253 if(!atm 254 || !stream 255 || items[0] >= darray_accel_struct_size_get(&atm->accel_structs) 256 || items[1] >= darray_accel_struct_size_get(&atm->accel_structs) 257 || items[0] > items[1]) { 258 res = RES_BAD_ARG; 259 goto error; 260 } 261 262 /* Start time recording */ 263 log_info(atm, "write VTK octrees\n"); 264 time_current(&t0); 265 266 octree_data_init(atm, &data); 267 FOR_EACH(i, items[0], items[1]+1) { 268 res = write_vtk_octree(atm, items[0]+i, &data, stream); 269 if(res != RES_OK) goto error; 270 } 271 octree_data_release(&data); 272 273 /* Log elapsed time */ 274 time_sub(&t0, time_current(&t1), &t0); 275 time_dump(&t0, TIME_ALL, NULL, buf, sizeof(buf)); 276 log_info(atm, "octrees written in %s\n", buf); 277 278 exit: 279 return res; 280 error: 281 goto exit; 282 }