sdis_heat_path_boundary_Xd_c.h (40039B)
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_green.h" 17 #include "sdis_heat_path_boundary_c.h" 18 #include "sdis_interface_c.h" 19 #include "sdis_log.h" 20 #include "sdis_medium_c.h" 21 #include "sdis_misc.h" 22 #include "sdis_scene_c.h" 23 24 #include <star/ssp.h> 25 26 #include "sdis_Xd_begin.h" 27 28 struct XD(find_reinjection_ray_args) { 29 const struct rwalk* rwalk; /* Current random walk state */ 30 float dir0[DIM]; /* Challenged ray direction */ 31 float dir1[DIM]; /* Challenged ray direction */ 32 double distance; /* Maximum reinjection distance */ 33 unsigned solid_enc_id; /* Enclosure id into which the reinjection occurs */ 34 35 /* Define if the random walk position can be moved or not to find a valid 36 * reinjection direction */ 37 int can_move; 38 }; 39 static const struct XD(find_reinjection_ray_args) 40 XD(FIND_REINJECTION_RAY_ARGS_NULL) = { NULL, {0}, {0}, 0, ENCLOSURE_ID_NULL, 0 }; 41 42 struct XD(reinjection_ray) { 43 double org[DIM]; /* Origin of the reinjection */ 44 float dir[DIM]; /* Direction of the reinjection */ 45 float dst; /* Reinjection distance along dir */ 46 struct sXd(hit) hit; /* Hit along the reinjection dir */ 47 48 /* Define whether or not the random walk was moved to find this reinjection 49 * ray */ 50 int position_was_moved; 51 }; 52 static const struct XD(reinjection_ray) 53 XD(REINJECTION_RAY_NULL) = { {0}, {0}, 0, SXD_HIT_NULL__, 0 }; 54 55 /******************************************************************************* 56 * Helper functions 57 ******************************************************************************/ 58 static INLINE res_T 59 XD(check_find_reinjection_ray_args) 60 (struct sdis_scene* scn, 61 const struct XD(find_reinjection_ray_args)* args) 62 { 63 const struct enclosure* enc = NULL; 64 struct sdis_medium* mdm = NULL; 65 res_T res = RES_OK; 66 ASSERT(scn); 67 68 /* Check pointers */ 69 if(!args || !args->rwalk) return RES_BAD_ARG; 70 71 /* Check distance */ 72 if(args->distance <= 0) return RES_BAD_ARG; 73 74 /* Check directions */ 75 if(!fX(is_normalized)(args->dir0) || !fX(is_normalized)(args->dir1)) { 76 return RES_BAD_ARG; 77 } 78 79 /* Check enclosure id */ 80 if(args->solid_enc_id == ENCLOSURE_ID_NULL) { 81 return RES_BAD_ARG; 82 } 83 enc = scene_get_enclosure(scn, args->solid_enc_id); 84 85 /* Check the enclosure */ 86 enc = scene_get_enclosure(scn, args->solid_enc_id); 87 if(enc->medium_id != MEDIUM_ID_MULTI) { 88 if((res = scene_get_enclosure_medium(scn, enc, &mdm)) != RES_OK) return res; 89 if(sdis_medium_get_type(mdm) != SDIS_SOLID) { 90 res = RES_BAD_ARG; 91 } 92 } 93 94 return RES_OK; 95 } 96 97 static INLINE res_T 98 XD(check_sample_reinjection_step_args) 99 (struct sdis_scene* scn, 100 const struct sample_reinjection_step_args* args) 101 { 102 const struct enclosure* enc = NULL; 103 struct sdis_medium* mdm = NULL; 104 res_T res = RES_OK; 105 ASSERT(scn); 106 107 /* Check pointers */ 108 if(!args || !args->rng || !args->rwalk) return RES_BAD_ARG; 109 110 /* Check distance */ 111 if(args->distance <= 0) return RES_BAD_ARG; 112 113 /* Check side */ 114 if((unsigned)args->side >= SDIS_SIDE_NULL__) return RES_BAD_ARG; 115 116 /* Check enclosure id */ 117 if(args->solid_enc_id == ENCLOSURE_ID_NULL) { 118 return RES_BAD_ARG; 119 } 120 121 /* Check the enclosure */ 122 enc = scene_get_enclosure(scn, args->solid_enc_id); 123 if(enc->medium_id != MEDIUM_ID_MULTI) { 124 if((res = scene_get_enclosure_medium(scn, enc, &mdm)) != RES_OK) return res; 125 if(sdis_medium_get_type(mdm) != SDIS_SOLID) { 126 return RES_BAD_ARG; 127 } 128 } 129 130 return RES_OK; 131 } 132 133 static INLINE res_T 134 XD(check_reinjection_step)(const struct reinjection_step* step) 135 { 136 /* Check pointer */ 137 if(!step) return RES_BAD_ARG; 138 139 /* Check direction */ 140 if(!fX(is_normalized)(step->direction)) return RES_BAD_ARG; 141 142 /* Check distance */ 143 if(step->distance <= 0) return RES_BAD_ARG; 144 145 return RES_OK; 146 } 147 148 static INLINE res_T 149 XD(check_solid_reinjection_args)(const struct solid_reinjection_args* args) 150 { 151 /* Check pointers */ 152 if(!args || !args->rng || !args->rwalk || !args->rwalk_ctx || !args->T) 153 return RES_BAD_ARG; 154 155 /* Check unit */ 156 if(args->fp_to_meter <= 0) return RES_BAD_ARG; 157 158 return XD(check_reinjection_step)(args->reinjection); 159 } 160 161 /* Check that the interface fragment is consistent with the current state of 162 * the random walk */ 163 static INLINE res_T 164 XD(check_rwalk_fragment_consistency) 165 (const struct rwalk* rwalk, 166 const struct sdis_interface_fragment* frag) 167 { 168 double N[DIM]; 169 double uv[2] = {0, 0}; 170 ASSERT(rwalk && frag); 171 172 /* Check intersection */ 173 if(SXD_HIT_NONE(&rwalk->XD(hit))) return RES_BAD_ARG; 174 175 /* Check positions */ 176 if(!dX(eq_eps)(rwalk->vtx.P, frag->P, 1.e-6)) return RES_BAD_ARG; 177 178 /* Check normals */ 179 dX(normalize)(N, dX_set_fX(N, rwalk->XD(hit).normal)); 180 if(!dX(eq_eps)(N, frag->Ng, 1.e-6)) return RES_BAD_ARG; 181 182 /* Check time */ 183 if(!eq_eps(rwalk->vtx.time, frag->time, 1.e-6) 184 && !(IS_INF(rwalk->vtx.time) && IS_INF(frag->time))) { 185 return RES_BAD_ARG; 186 } 187 188 /* Check parametric coordinates */ 189 #if (SDIS_XD_DIMENSION == 2) 190 uv[0] = rwalk->XD(hit).u; 191 #else 192 d2_set_f2(uv, rwalk->XD(hit).uv); 193 #endif 194 if(!d2_eq_eps(uv, frag->uv, 1.e-6)) return RES_BAD_ARG; 195 196 return RES_OK; 197 } 198 199 static void 200 XD(sample_reinjection_dir) 201 (const struct rwalk* rwalk, 202 struct ssp_rng* rng, 203 float dir[DIM]) 204 { 205 #if DIM == 2 206 /* The sampled directions is defined by rotating the normal around the Z axis 207 * of an angle of PI/4 or -PI/4. Let the rotation matrix defined as 208 * | cos(a) -sin(a) | 209 * | sin(a) cos(a) | 210 * with a = PI/4, dir = sqrt(2)/2 * | 1 -1 | . N 211 * | 1 1 | 212 * with a =-PI/4, dir = sqrt(2)/2 * | 1 1 | . N 213 * |-1 1 | 214 * Note that since the sampled direction is finally normalized, we can 215 * discard the sqrt(2)/2 constant. */ 216 const uint64_t r = ssp_rng_uniform_uint64(rng, 0, 1); 217 ASSERT(rwalk && rng && dir); 218 ASSERT(!SXD_HIT_NONE(&rwalk->XD(hit))); 219 ASSERT(rwalk->enc_id == ENCLOSURE_ID_NULL); 220 221 if(r) { 222 dir[0] = rwalk->XD(hit).normal[0] - rwalk->XD(hit).normal[1]; 223 dir[1] = rwalk->XD(hit).normal[0] + rwalk->XD(hit).normal[1]; 224 } else { 225 dir[0] = rwalk->XD(hit).normal[0] + rwalk->XD(hit).normal[1]; 226 dir[1] =-rwalk->XD(hit).normal[0] + rwalk->XD(hit).normal[1]; 227 } 228 f2_normalize(dir, dir); 229 #else 230 /* Sample a random direction around the normal whose cosine is 1/sqrt(3). To 231 * do so we sample a position onto a cone whose height is 1/sqrt(2) and the 232 * radius of its base is 1. */ 233 float frame[9]; 234 ASSERT(rwalk && rng && dir); 235 ASSERT(!SXD_HIT_NONE(&rwalk->XD(hit))); 236 ASSERT(rwalk->enc_id == ENCLOSURE_ID_NULL); 237 ASSERT(fX(is_normalized)(rwalk->XD(hit).normal)); 238 239 ssp_ran_circle_uniform_float(rng, dir, NULL); 240 dir[2] = (float)(1.0/sqrt(2)); 241 242 f33_basis(frame, rwalk->XD(hit).normal); 243 f33_mulf3(dir, frame, dir); 244 f3_normalize(dir, dir); 245 ASSERT(eq_epsf(f3_dot(dir, rwalk->XD(hit).normal), (float)(1.0/sqrt(3)), 1.e-4f)); 246 #endif 247 } 248 249 static res_T 250 XD(find_reinjection_ray) 251 (struct sdis_scene* scn, 252 const struct XD(find_reinjection_ray_args)* args, 253 struct XD(reinjection_ray)* ray) 254 { 255 /* Emperical scale factor applied to the challenged reinjection distance. If 256 * the distance to reinject is less than this adjusted value, the solver will 257 * try to discard the reinjection distance if possible in order to avoid 258 * numerical issues. */ 259 const float REINJECT_DST_MIN_SCALE = 0.125f; 260 261 /* # attempts to find a ray direction */ 262 int MAX_ATTEMPTS = 1; 263 264 /* Enclosures */ 265 unsigned enc_ids[2] = {ENCLOSURE_ID_NULL, ENCLOSURE_ID_NULL}; 266 unsigned enc0_id = ENCLOSURE_ID_NULL; 267 unsigned enc1_id = ENCLOSURE_ID_NULL; 268 269 struct hit_filter_data filter_data = HIT_FILTER_DATA_NULL; 270 struct sXd(hit) hit; 271 struct sXd(hit) hit0; 272 struct sXd(hit) hit1; 273 double tmp[DIM]; 274 double dst; 275 double dst0; 276 double dst1; 277 const float* dir; 278 float reinject_threshold; 279 double dst_adjusted; 280 float org[DIM]; 281 const float range[2] = {0, FLT_MAX}; 282 enum sdis_side side; 283 int iattempt = 0; 284 res_T res = RES_OK; 285 286 ASSERT(scn && args && ray); 287 ASSERT(XD(check_find_reinjection_ray_args)(scn, args) == RES_OK); 288 289 *ray = XD(REINJECTION_RAY_NULL); 290 MAX_ATTEMPTS = args->can_move ? 20 : 1; 291 292 dst_adjusted = args->distance * RAY_RANGE_MAX_SCALE; 293 reinject_threshold = (float)args->distance * REINJECT_DST_MIN_SCALE; 294 295 dX(set)(ray->org, args->rwalk->vtx.P); 296 297 do { 298 fX_set_dX(org, ray->org); 299 filter_data.XD(hit) = args->rwalk->XD(hit); 300 301 /* Limit the epsilon to 1.e-6, as Star-3D's single-precision floating-point 302 * representation will inevitably present numerical accuracy problems below 303 * this threshold. There's no point in going any lower */ 304 /*filter_data.epsilon = MMAX(args->distance * 0.01, 1e-6);*/ 305 filter_data.epsilon = args->distance * 0.01; 306 307 SXD(scene_view_trace_ray 308 (scn->sXd(view), org, args->dir0, range, &filter_data, &hit0)); 309 SXD(scene_view_trace_ray 310 (scn->sXd(view), org, args->dir1, range, &filter_data, &hit1)); 311 312 /* Retrieve the enclosure at the reinjection pos along dir0 */ 313 if(!SXD_HIT_NONE(&hit0)) { 314 scene_get_enclosure_ids(scn, hit0.prim.prim_id, enc_ids); 315 side = fX(dot)(args->dir0, hit0.normal) < 0 ? SDIS_FRONT : SDIS_BACK; 316 enc0_id = enc_ids[side]; 317 } else { 318 XD(move_pos)(dX(set)(tmp, ray->org), args->dir0, (float)args->distance); 319 res = scene_get_enclosure_id_in_closed_boundaries(scn, tmp, &enc0_id); 320 if(res == RES_BAD_OP) { enc0_id = ENCLOSURE_ID_NULL; res = RES_OK; } 321 if(res != RES_OK) goto error; 322 } 323 324 /* Retrieve the enclosure at the reinjection pos along dir1 */ 325 if(!SXD_HIT_NONE(&hit1)) { 326 scene_get_enclosure_ids(scn, hit1.prim.prim_id, enc_ids); 327 side = fX(dot)(args->dir1, hit1.normal) < 0 ? SDIS_FRONT : SDIS_BACK; 328 enc1_id = enc_ids[side]; 329 } else { 330 XD(move_pos)(dX(set)(tmp, ray->org), args->dir1, (float)args->distance); 331 res = scene_get_enclosure_id_in_closed_boundaries(scn, tmp, &enc1_id); 332 if(res == RES_BAD_OP) { enc1_id = ENCLOSURE_ID_NULL; res = RES_OK; } 333 if(res != RES_OK) goto error; 334 } 335 336 dst0 = dst1 = -1; 337 if(enc0_id == args->solid_enc_id) { /* Check reinjection consistency */ 338 if(hit0.distance <= dst_adjusted) { 339 dst0 = hit0.distance; 340 } else { 341 dst0 = args->distance; 342 hit0 = SXD_HIT_NULL; 343 } 344 } 345 if(enc1_id == args->solid_enc_id) { /* Check reinjection consistency */ 346 if(hit1.distance <= dst_adjusted) { 347 dst1 = hit1.distance; 348 } else { 349 dst1 = args->distance; 350 hit1 = SXD_HIT_NULL; 351 } 352 } 353 354 /* No valid reinjection. Maybe the random walk is near a sharp corner and 355 * thus the ray-tracing misses the enclosure geometry. Another possibility 356 * is that the random walk lies roughly on an edge. In this case, sampled 357 * reinjection dirs can intersect the primitive on the other side of the 358 * edge. Normally, this primitive should be filtered by the "hit_filter" 359 * function but this may be not the case due to a "threshold effect". In 360 * both situations, try to slightly move away from the primitive boundaries 361 * and retry to find a valid reinjection. */ 362 if(dst0 == -1 && dst1 == -1 363 && iattempt < MAX_ATTEMPTS - 1) { /* Is there still a trial to be done? */ 364 XD(move_away_primitive_boundaries) 365 (&args->rwalk->XD(hit), args->distance, ray->org); 366 ray->position_was_moved = 1; 367 } 368 } while(dst0 == -1 && dst1 == -1 && ++iattempt < MAX_ATTEMPTS); 369 370 if(dst0 == -1 && dst1 == -1) { /* No valid reinjection */ 371 log_err(scn->dev, "%s: no valid reinjection direction at {"FORMAT_VECX"}.\n", 372 FUNC_NAME, SPLITX(ray->org)); 373 res = RES_BAD_OP_IRRECOVERABLE; 374 goto error; 375 } 376 377 if(dst0 == -1) { 378 /* Invalid dir0 -> move along dir1 */ 379 dir = args->dir1; 380 dst = dst1; 381 hit = hit1; 382 } else if(dst1 == -1) { 383 /* Invalid dir1 -> move along dir0 */ 384 dir = args->dir0; 385 dst = dst0; 386 hit = hit0; 387 } else if(dst0 < reinject_threshold && dst1 < reinject_threshold) { 388 /* The displacement along dir0 and dir1 are both below the reinjection 389 * threshold that defines a distance under which the temperature gradients 390 * are ignored. Move along the direction that allows the maximum 391 * displacement. */ 392 if(dst0 > dst1) { 393 dir = args->dir0; 394 dst = dst0; 395 hit = hit0; 396 } else { 397 dir = args->dir1; 398 dst = dst1; 399 hit = hit1; 400 } 401 } else if(dst0 < reinject_threshold) { 402 /* Ingore dir0 that is bellow the reinject threshold */ 403 dir = args->dir1; 404 dst = dst1; 405 hit = hit1; 406 } else if(dst1 < reinject_threshold) { 407 /* Ingore dir1 that is bellow the reinject threshold */ 408 dir = args->dir0; 409 dst = dst0; 410 hit = hit0; 411 } else { 412 /* All reinjection directions are valid. Choose the first 1 that was 413 * randomly selected by the sample_reinjection_dir procedure and adjust 414 * the displacement distance. */ 415 dir = args->dir0; 416 417 /* Define the reinjection distance along dir0 and its corresponding hit */ 418 if(dst0 <= dst1) { 419 dst = dst0; 420 hit = hit0; 421 } else { 422 dst = dst1; 423 hit = SXD_HIT_NULL; 424 } 425 426 /* If the displacement distance is too close of a boundary, move to the 427 * boundary in order to avoid numerical uncertainty. */ 428 if(!SXD_HIT_NONE(&hit0) 429 && dst0 != dst 430 && eq_eps(dst0, dst, dst0*0.1)) { 431 dst = dst0; 432 hit = hit0; 433 } 434 } 435 436 /* Setup the ray */ 437 fX(set)(ray->dir, dir); 438 ray->dst = (float)dst; 439 ray->hit = hit; 440 441 exit: 442 return res; 443 error: 444 goto exit; 445 } 446 447 static res_T 448 XD(find_reinjection_ray_and_check_validity) 449 (struct sdis_scene* scn, 450 const struct XD(find_reinjection_ray_args)* args, 451 struct XD(reinjection_ray)* ray) 452 { 453 double pos[DIM]; 454 res_T res = RES_OK; 455 456 ASSERT(scn && args && ray); 457 ASSERT(XD(check_find_reinjection_ray_args)(scn, args) == RES_OK); 458 459 /* Select a reinjection direction */ 460 res = XD(find_reinjection_ray)(scn, args, ray); 461 if(res != RES_OK) goto error; 462 463 if(SXD_HIT_NONE(&ray->hit)) { 464 unsigned enc_id = ENCLOSURE_ID_NULL; 465 466 /* Obtain the enclosure in which the reinjection position lies */ 467 XD(move_pos)(dX(set)(pos, ray->org), ray->dir, (float)ray->dst); 468 res = scene_get_enclosure_id_in_closed_boundaries(scn, pos, &enc_id); 469 if(res != RES_OK) goto error; 470 471 /* Check enclosure consistency at the reinjection position */ 472 if(enc_id != args->solid_enc_id) { 473 res = RES_BAD_OP; 474 goto error; 475 } 476 } 477 478 exit: 479 return res; 480 error: 481 goto exit; 482 } 483 484 static res_T 485 XD(handle_volumic_power) 486 (struct sdis_medium* solid, 487 struct rwalk_context* rwalk_ctx, 488 struct rwalk* rwalk, 489 const double reinject_dst_m, 490 struct temperature* T) 491 { 492 double power; 493 double lambda; 494 double power_term; 495 size_t picard_order; 496 res_T res = RES_OK; 497 498 /* Check pre-conditions */ 499 ASSERT(solid && rwalk_ctx && rwalk && T && reinject_dst_m > 0); 500 501 /* Fetch the volumic power */ 502 power = solid_get_volumic_power(solid, &rwalk->vtx); 503 if(power == SDIS_VOLUMIC_POWER_NONE) goto exit; /* Do nothing */ 504 505 /* Currently, the power term can be correctly taken into account only when 506 * the radiative temperature is linearized, i.e. when the picard order is 507 * equal to 1 */ 508 picard_order = get_picard_order(rwalk_ctx); 509 if(picard_order > 1) { 510 log_err(solid->dev, 511 "%s: invalid not null volumic power '%g' kg/m^3. Could not manage a " 512 "volumic power when the picard order is not equal to 1; Picard order is " 513 "currently set to %lu.\n", 514 FUNC_NAME, power, (unsigned long)picard_order); 515 res = RES_BAD_ARG; 516 goto error; 517 } 518 519 /* Fetch the conductivity */ 520 lambda = solid_get_thermal_conductivity(solid, &rwalk->vtx); 521 522 /* Compute the power term and handle the volumic power */ 523 power_term = (reinject_dst_m * reinject_dst_m)/ (2.0 * DIM * lambda); 524 T->value += power * power_term; 525 526 /* Update the green path with the power term */ 527 if(rwalk_ctx->green_path) { 528 res = green_path_add_power_term 529 (rwalk_ctx->green_path, solid, &rwalk->vtx, power_term); 530 if(res != RES_OK) goto error; 531 } 532 533 exit: 534 return res; 535 error: 536 goto exit; 537 } 538 539 /******************************************************************************* 540 * Local functions 541 ******************************************************************************/ 542 res_T 543 XD(sample_reinjection_step_solid_fluid) 544 (struct sdis_scene* scn, 545 const struct sample_reinjection_step_args* args, 546 struct reinjection_step* step) 547 { 548 /* Input/output data of the function finding a valid reinjection ray */ 549 struct XD(find_reinjection_ray_args) find_reinject_ray_args = 550 XD(FIND_REINJECTION_RAY_ARGS_NULL); 551 struct XD(reinjection_ray) ray = XD(REINJECTION_RAY_NULL); 552 553 /* Enclosures */ 554 const struct enclosure* solid_enc = NULL; 555 unsigned enc_ids[2] = {ENCLOSURE_ID_NULL, ENCLOSURE_ID_NULL}; 556 557 /* In 2D it is useless to try to resample a reinjection direction since there 558 * is only one possible direction */ 559 const int MAX_ATTEMPTS = DIM == 2 ? 1 : 10; 560 561 /* Miscellaneous variables */ 562 float dir0[DIM]; /* Sampled direction */ 563 float dir1[DIM]; /* Sampled direction reflected */ 564 int iattempt = 0; /* #attempts to find a reinjection dir */ 565 res_T res = RES_OK; 566 567 /* Pre-conditions */ 568 ASSERT(scn && args && step); 569 ASSERT(XD(check_sample_reinjection_step_args)(scn, args) == RES_OK); 570 571 /* Initialise the reinjection step */ 572 *step = REINJECTION_STEP_NULL; 573 574 /* Check whether the solid enclosure is part of the system, or whether it is 575 * only there to describe boundary conditions. In the latter case, the 576 * enclosure has no geometrical existence and it is sufficient to return a 577 * valid reinjection step which will be used to select the next step. Note 578 * that if the trajectory passes through the solid enclosure, it will stop, 579 * i.e. the temperature of the solid should be fixed. If it doesn't, it's a 580 * error */ 581 scene_get_enclosure_ids(scn, args->rwalk->XD(hit).prim.prim_id, enc_ids); 582 solid_enc = scene_get_enclosure(scn, enc_ids[args->side]); 583 if(solid_enc->medium_id == MEDIUM_ID_MULTI) { 584 step->XD(hit) = SXD_HIT_NULL; 585 fX(normalize)(step->direction, args->rwalk->XD(hit).normal); 586 if(args->side == SDIS_BACK) fX(minus)(step->direction, step->direction); 587 step->distance = (float)args->distance; 588 goto exit; /* That's all folks! */ 589 } 590 591 iattempt = 0; 592 do { 593 /* Sample a reinjection direction */ 594 XD(sample_reinjection_dir)(args->rwalk, args->rng, dir0); 595 596 /* Reflect the sampled direction around the normal */ 597 XD(reflect)(dir1, dir0, args->rwalk->XD(hit).normal); 598 599 /* Flip the sampled directions if one wants to reinject to back side */ 600 if(args->side == SDIS_BACK) { 601 fX(minus)(dir0, dir0); 602 fX(minus)(dir1, dir1); 603 } 604 605 /* Find the reinjection step */ 606 find_reinject_ray_args.solid_enc_id = args->solid_enc_id; 607 find_reinject_ray_args.rwalk = args->rwalk; 608 find_reinject_ray_args.distance = args->distance; 609 find_reinject_ray_args.can_move = 1; 610 fX(set)(find_reinject_ray_args.dir0, dir0); 611 fX(set)(find_reinject_ray_args.dir1, dir1); 612 res = XD(find_reinjection_ray_and_check_validity) 613 (scn, &find_reinject_ray_args, &ray); 614 if(res == RES_BAD_OP) continue; /* Cannot find a valid reinjection ray. Retry */ 615 if(res != RES_OK) goto error; 616 617 } while(res != RES_OK && ++iattempt < MAX_ATTEMPTS); 618 619 /* Could not find a valid reinjecton step */ 620 if(iattempt >= MAX_ATTEMPTS) { 621 log_err(scn->dev, 622 "%s: could not find a valid reinjection step at `%g %g %g'.\n", 623 FUNC_NAME, SPLIT3(args->rwalk->vtx.P)); 624 res = RES_BAD_OP_IRRECOVERABLE; 625 goto error; 626 } 627 628 /* Setup the reinjection step */ 629 step->XD(hit) = ray.hit; 630 step->distance = ray.dst; 631 fX(set)(step->direction, ray.dir); 632 633 /* Update the random walk position if necessary */ 634 if(ray.position_was_moved) { 635 dX(set)(args->rwalk->vtx.P, ray.org); 636 } 637 638 /* Post-conditions */ 639 ASSERT(dX(eq)(args->rwalk->vtx.P, ray.org)); 640 ASSERT(XD(check_reinjection_step)(step) == RES_OK); 641 642 exit: 643 return res; 644 error: 645 goto exit; 646 } 647 648 res_T 649 XD(sample_reinjection_step_solid_solid) 650 (struct sdis_scene* scn, 651 const struct sample_reinjection_step_args* args_frt, 652 const struct sample_reinjection_step_args* args_bck, 653 struct reinjection_step* step_frt, 654 struct reinjection_step* step_bck) 655 { 656 /* Input/output data of the function finding a valid reinjection ray */ 657 struct XD(find_reinjection_ray_args) find_reinject_ray_frt_args = 658 XD(FIND_REINJECTION_RAY_ARGS_NULL); 659 struct XD(find_reinjection_ray_args) find_reinject_ray_bck_args = 660 XD(FIND_REINJECTION_RAY_ARGS_NULL); 661 struct XD(reinjection_ray) ray_frt = XD(REINJECTION_RAY_NULL); 662 struct XD(reinjection_ray) ray_bck = XD(REINJECTION_RAY_NULL); 663 664 /* Initial random walk position used as a backup */ 665 double rwalk_pos_backup[DIM]; 666 667 /* Variables shared by the two side */ 668 struct rwalk* rwalk = NULL; 669 struct ssp_rng* rng = NULL; 670 671 /* In 2D it is useless to try to resample a reinjection direction since there 672 * is only one possible direction */ 673 const int MAX_ATTEMPTS = DIM == 2 ? 1 : 10; 674 675 /* Enclosure */ 676 unsigned enc_ids[2]; 677 const struct enclosure* enc_frt = NULL; 678 const struct enclosure* enc_bck = NULL; 679 680 float dir_frt_samp[DIM]; /* Sampled direction */ 681 float dir_frt_refl[DIM]; /* Sampled direction reflected */ 682 float dir_bck_samp[DIM]; /* Negated sampled direction */ 683 float dir_bck_refl[DIM]; /* Negated sampled direction reflected */ 684 int multi_frt = 0; 685 int multi_bck = 0; 686 int iattempt = 0; /* #attempts to find a reinjection dir */ 687 res_T res = RES_OK; 688 689 /* Pre-conditions */ 690 ASSERT(scn && args_frt && args_bck && step_frt && step_bck); 691 ASSERT(XD(check_sample_reinjection_step_args)(scn, args_frt) == RES_OK); 692 ASSERT(XD(check_sample_reinjection_step_args)(scn, args_bck) == RES_OK); 693 ASSERT(args_frt->side == SDIS_FRONT); 694 ASSERT(args_bck->side == SDIS_BACK); 695 ASSERT(SXD_PRIMITIVE_EQ(&args_frt->rwalk->XD(hit).prim, &args_bck->rwalk->XD(hit).prim)); 696 697 rng = args_frt->rng; 698 rwalk = args_frt->rwalk; 699 ASSERT(args_bck->rng == rng); 700 ASSERT(args_bck->rwalk == rwalk); 701 702 /* Initialise the reinjection steps */ 703 *step_frt = REINJECTION_STEP_NULL; 704 *step_bck = REINJECTION_STEP_NULL; 705 706 /* Check whether the solid enclosure is part of the system, or whether it is 707 * only there to describe boundary conditions. In the latter case, the 708 * enclosure has no geometrical existence and it is sufficient to return a 709 * valid reinjection step which will be used to select the next step. Note 710 * that if the trajectory passes through the solid enclosure, it will stop, 711 * i.e. the temperature of the solid should be fixed. If it doesn't, it's an 712 * error */ 713 scene_get_enclosure_ids(scn, args_frt->rwalk->XD(hit).prim.prim_id, enc_ids); 714 enc_frt = scene_get_enclosure(scn, enc_ids[SDIS_FRONT]); 715 enc_bck = scene_get_enclosure(scn, enc_ids[SDIS_BACK]); 716 if(enc_frt->medium_id == MEDIUM_ID_MULTI) { 717 step_frt->XD(hit) = SXD_HIT_NULL; 718 fX(normalize)(step_frt->direction, args_frt->rwalk->XD(hit).normal); 719 step_frt->distance = (float)args_frt->distance; 720 multi_frt = 1; 721 } 722 if(enc_bck->medium_id == MEDIUM_ID_MULTI) { 723 step_bck->XD(hit) = SXD_HIT_NULL; 724 fX(normalize)(step_bck->direction, args_bck->rwalk->XD(hit).normal); 725 step_bck->distance = (float)args_bck->distance; 726 multi_bck = 1; 727 } 728 729 if(multi_frt && multi_bck) goto exit; /* That's all folks */ 730 731 dX(set)(rwalk_pos_backup, rwalk->vtx.P); 732 iattempt = 0; 733 do { 734 /* Restore random walk pos */ 735 if(iattempt != 0) dX(set)(rwalk->vtx.P, rwalk_pos_backup); 736 737 /* Sample a reinjection direction and reflect it around the normal. Then 738 * reflect them on the back side of the interface. */ 739 XD(sample_reinjection_dir)(rwalk, rng, dir_frt_samp); 740 XD(reflect)(dir_frt_refl, dir_frt_samp, rwalk->XD(hit).normal); 741 fX(minus)(dir_bck_samp, dir_frt_samp); 742 fX(minus)(dir_bck_refl, dir_frt_refl); 743 744 /* Reject the sampling of the re-injection step if it has already been 745 * defined, i.e. if the enclosure is a limit condition */ 746 if(!multi_frt) { 747 /* Find the reinjection ray for the front side */ 748 find_reinject_ray_frt_args.solid_enc_id = args_frt->solid_enc_id; 749 find_reinject_ray_frt_args.rwalk = args_frt->rwalk; 750 find_reinject_ray_frt_args.distance = args_frt->distance; 751 find_reinject_ray_frt_args.can_move = 1; 752 fX(set)(find_reinject_ray_frt_args.dir0, dir_frt_samp); 753 fX(set)(find_reinject_ray_frt_args.dir1, dir_frt_refl); 754 res = XD(find_reinjection_ray_and_check_validity) 755 (scn, &find_reinject_ray_frt_args, &ray_frt); 756 if(res == RES_BAD_OP) continue; 757 if(res != RES_OK) goto error; 758 759 /* Update the random walk position if necessary */ 760 if(ray_frt.position_was_moved) dX(set)(rwalk->vtx.P, ray_frt.org); 761 } 762 763 /* Reject the sampling of the re-injection step if it has already been 764 * defined, i.e. if the enclosure is a limit condition */ 765 if(!multi_bck) { 766 /* Select the reinjection direction and distance for the back side */ 767 find_reinject_ray_bck_args.solid_enc_id = args_bck->solid_enc_id; 768 find_reinject_ray_bck_args.rwalk = args_bck->rwalk; 769 find_reinject_ray_bck_args.distance = args_bck->distance; 770 find_reinject_ray_bck_args.can_move = 1; 771 fX(set)(find_reinject_ray_bck_args.dir0, dir_bck_samp); 772 fX(set)(find_reinject_ray_bck_args.dir1, dir_bck_refl); 773 res = XD(find_reinjection_ray_and_check_validity) 774 (scn, &find_reinject_ray_bck_args, &ray_bck); 775 if(res == RES_BAD_OP) continue; 776 if(res != RES_OK) goto error; 777 778 /* Update the random walk position if necessary */ 779 if(ray_bck.position_was_moved) dX(set)(rwalk->vtx.P, ray_bck.org); 780 781 /* If random walk was moved to find a valid rinjection ray on back side, 782 * one has to find a valid reinjection on front side from the new pos */ 783 if(ray_bck.position_was_moved) { 784 find_reinject_ray_frt_args.can_move = 0; 785 res = XD(find_reinjection_ray_and_check_validity) 786 (scn, &find_reinject_ray_frt_args, &ray_frt); 787 if(res == RES_BAD_OP) continue; 788 if(res != RES_OK) goto error; 789 790 /* Update the random walk position if necessary */ 791 if(ray_frt.position_was_moved) dX(set)(rwalk->vtx.P, ray_frt.org); 792 } 793 } 794 } while(res != RES_OK && ++iattempt < MAX_ATTEMPTS); 795 796 /* Could not find a valid reinjection */ 797 if(iattempt >= MAX_ATTEMPTS) { 798 dX(set)(rwalk->vtx.P, rwalk_pos_backup); /* Restore random walk pos */ 799 log_warn(scn->dev, 800 "%s: could not find a valid solid/solid reinjection at {%g, %g, %g}.\n", 801 FUNC_NAME, SPLIT3(rwalk->vtx.P)); 802 res = RES_BAD_OP_IRRECOVERABLE; 803 goto error; 804 } 805 806 /* Setup the front and back reinjection steps */ 807 if(!multi_frt) { 808 step_frt->XD(hit) = ray_frt.hit; 809 step_frt->distance = ray_frt.dst; 810 fX(set)(step_frt->direction, ray_frt.dir); 811 ASSERT(XD(check_reinjection_step)(step_frt) == RES_OK); /* Post-condition */ 812 } 813 if(!multi_bck) { 814 step_bck->XD(hit) = ray_bck.hit; 815 step_bck->distance = ray_bck.dst; 816 fX(set)(step_bck->direction, ray_bck.dir); 817 ASSERT(XD(check_reinjection_step)(step_bck) == RES_OK); /* Post-condition */ 818 } 819 820 exit: 821 return res; 822 error: 823 goto exit; 824 } 825 826 res_T 827 XD(solid_reinjection) 828 (struct sdis_scene* scn, 829 const unsigned solid_enc_id, 830 struct solid_reinjection_args* args) 831 { 832 /* Properties */ 833 struct solid_props props = SOLID_PROPS_NULL; 834 struct sdis_medium* solid = NULL; 835 const struct enclosure* enc = NULL; 836 837 double reinject_dst_m; /* Reinjection distance in meters */ 838 double mu; 839 res_T res = RES_OK; 840 ASSERT(XD(check_solid_reinjection_args)(args) == RES_OK); 841 ASSERT(solid_enc_id != ENCLOSURE_ID_NULL); 842 843 reinject_dst_m = args->reinjection->distance * args->fp_to_meter; 844 845 /* Get the enclosure medium properties */ 846 enc = scene_get_enclosure(scn, solid_enc_id); 847 res = scene_get_enclosure_medium(scn, enc, &solid); 848 if(res != RES_OK) goto error; 849 ASSERT(sdis_medium_get_type(solid) == SDIS_SOLID); 850 res = solid_get_properties(solid, &args->rwalk->vtx, &props); 851 if(res != RES_OK) goto error; 852 853 /* Manage the volumic power */ 854 res = XD(handle_volumic_power) 855 (solid, args->rwalk_ctx, args->rwalk, reinject_dst_m, args->T); 856 if(res != RES_OK) goto error; 857 858 /* Time rewind */ 859 args->rwalk->enc_id = solid_enc_id; /* Enclosure into which the time is rewind */ 860 mu = (2*DIM*props.lambda)/(props.rho*props.cp*reinject_dst_m*reinject_dst_m); 861 res = time_rewind 862 (scn, mu, props.t0, args->rng, args->rwalk, args->rwalk_ctx, args->T); 863 if(res != RES_OK) goto error; 864 865 /* Test if a limit condition was reached */ 866 if(args->T->done) goto exit; 867 868 /* Move the random walk to the reinjection position */ 869 XD(move_pos) 870 (args->rwalk->vtx.P, 871 args->reinjection->direction, 872 args->reinjection->distance); 873 874 /* The random walk is in the solid */ 875 if(args->reinjection->XD(hit).distance != args->reinjection->distance) { 876 args->T->func = XD(conductive_path); 877 args->rwalk->enc_id = solid_enc_id; 878 args->rwalk->XD(hit) = SXD_HIT_NULL; 879 args->rwalk->hit_side = SDIS_SIDE_NULL__; 880 881 /* The random walk is at a boundary */ 882 } else { 883 args->T->func = XD(boundary_path); 884 args->rwalk->enc_id = ENCLOSURE_ID_NULL; 885 args->rwalk->XD(hit) = args->reinjection->XD(hit); 886 if(fX(dot)(args->reinjection->XD(hit).normal, args->reinjection->direction) < 0) { 887 args->rwalk->hit_side = SDIS_FRONT; 888 } else { 889 args->rwalk->hit_side = SDIS_BACK; 890 } 891 } 892 893 /* Register the new vertex against the heat path */ 894 res = register_heat_vertex 895 (args->rwalk_ctx->heat_path, 896 &args->rwalk->vtx, 897 args->T->value, 898 SDIS_HEAT_VERTEX_CONDUCTION, 899 (int)args->rwalk_ctx->nbranchings); 900 if(res != RES_OK) goto error; 901 902 exit: 903 return res; 904 error: 905 goto exit; 906 } 907 908 res_T 909 XD(handle_net_flux) 910 (const struct sdis_scene* scn, 911 const struct handle_net_flux_args* args, 912 struct temperature* T) 913 { 914 double flux_term; 915 double phi; 916 res_T res = RES_OK; 917 CHK(scn && T); 918 CHK(args->interf && args->frag); 919 CHK(args->h_cond >= 0); 920 CHK(args->h_conv >= 0); 921 CHK(args->h_radi >= 0); 922 CHK(args->h_cond + args->h_conv + args->h_radi > 0); 923 924 phi = interface_side_get_flux(args->interf, args->frag); 925 if(phi == SDIS_FLUX_NONE) goto exit; /* No flux. Do nothing */ 926 927 if(args->picard_order > 1 && phi != 0) { 928 log_err(scn->dev, 929 "%s: invalid flux '%g' W/m^2. Could not manage a flux != 0 when the " 930 "picard order is not equal to 1; Picard order is currently set to %lu.\n", 931 FUNC_NAME, phi, (unsigned long)args->picard_order); 932 res = RES_BAD_ARG; 933 goto error; 934 } 935 936 flux_term = 1.0 / (args->h_cond + args->h_conv + args->h_radi); 937 T->value += phi * flux_term; 938 939 /* Register the net flux term */ 940 if(args->green_path) { 941 res = green_path_add_flux_term 942 (args->green_path, args->interf, args->frag, flux_term); 943 if(res != RES_OK) goto error; 944 } 945 946 exit: 947 return res; 948 error: 949 goto exit; 950 } 951 952 res_T 953 XD(check_Tref) 954 (const struct sdis_scene* scn, 955 const double pos[DIM], 956 const double Tref, 957 const char* func_name) 958 { 959 ASSERT(scn && pos && func_name); 960 961 #define CHECK_TBOUND(Bound, Name) { \ 962 if(SDIS_TEMPERATURE_IS_UNKNOWN(Bound)) { \ 963 log_err(scn->dev, \ 964 "%s: the "Name" temperature cannot be unknown " \ 965 "to sampling a radiative path.\n", \ 966 func_name); \ 967 return RES_BAD_OP_IRRECOVERABLE; \ 968 } \ 969 \ 970 if((Bound) < 0) { \ 971 log_err(scn->dev, \ 972 "%s: the "Name" temperature cannot be negative " \ 973 "to sample a radiative path -- T"Name" = %g K\n", \ 974 func_name, (Bound)); \ 975 return RES_BAD_OP_IRRECOVERABLE; \ 976 } \ 977 } (void) 0 978 CHECK_TBOUND(scn->tmin, "min"); 979 CHECK_TBOUND(scn->tmax, "max"); 980 #undef CHECK_TBOUND 981 982 if(scn->tmin > scn->tmax) { 983 log_err(scn->dev, 984 "%s: the temperature range cannot be degenerated to sample a radiative " 985 "path (Tmin = %g K; Tmax = %g K).\n", 986 func_name, scn->tmin, scn->tmax); 987 return RES_BAD_OP_IRRECOVERABLE; 988 } 989 990 if(SDIS_TEMPERATURE_IS_UNKNOWN(Tref)) { 991 log_err(scn->dev, 992 "%s: the reference temperature is unknown at ("FORMAT_VECX"). " 993 "Sampling a radiative path requires a valid reference temperature field.\n", 994 func_name, SPLITX(pos)); 995 return RES_BAD_OP_IRRECOVERABLE; 996 } 997 998 if(Tref < 0) { 999 log_err(scn->dev, 1000 "%s: the reference temperature is negative at ("FORMAT_VECX") and Tref = %g K. " 1001 "Sampling a radiative path requires a known, positive reference " 1002 "temperature field.\n", 1003 func_name, SPLITX(pos), Tref); 1004 return RES_BAD_OP_IRRECOVERABLE; 1005 } 1006 1007 if(Tref < scn->tmin || scn->tmax < Tref) { 1008 log_err(scn->dev, 1009 "%s: invalid reference temperature at ("FORMAT_VECX") and Tref=%g K. " 1010 "It must be included in the provided temperature range " 1011 "(Tmin = %g K; Tmax = %g K)\n", 1012 func_name, SPLITX(pos), Tref, scn->tmin, scn->tmax); 1013 return RES_BAD_OP_IRRECOVERABLE; 1014 } 1015 1016 return RES_OK; 1017 } 1018 1019 /* This function checks whether the random walk is on a boundary and, if so, 1020 * verifies that the temperature of the medium attached to the interface is 1021 * known. This medium can be different from the medium of the enclosure. Indeed, 1022 * the enclosure can contain several media used to set the temperatures of 1023 * several boundary conditions. Hence this function, which queries the medium on 1024 * the path coming from a boundary */ 1025 res_T 1026 XD(query_medium_temperature_from_boundary) 1027 (struct sdis_scene* scn, 1028 struct rwalk_context* ctx, 1029 struct rwalk* rwalk, 1030 struct temperature* T) 1031 { 1032 struct sdis_interface* interf = NULL; 1033 struct sdis_medium* mdm = NULL; 1034 double temperature = SDIS_TEMPERATURE_NONE; 1035 res_T res = RES_OK; 1036 ASSERT(scn && ctx && rwalk && T); 1037 1038 /* Not at an interface */ 1039 if(SXD_HIT_NONE(&rwalk->XD(hit))) return RES_OK; /* Nothing to do */ 1040 1041 interf = scene_get_interface(scn, rwalk->XD(hit).prim.prim_id); 1042 mdm = rwalk->hit_side==SDIS_FRONT ? interf->medium_front: interf->medium_back; 1043 1044 temperature = medium_get_temperature(mdm, &rwalk->vtx); 1045 1046 /* Check if the temperature is known */ 1047 if(SDIS_TEMPERATURE_IS_UNKNOWN(temperature)) goto exit; 1048 1049 T->value += temperature; 1050 T->done = 1; 1051 1052 if(ctx->green_path) { 1053 res = green_path_set_limit_vertex 1054 (ctx->green_path, mdm, &rwalk->vtx, rwalk->elapsed_time); 1055 if(res != RES_OK) goto error; 1056 } 1057 1058 if(ctx->heat_path) { 1059 heat_path_get_last_vertex(ctx->heat_path)->weight = T->value; 1060 } 1061 1062 exit: 1063 return res; 1064 error: 1065 goto exit; 1066 } 1067 1068 #if DIM == 2 1069 void 1070 XD(move_away_primitive_boundaries) 1071 (const struct sXd(hit)* hit, 1072 const double delta, 1073 double position[DIM]) /* Position to move */ 1074 { 1075 struct sXd(attrib) attr; 1076 float pos[DIM]; 1077 float dir[DIM]; 1078 float len; 1079 const float st = 0.5f; 1080 ASSERT(!SXD_HIT_NONE(hit) && delta > 0); 1081 1082 SXD(primitive_get_attrib(&hit->prim, SXD_POSITION, st, &attr)); 1083 1084 fX_set_dX(pos, position); 1085 fX(sub)(dir, attr.value, pos); 1086 len = fX(normalize)(dir, dir); 1087 len = MMIN(len, (float)(delta*0.1)); 1088 1089 XD(move_pos)(position, dir, len); 1090 } 1091 #else 1092 /* Move the submitted position away from the primitive boundaries to avoid 1093 * numerical issues leading to inconsistent random walks. */ 1094 void 1095 XD(move_away_primitive_boundaries) 1096 (const struct sXd(hit)* hit, 1097 const double delta, 1098 double position[DIM]) 1099 { 1100 struct s3d_attrib v0, v1, v2; /* Triangle vertices */ 1101 float E[3][4]; /* 3D edge equations */ 1102 float dst[3]; /* Distance from current position to edge equation */ 1103 float N[3]; /* Triangle normal */ 1104 float P[3]; /* Random walk position */ 1105 float tmp[3]; 1106 float min_dst, max_dst; 1107 float len; 1108 int imax = 0; 1109 int imin = 0; 1110 int imid = 0; 1111 int i; 1112 ASSERT(delta > 0 && !S3D_HIT_NONE(hit)); 1113 1114 fX_set_dX(P, position); 1115 1116 /* Fetch triangle vertices */ 1117 S3D(triangle_get_vertex_attrib(&hit->prim, 0, S3D_POSITION, &v0)); 1118 S3D(triangle_get_vertex_attrib(&hit->prim, 1, S3D_POSITION, &v1)); 1119 S3D(triangle_get_vertex_attrib(&hit->prim, 2, S3D_POSITION, &v2)); 1120 1121 /* Compute the edge vector */ 1122 f3_sub(E[0], v1.value, v0.value); 1123 f3_sub(E[1], v2.value, v1.value); 1124 f3_sub(E[2], v0.value, v2.value); 1125 1126 /* Compute the triangle normal */ 1127 f3_cross(N, E[1], E[0]); 1128 1129 /* Compute the 3D edge equation */ 1130 f3_normalize(E[0], f3_cross(E[0], E[0], N)); 1131 f3_normalize(E[1], f3_cross(E[1], E[1], N)); 1132 f3_normalize(E[2], f3_cross(E[2], E[2], N)); 1133 E[0][3] = -f3_dot(E[0], v0.value); 1134 E[1][3] = -f3_dot(E[1], v1.value); 1135 E[2][3] = -f3_dot(E[2], v2.value); 1136 1137 /* Compute the distance from current position to the edges */ 1138 dst[0] = f3_dot(E[0], P) + E[0][3]; 1139 dst[1] = f3_dot(E[1], P) + E[1][3]; 1140 dst[2] = f3_dot(E[2], P) + E[2][3]; 1141 1142 /* Retrieve the min and max distance from random walk position to triangle 1143 * edges */ 1144 min_dst = MMIN(MMIN(dst[0], dst[1]), dst[2]); 1145 max_dst = MMAX(MMAX(dst[0], dst[1]), dst[2]); 1146 1147 /* Sort the edges with respect to their distance to the random walk position */ 1148 FOR_EACH(i, 0, 3) { 1149 if(dst[i] == min_dst) { 1150 imin = i; 1151 } else if(dst[i] == max_dst) { 1152 imax = i; 1153 } else { 1154 imid = i; 1155 } 1156 } 1157 (void)imax; 1158 1159 if(eq_eps(dst[imin], 0, delta*1e-3) && eq_eps(dst[imid], 0, delta*1e-3)) { 1160 /* The random position is in a corner, meaning that its distance to the two 1161 * nearest edges is approximately equal to 0. Move it toward the farthest 1162 * edge along its normal to avoid moving too little. */ 1163 len = MMIN(dst[imax]*0.5f, (float)delta*0.1f); 1164 XD(move_pos)(position, f3_minus(tmp, E[imax]), len); 1165 1166 } else { 1167 /* Compute the distance `dst' from the current position to the edges to move 1168 * to, along the normal of the edge from which the random walk is the nearest 1169 * 1170 * +. cos(a) = d / dst => dst = d / cos_a 1171 * / `*. 1172 * / | `*. 1173 * / dst| a /`*. 1174 * / | / `*. 1175 * / | / d `*. 1176 * / |/ `*. 1177 * +---------o----------------+ */ 1178 const float cos_a1 = f3_dot(E[imin], f3_minus(tmp, E[imid])); 1179 const float cos_a2 = f3_dot(E[imin], f3_minus(tmp, E[imax])); 1180 dst[imid] = cos_a1 > 0 ? dst[imid] / cos_a1 : FLT_MAX; 1181 dst[imax] = cos_a2 > 0 ? dst[imax] / cos_a2 : FLT_MAX; 1182 len = MMIN(dst[imid], dst[imax]); 1183 ASSERT(len != FLT_MAX); 1184 1185 /* Define the displacement distance as the minimum between 10 percent of 1186 * delta and len / 2. */ 1187 len = MMIN(len*0.5f, (float)(delta*0.1)); 1188 XD(move_pos)(position, E[imin], len); 1189 } 1190 1191 } 1192 #endif 1193 1194 #include "sdis_Xd_end.h"