stardis-output.c (76782B)
1 /* Copyright (C) 2018-2025 |Méso|Star> (contact@meso-star.com) 2 * 3 * This program is free software: you can redistribute it and/or modify 4 * it under the terms of the GNU General Public License as published by 5 * the Free Software Foundation, either version 3 of the License, or 6 * (at your option) any later version. 7 * 8 * This program is distributed in the hope that it will be useful, 9 * but WITHOUT ANY WARRANTY; without even the implied warranty of 10 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 11 * GNU General Public License for more details. 12 * 13 * You should have received a copy of the GNU General Public License 14 * along with this program. If not, see <http://www.gnu.org/licenses/>. */ 15 16 #define _POSIX_C_SOURCE 200112L /* snprintf */ 17 #include "stardis-output.h" 18 #include "stardis-args.h" 19 #include "stardis-compute.h" 20 #include "stardis-fluid.h" 21 #include "stardis-solid.h" 22 #include "stardis-hbound.h" 23 #include "stardis-tbound.h" 24 #include "stardis-fbound.h" 25 #include "stardis-ssconnect.h" 26 #include "stardis-sfconnect.h" 27 #include "stardis-intface.h" 28 #include "stardis-app.h" 29 #include "stardis-green-types.h" 30 31 #include <sdis.h> 32 33 #include<star/ssp.h> 34 #include<star/senc3d.h> 35 #include<star/sg3d.h> 36 37 #include <rsys/math.h> 38 #include <rsys/mem_allocator.h> 39 #include <rsys/dynamic_array_uint.h> 40 #include <rsys/hash_table.h> 41 #include <rsys/logger.h> 42 #include <rsys/str.h> 43 #include <rsys/clock_time.h> 44 45 #include <limits.h> 46 #include <string.h> 47 #include <stdio.h> 48 49 #define HTABLE_NAME weigth 50 #define HTABLE_DATA double 51 #define HTABLE_KEY unsigned 52 #include <rsys/hash_table.h> 53 54 struct w_ctx { 55 struct mem_allocator* alloc; 56 const struct darray_descriptions* desc; 57 struct htable_weigth pw; 58 struct htable_weigth flux; 59 FILE* stream; 60 }; 61 62 struct e_ctx { 63 const struct darray_descriptions* desc; 64 FILE* stream; 65 }; 66 67 /******************************************************************************* 68 * Local Type; for documentation purpose 69 * These values are used in dumps 70 ******************************************************************************/ 71 enum enclosure_errors_t { 72 NO_ENCLOSURE_ERROR = BIT(0), 73 ENCLOSURE_WITH_N_MEDIA = BIT(1), 74 ENCLOSURE_WITH_UNDEF_MEDIUM = BIT(2) 75 }; 76 77 static res_T 78 copy_desc_to_green_desc 79 (struct green_description* gdesc, 80 const struct darray_descriptions* descriptions, 81 const size_t idx) 82 { 83 size_t sz; 84 const struct description* desc; 85 ASSERT(gdesc && descriptions); 86 sz = darray_descriptions_size_get(descriptions); 87 CHK(idx < sz); 88 desc = darray_descriptions_cdata_get(descriptions) + idx; 89 switch(desc->type) { 90 case DESC_MAT_SOLID: 91 gdesc->type = GREEN_MAT_SOLID; 92 strcpy(gdesc->d.solid.name, str_cget(&desc->d.solid->name)); 93 gdesc->d.solid.conductivity = desc->d.solid->lambda; 94 gdesc->d.solid.volumic_mass = desc->d.solid->rho; 95 gdesc->d.solid.calorific_capacity = desc->d.solid->cp; 96 gdesc->d.solid.volumic_power = desc->d.solid->vpower; 97 gdesc->d.solid.initial_temperature = desc->d.solid->tinit; 98 gdesc->d.solid.imposed_temperature = desc->d.solid->imposed_temperature; 99 break; 100 case DESC_MAT_FLUID: 101 gdesc->type = GREEN_MAT_FLUID; 102 strcpy(gdesc->d.fluid.name, str_cget(&desc->d.fluid->name)); 103 gdesc->d.fluid.volumic_mass = desc->d.fluid->rho; 104 gdesc->d.fluid.calorific_capacity = desc->d.fluid->cp; 105 gdesc->d.fluid.initial_temperature = desc->d.fluid->tinit; 106 gdesc->d.fluid.imposed_temperature = desc->d.fluid->imposed_temperature; 107 break; 108 case DESC_BOUND_H_FOR_FLUID: 109 case DESC_BOUND_H_FOR_SOLID: 110 gdesc->type = GREEN_BOUND_H; 111 strcpy(gdesc->d.h_boundary.name, str_cget(&desc->d.h_boundary->name)); 112 gdesc->d.h_boundary.reference_temperature 113 = desc->d.h_boundary->ref_temperature; 114 gdesc->d.h_boundary.emissivity = desc->d.h_boundary->emissivity; 115 gdesc->d.h_boundary.specular_fraction 116 = desc->d.h_boundary->specular_fraction; 117 gdesc->d.h_boundary.convection_coefficient = desc->d.h_boundary->hc; 118 gdesc->d.h_boundary.imposed_temperature 119 = desc->d.h_boundary->imposed_temperature; 120 break; 121 case DESC_BOUND_T_FOR_SOLID: 122 gdesc->type = GREEN_BOUND_T; 123 strcpy(gdesc->d.t_boundary.name, str_cget(&desc->d.t_boundary->name)); 124 gdesc->d.t_boundary.imposed_temperature 125 = desc->d.t_boundary->imposed_temperature; 126 break; 127 case DESC_BOUND_F_FOR_SOLID: 128 gdesc->type = GREEN_BOUND_F; 129 strcpy(gdesc->d.f_boundary.name, str_cget(&desc->d.f_boundary->name)); 130 gdesc->d.f_boundary.imposed_flux 131 = desc->d.f_boundary->imposed_flux; 132 break; 133 case DESC_SOLID_FLUID_CONNECT: 134 gdesc->type = GREEN_SOLID_FLUID_CONNECT; 135 strcpy(gdesc->d.sf_connect.name, str_cget(&desc->d.sf_connect->name)); 136 gdesc->d.sf_connect.reference_temperature 137 = desc->d.sf_connect->ref_temperature; 138 gdesc->d.sf_connect.emissivity = desc->d.sf_connect->emissivity; 139 gdesc->d.sf_connect.specular_fraction 140 = desc->d.sf_connect->specular_fraction; 141 gdesc->d.sf_connect.convection_coefficient = desc->d.sf_connect->hc; 142 break; 143 case DESC_SOLID_SOLID_CONNECT: 144 gdesc->type = GREEN_SOLID_SOLID_CONNECT; 145 strcpy(gdesc->d.ss_connect.name, str_cget(&desc->d.ss_connect->name)); 146 gdesc->d.ss_connect.thermal_contact_resistance = desc->d.ss_connect->tcr; 147 break; 148 default: return RES_BAD_ARG; 149 } 150 return RES_OK; 151 } 152 153 /******************************************************************************* 154 * Local Functions 155 ******************************************************************************/ 156 157 static res_T 158 merge_flux_terms 159 (struct sdis_interface* interf, 160 const enum sdis_side side, 161 const double flux_term, 162 void* ctx) 163 { 164 res_T res = RES_OK; 165 struct sdis_data* data = NULL; 166 struct intface* d__; 167 unsigned desc_id; 168 struct w_ctx* w_ctx = ctx; 169 const struct description* descs; 170 171 ASSERT(interf && w_ctx); 172 (void)side; 173 174 data = sdis_interface_get_data(interf); 175 d__ = sdis_data_get(data); 176 desc_id = d__->desc_id; 177 CHK(desc_id < darray_descriptions_size_get(w_ctx->desc)); 178 descs = darray_descriptions_cdata_get(w_ctx->desc); 179 180 switch (descs[desc_id].type) { 181 case DESC_BOUND_T_FOR_SOLID: 182 case DESC_BOUND_H_FOR_SOLID: 183 case DESC_BOUND_H_FOR_FLUID: 184 FATAL("Cannot have a flux term here.\n"); break; 185 case DESC_BOUND_F_FOR_SOLID: { 186 double* w; 187 w = htable_weigth_find(&w_ctx->flux, &desc_id); 188 if(w) *w += flux_term; 189 else ERR(htable_weigth_set(&w_ctx->flux, &desc_id, &flux_term)); 190 break; 191 } 192 default: FATAL("Unreachable code.\n"); break; 193 } 194 end: 195 return res; 196 error: 197 goto end; 198 } 199 200 static res_T 201 merge_power_terms 202 (struct sdis_medium* mdm, 203 const double power_term, 204 void* ctx) 205 { 206 res_T res = RES_OK; 207 struct sdis_data* data = NULL; 208 enum sdis_medium_type type; 209 struct w_ctx* w_ctx = ctx; 210 size_t sz; 211 212 ASSERT(mdm && w_ctx); 213 214 sz = darray_descriptions_size_get(w_ctx->desc); 215 data = sdis_medium_get_data(mdm); 216 type = sdis_medium_get_type(mdm); 217 218 switch (type) { 219 case SDIS_FLUID: { 220 /* Could be OK, but unimplemented in stardis */ 221 FATAL("Unexpected power term in fluid"); 222 } 223 case SDIS_SOLID: { 224 struct solid** psolid = sdis_data_get(data); 225 double* w; 226 unsigned id = (*psolid)->desc_id; 227 CHK(id < sz); 228 w = htable_weigth_find(&w_ctx->pw, &id); 229 if(w) *w += power_term; 230 else ERR(htable_weigth_set(&w_ctx->pw, &id, &power_term)); 231 break; 232 } 233 default: FATAL("Unreachable code.\n"); break; 234 } 235 end: 236 return res; 237 error: 238 goto end; 239 } 240 241 /******************************************************************************* 242 * Public Functions 243 ******************************************************************************/ 244 245 res_T 246 dump_path 247 (const struct sdis_heat_path* path, 248 void* context) 249 { 250 res_T res = RES_OK; 251 struct dump_path_context* dump_ctx = context; 252 FILE* stream = NULL; 253 char* name = NULL; 254 enum sdis_heat_path_flag status = SDIS_HEAT_PATH_NONE; 255 size_t vcount_, scount_, offset, name_sz, istrip; 256 unsigned long scount, vcount, strip_1, type_changes, *strip_type_changes = NULL; 257 258 ASSERT(path && dump_ctx 259 && dump_ctx->stardis 260 && !str_is_empty(&dump_ctx->stardis->paths_filename)); 261 262 ERR(sdis_heat_path_get_status(path, &status)); 263 264 name_sz = 20 + str_len(&dump_ctx->stardis->paths_filename); 265 name = MEM_CALLOC(dump_ctx->stardis->allocator, name_sz, sizeof(*name)); 266 if(!name) { 267 res = RES_MEM_ERR; 268 goto error; 269 } 270 snprintf(name, name_sz, "%s%08lu%s.vtk", 271 str_cget(&dump_ctx->stardis->paths_filename), dump_ctx->rank++, 272 (status == SDIS_HEAT_PATH_FAILURE ? "_err" : "")); 273 274 stream = fopen(name, "w"); 275 if(!stream) { 276 logger_print(dump_ctx->stardis->logger, LOG_ERROR, 277 "cannot open file '%s' for writing.\n", name); 278 res = RES_IO_ERR; 279 goto error; 280 } 281 282 /* Get counts */ 283 ERR(sdis_heat_path_get_line_strips_count(path, &scount_)); 284 if(scount_ > ULONG_MAX) goto abort; 285 scount = (unsigned long)scount_; 286 vcount_ = 0; 287 strip_1 = 0; 288 type_changes = 0; 289 strip_type_changes = MEM_CALLOC(dump_ctx->stardis->allocator, scount, 290 sizeof(*strip_type_changes)); 291 if(!strip_type_changes) { 292 res = RES_MEM_ERR; 293 goto error; 294 } 295 296 FOR_EACH(istrip, 0, scount_) { 297 size_t ivert, nverts; 298 enum sdis_heat_vertex_type prev_type = SDIS_HEAT_VERTEX_CONDUCTION; 299 ERR(sdis_heat_path_line_strip_get_vertices_count(path, istrip, &nverts)); 300 if(nverts == 0 || nverts > ULONG_MAX) goto abort; 301 if(nverts == 1) strip_1++; 302 FOR_EACH(ivert, 0, nverts) { 303 struct sdis_heat_vertex vtx; 304 ERR(sdis_heat_path_line_strip_get_vertex(path, istrip, ivert, &vtx)); 305 /* Count changes in vertex type along strips 306 * As we use vertex type instead of segment type, the vertices where type 307 * change are duplicated (with the only difference of their types) so that 308 * segments where type change are zero-length. 309 * This way Paraview displays segments with no misleading color gradient */ 310 if(ivert != 0 && vtx.type != prev_type) { 311 type_changes++; 312 strip_type_changes[istrip]++; 313 } 314 prev_type = vtx.type; 315 } 316 vcount_+= nverts; 317 } 318 if(vcount_ > ULONG_MAX) goto abort; 319 vcount = (unsigned long)vcount_; 320 /* Header */ 321 fprintf(stream, "# vtk DataFile Version 2.0\n"); 322 fprintf(stream, "Heat path\n"); 323 fprintf(stream, "ASCII\n"); 324 fprintf(stream, "DATASET POLYDATA\n"); 325 /* Write path positions */ 326 fprintf(stream, "POINTS %lu double\n", vcount + type_changes); 327 FOR_EACH(istrip, 0, scount_) { 328 size_t ivert, nverts; 329 enum sdis_heat_vertex_type prev_type = SDIS_HEAT_VERTEX_CONDUCTION; 330 ERR(sdis_heat_path_line_strip_get_vertices_count(path, istrip, &nverts)); 331 FOR_EACH(ivert, 0, nverts) { 332 struct sdis_heat_vertex vtx; 333 ERR(sdis_heat_path_line_strip_get_vertex(path, istrip, ivert, &vtx)); 334 if(ivert != 0 && vtx.type != prev_type) { 335 /* duplicate the previous vertex */ 336 struct sdis_heat_vertex v; 337 ERR(sdis_heat_path_line_strip_get_vertex(path, istrip, ivert-1, &v)); 338 fprintf(stream, "%g %g %g\n", SPLIT3(v.P)); 339 } 340 fprintf(stream, "%g %g %g\n", SPLIT3(vtx.P)); 341 prev_type = vtx.type; 342 } 343 } 344 /* Write the strips of the path 345 * Workaround a Paraview crash by creating 2-vertices-long paths from 346 * single-vertex paths */ 347 fprintf(stream, "LINES %lu %lu\n", 348 scount, scount + vcount + strip_1 + type_changes); 349 offset = 0; 350 FOR_EACH(istrip, 0, scount_) { 351 size_t nverts; 352 ERR(sdis_heat_path_line_strip_get_vertices_count(path, istrip, &nverts)); 353 if(nverts == 1) { 354 fprintf(stream, "2 %lu %lu\n", (unsigned long)offset, (unsigned long)offset); 355 } else { 356 size_t ivert; 357 fprintf(stream, "%lu", (unsigned long)nverts + strip_type_changes[istrip]); 358 FOR_EACH(ivert, 0, nverts + strip_type_changes[istrip]) { 359 if(ivert + offset > ULONG_MAX) goto abort; 360 fprintf(stream, " %lu", (unsigned long)(ivert + offset)); 361 } 362 fprintf(stream, "\n"); 363 } 364 offset += nverts + strip_type_changes[istrip]; 365 } 366 367 /* Write path status on strips */ 368 fprintf(stream, "CELL_DATA %lu\n", scount); 369 fprintf(stream, "SCALARS Path_Failure unsigned_char 1\n"); 370 fprintf(stream, "LOOKUP_TABLE default\n"); 371 FOR_EACH(istrip, 0, scount_) { 372 switch (status) { 373 case SDIS_HEAT_PATH_SUCCESS: fprintf(stream, "0\n"); break; 374 case SDIS_HEAT_PATH_FAILURE: fprintf(stream, "1\n"); break; 375 default: FATAL("Unreachable code.\n"); break; 376 } 377 } 378 fprintf(stream, "POINT_DATA %lu\n", vcount + type_changes); 379 /* Write the type of the random walk vertices */ 380 fprintf(stream, "SCALARS Segment_Type unsigned_char 1\n"); 381 fprintf(stream, "LOOKUP_TABLE default\n"); 382 FOR_EACH(istrip, 0, scount_) { 383 size_t ivert, nverts; 384 enum sdis_heat_vertex_type prev_type = SDIS_HEAT_VERTEX_CONDUCTION; 385 ERR(sdis_heat_path_line_strip_get_vertices_count(path, istrip, &nverts)); 386 FOR_EACH(ivert, 0, nverts) { 387 struct sdis_heat_vertex vtx; 388 int t; 389 ERR(sdis_heat_path_line_strip_get_vertex(path, istrip, ivert, &vtx)); 390 if((size_t)vtx.type > UCHAR_MAX) goto abort; 391 switch(vtx.type) { 392 case SDIS_HEAT_VERTEX_CONDUCTION: t = 0; break; 393 case SDIS_HEAT_VERTEX_CONVECTION: t = 1; break; 394 case SDIS_HEAT_VERTEX_RADIATIVE: t = 2; break; 395 default: FATAL("Unreachable code.\n"); break; 396 } 397 fprintf(stream, "%d\n", t); 398 if(ivert != 0 && vtx.type != prev_type) { 399 /* duplicate the previous vertex, except its type which is defined as 400 * the type of the current vertex */ 401 fprintf(stream, "%d\n", t); 402 } 403 prev_type = vtx.type; 404 } 405 } 406 /* Write the weights of the random walk vertices */ 407 fprintf(stream, "SCALARS Weight double 1\n"); 408 fprintf(stream, "LOOKUP_TABLE default\n"); 409 FOR_EACH(istrip, 0, scount_) { 410 size_t ivert, nverts; 411 enum sdis_heat_vertex_type prev_type = SDIS_HEAT_VERTEX_CONDUCTION; 412 double prev_w = 0; 413 ERR(sdis_heat_path_line_strip_get_vertices_count(path, istrip, &nverts)); 414 FOR_EACH(ivert, 0, nverts) { 415 struct sdis_heat_vertex vtx; 416 ERR(sdis_heat_path_line_strip_get_vertex(path, istrip, ivert, &vtx)); 417 fprintf(stream, "%g\n", vtx.weight); 418 if(ivert != 0 && vtx.type != prev_type) { 419 /* duplicate the previous vertex */ 420 fprintf(stream, "%g\n", prev_w); 421 } 422 prev_type = vtx.type; 423 prev_w = vtx.weight; 424 } 425 } 426 /* Write the branch_id of the random walk vertices */ 427 fprintf(stream, "SCALARS Branch_id int 1\n"); 428 fprintf(stream, "LOOKUP_TABLE default\n"); 429 FOR_EACH(istrip, 0, scount_) { 430 size_t ivert, nverts; 431 enum sdis_heat_vertex_type prev_type = SDIS_HEAT_VERTEX_CONDUCTION; 432 int prev_id = 0; 433 ERR(sdis_heat_path_line_strip_get_vertices_count(path, istrip, &nverts)); 434 FOR_EACH(ivert, 0, nverts) { 435 struct sdis_heat_vertex vtx; 436 ERR(sdis_heat_path_line_strip_get_vertex(path, istrip, ivert, &vtx)); 437 fprintf(stream, "%d\n", vtx.branch_id); 438 if(ivert != 0 && vtx.type != prev_type) { 439 /* duplicate the previous vertex */ 440 fprintf(stream, "%d\n", prev_id); 441 } 442 prev_type = vtx.type; 443 prev_id = vtx.branch_id; 444 } 445 } 446 /* If computation time is not INF, 447 * write the time of the random walk vertices */ 448 FOR_EACH(istrip, 0, scount_) { 449 size_t ivert, nverts; 450 enum sdis_heat_vertex_type prev_type = SDIS_HEAT_VERTEX_CONDUCTION; 451 double prev_time = 0; 452 ERR(sdis_heat_path_line_strip_get_vertices_count(path, istrip, &nverts)); 453 FOR_EACH(ivert, 0, nverts) { 454 struct sdis_heat_vertex vtx; 455 ERR(sdis_heat_path_line_strip_get_vertex(path, istrip, ivert, &vtx)); 456 if(IS_INF(vtx.time)) { 457 if(istrip || ivert) goto abort; 458 goto end_prt_time; 459 } 460 if(istrip == 0 && ivert == 0) { 461 fprintf(stream, "SCALARS Time double 1\n"); 462 fprintf(stream, "LOOKUP_TABLE default\n"); 463 } 464 fprintf(stream, "%g\n", vtx.time); 465 if(ivert != 0 && vtx.type != prev_type) { 466 /* duplicate the previous vertex */ 467 fprintf(stream, "%g\n", prev_time); 468 } 469 prev_type = vtx.type; 470 prev_time = vtx.time; 471 } 472 } 473 end_prt_time: 474 475 end: 476 MEM_RM(dump_ctx->stardis->allocator, name); 477 MEM_RM(dump_ctx->stardis->allocator, strip_type_changes); 478 if(stream) fclose(stream); 479 return res; 480 error: 481 goto end; 482 abort: 483 res = RES_BAD_ARG; 484 goto error; 485 } 486 487 res_T 488 dump_vtk_image 489 (const struct sdis_estimator_buffer* buf, 490 FILE* stream) 491 { 492 res_T res = RES_OK; 493 size_t def[2]; 494 unsigned long definition[2]; 495 double* temps = NULL; 496 size_t ix, iy; 497 498 ASSERT(buf && stream); 499 ERR(sdis_estimator_buffer_get_definition(buf, def)); 500 if(def[0] == 0 || def[1] == 0 || def[0] * def[1] > ULONG_MAX) goto abort; 501 definition[0] = (unsigned long)def[0]; 502 definition[1] = (unsigned long)def[1]; 503 504 temps = mem_alloc(definition[0] * definition[1] * sizeof(double)); 505 if(temps == NULL) { 506 res = RES_MEM_ERR; 507 goto error; 508 } 509 510 /* Compute the per pixel temperature */ 511 fprintf(stream, "# vtk DataFile Version 2.0\n"); 512 fprintf(stream, "Infrared Image\n"); 513 fprintf(stream, "ASCII\n"); 514 fprintf(stream, "DATASET STRUCTURED_POINTS\n"); 515 fprintf(stream, "DIMENSIONS %lu %lu 1\n", definition[0], definition[1]); 516 fprintf(stream, "ORIGIN 0 0 0\n"); 517 fprintf(stream, "SPACING 1 1 1\n"); 518 fprintf(stream, "POINT_DATA %lu\n", definition[0] * definition[1]); 519 fprintf(stream, "SCALARS temperature_estimate float 1\n"); 520 fprintf(stream, "LOOKUP_TABLE default\n"); 521 FOR_EACH(iy, 0, definition[1]) { 522 FOR_EACH(ix, 0, definition[0]) { 523 const struct sdis_estimator* estimator; 524 struct sdis_mc T; 525 ERR(sdis_estimator_buffer_at(buf, ix, iy, &estimator)); 526 ERR(sdis_estimator_get_temperature(estimator, &T)); 527 fprintf(stream, "%f\n", T.E); 528 } 529 } 530 fprintf(stream, "SCALARS temperature_std_dev float 1\n"); 531 fprintf(stream, "LOOKUP_TABLE default\n"); 532 FOR_EACH(iy, 0, definition[1]) { 533 FOR_EACH(ix, 0, definition[0]) { 534 const struct sdis_estimator* estimator; 535 struct sdis_mc T; 536 ERR(sdis_estimator_buffer_at(buf, ix, iy, &estimator)); 537 ERR(sdis_estimator_get_temperature(estimator, &T)); 538 fprintf(stream, "%f\n", T.SE); 539 } 540 } 541 fprintf(stream, "SCALARS computation_time float 1\n"); 542 fprintf(stream, "LOOKUP_TABLE default\n"); 543 FOR_EACH(iy, 0, definition[1]) { 544 FOR_EACH(ix, 0, definition[0]) { 545 const struct sdis_estimator* estimator; 546 struct sdis_mc time; 547 ERR(sdis_estimator_buffer_at(buf, ix, iy, &estimator)); 548 ERR(sdis_estimator_get_realisation_time(estimator, &time)); 549 fprintf(stream, "%f\n", time.E); 550 } 551 } 552 fprintf(stream, "SCALARS computation_time_std_dev float 1\n"); 553 fprintf(stream, "LOOKUP_TABLE default\n"); 554 FOR_EACH(iy, 0, definition[1]) { 555 FOR_EACH(ix, 0, definition[0]) { 556 const struct sdis_estimator* estimator; 557 struct sdis_mc time; 558 ERR(sdis_estimator_buffer_at(buf, ix, iy, &estimator)); 559 ERR(sdis_estimator_get_realisation_time(estimator, &time)); 560 fprintf(stream, "%f\n", time.SE); 561 } 562 } 563 fprintf(stream, "SCALARS failures_count unsigned_long_long 1\n"); 564 fprintf(stream, "LOOKUP_TABLE default\n"); 565 FOR_EACH(iy, 0, definition[1]) { 566 FOR_EACH(ix, 0, definition[0]) { 567 const struct sdis_estimator* estimator; 568 size_t nfails; 569 ERR(sdis_estimator_buffer_at(buf, ix, iy, &estimator)); 570 ERR(sdis_estimator_get_failure_count(estimator, &nfails)); 571 if(nfails > ULONG_MAX) goto abort; 572 fprintf(stream, "%lu\n", (unsigned long)nfails); 573 } 574 } 575 mem_rm(temps); 576 577 end: 578 return res; 579 error: 580 goto end; 581 abort: 582 res = RES_BAD_ARG; 583 goto error; 584 } 585 586 res_T 587 dump_ht_image 588 (const struct sdis_estimator_buffer* buf, 589 FILE* stream) 590 { 591 res_T res = RES_OK; 592 size_t def[2]; 593 unsigned long definition[2]; 594 size_t ix, iy; 595 596 ASSERT(buf && stream); 597 ERR(sdis_estimator_buffer_get_definition(buf, def)); 598 if(def[0] > ULONG_MAX || def[1] > ULONG_MAX) goto abort; 599 definition[0] = (unsigned long)def[0]; 600 definition[1] = (unsigned long)def[1]; 601 602 fprintf(stream, "%lu %lu\n", definition[0], definition[1]); 603 FOR_EACH(iy, 0, definition[1]) { 604 FOR_EACH(ix, 0, definition[0]) { 605 const struct sdis_estimator* estimator; 606 struct sdis_mc T; 607 struct sdis_mc time; 608 ERR(sdis_estimator_buffer_at(buf, ix, iy, &estimator)); 609 ERR(sdis_estimator_get_realisation_time(estimator, &time)); 610 ERR(sdis_estimator_get_temperature(estimator, &T)); 611 fprintf(stream, "%f %f 0 0 0 0 %f %f\n", 612 T.E, T.SE, time.E, time.SE); 613 } 614 }; 615 616 end: 617 return res; 618 error: 619 goto end; 620 abort: 621 res = RES_BAD_ARG; 622 goto error; 623 } 624 625 #define FW(Ptr, Count) \ 626 if((Count) != fwrite((Ptr), sizeof(*(Ptr)), (Count), stream)) { \ 627 res = RES_IO_ERR; \ 628 goto error; \ 629 } 630 631 static FINLINE double 632 medium_get_t0 633 (struct sdis_medium* medium) 634 { 635 struct sdis_data* data = NULL; 636 enum sdis_medium_type type; 637 ASSERT(medium); 638 type = sdis_medium_get_type(medium); 639 data = sdis_medium_get_data(medium); 640 if(type == SDIS_FLUID) { 641 const struct fluid* fluid_props = *((const struct fluid**)sdis_data_cget(data)); 642 return fluid_props->t0; 643 } else { 644 const struct solid* solid_props = *((const struct solid**)sdis_data_cget(data)); 645 ASSERT(type == SDIS_SOLID); 646 return solid_props->t0; 647 } 648 } 649 650 static res_T 651 dump_sample_end(struct sdis_green_path* path, void* ctx) 652 { 653 /* Stardis */ 654 struct sdis_green_path_end end = SDIS_GREEN_PATH_END_NULL; 655 struct sdis_data* data = NULL; 656 enum sdis_medium_type type; 657 658 /* Stream */ 659 struct e_ctx* e_ctx = ctx; 660 FILE* stream; 661 662 /* Miscellaneous */ 663 const struct description* descs; 664 double* pos; 665 double elapsed; 666 size_t sz; 667 unsigned trad_id; 668 unsigned id; 669 res_T res = RES_OK; 670 671 CHK(path && ctx); 672 673 stream = e_ctx->stream; 674 ERR(sdis_green_path_get_elapsed_time(path, &elapsed)); 675 ERR(sdis_green_path_get_end(path, &end)); 676 677 sz = darray_descriptions_size_get(e_ctx->desc); 678 if(sz > UINT_MAX) { res = RES_BAD_ARG; goto error; } 679 trad_id = (unsigned)sz; 680 681 descs = darray_descriptions_cdata_get(e_ctx->desc); 682 switch(end.type) { 683 case SDIS_GREEN_PATH_END_AT_RADIATIVE_ENV: 684 /* End, End ID, X, Y, Z, Elapsed time */ 685 fprintf(stream, "TRAD, %u, 0, 0, 0, %g\n", trad_id, elapsed); 686 break; 687 case SDIS_GREEN_PATH_END_AT_INTERFACE: { 688 struct intface* d__; 689 data = sdis_interface_get_data(end.data.itfrag.intface); 690 pos = end.data.itfrag.fragment.P; 691 d__ = sdis_data_get(data); 692 id = d__->desc_id; 693 CHK(DESC_IS_T(descs+id) || DESC_IS_H(descs+id)); 694 /* End, End ID, X, Y, Z, Elapsed time */ 695 fprintf(stream, "%s, %u, %g, %g, %g, %g\n", 696 str_cget(get_description_name(descs + id)), id, SPLIT3(pos), elapsed); 697 break; 698 } 699 case SDIS_GREEN_PATH_END_IN_VOLUME: 700 type = sdis_medium_get_type(end.data.mdmvert.medium); 701 data = sdis_medium_get_data(end.data.mdmvert.medium); 702 pos = end.data.mdmvert.vertex.P; 703 if(end.data.mdmvert.vertex.P[0] == INF) { 704 /* Radiative output (Trad) */ 705 id = trad_id; 706 } 707 else if(type == SDIS_FLUID) { 708 struct fluid** pfluid = sdis_data_get(data); 709 id = (*pfluid)->desc_id; 710 } else { 711 struct solid** psolid = sdis_data_get(data); 712 ASSERT(type == SDIS_SOLID); 713 ASSERT(!(*psolid)->is_outside); /* FIXME: what if in external solid? */ 714 id = (*psolid)->desc_id; 715 } 716 /* End, End ID, X, Y, Z, Elapsed time */ 717 fprintf(stream, "%s, %u, %g, %g, %g, %g\n", 718 str_cget(get_description_name(descs + id)), id, SPLIT3(pos), elapsed); 719 break; 720 default: FATAL("Unreachable code.\n"); break; 721 } 722 723 end: 724 return res; 725 error: 726 goto end; 727 } 728 729 static res_T 730 dump_sample 731 (struct sdis_green_path* path, 732 void* ctx) 733 { 734 /* Stardis variables */ 735 struct sdis_green_path_end path_end = SDIS_GREEN_PATH_END_NULL; 736 struct sdis_data* data = NULL; 737 enum sdis_medium_type type; 738 739 /* Miscellaneous variables */ 740 struct w_ctx* w_ctx = ctx; 741 FILE* stream; 742 const struct description* descs; 743 struct htable_weigth_iterator it, end; 744 struct green_sample_header header; 745 unsigned* ids = NULL; 746 double* weights = NULL; 747 double t0; 748 size_t sz, i; 749 unsigned trad_id; 750 res_T res = RES_OK; 751 752 CHK(path && ctx); 753 754 stream = w_ctx->stream; 755 ERR(sdis_green_path_get_end(path, &path_end)); 756 sz = darray_descriptions_size_get(w_ctx->desc); 757 if(sz > UINT_MAX) { res = RES_BAD_ARG; goto error; } 758 trad_id = (unsigned)sz; 759 760 /* For each path, dump: 761 * # end_id #power_terms #flux_terms 762 * power_id_1 ... power_id_n flux_id_1 ... flux_id_n 763 * power_factor_1 ... power_factor_n flux_factor_1 ... flux_factor_n 764 */ 765 766 descs = darray_descriptions_cdata_get(w_ctx->desc); 767 switch(path_end.type) { 768 case SDIS_GREEN_PATH_END_AT_RADIATIVE_ENV: 769 header.at_initial = 0; 770 header.sample_end_description_id = trad_id; 771 break; 772 case SDIS_GREEN_PATH_END_AT_INTERFACE: { 773 struct intface* d__; 774 unsigned desc_id; 775 data = sdis_interface_get_data(path_end.data.itfrag.intface); 776 d__ = sdis_data_get(data); 777 desc_id = d__->desc_id; 778 CHK(DESC_IS_T(descs+desc_id) || DESC_IS_H(descs+desc_id)); 779 header.sample_end_description_id = desc_id; 780 header.at_initial = 0; 781 break; 782 } 783 case SDIS_GREEN_PATH_END_IN_VOLUME: 784 type = sdis_medium_get_type(path_end.data.mdmvert.medium); 785 data = sdis_medium_get_data(path_end.data.mdmvert.medium); 786 t0 = medium_get_t0(path_end.data.mdmvert.medium); 787 header.at_initial = (path_end.data.mdmvert.vertex.time <= t0); 788 if(path_end.data.mdmvert.vertex.P[0] == INF) { 789 /* Radiative output (Trad) */ 790 header.sample_end_description_id = trad_id; 791 } 792 else if(type == SDIS_FLUID) { 793 struct fluid** pfluid = sdis_data_get(data); 794 header.sample_end_description_id = (*pfluid)->desc_id; 795 } else { 796 struct solid** psolid = sdis_data_get(data); 797 ASSERT(type == SDIS_SOLID); 798 ASSERT(!(*psolid)->is_outside); /* FIXME: what if in external solid? */ 799 header.sample_end_description_id = (*psolid)->desc_id; 800 } 801 break; 802 default: FATAL("Unreachable code.\n"); break; 803 } 804 805 /* Merge power and flux terms */ 806 htable_weigth_clear(&w_ctx->pw); 807 htable_weigth_clear(&w_ctx->flux); 808 ERR(sdis_green_path_for_each_power_term(path, merge_power_terms, w_ctx)); 809 ERR(sdis_green_path_for_each_flux_term(path, merge_flux_terms, w_ctx)); 810 sz = htable_weigth_size_get(&w_ctx->pw); 811 if(sz > UINT_MAX) { res = RES_BAD_ARG; goto error; } 812 header.pw_count = (unsigned)sz; 813 sz = htable_weigth_size_get(&w_ctx->flux); 814 if(sz > UINT_MAX) { res = RES_BAD_ARG; goto error; } 815 header.fx_count = (unsigned)sz; 816 817 /* Write path's header */ 818 FW(&header, 1); 819 820 /* Allocate buffers */ 821 sz = header.pw_count + header.fx_count; 822 ids = MEM_CALLOC(w_ctx->alloc, sz, sizeof(*ids)); 823 weights = MEM_CALLOC(w_ctx->alloc, sz, sizeof(*weights)); 824 if(!ids || !weights) { 825 res = RES_MEM_ERR; 826 goto error; 827 } 828 829 /* Write terms */ 830 htable_weigth_begin(&w_ctx->pw, &it); 831 htable_weigth_end(&w_ctx->pw, &end); 832 i = 0; 833 while(!htable_weigth_iterator_eq(&it, &end)) { 834 double* w = htable_weigth_iterator_data_get(&it); 835 unsigned* k = htable_weigth_iterator_key_get(&it); 836 CHK(*k <= trad_id); 837 ids[i] = *k; 838 weights[i] = *w; 839 htable_weigth_iterator_next(&it); 840 i++; 841 } 842 CHK(i == header.pw_count); 843 844 htable_weigth_begin(&w_ctx->flux, &it); 845 htable_weigth_end(&w_ctx->flux, &end); 846 while (!htable_weigth_iterator_eq(&it, &end)) { 847 double* w = htable_weigth_iterator_data_get(&it); 848 unsigned* k = htable_weigth_iterator_key_get(&it); 849 CHK(*k <= trad_id); 850 ids[i] = *k; 851 weights[i] = *w; 852 htable_weigth_iterator_next(&it); 853 i++; 854 } 855 CHK(i == header.pw_count + header.fx_count); 856 857 FW(ids, sz); 858 FW(weights, sz); 859 860 end: 861 MEM_RM(w_ctx->alloc, ids); 862 MEM_RM(w_ctx->alloc, weights); 863 return res; 864 error: 865 goto end; 866 } 867 868 res_T 869 dump_green_bin 870 (struct sdis_green_function* green, 871 const struct stardis* stardis, 872 FILE* stream) 873 { 874 /* The following type must be identical to its stardis-green counterpart! */ 875 struct green_file_header header; 876 const struct radiative_env_const* radenv_const = NULL; 877 struct w_ctx w_ctx; 878 size_t sz, i; 879 int table_initialized = 0; 880 res_T res = RES_OK; 881 882 ASSERT(green && stardis && stream); 883 884 /* Stardis can produce the green function on systems 885 * with constant properties only */ 886 ASSERT(stardis->radenv.type == RADIATIVE_ENV_CONST); 887 radenv_const = &stardis->radenv.data.cst; 888 889 /* Init header */ 890 strcpy(header.green_string, BIN_FILE_IDENT_STRING); 891 header.file_format_version = GREEN_FILE_FORMAT_VERSION; 892 header.solid_count = stardis->counts.smed_count; 893 header.fluid_count = stardis->counts.fmed_count; 894 header.tbound_count = stardis->counts.tbound_count; 895 header.hbound_count = stardis->counts.hbound_count; 896 header.fbound_count = stardis->counts.fbound_count; 897 header.sfconnect_count = stardis->counts.sfconnect_count; 898 header.ssconnect_count = stardis->counts.ssconnect_count; 899 ERR(sdis_green_function_get_paths_count(green, &header.ok_count)); 900 ERR(sdis_green_function_get_invalid_paths_count(green, &header.failed_count)); 901 sz = darray_descriptions_size_get(&stardis->descriptions); 902 if(sz > UINT_MAX) goto abort; 903 ASSERT(sz == 904 (stardis->counts.smed_count + stardis->counts.fmed_count 905 + stardis->counts.tbound_count + stardis->counts.hbound_count 906 + stardis->counts.fbound_count + stardis->counts.sfconnect_count 907 + stardis->counts.ssconnect_count)); 908 header.description_count = (unsigned)sz; 909 header.ambient_radiative_temperature = radenv_const->temperature; 910 header.ambient_radiative_temperature_reference = 911 radenv_const->reference_temperature; 912 d2_set(header.time_range, stardis->time_range); 913 914 /* Write header */ 915 FW(&header, 1); 916 917 /* Write descriptions*/ 918 for(i = 0; i < sz; i++) { 919 struct green_description desc; 920 ERR(copy_desc_to_green_desc(&desc, &stardis->descriptions, i)); 921 FW(&desc, 1); 922 } 923 924 w_ctx.alloc = stardis->allocator; 925 w_ctx.desc = &stardis->descriptions; 926 htable_weigth_init(stardis->allocator, &w_ctx.pw); 927 htable_weigth_init(stardis->allocator, &w_ctx.flux); 928 w_ctx.stream = stream; 929 table_initialized = 1; 930 931 /* Write samples */ 932 ERR(sdis_green_function_for_each_path(green, dump_sample, &w_ctx)); 933 934 end: 935 if(table_initialized) htable_weigth_release(&w_ctx.pw); 936 if(table_initialized) htable_weigth_release(&w_ctx.flux); 937 return res; 938 error: 939 goto end; 940 abort: 941 res = RES_BAD_ARG; 942 goto error; 943 } 944 945 res_T 946 dump_paths_end 947 (struct sdis_green_function* green, 948 const struct stardis* stardis, 949 FILE* stream) 950 { 951 res_T res = RES_OK; 952 struct e_ctx e_ctx = { 0 }; 953 954 ASSERT(green && stardis && stream); 955 956 e_ctx.desc = &stardis->descriptions; 957 e_ctx.stream = stream; 958 959 fprintf(stream, "\"End\", \"End ID\", \"X\", \"Y\", \"Z\", \"Elapsed time\"\n"); 960 ERR(sdis_green_function_for_each_path(green, dump_sample_end, &e_ctx)); 961 962 end: 963 return res; 964 error: 965 goto end; 966 } 967 968 res_T 969 print_sample 970 (struct sdis_green_path* path, 971 void* ctx) 972 { 973 res_T res = RES_OK; 974 struct sdis_green_path_end path_end = SDIS_GREEN_PATH_END_NULL; 975 struct sdis_data* data = NULL; 976 enum sdis_medium_type type; 977 struct htable_weigth_iterator it, end; 978 unsigned desc_id; 979 size_t pw_count, fx_count; 980 struct w_ctx* w_ctx = ctx; 981 const struct description* descs; 982 CHK(path && ctx); 983 984 ERR(sdis_green_path_get_end(path, &path_end)); 985 986 /* For each path, prints: 987 * # end #power_terms #flux_terms power_term_1 ... power_term_n flux_term_1 ... flux_term_n 988 * with: 989 * - end = end_type end_id; end_type = T | H | R | F | S 990 * - power_term_i = power_type_i power_id_i factor_i 991 * - flux_term_i = flux_id_i factor_i 992 */ 993 994 descs = darray_descriptions_cdata_get(w_ctx->desc); 995 switch (path_end.type) { 996 case SDIS_GREEN_PATH_END_AT_INTERFACE: { 997 struct intface* d__; 998 data = sdis_interface_get_data(path_end.data.itfrag.intface); 999 d__ = sdis_data_get(data); 1000 desc_id = d__->desc_id; 1001 switch (descs[desc_id].type) { 1002 case DESC_BOUND_T_FOR_SOLID: 1003 fprintf(w_ctx->stream, "T\t%u", desc_id); 1004 break; 1005 case DESC_BOUND_H_FOR_SOLID: 1006 case DESC_BOUND_H_FOR_FLUID: 1007 fprintf(w_ctx->stream, "H\t%u", desc_id); 1008 break; 1009 case DESC_BOUND_F_FOR_SOLID: 1010 FATAL("Heat path cannot end at a flux boundary.\n"); break; 1011 break; 1012 default: FATAL("Unreachable code.\n"); break; 1013 } 1014 break; 1015 } 1016 case SDIS_GREEN_PATH_END_IN_VOLUME: 1017 type = sdis_medium_get_type(path_end.data.mdmvert.medium); 1018 data = sdis_medium_get_data(path_end.data.mdmvert.medium); 1019 if(path_end.data.mdmvert.vertex.P[0] == INF) { 1020 /* Radiative output (Trad)*/ 1021 size_t sz = darray_descriptions_size_get(w_ctx->desc); 1022 if(sz > UINT_MAX) goto abort; 1023 fprintf(w_ctx->stream, "R\t%u", (unsigned)sz); 1024 } 1025 else if(type == SDIS_FLUID) { 1026 struct fluid** pfluid = sdis_data_get(data); 1027 desc_id = (*pfluid)->desc_id; 1028 if((*pfluid)->is_outside) 1029 /* If outside the model and in a fluid with known temperature, 1030 * its a fluid attached to a H boundary */ 1031 fprintf(w_ctx->stream, "H\t%u", desc_id); 1032 /* In a standard fluid with known temperature */ 1033 else fprintf(w_ctx->stream, "F\t%u", desc_id); 1034 } else { 1035 struct solid** psolid = sdis_data_get(data); 1036 ASSERT(type == SDIS_SOLID); 1037 ASSERT(!(*psolid)->is_outside); /* FIXME: what if in external solid? */ 1038 desc_id = (*psolid)->desc_id; 1039 fprintf(w_ctx->stream, "S\t%u", desc_id); 1040 } 1041 break; 1042 default: FATAL("Unreachable code.\n"); break; 1043 } 1044 1045 ERR(sdis_green_function_get_power_terms_count(path, &pw_count)); 1046 1047 htable_weigth_clear(&w_ctx->pw); 1048 htable_weigth_clear(&w_ctx->flux); 1049 ERR(sdis_green_path_for_each_power_term(path, merge_power_terms, w_ctx)); 1050 ERR(sdis_green_path_for_each_flux_term(path, merge_flux_terms, w_ctx)); 1051 fx_count = htable_weigth_size_get(&w_ctx->flux); 1052 1053 if(pw_count > ULONG_MAX || fx_count > ULONG_MAX) goto abort; 1054 fprintf(w_ctx->stream, "\t%lu\t%lu", 1055 (unsigned long)pw_count, (unsigned long)fx_count); 1056 1057 htable_weigth_begin(&w_ctx->pw, &it); 1058 htable_weigth_end(&w_ctx->pw, &end); 1059 while(!htable_weigth_iterator_eq(&it, &end)) { 1060 double* w = htable_weigth_iterator_data_get(&it); 1061 unsigned* k = htable_weigth_iterator_key_get(&it); 1062 fprintf(w_ctx->stream, "\t%u\t%g", *k, *w); 1063 htable_weigth_iterator_next(&it); 1064 } 1065 1066 htable_weigth_begin(&w_ctx->flux, &it); 1067 htable_weigth_end(&w_ctx->flux, &end); 1068 while (!htable_weigth_iterator_eq(&it, &end)) { 1069 double* w = htable_weigth_iterator_data_get(&it); 1070 unsigned* k = htable_weigth_iterator_key_get(&it); 1071 fprintf(w_ctx->stream, "\t%u\t%g", *k, *w); 1072 htable_weigth_iterator_next(&it); 1073 } 1074 fprintf(w_ctx->stream, "\n"); 1075 1076 end: 1077 return res; 1078 error: 1079 goto end; 1080 abort: 1081 res = RES_BAD_ARG; 1082 goto error; 1083 } 1084 1085 res_T 1086 dump_green_ascii 1087 (struct sdis_green_function* green, 1088 const struct stardis* stardis, 1089 FILE* stream) 1090 { 1091 res_T res = RES_OK; 1092 const struct radiative_env_const* radenv_const = NULL; 1093 unsigned ok_count, failed_count; 1094 size_t sz; 1095 struct w_ctx w_ctx; 1096 int table_initialized = 0; 1097 unsigned i, szd; 1098 const struct description* descs; 1099 1100 ASSERT(green && stardis && stream); 1101 1102 /* Stardis can produce the green function on systems 1103 * with constant properties only */ 1104 ASSERT(stardis->radenv.type == RADIATIVE_ENV_CONST); 1105 radenv_const = &stardis->radenv.data.cst; 1106 1107 ERR(sdis_green_function_get_paths_count(green, &sz)); 1108 if(sz > UINT_MAX) goto abort; 1109 ok_count = (unsigned)sz; 1110 ERR(sdis_green_function_get_invalid_paths_count(green, &sz)); 1111 if(sz > UINT_MAX) goto abort; 1112 failed_count = (unsigned)sz; 1113 sz = darray_descriptions_size_get(&stardis->descriptions); 1114 if(sz > UINT_MAX) goto abort; 1115 szd = (unsigned)sz; 1116 descs = darray_descriptions_cdata_get(&stardis->descriptions); 1117 1118 /* Output counts */ 1119 fprintf(stream, "---BEGIN GREEN---\n"); 1120 fprintf(stream, "# time range\n"); 1121 fprintf(stream, "%g %g\n", 1122 SPLIT2(stardis->time_range)); 1123 fprintf(stream, 1124 "# #solids #fluids #t_boundaries #h_boundaries #f_boundaries #ok #failures\n"); 1125 fprintf(stream, "%u %u %u %u %u %u %u\n", 1126 stardis->counts.smed_count, stardis->counts.fmed_count, 1127 stardis->counts.tbound_count, stardis->counts.hbound_count, 1128 stardis->counts.fbound_count, ok_count, failed_count); 1129 1130 /* List Media */ 1131 if(stardis->counts.smed_count) { 1132 fprintf(stream, "# Solids\n"); 1133 fprintf(stream, "# ID Name lambda rho cp power initial_temp imposed_temp\n"); 1134 FOR_EACH(i, 0, szd) { 1135 const struct description* desc = descs + i; 1136 const struct solid* sl; 1137 if(desc->type != DESC_MAT_SOLID) continue; 1138 sl = desc->d.solid; 1139 fprintf(stream, "%u\t%s\t%g\t%g\t%g\t%g", 1140 i, str_cget(&sl->name), sl->lambda, sl->rho, sl->cp, sl->vpower); 1141 if(SDIS_TEMPERATURE_IS_KNOWN(sl->tinit)) { 1142 fprintf(stream, "\t%g", sl->tinit); 1143 } else { 1144 fprintf(stream, "\tNONE"); 1145 } 1146 if(SDIS_TEMPERATURE_IS_KNOWN(sl->imposed_temperature)) { 1147 fprintf(stream, "\t%g\n", sl->imposed_temperature); 1148 } else { 1149 fprintf(stream, "\tNONE\n"); 1150 } 1151 } 1152 } 1153 if(stardis->counts.fmed_count) { 1154 fprintf(stream, "# Fluids\n"); 1155 fprintf(stream, "# ID Name rho cp initial_temp imposed_temp\n"); 1156 FOR_EACH(i, 0, szd) { 1157 const struct description* desc = descs + i; 1158 const struct fluid* fl; 1159 if(desc->type != DESC_MAT_FLUID) continue; 1160 fl = desc->d.fluid; 1161 if(SDIS_TEMPERATURE_IS_KNOWN(fl->imposed_temperature)) { 1162 fprintf(stream, "%u\t%s\t%g\t%g", 1163 i, str_cget(&fl->name), fl->rho, fl->cp); 1164 } else { 1165 fprintf(stream, "%u\t%s\t%g\t%g", 1166 i, str_cget(&fl->name), fl->rho, fl->cp); 1167 } 1168 if(SDIS_TEMPERATURE_IS_KNOWN(fl->tinit)) { 1169 fprintf(stream, "\t%g", fl->tinit); 1170 } else { 1171 fprintf(stream, "\tNONE"); 1172 } 1173 if(SDIS_TEMPERATURE_IS_KNOWN(fl->imposed_temperature)) { 1174 fprintf(stream, "\t%g\n", fl->imposed_temperature); 1175 } else { 1176 fprintf(stream, "\tNONE\n"); 1177 } 1178 } 1179 } 1180 1181 /* List Boundaries */ 1182 if(stardis->counts.tbound_count) { 1183 fprintf(stream, "# T Boundaries\n"); 1184 fprintf(stream, "# ID Name temperature\n"); 1185 FOR_EACH(i, 0, szd) { 1186 const struct description* desc = descs + i; 1187 const struct t_boundary* bd; 1188 bd = desc->d.t_boundary; 1189 fprintf(stream, "%u\t%s\t%g\n", 1190 i, str_cget(&bd->name), bd->imposed_temperature); 1191 } 1192 } 1193 if(stardis->counts.hbound_count) { 1194 fprintf(stream, "# H Boundaries\n"); 1195 fprintf(stream, "# ID Name ref_temperature emissivity specular_fraction hc T_env\n"); 1196 FOR_EACH(i, 0, szd) { 1197 const struct description* desc = descs + i; 1198 const struct h_boundary* bd; 1199 if(desc->type != DESC_BOUND_H_FOR_SOLID 1200 && desc->type != DESC_BOUND_H_FOR_FLUID) continue; 1201 bd = desc->d.h_boundary; 1202 fprintf(stream, "%u\t%s\t%g\t%g\t%g\t%g\t%g\n", 1203 i, str_cget(&bd->name), bd->ref_temperature, bd->emissivity, 1204 bd->specular_fraction, bd->hc, bd->imposed_temperature); 1205 } 1206 } 1207 if(stardis->counts.fbound_count) { 1208 fprintf(stream, "# F Boundaries\n"); 1209 fprintf(stream, "# ID Name flux\n"); 1210 FOR_EACH(i, 0, szd) { 1211 const struct description* desc = descs + i; 1212 const struct f_boundary* bd; 1213 if(desc->type != DESC_BOUND_F_FOR_SOLID) continue; 1214 bd = desc->d.f_boundary; 1215 fprintf(stream, "%u\t%s\t%g\n", 1216 i, str_cget(&bd->name), bd->imposed_flux); 1217 } 1218 } 1219 1220 /* Radiative Temperature */ 1221 fprintf(stream, "# Radiative Temperatures\n"); 1222 fprintf(stream, "# ID Trad Trad_Ref\n"); 1223 fprintf(stream, "%u\t%g\t%g\n", 1224 szd, radenv_const->temperature, radenv_const->reference_temperature); 1225 1226 fprintf(stream, "# Samples\n"); 1227 fprintf(stream, 1228 "# end #power_terms #flux_terms power_term_1 ... power_term_n flux_term_1 ... flux_term_n\n"); 1229 fprintf(stream, "# end = end_type end_id; end_type = T | H | R | F | S\n"); 1230 fprintf(stream, "# power_term_i = power_id_i factor_i\n"); 1231 fprintf(stream, "# flux_term_i = flux_id_i factor_i\n"); 1232 1233 w_ctx.alloc = stardis->allocator; 1234 w_ctx.desc = &stardis->descriptions; 1235 htable_weigth_init(stardis->allocator, &w_ctx.pw); 1236 htable_weigth_init(stardis->allocator, &w_ctx.flux); 1237 w_ctx.stream = stream; 1238 table_initialized = 1; 1239 1240 ERR(sdis_green_function_for_each_path(green, print_sample, &w_ctx)); 1241 1242 fprintf(stream, "---END GREEN---\n"); 1243 1244 end: 1245 if(table_initialized) htable_weigth_release(&w_ctx.pw); 1246 if(table_initialized) htable_weigth_release(&w_ctx.flux); 1247 return res; 1248 error: 1249 goto end; 1250 abort: 1251 res = RES_BAD_ARG; 1252 goto error; 1253 } 1254 1255 res_T 1256 dump_boundaries_at_the_end_of_vtk 1257 (const struct stardis* stardis, 1258 FILE* stream) 1259 { 1260 res_T res = RES_OK; 1261 const struct description* descriptions; 1262 unsigned tsz, t; 1263 ASSERT(stardis && stream); 1264 1265 ERR(sg3d_geometry_get_unique_triangles_count(stardis->geometry.sg3d, &tsz)); 1266 descriptions = darray_descriptions_cdata_get(&stardis->descriptions); 1267 1268 fprintf(stream, "SCALARS Boundaries unsigned_int 1\n"); 1269 fprintf(stream, "LOOKUP_TABLE default\n"); 1270 FOR_EACH(t, 0, tsz) { 1271 unsigned descr[SG3D_PROP_TYPES_COUNT__]; 1272 ERR(sg3d_geometry_get_unique_triangle_properties(stardis->geometry.sg3d, t, 1273 descr)); 1274 if(descr[SG3D_INTFACE] != SG3D_UNSPECIFIED_PROPERTY 1275 && DESC_IS_BOUNDARY(descriptions+descr[SG3D_INTFACE])) 1276 /* Descriptions are numbered from 1 in the log (so the 1+ below) */ 1277 fprintf(stream, "%u\n", 1 + descr[SG3D_INTFACE]); 1278 else fprintf(stream, "%u\n", SG3D_UNSPECIFIED_PROPERTY); 1279 } 1280 1281 exit: 1282 return res; 1283 error: 1284 goto exit; 1285 } 1286 1287 res_T 1288 dump_enclosure_related_stuff_at_the_end_of_vtk 1289 (struct stardis* stardis, 1290 FILE* stream) 1291 { 1292 res_T res = RES_OK; 1293 unsigned* trgs = NULL; 1294 struct senc3d_enclosure* enc = NULL; 1295 unsigned tsz, e, s, t, scount, ecount, ocount; 1296 int* enc_status = NULL; 1297 const struct description* descriptions; 1298 int undef_count = 0, multi_count = 0; 1299 struct str msg; 1300 ASSERT(stardis && stream); 1301 1302 str_init(stardis->allocator, &msg); 1303 descriptions = darray_descriptions_cdata_get(&stardis->descriptions); 1304 ERR(sg3d_geometry_get_unique_triangles_count(stardis->geometry.sg3d, &tsz)); 1305 trgs = MEM_CALLOC(stardis->allocator, tsz, sizeof(*trgs)); 1306 if(!trgs) { 1307 res = RES_MEM_ERR; 1308 goto error; 1309 } 1310 1311 /* If enclosure where not extracted, dump only errors */ 1312 ERR(senc3d_scene_get_overlapping_triangles_count(stardis->senc3d_scn, &ocount)); 1313 if(ocount) { 1314 FOR_EACH(t, 0, tsz) trgs[t] = 0; 1315 FOR_EACH(t, 0, ocount) { 1316 unsigned trid; 1317 ERR(senc3d_scene_get_overlapping_triangle(stardis->senc3d_scn, t, &trid)); 1318 trgs[trid] = 1; 1319 } 1320 fprintf(stream, "SCALARS Overlapping_triangles unsigned_int 1\n"); 1321 fprintf(stream, "LOOKUP_TABLE default\n"); 1322 FOR_EACH(t, 0, tsz) fprintf(stream, "%u\n", trgs[t]); 1323 goto exit; 1324 } 1325 1326 /* Keep the segments involved in holes (not the vertices) */ 1327 ERR(senc3d_scene_get_frontier_segments_count(stardis->senc3d_scn, &scount)); 1328 if(scount) { 1329 /* Room to store frontier triangles */ 1330 FOR_EACH(s, 0, scount) { 1331 unsigned vrtc[2], trid; 1332 ERR(senc3d_scene_get_frontier_segment(stardis->senc3d_scn, s, vrtc, &trid)); 1333 trgs[trid] = 1; 1334 } 1335 logger_print(stardis->logger, LOG_WARNING, "Model contains hole(s).\n"); 1336 fprintf(stream, "SCALARS Hole_frontiers unsigned_int 1\n"); 1337 fprintf(stream, "LOOKUP_TABLE default\n"); 1338 FOR_EACH(t, 0, tsz) fprintf(stream, "%u\n", trgs[t]); 1339 } 1340 1341 /* Dump enclosure information */ 1342 ERR(senc3d_scene_get_enclosure_count(stardis->senc3d_scn, &ecount)); 1343 enc_status = MEM_CALLOC(stardis->allocator, ecount, sizeof(*enc_status)); 1344 if(!enc_status) { 1345 res = RES_MEM_ERR; 1346 goto error; 1347 } 1348 ASSERT(stardis->undefined_medium_behind_boundary_id != SENC3D_UNSPECIFIED_MEDIUM); 1349 FOR_EACH(e, 0, ecount) { 1350 struct senc3d_enclosure_header header; 1351 unsigned tid, med, enc_fst_med = SENC3D_UNSPECIFIED_MEDIUM; 1352 int is_fst_med = 1, is_err_cs = 0; 1353 1354 enc_status[e] = NO_ENCLOSURE_ERROR; 1355 ERR(senc3d_scene_get_enclosure(stardis->senc3d_scn, e, &enc)); 1356 ERR(senc3d_enclosure_get_header(enc, &header)); 1357 1358 FOR_EACH(t, 0, header.unique_primitives_count) { 1359 unsigned prop[SG3D_PROP_TYPES_COUNT__]; 1360 enum senc3d_side side; 1361 size_t j; 1362 ERR(senc3d_enclosure_get_triangle_id(enc, t, &tid, &side)); 1363 FOR_EACH(j, 0, darray_uint_size_get(&stardis->compute_surface.err_triangles)) 1364 { 1365 unsigned prim 1366 = darray_uint_cdata_get(&stardis->compute_surface.err_triangles)[j]; 1367 if(prim == tid) { 1368 is_err_cs = 1; 1369 break; 1370 } 1371 } 1372 if(is_err_cs) 1373 /* Don't flag an enclosure invalid because of a triangle that is 1374 * considered not member of it */ 1375 continue; 1376 1377 ERR(sg3d_geometry_get_unique_triangle_properties(stardis->geometry.sg3d, 1378 tid, prop)); 1379 1380 if(prop[side] != SG3D_UNSPECIFIED_PROPERTY) { 1381 ASSERT(prop[side] < darray_descriptions_size_get(&stardis->descriptions)); 1382 description_get_medium_id(descriptions + prop[side], &med); 1383 } else { 1384 /* If unspecified behind a boundary, use a specific ID to avoid to flag 1385 * using different boundaries for a given enclosure as invalid */ 1386 int properties_conflict_status; 1387 validate_properties(t, prop, stardis, &properties_conflict_status); 1388 if(properties_conflict_status == NO_PROPERTY_CONFLICT) 1389 med = stardis->undefined_medium_behind_boundary_id; 1390 else { 1391 if(!(enc_status[e] & ENCLOSURE_WITH_UNDEF_MEDIUM)) undef_count++; 1392 enc_status[e] |= ENCLOSURE_WITH_UNDEF_MEDIUM; 1393 continue; /* Don't flag N_MEDIA at the same time */ 1394 } 1395 } 1396 if(is_fst_med) { 1397 is_fst_med = 0; 1398 enc_fst_med = med; 1399 } else { 1400 if(enc_fst_med != med) { 1401 /* The external (infinite) enclosure can have multiple media */ 1402 if(!header.is_infinite && !(enc_status[e] & ENCLOSURE_WITH_N_MEDIA)) 1403 multi_count++; 1404 enc_status[e] |= ENCLOSURE_WITH_N_MEDIA; 1405 } 1406 } 1407 } 1408 ERR(senc3d_enclosure_ref_put(enc)); 1409 } 1410 if(multi_count) { 1411 int fst = 1; 1412 str_printf(&msg, 1413 "Found %d enclosure(s) with more than 1 medium:", multi_count); 1414 FOR_EACH(e, 0, ecount) { 1415 if(!(enc_status[e] & ENCLOSURE_WITH_N_MEDIA)) continue; 1416 str_append_printf(&msg, (fst ? " %u" : ", %u"), e); 1417 fst = 0; 1418 } 1419 logger_print(stardis->logger, LOG_OUTPUT, "%s.\n", str_cget(&msg)); 1420 } 1421 if(undef_count) { 1422 int fst = 1; 1423 str_printf(&msg, 1424 "Found %d enclosure(s) with undefined medium:", undef_count); 1425 FOR_EACH(e, 0, ecount) { 1426 if(!(enc_status[e] & ENCLOSURE_WITH_UNDEF_MEDIUM)) continue; 1427 str_append_printf(&msg, (fst ? " %u" : ", %u"), e); 1428 fst = 0; 1429 } 1430 logger_print(stardis->logger, LOG_OUTPUT, "%s.\n", str_cget(&msg)); 1431 } 1432 fprintf(stream, "FIELD EnclosuresData 2\n"); 1433 fprintf(stream, "Enclosures %d %d unsigned_char\n", ecount, tsz); 1434 FOR_EACH(t, 0, tsz) { 1435 unsigned encs[2], is_err_cs = 0; 1436 size_t j; 1437 FOR_EACH(j, 0, darray_uint_size_get(&stardis->compute_surface.err_triangles)) 1438 { 1439 unsigned prim 1440 = darray_uint_cdata_get(&stardis->compute_surface.err_triangles)[j]; 1441 if(prim == t) { 1442 is_err_cs = 1; 1443 break; 1444 } 1445 } 1446 if(is_err_cs) { 1447 /* Triangles in compute surface with error are considered member of no 1448 * enclosure */ 1449 FOR_EACH(e, 1, ecount) fprintf(stream, "0 "); 1450 fprintf(stream, "0\n"); 1451 continue; 1452 } 1453 ERR(senc3d_scene_get_triangle_enclosures(stardis->senc3d_scn, t, encs)); 1454 FOR_EACH(e, 0, ecount) { 1455 unsigned c = (e == encs[SENC3D_FRONT] || e == encs[SENC3D_BACK]) 1456 ? (unsigned char)enc_status[e] : 0; 1457 if(e == ecount - 1) 1458 fprintf(stream, "%u\n", c); 1459 else fprintf(stream, "%u ", c); 1460 } 1461 } 1462 1463 #define ENC_NOT_MEMBER SENC3D_UNSPECIFIED_MEDIUM 1464 #define ENC_MEMBER_2_DISTINT_MEDIA (ENC_NOT_MEMBER - 1) 1465 #define ENC_MEMBER_NO_MEDIUM (ENC_NOT_MEMBER - 2) 1466 1467 fprintf(stream, "Enclosures_internal_media %d %d unsigned_int\n", ecount, tsz); 1468 FOR_EACH(t, 0, tsz) { 1469 unsigned descr[SG3D_PROP_TYPES_COUNT__]; 1470 unsigned encs[2]; 1471 unsigned is_err_cs = 0; 1472 size_t j; 1473 ERR(senc3d_scene_get_triangle_enclosures(stardis->senc3d_scn, t, encs)); 1474 ERR(sg3d_geometry_get_unique_triangle_properties(stardis->geometry.sg3d, t, 1475 descr)); 1476 1477 /* Special value for triangles in compute surface with error */ 1478 FOR_EACH(j, 0, darray_uint_size_get(&stardis->compute_surface.err_triangles)) 1479 { 1480 unsigned prim 1481 = darray_uint_cdata_get(&stardis->compute_surface.err_triangles)[j]; 1482 if(prim == t) { 1483 is_err_cs = 1; 1484 break; 1485 } 1486 } 1487 if(is_err_cs) { 1488 /* Triangles in compute surface with error are considered member of no 1489 * enclosure */ 1490 FOR_EACH(e, 1, ecount) fprintf(stream, "%u ", ENC_NOT_MEMBER); 1491 fprintf(stream, "%u\n", ENC_NOT_MEMBER); 1492 continue; 1493 } 1494 FOR_EACH(e, 0, ecount) { 1495 unsigned mid; 1496 if(e == encs[SENC3D_FRONT] && e == encs[SENC3D_BACK]) { 1497 /* Both sides of this triangle are in enclosure #e */ 1498 unsigned fmid, bmid; 1499 if(descr[SG3D_FRONT] != SG3D_UNSPECIFIED_PROPERTY) 1500 description_get_medium_id(descriptions + descr[SG3D_FRONT], &fmid); 1501 else if(descr[SG3D_INTFACE] != SG3D_UNSPECIFIED_PROPERTY 1502 && DESC_IS_BOUNDARY(descriptions+descr[SG3D_INTFACE])) 1503 { 1504 description_get_medium_id(descriptions + descr[SG3D_INTFACE], &fmid); 1505 } 1506 else fmid = ENC_MEMBER_NO_MEDIUM; 1507 if(descr[SENC3D_BACK] != SG3D_UNSPECIFIED_PROPERTY) 1508 description_get_medium_id(descriptions + descr[SENC3D_BACK], &bmid); 1509 else if(descr[SG3D_INTFACE] != SG3D_UNSPECIFIED_PROPERTY 1510 && DESC_IS_BOUNDARY(descriptions+descr[SG3D_INTFACE])) 1511 { 1512 description_get_medium_id(descriptions + descr[SG3D_INTFACE], &bmid); 1513 } 1514 else bmid = ENC_MEMBER_NO_MEDIUM; 1515 mid = (fmid == bmid) ? fmid : ENC_MEMBER_2_DISTINT_MEDIA; 1516 } 1517 else if(e == encs[SENC3D_FRONT]) { 1518 /* Member of enclosure #e (front side only) */ 1519 if(descr[SG3D_FRONT] != SG3D_UNSPECIFIED_PROPERTY) 1520 description_get_medium_id(descriptions + descr[SG3D_FRONT], &mid); 1521 else if(descr[SG3D_INTFACE] != SG3D_UNSPECIFIED_PROPERTY 1522 && DESC_IS_BOUNDARY(descriptions+descr[SG3D_INTFACE])) 1523 { 1524 description_get_medium_id(descriptions + descr[SG3D_INTFACE], &mid); 1525 } 1526 else mid = ENC_MEMBER_NO_MEDIUM; 1527 } 1528 else if(e == encs[SENC3D_BACK]) { 1529 /* Member of enclosure #e (back side only) */ 1530 if(descr[SENC3D_BACK] != SG3D_UNSPECIFIED_PROPERTY) 1531 description_get_medium_id(descriptions + descr[SENC3D_BACK], &mid); 1532 else if(descr[SG3D_INTFACE] != SG3D_UNSPECIFIED_PROPERTY 1533 && DESC_IS_BOUNDARY(descriptions+descr[SG3D_INTFACE])) 1534 { 1535 description_get_medium_id(descriptions + descr[SG3D_INTFACE], &mid); 1536 } 1537 else mid = ENC_MEMBER_NO_MEDIUM; 1538 } else { 1539 /* Not member of enclosure #e */ 1540 mid = ENC_NOT_MEMBER; 1541 } 1542 if(e == ecount - 1) 1543 fprintf(stream, "%u\n", mid); 1544 else fprintf(stream, "%u ", mid); 1545 } 1546 } 1547 1548 #undef ENC_NOT_MEMBER 1549 #undef ENC_ERR_COMPUTE_SURFACE 1550 #undef ENC_MEMBER_2_DISTINT_MEDIA 1551 #undef ENC_MEMBER_NO_MEDIUM 1552 1553 exit: 1554 str_release(&msg); 1555 MEM_RM(stardis->allocator, trgs); 1556 MEM_RM(stardis->allocator, enc_status); 1557 return res; 1558 error: 1559 goto exit; 1560 } 1561 1562 res_T 1563 print_computation_time 1564 (struct sdis_mc* time_per_realisation, 1565 struct stardis* stardis, 1566 struct time* start, 1567 struct time* compute_start, 1568 struct time* compute_end, 1569 struct time* output_end) 1570 { 1571 struct time tmp; 1572 char buf[128]; 1573 const int flag = TIME_MSEC | TIME_SEC | TIME_MIN | TIME_HOUR | TIME_DAY; 1574 1575 ASSERT(stardis && start && compute_start && compute_end); 1576 1577 /* Only master prints or reads estimators */ 1578 ASSERT(!stardis->mpi_initialized || stardis->mpi_rank == 0); 1579 1580 time_sub(&tmp, compute_start, start); 1581 time_dump(&tmp, flag, NULL, buf, sizeof(buf)); 1582 logger_print(stardis->logger, LOG_OUTPUT, 1583 "Initialisation time = %s\n", buf); 1584 time_sub(&tmp, compute_end, compute_start); 1585 time_dump(&tmp, flag, NULL, buf, sizeof(buf)); 1586 logger_print(stardis->logger, LOG_OUTPUT, 1587 "Computation time = %s\n", buf); 1588 if(output_end) { 1589 time_sub(&tmp, output_end, compute_end); 1590 time_dump(&tmp, flag, NULL, buf, sizeof(buf)); 1591 logger_print(stardis->logger, LOG_OUTPUT, 1592 "Result output time = %s\n", buf); 1593 } 1594 1595 if(time_per_realisation) { 1596 logger_print(stardis->logger, LOG_OUTPUT, 1597 "Time per realisation (in usec) = %g +/- %g\n", 1598 time_per_realisation->E, time_per_realisation->SE); 1599 } 1600 1601 return RES_OK; 1602 } 1603 1604 res_T 1605 print_single_MC_result 1606 (struct sdis_estimator* estimator, 1607 struct stardis* stardis, 1608 FILE* stream) 1609 { 1610 res_T res = RES_OK; 1611 struct sdis_mc result; 1612 size_t nfailures_; 1613 unsigned long nfailures, nsamples; 1614 1615 ASSERT(estimator && stardis && stream); 1616 1617 /* Only master prints or reads estimators */ 1618 ASSERT(!stardis->mpi_initialized || stardis->mpi_rank == 0); 1619 1620 /* Fetch the estimation data */ 1621 ERR(sdis_estimator_get_temperature(estimator, &result)); 1622 ERR(sdis_estimator_get_failure_count(estimator, &nfailures_)); 1623 if(nfailures_ > ULONG_MAX || stardis->samples > ULONG_MAX) goto abort; 1624 nfailures = (unsigned long)nfailures_; 1625 nsamples = (unsigned long)stardis->samples; 1626 if(nfailures == nsamples) { 1627 logger_print(stardis->logger, LOG_ERROR, 1628 "All the %lu samples failed. No result to display.\n", nsamples); 1629 res = RES_BAD_OP; 1630 goto error; 1631 } 1632 1633 /* Print the results */ 1634 switch (stardis->mode & COMPUTE_MODES) { 1635 case MODE_COMPUTE_PROBE_TEMP_ON_VOL: 1636 if(stardis->mode & MODE_EXTENDED_RESULTS) { 1637 if(stardis->time_range[0] == stardis->time_range[1]) 1638 fprintf(stream, "Temperature at [%g, %g, %g] at t=%g = %g K +/- %g\n", 1639 SPLIT3(stardis->probe), stardis->time_range[0], 1640 result.E, /* Expected value */ 1641 result.SE); /* Standard error */ 1642 else 1643 fprintf(stream, 1644 "Temperature at [%g, %g, %g] with t in [%g %g] = %g K +/- %g\n", 1645 SPLIT3(stardis->probe), SPLIT2(stardis->time_range), 1646 result.E, /* Expected value */ 1647 result.SE); /* Standard error */ 1648 } 1649 else fprintf(stream, "%g %g %lu %lu\n", 1650 result.E, result.SE, nfailures, nsamples); 1651 break; 1652 case MODE_COMPUTE_PROBE_TEMP_ON_SURF: 1653 case MODE_COMPUTE_LIST_PROBE_TEMP_ON_SURF: 1654 { 1655 const struct stardis_probe_boundary* probe = NULL; 1656 probe = darray_probe_boundary_cdata_get(&stardis->probe_boundary_list); 1657 ERR(print_single_MC_result_probe_boundary 1658 (stardis, probe, estimator, stream)); 1659 break; 1660 } 1661 case MODE_COMPUTE_TEMP_MEAN_IN_MEDIUM: 1662 if(stardis->mode & MODE_EXTENDED_RESULTS) { 1663 if(stardis->time_range[0] == stardis->time_range[1]) 1664 fprintf(stream, "Temperature in medium '%s' at t=%g = %g K +/- %g\n", 1665 str_cget(&stardis->solve_name), stardis->time_range[0], 1666 result.E, /* Expected value */ 1667 result.SE); /* Standard error */ 1668 else 1669 fprintf(stream, 1670 "Temperature in medium '%s' with t in [%g %g] = %g K +/- %g\n", 1671 str_cget(&stardis->solve_name), SPLIT2(stardis->time_range), 1672 result.E, /* Expected value */ 1673 result.SE); /* Standard error */ 1674 } 1675 else fprintf(stream, "%g %g %lu %lu\n", 1676 result.E, result.SE, nfailures, nsamples); 1677 break; 1678 case MODE_COMPUTE_TEMP_MEAN_ON_SURF: 1679 if(stardis->mode & MODE_EXTENDED_RESULTS) { 1680 if(stardis->time_range[0] == stardis->time_range[1]) 1681 fprintf(stream, "Temperature at boundary '%s' at t=%g = %g K +/- %g\n", 1682 str_cget(&stardis->solve_name), stardis->time_range[0], 1683 result.E, /* Expected value */ 1684 result.SE); /* Standard error */ 1685 else 1686 fprintf(stream, 1687 "Temperature at boundary '%s' with t in [%g %g] = %g K +/- %g\n", 1688 str_cget(&stardis->solve_name), SPLIT2(stardis->time_range), 1689 result.E, /* Expected value */ 1690 result.SE); /* Standard error */ 1691 } 1692 else fprintf(stream, "%g %g %lu %lu\n", 1693 result.E, result.SE, nfailures, nsamples); 1694 break; 1695 1696 case MODE_COMPUTE_PROBE_FLUX_DNSTY_ON_SURF: 1697 case MODE_COMPUTE_LIST_PROBE_FLUX_DNSTY_ON_SURF: 1698 { 1699 enum sdis_estimator_type type; 1700 ERR(sdis_estimator_get_type(estimator, &type)); 1701 ASSERT(type == SDIS_ESTIMATOR_FLUX); 1702 1703 if(stardis->mode & MODE_EXTENDED_RESULTS) { 1704 if(stardis->time_range[0] == stardis->time_range[1]) 1705 fprintf(stream, "Temperature at [%g, %g, %g] at t=%g = %g K +/- %g\n", 1706 SPLIT3(stardis->probe), stardis->time_range[0], 1707 result.E, /* Expected value */ 1708 result.SE); /* Standard error */ 1709 else 1710 fprintf(stream, 1711 "Temperature at [%g, %g, %g] with t in [%g %g] = %g K +/- %g\n", 1712 SPLIT3(stardis->probe), SPLIT2(stardis->time_range), 1713 result.E, /* Expected value */ 1714 result.SE); /* Standard error */ 1715 ERR(sdis_estimator_get_convective_flux(estimator, &result)); 1716 if(stardis->time_range[0] == stardis->time_range[1]) 1717 fprintf(stream, "Convective flux density at [%g, %g, %g] at t=%g = %g W/m² +/- %g\n", 1718 SPLIT3(stardis->probe), stardis->time_range[0], 1719 result.E, /* Expected value */ 1720 result.SE); /* Standard error */ 1721 else 1722 fprintf(stream, 1723 "Convective flux density at [%g, %g, %g] with t in [%g %g] = %g W/m² +/- %g\n", 1724 SPLIT3(stardis->probe), SPLIT2(stardis->time_range), 1725 result.E, /* Expected value */ 1726 result.SE); /* Standard error */ 1727 ERR(sdis_estimator_get_radiative_flux(estimator, &result)); 1728 if(stardis->time_range[0] == stardis->time_range[1]) 1729 fprintf(stream, "Radiative flux density at [%g, %g, %g] at t=%g = %g W/m² +/- %g\n", 1730 SPLIT3(stardis->probe), stardis->time_range[0], 1731 result.E, /* Expected value */ 1732 result.SE); /* Standard error */ 1733 else 1734 fprintf(stream, 1735 "Radiative flux density at [%g, %g, %g] with t in [%g %g] = %g W/m² +/- %g\n", 1736 SPLIT3(stardis->probe), SPLIT2(stardis->time_range), 1737 result.E, /* Expected value */ 1738 result.SE); /* Standard error */ 1739 ERR(sdis_estimator_get_imposed_flux(estimator, &result)); 1740 if(stardis->time_range[0] == stardis->time_range[1]) 1741 fprintf(stream, "Imposed flux density at [%g, %g, %g] at t=%g = %g W/m² +/- %g\n", 1742 SPLIT3(stardis->probe), stardis->time_range[0], 1743 result.E, /* Expected value */ 1744 result.SE); /* Standard error */ 1745 else 1746 fprintf(stream, 1747 "Imposed flux density at [%g, %g, %g] with t in [%g %g] = %g W/m² +/- %g\n", 1748 SPLIT3(stardis->probe), SPLIT2(stardis->time_range), 1749 result.E, /* Expected value */ 1750 result.SE); /* Standard error */ 1751 ERR(sdis_estimator_get_total_flux(estimator, &result)); 1752 if(stardis->time_range[0] == stardis->time_range[1]) 1753 fprintf(stream, "Total flux density at [%g, %g, %g] at t=%g = %g W/m² +/- %g\n", 1754 SPLIT3(stardis->probe), stardis->time_range[0], 1755 result.E, /* Expected value */ 1756 result.SE); /* Standard error */ 1757 else 1758 fprintf(stream, 1759 "Total flux density at [%g, %g, %g] with t in [%g %g] = %g W/m² +/- %g\n", 1760 SPLIT3(stardis->probe), SPLIT2(stardis->time_range), 1761 result.E, /* Expected value */ 1762 result.SE); /* Standard error */ 1763 } else { 1764 fprintf(stream, "%g %g ", result.E, result.SE); /* Temperature */ 1765 ERR(sdis_estimator_get_convective_flux(estimator, &result)); 1766 fprintf(stream, "%g %g ", 1767 result.E, 1768 result.SE); 1769 ERR(sdis_estimator_get_radiative_flux(estimator, &result)); 1770 fprintf(stream, "%g %g ", 1771 result.E, 1772 result.SE); 1773 ERR(sdis_estimator_get_imposed_flux(estimator, &result)); 1774 fprintf(stream, "%g %g ", 1775 result.E, 1776 result.SE); 1777 ERR(sdis_estimator_get_total_flux(estimator, &result)); 1778 fprintf(stream, "%g %g ", 1779 result.E, 1780 result.SE); 1781 fprintf(stream, "%lu %lu\n", nfailures, nsamples); 1782 } 1783 break; 1784 } 1785 1786 case MODE_COMPUTE_FLUX_THROUGH_SURF: 1787 { 1788 enum sdis_estimator_type type; 1789 ERR(sdis_estimator_get_type(estimator, &type)); 1790 ASSERT(type == SDIS_ESTIMATOR_FLUX); 1791 1792 if(stardis->mode & MODE_EXTENDED_RESULTS) { 1793 if(stardis->time_range[0] == stardis->time_range[1]) 1794 fprintf(stream, "Temperature at boundary '%s' at t=%g = %g K +/- %g\n", 1795 str_cget(&stardis->solve_name), stardis->time_range[0], 1796 result.E, /* Expected value */ 1797 result.SE); /* Standard error */ 1798 else 1799 fprintf(stream, 1800 "Temperature at boundary '%s' with t in [%g %g] = %g K +/- %g\n", 1801 str_cget(&stardis->solve_name), SPLIT2(stardis->time_range), 1802 result.E, /* Expected value */ 1803 result.SE); /* Standard error */ 1804 ERR(sdis_estimator_get_convective_flux(estimator, &result)); 1805 if(stardis->time_range[0] == stardis->time_range[1]) 1806 fprintf(stream, "Convective flux at boundary '%s' at t=%g = %g W +/- %g\n", 1807 str_cget(&stardis->solve_name), stardis->time_range[0], 1808 stardis->compute_surface.area * result.E, /* Expected value */ 1809 stardis->compute_surface.area * result.SE); /* Standard error */ 1810 else 1811 fprintf(stream, 1812 "Convective flux at boundary '%s' with t in [%g %g] = %g W +/- %g\n", 1813 str_cget(&stardis->solve_name), SPLIT2(stardis->time_range), 1814 stardis->compute_surface.area * result.E, /* Expected value */ 1815 stardis->compute_surface.area * result.SE); /* Standard error */ 1816 ERR(sdis_estimator_get_radiative_flux(estimator, &result)); 1817 if(stardis->time_range[0] == stardis->time_range[1]) 1818 fprintf(stream, "Radiative flux at boundary '%s' at t=%g = %g W +/- %g\n", 1819 str_cget(&stardis->solve_name), stardis->time_range[0], 1820 stardis->compute_surface.area * result.E, /* Expected value */ 1821 stardis->compute_surface.area * result.SE); /* Standard error */ 1822 else 1823 fprintf(stream, 1824 "Radiative flux at boundary '%s' with t in [%g %g] = %g W +/- %g\n", 1825 str_cget(&stardis->solve_name), SPLIT2(stardis->time_range), 1826 stardis->compute_surface.area * result.E, /* Expected value */ 1827 stardis->compute_surface.area * result.SE); /* Standard error */ 1828 ERR(sdis_estimator_get_imposed_flux(estimator, &result)); 1829 if(stardis->time_range[0] == stardis->time_range[1]) 1830 fprintf(stream, "Imposed flux at boundary '%s' at t=%g = %g W +/- %g\n", 1831 str_cget(&stardis->solve_name), stardis->time_range[0], 1832 stardis->compute_surface.area * result.E, /* Expected value */ 1833 stardis->compute_surface.area * result.SE); /* Standard error */ 1834 else 1835 fprintf(stream, 1836 "Imposed flux at boundary '%s' with t in [%g %g] = %g W +/- %g\n", 1837 str_cget(&stardis->solve_name), SPLIT2(stardis->time_range), 1838 stardis->compute_surface.area * result.E, /* Expected value */ 1839 stardis->compute_surface.area * result.SE); /* Standard error */ 1840 ERR(sdis_estimator_get_total_flux(estimator, &result)); 1841 if(stardis->time_range[0] == stardis->time_range[1]) 1842 fprintf(stream, "Total flux at boundary '%s' at t=%g = %g W +/- %g\n", 1843 str_cget(&stardis->solve_name), stardis->time_range[0], 1844 stardis->compute_surface.area * result.E, /* Expected value */ 1845 stardis->compute_surface.area * result.SE); /* Standard error */ 1846 else 1847 fprintf(stream, 1848 "Total flux at boundary '%s' with t in [%g %g] = %g W +/- %g\n", 1849 str_cget(&stardis->solve_name), SPLIT2(stardis->time_range), 1850 stardis->compute_surface.area * result.E, /* Expected value */ 1851 stardis->compute_surface.area * result.SE); /* Standard error */ 1852 } else { 1853 fprintf(stream, "%g %g ", result.E, result.SE); /* Temperature */ 1854 ERR(sdis_estimator_get_convective_flux(estimator, &result)); 1855 fprintf(stream, "%g %g ", 1856 stardis->compute_surface.area * result.E, 1857 stardis->compute_surface.area * result.SE); 1858 ERR(sdis_estimator_get_radiative_flux(estimator, &result)); 1859 fprintf(stream, "%g %g ", 1860 stardis->compute_surface.area * result.E, 1861 stardis->compute_surface.area * result.SE); 1862 ERR(sdis_estimator_get_imposed_flux(estimator, &result)); 1863 fprintf(stream, "%g %g ", 1864 stardis->compute_surface.area * result.E, 1865 stardis->compute_surface.area * result.SE); 1866 ERR(sdis_estimator_get_total_flux(estimator, &result)); 1867 fprintf(stream, "%g %g ", 1868 stardis->compute_surface.area * result.E, 1869 stardis->compute_surface.area * result.SE); 1870 fprintf(stream, "%lu %lu\n", nfailures, nsamples); 1871 } 1872 break; 1873 } 1874 default: FATAL("Invalid mode\n."); 1875 } 1876 if(stardis->mode & MODE_EXTENDED_RESULTS) 1877 fprintf(stream, "#failures: %lu/%lu\n", nfailures, nsamples); 1878 if(nfailures) 1879 logger_print(stardis->logger, LOG_ERROR, 1880 "#failures: %lu/%lu\n", nfailures, nsamples); 1881 1882 end: 1883 return res; 1884 error: 1885 goto end; 1886 abort: 1887 res = RES_BAD_ARG; 1888 goto error; 1889 } 1890 1891 res_T 1892 print_single_MC_result_probe_boundary 1893 (struct stardis* stardis, 1894 const struct stardis_probe_boundary* probe, 1895 const struct sdis_estimator* estimator, 1896 FILE* stream) 1897 { 1898 struct sdis_mc result = SDIS_MC_NULL; 1899 size_t nfailures = 0; 1900 size_t nsamples = 0; 1901 res_T res = RES_OK; 1902 1903 ASSERT(stardis && probe && estimator); 1904 ASSERT((stardis->mode & MODE_COMPUTE_PROBE_TEMP_ON_SURF) 1905 || (stardis->mode & MODE_COMPUTE_LIST_PROBE_TEMP_ON_SURF)); 1906 1907 /* Only master prints or reads estimators */ 1908 ASSERT(!stardis->mpi_initialized || stardis->mpi_rank == 0); 1909 1910 /* Fetch the estimation data */ 1911 ERR(sdis_estimator_get_temperature(estimator, &result)); 1912 ERR(sdis_estimator_get_failure_count(estimator, &nfailures)); 1913 nsamples = stardis->samples; 1914 1915 if(nfailures == nsamples) { 1916 logger_print(stardis->logger, LOG_ERROR, 1917 "All the %lu samples failed. No result to display.\n", nsamples); 1918 res = RES_BAD_OP; 1919 goto error; 1920 } 1921 1922 /* Raw output */ 1923 if((stardis->mode & MODE_EXTENDED_RESULTS) == 0) { 1924 fprintf(stream, "%g %g %lu %lu\n", 1925 result.E, result.SE, (unsigned long)nfailures, (unsigned long)nsamples); 1926 1927 /* Extended output */ 1928 } else if(stardis->time_range[0] == stardis->time_range[1]) { 1929 fprintf(stream, 1930 "Boundary temperature at [%g, %g, %g] at t=%g = %g K +/- %g\n", 1931 SPLIT3(probe->position), probe->time[0], result.E, result.SE); 1932 1933 /* Extended output with time range */ 1934 } else { 1935 fprintf(stream, 1936 "Boundary temperature at [%g, %g, %g] with t in [%g %g] = %g K +/- %g\n", 1937 SPLIT3(probe->position), SPLIT2(probe->time), result.E, result.SE); 1938 } 1939 1940 exit: 1941 return res; 1942 error: 1943 goto exit; 1944 } 1945 1946 res_T 1947 dump_map 1948 (const struct stardis* stardis, 1949 const struct darray_estimators* estimators, 1950 FILE* stream) 1951 { 1952 res_T res = RES_OK; 1953 unsigned i, vcount, tcount, last_v = 0; 1954 const size_t* idx; 1955 size_t sz; 1956 unsigned szp; 1957 struct sdis_estimator* const* est; 1958 1959 ASSERT(stardis && estimators && stream); 1960 1961 /* Only master prints or reads estimators */ 1962 ASSERT(!stardis->mpi_initialized || stardis->mpi_rank == 0); 1963 1964 est = darray_estimators_cdata_get(estimators); 1965 idx = darray_size_t_cdata_get(&stardis->compute_surface.primitives); 1966 sz = darray_size_t_size_get(&stardis->compute_surface.primitives); 1967 if(sz > UINT_MAX) goto abort; 1968 szp = (unsigned)sz; 1969 SG3D(geometry_get_unique_vertices_count(stardis->geometry.sg3d, &vcount)); 1970 SG3D(geometry_get_unique_triangles_count(stardis->geometry.sg3d, &tcount)); 1971 1972 /* Find last used vertex */ 1973 for(i = 0; i < szp; ++i) { 1974 unsigned t; 1975 unsigned indices[3]; 1976 if(idx[i] > UINT_MAX) goto abort; 1977 t = (unsigned)idx[i]; 1978 SG3D(geometry_get_unique_triangle_vertices(stardis->geometry.sg3d, t, 1979 indices)); 1980 last_v = MMAX(MMAX(last_v, indices[0]), MMAX(indices[1], indices[2])); 1981 } 1982 1983 /* Dump vertices up to last_v, even unused ones, to avoid reindexing */ 1984 fprintf(stream, "# vtk DataFile Version 2.0\n"); 1985 fprintf(stream, "Temperature Map\n"); 1986 fprintf(stream, "ASCII\n"); 1987 fprintf(stream, "DATASET POLYDATA\n"); 1988 fprintf(stream, "POINTS %u float\n\n", last_v + 1); 1989 for(i = 0; i <= last_v; ++i) { 1990 double coord[3]; 1991 SG3D(geometry_get_unique_vertex(stardis->geometry.sg3d, i, coord)); 1992 fprintf(stream, "%f %f %f\n", SPLIT3(coord)); 1993 } 1994 /* Dump only primitives in boundary */ 1995 fprintf(stream, "\nPOLYGONS %u %u\n", szp, 4 * szp); 1996 for(i = 0; i < szp; ++i) { 1997 unsigned t; 1998 unsigned indices[3]; 1999 if(idx[i] > UINT_MAX) goto abort; 2000 t = (unsigned)idx[i]; 2001 SG3D(geometry_get_unique_triangle_vertices(stardis->geometry.sg3d, t, 2002 indices)); 2003 fprintf(stream, "3 %u %u %u\n", SPLIT3(indices)); 2004 } 2005 fprintf(stream, "\nCELL_DATA %u\n", szp); 2006 fprintf(stream, "SCALARS temperature_estimate float 1\n"); 2007 fprintf(stream, "LOOKUP_TABLE default\n"); 2008 for(i = 0; i < szp; ++i) { 2009 struct sdis_mc T; 2010 SDIS(estimator_get_temperature(est[i], &T)); 2011 fprintf(stream, "%f\n", T.E); 2012 } 2013 fprintf(stream, "SCALARS temperature_std_dev float 1\n"); 2014 fprintf(stream, "LOOKUP_TABLE default\n"); 2015 for(i = 0; i < szp; ++i) { 2016 struct sdis_mc T; 2017 SDIS(estimator_get_temperature(est[i], &T)); 2018 fprintf(stream, "%f\n", T.SE); 2019 } 2020 fprintf(stream, "SCALARS failures_count unsigned_long_long 1\n"); 2021 fprintf(stream, "LOOKUP_TABLE default\n"); 2022 for(i = 0; i < szp; ++i) { 2023 size_t nfails; 2024 SDIS(estimator_get_failure_count(est[i], &nfails)); 2025 if(nfails > UINT_MAX) goto abort; 2026 fprintf(stream, "%u\n", (unsigned)nfails); 2027 } 2028 fprintf(stream, "SCALARS computation_time_estimate float 1\n"); 2029 fprintf(stream, "LOOKUP_TABLE default\n"); 2030 for(i = 0; i < szp; ++i) { 2031 struct sdis_mc time; 2032 SDIS(estimator_get_realisation_time(est[i], &time)); 2033 fprintf(stream, "%f\n", time.E); 2034 } 2035 fprintf(stream, "SCALARS computation_time_std_dev float 1\n"); 2036 fprintf(stream, "LOOKUP_TABLE default\n"); 2037 for(i = 0; i < szp; ++i) { 2038 struct sdis_mc time; 2039 SDIS(estimator_get_realisation_time(est[i], &time)); 2040 fprintf(stream, "%f\n", time.SE); 2041 } 2042 end: 2043 return res; 2044 error: 2045 goto end; 2046 abort: 2047 res = RES_BAD_ARG; 2048 goto error; 2049 } 2050 2051 res_T 2052 dump_compute_region_at_the_end_of_vtk 2053 (struct stardis* stardis, 2054 FILE* stream) 2055 { 2056 res_T res = RES_OK; 2057 unsigned char* v = NULL; 2058 unsigned tsz, i; 2059 size_t j, psz; 2060 ASSERT(stardis && stream); 2061 2062 /* Only master prints or reads estimators */ 2063 ASSERT(!stardis->mpi_initialized || stardis->mpi_rank == 0); 2064 2065 psz = darray_size_t_size_get(&stardis->compute_surface.primitives); 2066 ASSERT(psz == darray_sides_size_get(&stardis->compute_surface.sides)); 2067 2068 ERR(sg3d_geometry_get_unique_triangles_count(stardis->geometry.sg3d, &tsz)); 2069 /* For triangles not in compute region v==0 */ 2070 v = MEM_CALLOC(stardis->allocator, tsz, sizeof(*v)); 2071 if(!v) { 2072 res = RES_MEM_ERR; 2073 goto error; 2074 } 2075 2076 if(stardis->mode & SURFACE_COMPUTE_MODES) { 2077 /* For triangles in compute surface, v==1 if FRONT or v==2 for BACK */ 2078 FOR_EACH(j, 0, psz) { 2079 size_t prim 2080 = darray_size_t_cdata_get(&stardis->compute_surface.primitives)[j]; 2081 enum sdis_side side 2082 = darray_sides_cdata_get(&stardis->compute_surface.sides)[j]; 2083 ASSERT(prim <= tsz); 2084 v[(unsigned)prim] = 2085 (unsigned char)(v[(unsigned)prim] | (side == SDIS_FRONT ? 1 : 2)); 2086 } 2087 2088 /* For triangles in compute surface with error v==MAX */ 2089 FOR_EACH(j, 0, darray_uint_size_get(&stardis->compute_surface.err_triangles)) 2090 { 2091 unsigned prim 2092 = darray_uint_cdata_get(&stardis->compute_surface.err_triangles)[j]; 2093 ASSERT(prim <= tsz); 2094 v[(unsigned)prim] = UCHAR_MAX; 2095 } 2096 } else { 2097 unsigned descr[SG3D_PROP_TYPES_COUNT__]; 2098 struct sdis_medium* medium; 2099 const struct description* descriptions; 2100 unsigned medium_id; 2101 ASSERT(stardis->mode & MODE_COMPUTE_TEMP_MEAN_IN_MEDIUM); 2102 medium = find_medium_by_name 2103 (stardis, str_cget(&stardis->solve_name), &medium_id); 2104 ASSERT(medium != NULL); (void)medium; 2105 descriptions = darray_descriptions_cdata_get(&stardis->descriptions); 2106 FOR_EACH(i, 0, tsz) { 2107 unsigned f_mid, b_mid; 2108 /* Get the description IDs for this triangle */ 2109 ERR(sg3d_geometry_get_unique_triangle_properties(stardis->geometry.sg3d, 2110 i, descr)); 2111 /* For triangles in compute volume, 2112 * v==1 if FRONT, v==2 for BACK */ 2113 if(descr[SG3D_FRONT] == SG3D_UNSPECIFIED_PROPERTY) 2114 f_mid = UINT_MAX; 2115 else description_get_medium_id(descriptions + descr[SG3D_FRONT], &f_mid); 2116 if(descr[SG3D_BACK] == SG3D_UNSPECIFIED_PROPERTY) 2117 b_mid = UINT_MAX; 2118 else description_get_medium_id(descriptions + descr[SG3D_BACK], &b_mid); 2119 if(f_mid == medium_id && b_mid == medium_id) 2120 ; /* Keep v==0, not really a boundary */ 2121 else if(f_mid == medium_id) 2122 v[i] = 1; 2123 else if(b_mid == medium_id) 2124 v[i] = 2; 2125 } 2126 } 2127 2128 fprintf(stream, "SCALARS Compute_region unsigned_int 1\n"); 2129 fprintf(stream, "LOOKUP_TABLE default\n"); 2130 FOR_EACH(i, 0, tsz) 2131 fprintf(stream, "%u\n", v[i] == UCHAR_MAX ? UINT_MAX : v[i]); 2132 2133 exit: 2134 MEM_RM(stardis->allocator, v); 2135 return res; 2136 error: 2137 goto exit; 2138 } 2139 2140 res_T 2141 dump_model_as_c_chunks 2142 (struct stardis* stardis, 2143 FILE* stream) 2144 { 2145 res_T res = RES_OK; 2146 const char* prefix; 2147 unsigned n, vcount, tcount; 2148 2149 ASSERT(stardis && stream); 2150 2151 /* Only master prints or reads estimators */ 2152 ASSERT(!stardis->mpi_initialized || stardis->mpi_rank == 0); 2153 2154 prefix = str_cget(&stardis->chunks_prefix); 2155 ERR(sg3d_geometry_get_unique_vertices_count(stardis->geometry.sg3d, &vcount)); 2156 ERR(sg3d_geometry_get_unique_triangles_count(stardis->geometry.sg3d, &tcount)); 2157 2158 fprintf(stream, "#define %s_UNSPECIFIED_PROPERTY %u\n\n", 2159 prefix, SG3D_UNSPECIFIED_PROPERTY); 2160 2161 fprintf(stream, "static const unsigned\n"); 2162 fprintf(stream, "%s_vertices_count = %u;\n\n", prefix, vcount); 2163 2164 fprintf(stream, "static const unsigned\n"); 2165 fprintf(stream, "%s_triangles_count = %u;\n\n", prefix, tcount); 2166 2167 fprintf(stream, "static const double\n"); 2168 fprintf(stream, "%s_vertices[%u][3] = {\n", prefix, vcount); 2169 for(n = 0; n < vcount; n++) { 2170 double vertex[3]; 2171 ERR(sg3d_geometry_get_unique_vertex(stardis->geometry.sg3d, n, 2172 vertex)); 2173 fprintf(stream, " { %g, %g, %g }%c\n", 2174 SPLIT3(vertex), (n == vcount - 1 ? ' ' : ',')); 2175 } 2176 fprintf(stream, "};\n\n"); 2177 2178 fprintf(stream, "static const unsigned\n"); 2179 fprintf(stream, "%s_triangles[%u][3] = {\n", prefix, tcount); 2180 for(n = 0; n < tcount; n++) { 2181 unsigned triangle[3]; 2182 ERR(sg3d_geometry_get_unique_triangle_vertices(stardis->geometry.sg3d, n, 2183 triangle)); 2184 fprintf(stream, " { %u, %u, %u }%c\n", 2185 SPLIT3(triangle), (n == tcount - 1 ? ' ' : ',')); 2186 } 2187 fprintf(stream, "};\n\n"); 2188 2189 fprintf(stream, "static const unsigned\n"); 2190 fprintf(stream, "%s_properties[%u][3] = {\n", prefix, tcount); 2191 for(n = 0; n < tcount; n++) { 2192 unsigned properties[SG3D_PROP_TYPES_COUNT__]; 2193 ERR(sg3d_geometry_get_unique_triangle_properties(stardis->geometry.sg3d, n, 2194 properties)); 2195 if(properties[0] == SG3D_UNSPECIFIED_PROPERTY) 2196 fprintf(stream, " { %s_UNSPECIFIED_PROPERTY, ", prefix); 2197 else fprintf(stream, " { %u, ", properties[0]); 2198 if(properties[1] == SG3D_UNSPECIFIED_PROPERTY) 2199 fprintf(stream, "%s_UNSPECIFIED_PROPERTY, ", prefix); 2200 else fprintf(stream, "%u, ", properties[1]); 2201 if(properties[2] == SG3D_UNSPECIFIED_PROPERTY) 2202 fprintf(stream, "%s_UNSPECIFIED_PROPERTY }", prefix); 2203 else fprintf(stream, "%u }", properties[2]); 2204 if(n == tcount - 1) fprintf(stream, "\n"); else fprintf(stream, ",\n"); 2205 } 2206 fprintf(stream, "};\n\n"); 2207 2208 exit: 2209 return res; 2210 error: 2211 goto exit; 2212 } 2213 2214 res_T 2215 write_random_generator_state 2216 (struct sdis_estimator* estimator, 2217 FILE* stream) 2218 { 2219 res_T res; 2220 struct ssp_rng* state; 2221 res = sdis_estimator_get_rng_state(estimator, &state); 2222 if(res != RES_OK) return res; 2223 return ssp_rng_write(state, stream); 2224 } 2225 2226 res_T 2227 read_random_generator_state 2228 (struct ssp_rng* state, 2229 FILE* stream) 2230 { 2231 return ssp_rng_read(state, stream); 2232 }