test_sdis_solve_camera.c (27908B)
1 /* Copyright (C) 2016-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 #include "sdis.h" 17 #include "test_sdis_utils.h" 18 19 #include <star/ssp.h> 20 #include <star/s3dut.h> 21 22 #include <rsys/algorithm.h> 23 #include <rsys/image.h> 24 #include <rsys/double3.h> 25 #include <rsys/double33.h> 26 #include <rsys/math.h> 27 #include <rsys/stretchy_array.h> 28 29 #include <string.h> 30 31 #define IMG_WIDTH 157 32 #define IMG_HEIGHT 53 33 #define SPP 32 /* #Samples per pixel, i.e. #realisations per pixel */ 34 35 /* 36 * The scene consists of a solid cube whose temperature is unknown. The 37 * emissivity of the cube is 1 and the convection coefficient with the 38 * surrounding fluid at 300 K is 0.1. At the center of the cube is a spherical 39 * cavity of fluid with a temperature of 350 K. The convection coefficient 40 * between the solid and the cavity is 1, and the emissivity of this interface 41 * is zero. The ambient radiative temperature of the system is 300 K. 42 * 43 * Finally, a parallelepiped below the cube symbolizes the ground. The 44 * temperature of its Robin condition is 280 K. This geometry verifies that a 45 * camera can draw a scene in an enclosure containing several media, such as 46 * those used to define several boundary conditions. 47 * 48 * In this test, we calculate the radiative temperature that reaches a camera 49 * looking at the cube and produce an image of the result written to the 50 * standard output in htrdr-image(5) format. 51 * 52 * 53 * +----------------+ 54 * /' # # /| 55 * +----*--------*--+ | __\ 300 K 56 * | ' # # | | / / 57 * | ' # 350K # | | \__/ 58 * | ' # # | | 59 * | +.....#..#.....|.+ 60 * |/ |/ 61 * +----------------+ 280K 62 * +---------------------------------__\------+ 63 * / / / /| 64 * / \__/ / + 65 * +------------------------------------------+ / 66 * | |/ 67 * +------------------------------------------+ 68 */ 69 70 /******************************************************************************* 71 * Geometry 72 ******************************************************************************/ 73 struct map_interf { 74 size_t key; 75 struct sdis_interface* interf; 76 }; 77 78 static int 79 cmp_map_inter(const void* key, const void* elmt) 80 { 81 const size_t* iprim = key; 82 const struct map_interf* interf = elmt; 83 if(*iprim < interf->key) return -1; 84 else if(*iprim > interf->key) return 1; 85 else return 0; 86 } 87 88 struct geometry { 89 double* positions; 90 size_t* indices; 91 struct map_interf* interfaces; 92 }; 93 static const struct geometry GEOMETRY_NULL = {NULL, NULL, NULL}; 94 95 static void 96 geometry_add_shape 97 (struct geometry* geom, 98 const double* positions, 99 const size_t nverts, 100 const size_t* indices, 101 const size_t nprims, 102 const double transform[12], /* May be NULL <=> no transformation */ 103 struct sdis_interface* interf) 104 { 105 struct map_interf* geom_interf = NULL; 106 size_t nverts_prev = 0; 107 size_t i; 108 109 CHK(geom != NULL); 110 CHK(positions != NULL); 111 CHK(indices != NULL); 112 CHK(nverts != 0); 113 CHK(nprims != 0); 114 CHK(interf != NULL); 115 116 /* Save the previous number of vertices/primitives of the geometry */ 117 nverts_prev = sa_size(geom->positions) / 3; 118 119 /* Add the vertices */ 120 FOR_EACH(i, 0, nverts) { 121 double* pos = sa_add(geom->positions, 3); 122 d3_set(pos, positions + i*3); 123 if(transform) { 124 d33_muld3(pos, transform, pos); 125 d3_add(pos, transform+9, pos); 126 } 127 } 128 129 /* Add the indices */ 130 FOR_EACH(i, 0, nprims) { 131 sa_push(geom->indices, indices[i*3+0] + nverts_prev); 132 sa_push(geom->indices, indices[i*3+1] + nverts_prev); 133 sa_push(geom->indices, indices[i*3+2] + nverts_prev); 134 } 135 136 geom_interf = sa_add(geom->interfaces, 1); 137 geom_interf->key = sa_size(geom->indices) / 3 - 1; 138 geom_interf->interf = interf; 139 } 140 141 static void 142 geometry_release(struct geometry* geom) 143 { 144 CHK(geom != NULL); 145 sa_release(geom->positions); 146 sa_release(geom->indices); 147 sa_release(geom->interfaces); 148 *geom = GEOMETRY_NULL; 149 } 150 151 static void 152 geometry_get_indices(const size_t itri, size_t ids[3], void* ctx) 153 { 154 struct geometry* geom = ctx; 155 156 CHK(ids != NULL); 157 CHK(geom != NULL); 158 CHK(itri < sa_size(geom->indices)/3/*#indices per triangle*/); 159 160 /* Fetch the indices */ 161 ids[0] = geom->indices[itri*3+0]; 162 ids[1] = geom->indices[itri*3+1]; 163 ids[2] = geom->indices[itri*3+2]; 164 } 165 166 static void 167 geometry_get_position(const size_t ivert, double pos[3], void* ctx) 168 { 169 struct geometry* geom = ctx; 170 171 CHK(pos != NULL); 172 CHK(geom != NULL); 173 CHK(ivert < sa_size(geom->positions)/3/*#coords per triangle*/); 174 175 /* Fetch the vertices */ 176 pos[0] = geom->positions[ivert*3+0]; 177 pos[1] = geom->positions[ivert*3+1]; 178 pos[2] = geom->positions[ivert*3+2]; 179 } 180 181 static void 182 geometry_get_interface 183 (const size_t itri, 184 struct sdis_interface** bound, 185 void* ctx) 186 { 187 struct geometry* geom = ctx; 188 struct map_interf* interf = NULL; 189 190 CHK(bound != NULL); 191 CHK(geom != NULL); 192 CHK(itri < sa_size(geom->indices)/3/*#indices per triangle*/); 193 194 /* Find the interface of the triangle */ 195 interf = search_lower_bound(&itri, geom->interfaces, 196 sa_size(geom->interfaces), sizeof(struct map_interf), cmp_map_inter); 197 198 CHK(interf != NULL); 199 CHK(interf->key >= itri); 200 *bound = interf->interf; 201 } 202 203 static void 204 add_cube(struct geometry* geom, struct sdis_interface* interf) 205 { 206 struct s3dut_mesh_data msh_data; 207 struct s3dut_mesh* msh = NULL; 208 CHK(geom); 209 210 OK(s3dut_create_cuboid(NULL, 2, 2, 2, &msh)); 211 OK(s3dut_mesh_get_data(msh, &msh_data)); 212 geometry_add_shape(geom, msh_data.positions, msh_data.nvertices, 213 msh_data.indices, msh_data.nprimitives, NULL, interf); 214 OK(s3dut_mesh_ref_put(msh)); 215 } 216 217 static void 218 add_sphere(struct geometry* geom, struct sdis_interface* interf) 219 { 220 struct s3dut_mesh_data msh_data; 221 struct s3dut_mesh* msh = NULL; 222 CHK(geom); 223 224 OK(s3dut_create_sphere(NULL, 0.5, 32, 16, &msh)); 225 OK(s3dut_mesh_get_data(msh, &msh_data)); 226 geometry_add_shape(geom, msh_data.positions, msh_data.nvertices, 227 msh_data.indices, msh_data.nprimitives, NULL, interf); 228 OK(s3dut_mesh_ref_put(msh)); 229 } 230 231 static void 232 add_ground(struct geometry* geom, struct sdis_interface* interf) 233 { 234 struct s3dut_mesh_data msh_data; 235 struct s3dut_mesh* msh = NULL; 236 const double transform[12] = {1,0,0, 0,1,0, 0,0,1, 0,0,-4}; 237 CHK(geom); 238 239 OK(s3dut_create_cuboid(NULL, 10, 10, 2, &msh)); 240 OK(s3dut_mesh_get_data(msh, &msh_data)); 241 geometry_add_shape(geom, msh_data.positions, msh_data.nvertices, 242 msh_data.indices, msh_data.nprimitives, transform, interf); 243 OK(s3dut_mesh_ref_put(msh)); 244 } 245 246 /******************************************************************************* 247 * Media 248 ******************************************************************************/ 249 struct fluid { 250 double cp; 251 double rho; 252 double temperature; 253 }; 254 static const struct fluid FLUID_NULL = {0, 0, SDIS_TEMPERATURE_NONE}; 255 256 struct solid { 257 double cp; 258 double lambda; 259 double rho; 260 double delta; 261 double temperature; 262 }; 263 static const struct solid SOLID_NULL = {0, 0, 0, 0, SDIS_TEMPERATURE_NONE}; 264 265 #define MEDIUM_GETTER(Type, Prop) \ 266 static double \ 267 Type##_get_##Prop \ 268 (const struct sdis_rwalk_vertex* vtx, \ 269 struct sdis_data* data) \ 270 { \ 271 CHK(data && vtx); \ 272 return ((const struct Type*)sdis_data_cget(data))->Prop; \ 273 } 274 /* Fluid getters */ 275 MEDIUM_GETTER(fluid, cp) 276 MEDIUM_GETTER(fluid, rho) 277 MEDIUM_GETTER(fluid, temperature) 278 /* Solid getters */ 279 MEDIUM_GETTER(solid, cp) 280 MEDIUM_GETTER(solid, lambda) 281 MEDIUM_GETTER(solid, rho) 282 MEDIUM_GETTER(solid, delta) 283 MEDIUM_GETTER(solid, temperature) 284 #undef MEDIUM_GETTER 285 286 static struct sdis_medium* 287 create_fluid 288 (struct sdis_device* dev, 289 const struct fluid* param) 290 { 291 struct sdis_fluid_shader fluid_shader = DUMMY_FLUID_SHADER; 292 struct sdis_data* data = NULL; 293 struct sdis_medium* fluid = NULL; 294 295 CHK(param != NULL); 296 297 /* Copy the fluid parameters into the Stardis memory space */ 298 OK(sdis_data_create 299 (dev, sizeof(struct fluid), ALIGNOF(struct fluid), NULL, &data)); 300 memcpy(sdis_data_get(data), param, sizeof(struct fluid)); 301 302 /* Setup the fluid shader */ 303 fluid_shader.calorific_capacity = fluid_get_cp; 304 fluid_shader.volumic_mass = fluid_get_rho; 305 fluid_shader.temperature = fluid_get_temperature; 306 307 /* Create the fluid medium */ 308 OK(sdis_fluid_create(dev, &fluid_shader, data, &fluid)); 309 OK(sdis_data_ref_put(data)); 310 311 return fluid; 312 } 313 314 static struct sdis_medium* 315 create_solid 316 (struct sdis_device* dev, 317 const struct solid* param) 318 { 319 struct sdis_solid_shader solid_shader = DUMMY_SOLID_SHADER; 320 struct sdis_data* data = NULL; 321 struct sdis_medium* solid = NULL; 322 323 CHK(param != NULL); 324 325 /* Copy the solid parameters into the Stardis memory space */ 326 OK(sdis_data_create 327 (dev, sizeof(struct solid), ALIGNOF(struct solid), NULL, &data)); 328 memcpy(sdis_data_get(data), param, sizeof(struct solid)); 329 330 /* Setup the solid shader */ 331 solid_shader.calorific_capacity = solid_get_cp; 332 solid_shader.thermal_conductivity = solid_get_lambda; 333 solid_shader.volumic_mass = solid_get_rho; 334 solid_shader.delta = solid_get_delta; 335 solid_shader.temperature = solid_get_temperature; 336 337 /* Create the solid medium */ 338 OK(sdis_solid_create(dev, &solid_shader, data, &solid)); 339 OK(sdis_data_ref_put(data)); 340 341 return solid; 342 } 343 344 /******************************************************************************* 345 * Interfaces 346 ******************************************************************************/ 347 struct interf { 348 double hc; 349 double epsilon; 350 double specular_fraction; 351 double temperature; 352 double reference_temperature; 353 }; 354 static const struct interf INTERF_NULL = { 355 0, 0, 0, SDIS_TEMPERATURE_NONE, SDIS_TEMPERATURE_NONE 356 }; 357 358 #define INTERFACE_GETTER(Prop) \ 359 static double \ 360 interface_get_##Prop \ 361 (const struct sdis_interface_fragment* frag, \ 362 struct sdis_data* data) \ 363 { \ 364 CHK(data && frag); \ 365 return ((const struct interf*)sdis_data_cget(data))->Prop; \ 366 } 367 INTERFACE_GETTER(hc) 368 INTERFACE_GETTER(temperature) 369 INTERFACE_GETTER(reference_temperature) 370 #undef INTERFACE_GETTER 371 372 #define INTERFACE_GETTER(Prop) \ 373 static double \ 374 interface_get_##Prop \ 375 (const struct sdis_interface_fragment* frag, \ 376 const unsigned source_id, \ 377 struct sdis_data* data) \ 378 { \ 379 CHK(data && frag); \ 380 (void)source_id; \ 381 return ((const struct interf*)sdis_data_cget(data))->Prop; \ 382 } 383 INTERFACE_GETTER(epsilon) 384 INTERFACE_GETTER(specular_fraction) 385 #undef INTERFACE_GETTER 386 387 static struct sdis_interface* 388 create_interface 389 (struct sdis_device* dev, 390 struct sdis_medium* mdm_front, 391 struct sdis_medium* mdm_back, 392 const struct interf* param) 393 { 394 struct sdis_interface_shader interface_shader = SDIS_INTERFACE_SHADER_NULL; 395 struct sdis_data* data = NULL; 396 struct sdis_interface* interf = NULL; 397 398 CHK(mdm_front != NULL); 399 CHK(mdm_back != NULL); 400 401 /* Copy the interface parameters into the Stardis memory space */ 402 OK(sdis_data_create 403 (dev, sizeof(struct interf), ALIGNOF(struct interf), NULL, &data)); 404 memcpy(sdis_data_get(data), param, sizeof(struct interf)); 405 406 /* Setup the interface shader */ 407 interface_shader.convection_coef = interface_get_hc; 408 interface_shader.front.temperature = interface_get_temperature; 409 interface_shader.back.temperature = interface_get_temperature; 410 if(sdis_medium_get_type(mdm_front) == SDIS_FLUID) { 411 interface_shader.front.emissivity = interface_get_epsilon; 412 interface_shader.front.specular_fraction = interface_get_specular_fraction; 413 interface_shader.front.reference_temperature = 414 interface_get_reference_temperature; 415 } 416 if(sdis_medium_get_type(mdm_back) == SDIS_FLUID) { 417 interface_shader.back.emissivity = interface_get_epsilon; 418 interface_shader.back.specular_fraction = interface_get_specular_fraction; 419 interface_shader.back.reference_temperature = 420 interface_get_reference_temperature; 421 } 422 /* Create the interface */ 423 OK(sdis_interface_create 424 (dev, mdm_front, mdm_back, &interface_shader, data, &interf)); 425 OK(sdis_data_ref_put(data)); 426 427 return interf; 428 } 429 430 /******************************************************************************* 431 * Radiative environment 432 ******************************************************************************/ 433 struct radenv { 434 double temperature; /* [K] */ 435 double reference; /* [K] */ 436 }; 437 438 static double 439 radenv_get_temperature 440 (const struct sdis_radiative_ray* ray, 441 struct sdis_data* data) 442 { 443 (void)ray; 444 return ((const struct radenv*)sdis_data_cget(data))->temperature; 445 } 446 447 static double 448 radenv_get_reference_temperature 449 (const struct sdis_radiative_ray* ray, 450 struct sdis_data* data) 451 { 452 (void)ray; 453 return ((const struct radenv*)sdis_data_cget(data))->reference; 454 } 455 456 static struct sdis_radiative_env* 457 create_radenv 458 (struct sdis_device* dev, 459 const double temperature, 460 const double reference) 461 { 462 struct sdis_radiative_env_shader shader = SDIS_RADIATIVE_ENV_SHADER_NULL; 463 struct sdis_radiative_env* radenv = NULL; 464 struct sdis_data* data = NULL; 465 struct radenv* env = NULL; 466 467 OK(sdis_data_create(dev, sizeof(struct radenv), ALIGNOF(radenv), NULL, &data)); 468 env = sdis_data_get(data); 469 env->temperature = temperature; 470 env->reference = reference; 471 472 shader.temperature = radenv_get_temperature; 473 shader.reference_temperature = radenv_get_reference_temperature; 474 OK(sdis_radiative_env_create(dev, &shader, data, &radenv)); 475 OK(sdis_data_ref_put(data)); 476 return radenv; 477 } 478 479 /******************************************************************************* 480 * Scene 481 ******************************************************************************/ 482 static struct sdis_scene* 483 create_scene 484 (struct sdis_device* dev, 485 struct sdis_radiative_env* radenv, 486 struct geometry* geom) 487 { 488 struct sdis_scene_create_args scn_args = SDIS_SCENE_CREATE_ARGS_DEFAULT; 489 struct sdis_scene* scn = NULL; 490 size_t ntris = 0; 491 size_t npos = 0; 492 CHK(geom); 493 494 ntris = sa_size(geom->indices) / 3; /* #primitives */ 495 npos = sa_size(geom->positions) / 3; /* #positions */ 496 497 scn_args.get_indices = geometry_get_indices; 498 scn_args.get_interface = geometry_get_interface; 499 scn_args.get_position = geometry_get_position; 500 scn_args.nprimitives = ntris; 501 scn_args.nvertices = npos; 502 scn_args.t_range[0] = 300; 503 scn_args.t_range[1] = 350; 504 scn_args.radenv = radenv; 505 scn_args.context = geom; 506 507 OK(sdis_scene_create(dev, &scn_args, &scn)); 508 return scn; 509 } 510 511 /******************************************************************************* 512 * Rendering point of view 513 ******************************************************************************/ 514 static struct sdis_camera* 515 create_camera(struct sdis_device* dev) 516 { 517 const double pos[3] = {3, 3, 3}; 518 const double tgt[3] = {0, 0, 0}; 519 const double up[3] = {0, 0, 1}; 520 struct sdis_camera* cam = NULL; 521 522 OK(sdis_camera_create(dev, &cam)); 523 OK(sdis_camera_set_proj_ratio(cam, (double)IMG_WIDTH/(double)IMG_HEIGHT)); 524 OK(sdis_camera_set_fov(cam, MDEG2RAD(70))); 525 OK(sdis_camera_look_at(cam, pos, tgt, up)); 526 return cam; 527 } 528 529 /******************************************************************************* 530 * Draw the scene 531 ******************************************************************************/ 532 static void 533 check_pixel(const struct sdis_estimator* pixel) 534 { 535 struct sdis_mc T = SDIS_MC_NULL; 536 struct sdis_mc time = SDIS_MC_NULL; 537 size_t nfails = 0; 538 size_t nreals = 0; 539 540 OK(sdis_estimator_get_realisation_count(pixel, &nreals)); 541 OK(sdis_estimator_get_failure_count(pixel, &nfails)); 542 OK(sdis_estimator_get_temperature(pixel, &T)); 543 OK(sdis_estimator_get_realisation_time(pixel, &time)); 544 545 CHK(nreals + nfails == SPP); 546 CHK(T.E > 0); 547 CHK(time.E > 0); 548 } 549 550 static void 551 check_image(const struct sdis_estimator_buffer* img) 552 { 553 struct sdis_mc T = SDIS_MC_NULL; 554 struct sdis_mc time = SDIS_MC_NULL; 555 const struct sdis_estimator* pixel = NULL; 556 struct ssp_rng* rng_state = NULL; 557 size_t definition[2]; 558 size_t nreals, nfails; 559 size_t x, y; 560 561 BA(sdis_estimator_buffer_get_realisation_count(NULL, &nreals)); 562 BA(sdis_estimator_buffer_get_realisation_count(img, NULL)); 563 OK(sdis_estimator_buffer_get_realisation_count(img, &nreals)); 564 565 BA(sdis_estimator_buffer_get_failure_count(NULL, &nfails)); 566 BA(sdis_estimator_buffer_get_failure_count(img, NULL)); 567 OK(sdis_estimator_buffer_get_failure_count(img, &nfails)); 568 569 BA(sdis_estimator_buffer_get_temperature(NULL, &T)); 570 BA(sdis_estimator_buffer_get_temperature(img, NULL)); 571 OK(sdis_estimator_buffer_get_temperature(img, &T)); 572 573 BA(sdis_estimator_buffer_get_realisation_time(NULL, &time)); 574 BA(sdis_estimator_buffer_get_realisation_time(img, NULL)); 575 OK(sdis_estimator_buffer_get_realisation_time(img, &time)); 576 577 BA(sdis_estimator_buffer_get_rng_state(NULL, &rng_state)); 578 BA(sdis_estimator_buffer_get_rng_state(img, NULL)); 579 OK(sdis_estimator_buffer_get_rng_state(img, &rng_state)); 580 581 CHK(nreals + nfails == IMG_WIDTH*IMG_HEIGHT*SPP); 582 583 fprintf(stderr, "Overall temperature ~ %g +/- %g\n", T.E, T.SE); 584 fprintf(stderr, "Time per realisation (in usec) ~ %g +/- %g\n", time.E, time.SE); 585 fprintf(stderr, "#failures = %lu/%lu\n", 586 (unsigned long)nfails, (unsigned long)(IMG_WIDTH*IMG_HEIGHT*SPP)); 587 588 BA(sdis_estimator_buffer_get_definition(NULL, definition)); 589 BA(sdis_estimator_buffer_get_definition(img, NULL)); 590 OK(sdis_estimator_buffer_get_definition(img, definition)); 591 CHK(definition[0] == IMG_WIDTH); 592 CHK(definition[1] == IMG_HEIGHT); 593 594 BA(sdis_estimator_buffer_at(NULL, 0, 0, &pixel)); 595 BA(sdis_estimator_buffer_at(img, IMG_WIDTH+1,0 , &pixel)); 596 BA(sdis_estimator_buffer_at(img, 0, IMG_HEIGHT+1, &pixel)); 597 BA(sdis_estimator_buffer_at(img, 0, 0, NULL)); 598 599 /* Pixels ordered by row */ 600 FOR_EACH(y, 0, IMG_HEIGHT) { 601 FOR_EACH(x, 0, IMG_WIDTH) { 602 OK(sdis_estimator_buffer_at(img, x, y, &pixel)); 603 check_pixel(pixel); 604 } 605 } 606 } 607 608 /* Check that the images, although compatible from an estimation point of view, 609 * are not the same. This should be the case if different random sequences are 610 * used for each image */ 611 static void 612 check_image_difference 613 (const struct sdis_estimator_buffer* img0, 614 const struct sdis_estimator_buffer* img1) 615 { 616 struct sdis_mc T0 = SDIS_MC_NULL; 617 struct sdis_mc T1 = SDIS_MC_NULL; 618 size_t definition0[2]; 619 size_t definition1[2]; 620 size_t nreals0, nfails0; 621 size_t nreals1, nfails1; 622 623 OK(sdis_estimator_buffer_get_realisation_count(img0, &nreals0)); 624 OK(sdis_estimator_buffer_get_realisation_count(img1, &nreals1)); 625 626 OK(sdis_estimator_buffer_get_failure_count(img0, &nfails0)); 627 OK(sdis_estimator_buffer_get_failure_count(img1, &nfails1)); 628 CHK(nreals0 + nfails0 == nreals1 + nfails1); 629 630 OK(sdis_estimator_buffer_get_definition(img0, definition0)); 631 OK(sdis_estimator_buffer_get_definition(img1, definition1)); 632 CHK(definition0[0] == definition1[0]); 633 CHK(definition0[1] == definition1[1]); 634 635 OK(sdis_estimator_buffer_get_temperature(img0, &T0)); 636 OK(sdis_estimator_buffer_get_temperature(img1, &T1)); 637 CHK(T0.E != T1.E || (T0.SE == 0 && T1.SE == 0)); 638 CHK(T0.E + 3*T0.SE >= T1.E - 3*T1.SE); 639 CHK(T0.E - 3*T0.SE <= T1.E + 3*T1.SE); 640 } 641 642 /* Write an image in htrdr-image(5) format */ 643 static void 644 write_image(FILE* stream, struct sdis_estimator_buffer* img) 645 { 646 size_t x, y; 647 648 /* Header */ 649 fprintf(stream, "%d %d\n", IMG_WIDTH, IMG_HEIGHT); 650 651 /* Pixels ordered by row */ 652 FOR_EACH(y, 0, IMG_HEIGHT) { 653 FOR_EACH(x, 0, IMG_WIDTH) { 654 struct sdis_mc T = SDIS_MC_NULL; 655 const struct sdis_estimator* pixel = NULL; 656 657 OK(sdis_estimator_buffer_at(img, x, y, &pixel)); 658 OK(sdis_estimator_get_temperature(pixel, &T)); 659 fprintf(stream, "%g %g 0 0 0 0 0 0\n", T.E, T.SE); 660 } 661 } 662 } 663 664 static void 665 invalid_draw(struct sdis_scene* scn, struct sdis_camera* cam) 666 { 667 struct sdis_solve_camera_args args = SDIS_SOLVE_CAMERA_ARGS_DEFAULT; 668 struct sdis_estimator_buffer* img = NULL; 669 670 CHK(cam && scn); 671 672 args.cam = cam; 673 args.time_range[0] = INF; 674 args.time_range[1] = INF; 675 args.image_definition[0] = IMG_WIDTH; 676 args.image_definition[1] = IMG_HEIGHT; 677 args.spp = SPP; 678 679 BA(sdis_solve_camera(NULL, &args, &img)); 680 BA(sdis_solve_camera(scn, NULL, &img)); 681 BA(sdis_solve_camera(scn, &args, NULL)); 682 683 /* Invalid camera */ 684 args.cam = NULL; 685 BA(sdis_solve_camera(scn, &args, &img)); 686 args.cam = cam; 687 688 /* Invald time range */ 689 args.time_range[0] = args.time_range[1] = -1; 690 BA(sdis_solve_camera(scn, &args, &img)); 691 args.time_range[0] = 1; 692 BA(sdis_solve_camera(scn, &args, &img)); 693 args.time_range[1] = 0; 694 BA(sdis_solve_camera(scn, &args, &img)); 695 args.time_range[0] = args.time_range[1] = INF; 696 697 /* Invalid image definition */ 698 args.image_definition[0] = 0; 699 BA(sdis_solve_camera(scn, &args, &img)); 700 args.image_definition[0] = IMG_WIDTH; 701 args.image_definition[1] = 0; 702 BA(sdis_solve_camera(scn, &args, &img)); 703 args.image_definition[1] = IMG_HEIGHT; 704 705 /* Invalid number of samples per pixel */ 706 args.spp = 0; 707 BA(sdis_solve_camera(scn, &args, &img)); 708 args.spp = SPP; 709 } 710 711 static void 712 draw 713 (struct sdis_scene* scn, 714 struct sdis_camera* cam, 715 const int is_master_process) 716 { 717 struct sdis_solve_camera_args args = SDIS_SOLVE_CAMERA_ARGS_DEFAULT; 718 struct sdis_estimator_buffer* img0 = NULL; 719 struct sdis_estimator_buffer* img1 = NULL; 720 struct ssp_rng* rng = NULL; 721 722 args.cam = cam; 723 args.time_range[0] = INF; 724 args.time_range[1] = INF; 725 args.image_definition[0] = IMG_WIDTH; 726 args.image_definition[1] = IMG_HEIGHT; 727 args.spp = SPP; 728 729 OK(sdis_solve_camera(scn, &args, &img0)); 730 731 if(!is_master_process) { 732 CHK(img0 == NULL); 733 } else { 734 check_image(img0); 735 write_image(stdout, img0); 736 } 737 738 /* Provide a RNG state and check that the results, although compatible, are 739 * not the same */ 740 OK(ssp_rng_create(NULL, SSP_RNG_THREEFRY, &rng)); 741 OK(ssp_rng_discard(rng, 3141592653589)); /* Move the RNG state */ 742 args.rng_state = rng; 743 args.rng_type = SSP_RNG_TYPE_NULL; 744 OK(sdis_solve_camera(scn, &args, &img1)); 745 if(is_master_process) { 746 check_image_difference(img0, img1); 747 OK(sdis_estimator_buffer_ref_put(img1)); 748 } 749 OK(ssp_rng_ref_put(rng)); 750 751 /* Change the RNG type and check that the results, although compatible, are 752 * not the same */ 753 args.rng_state = NULL; 754 args.rng_type = SDIS_SOLVE_CAMERA_ARGS_DEFAULT.rng_type == SSP_RNG_THREEFRY 755 ? SSP_RNG_MT19937_64 : SSP_RNG_THREEFRY; 756 OK(sdis_solve_camera(scn, &args, &img1)); 757 if(is_master_process) { 758 check_image_difference(img0, img1); 759 OK(sdis_estimator_buffer_ref_put(img1)); 760 } 761 762 if(is_master_process) OK(sdis_estimator_buffer_ref_put(img0)); 763 } 764 765 /******************************************************************************* 766 * Test 767 ******************************************************************************/ 768 int 769 main(int argc, char** argv) 770 { 771 struct geometry geom = GEOMETRY_NULL; 772 struct sdis_camera* cam = NULL; 773 struct sdis_device* dev = NULL; 774 struct sdis_medium* solid = NULL; 775 struct sdis_medium* fluid0 = NULL; 776 struct sdis_medium* fluid1 = NULL; 777 struct sdis_medium* fluid2 = NULL; 778 struct sdis_interface* interf0 = NULL; 779 struct sdis_interface* interf1 = NULL; 780 struct sdis_interface* interf2 = NULL; 781 struct sdis_radiative_env* radenv = NULL; 782 struct sdis_scene* scn = NULL; 783 struct fluid fluid_args = FLUID_NULL; 784 struct solid solid_args = SOLID_NULL; 785 struct interf interface_args = INTERF_NULL; 786 int is_master_process; 787 (void)argc, (void)argv; 788 789 create_default_device(&argc, &argv, &is_master_process, &dev); 790 791 radenv = create_radenv(dev, 300, 300); 792 793 /* Create the fluid0 */ 794 fluid_args.temperature = 350; 795 fluid_args.rho = 0; 796 fluid_args.cp = 0; 797 fluid0 = create_fluid(dev, &fluid_args); 798 799 /* Create the fluid1 */ 800 fluid_args.temperature = 300; 801 fluid_args.rho = 0; 802 fluid_args.cp = 0; 803 fluid1 = create_fluid(dev, &fluid_args); 804 805 /* Create the fluid2 */ 806 fluid_args.temperature = 280; 807 fluid_args.rho = 0; 808 fluid_args.cp = 0; 809 fluid2 = create_fluid(dev, &fluid_args); 810 811 /* Create the solid medium */ 812 solid_args.cp = 1.0; 813 solid_args.lambda = 0.1; 814 solid_args.rho = 1.0; 815 solid_args.delta = 1.0/20.0; 816 solid_args.temperature = SDIS_TEMPERATURE_NONE; 817 solid = create_solid(dev, &solid_args); 818 819 /* Create the fluid0/solid interface */ 820 interface_args.hc = 1; 821 interface_args.epsilon = 0; 822 interface_args.specular_fraction = 0; 823 interface_args.temperature = SDIS_TEMPERATURE_NONE; 824 interf0 = create_interface(dev, solid, fluid0, &interface_args); 825 826 /* Create the fluid1/solid interface */ 827 interface_args.hc = 0.1; 828 interface_args.epsilon = 1; 829 interface_args.specular_fraction = 1; 830 interface_args.temperature = SDIS_TEMPERATURE_NONE; 831 interface_args.reference_temperature = 300; 832 interf1 = create_interface(dev, fluid1, solid, &interface_args); 833 834 /* Create the fluid2/ground interface */ 835 interface_args.hc = 0.2; 836 interface_args.epsilon = 1; 837 interface_args.specular_fraction = 1; 838 interface_args.temperature = SDIS_TEMPERATURE_NONE; 839 interface_args.reference_temperature = 300; 840 interf2 = create_interface(dev, fluid2, solid, &interface_args); 841 842 add_cube(&geom, interf1); 843 add_sphere(&geom, interf0); 844 add_ground(&geom, interf2); 845 846 #if 0 847 dump_mesh(stdout, geom.positions, sa_size(geom.positions)/3, geom.indices, 848 sa_size(geom.indices)/3); 849 #endif 850 851 scn = create_scene(dev, radenv, &geom); 852 cam = create_camera(dev); 853 854 invalid_draw(scn, cam); 855 draw(scn, cam, is_master_process); 856 857 /* Release memory */ 858 OK(sdis_medium_ref_put(solid)); 859 OK(sdis_medium_ref_put(fluid0)); 860 OK(sdis_medium_ref_put(fluid1)); 861 OK(sdis_medium_ref_put(fluid2)); 862 OK(sdis_radiative_env_ref_put(radenv)); 863 OK(sdis_scene_ref_put(scn)); 864 OK(sdis_camera_ref_put(cam)); 865 OK(sdis_interface_ref_put(interf0)); 866 OK(sdis_interface_ref_put(interf1)); 867 OK(sdis_interface_ref_put(interf2)); 868 free_default_device(dev); 869 geometry_release(&geom); 870 871 CHK(mem_allocated_size() == 0); 872 return 0; 873 }