commit bbc1fb6f9310b11d6055ddf07b24979fb7a66b8d
parent a55f64aa902c4c7383cbcc0046bef2d6749e3745
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Wed, 21 Feb 2018 20:46:30 +0100
Rename files
Diffstat:
5 files changed, 1242 insertions(+), 1241 deletions(-)
diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt
@@ -52,7 +52,7 @@ set(SDIS_FILES_SRC
sdis_interface.c
sdis_medium.c
sdis_scene.c
- sdis_solve_probe.c)
+ sdis_solve.c)
set(SDIS_FILES_INC_API
sdis.h)
@@ -63,7 +63,8 @@ set(SDIS_FILES_INC
sdis_estimator_c.h
sdis_interface_c.h
sdis_medium_c.h
- sdis_scene_c.h)
+ sdis_scene_c.h
+ sdis_solve_Xd.h)
set(SDIS_FILES_DOC COPYING README.md)
diff --git a/src/sdis_solve.c b/src/sdis_solve.c
@@ -0,0 +1,413 @@
+/* Copyright (C) |Meso|Star> 2016-2018 (contact@meso-star.com)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. */
+
+#include "sdis.h"
+#include "sdis_camera.h"
+#include "sdis_device_c.h"
+#include "sdis_estimator_c.h"
+#include "sdis_solve_Xd.h"
+
+/* Generate the 2D solver */
+#define SDIS_SOLVE_DIMENSION 2
+#include "sdis_solve_Xd.h"
+
+/* Generate the 3D solver */
+#define SDIS_SOLVE_DIMENSION 3
+#include "sdis_solve_Xd.h"
+
+#include <star/ssp.h>
+#include <omp.h>
+
+/*******************************************************************************
+ * Helper functions
+ ******************************************************************************/
+static FINLINE uint16_t
+morton2D_decode(const uint32_t u32)
+{
+ uint32_t x = u32 & 0x55555555;
+ x = (x | (x >> 1)) & 0x33333333;
+ x = (x | (x >> 2)) & 0x0F0F0F0F;
+ x = (x | (x >> 4)) & 0x00FF00FF;
+ x = (x | (x >> 8)) & 0x0000FFFF;
+ return (uint16_t)x;
+}
+
+static res_T
+solve_pixel
+ (struct sdis_scene* scn,
+ struct ssp_rng* rng,
+ const struct sdis_medium* mdm,
+ const struct sdis_camera* cam,
+ const double time, /* Observation time */
+ const double fp_to_meter, /* Scale from floating point units to meters */
+ const double Tarad, /* In Kelvin */
+ const double Tref, /* In Kelvin */
+ const size_t ipix[2], /* Pixel coordinate in the image plane */
+ const size_t nrealisations,
+ const double pix_sz[2], /* Pixel size in the normalized image plane */
+ struct sdis_mc* estimation)
+{
+ double weight = 0;
+ double sqr_weight = 0;
+ size_t N = 0; /* #realisations that do not fail */
+ size_t irealisation;
+ res_T res = RES_OK;
+ ASSERT(scn && mdm && rng && cam && ipix && nrealisations);
+ ASSERT(pix_sz && pix_sz[0] > 0 && pix_sz[1] > 0);
+
+ FOR_EACH(irealisation, 0, nrealisations) {
+ double samp[2]; /* Pixel sample */
+ double ray_pos[3];
+ double ray_dir[3];
+ double w = 0;
+
+ /* Generate a sample into the pixel to estimate */
+ samp[0] = (double)ipix[0] + ssp_rng_uniform_double(rng, 0, pix_sz[0]);
+ samp[1] = (double)ipix[1] + ssp_rng_uniform_double(rng, 0, pix_sz[1]);
+
+ /* Generate a ray starting from the camera position and passing through
+ * pixel sample */
+ camera_ray(cam, samp, ray_pos, ray_dir);
+
+ /* Launch the realisation */
+ res = ray_realisation_3d(scn, rng, mdm, ray_pos, ray_dir,
+ time, fp_to_meter, Tarad, Tref, &w);
+ if(res != RES_OK) goto error;
+
+ if(res == RES_OK) {
+ weight += w;
+ sqr_weight += w*w;
+ ++N;
+ } else if(res != RES_BAD_OP) {
+ goto error;
+ }
+ }
+
+ if(N == 0) {
+ estimation->E = NaN;
+ estimation->SE = NaN;
+ estimation->V = NaN;
+ } else {
+ estimation->E = weight/(double)N;
+ estimation->V = sqr_weight/(double)N - estimation->E*estimation->E;
+ estimation->SE = sqrt(estimation->V/(double)N);
+ }
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
+static res_T
+solve_tile
+ (struct sdis_scene* scn,
+ struct ssp_rng* rng,
+ const struct sdis_medium* mdm,
+ const struct sdis_camera* cam,
+ const double time,
+ const double fp_to_meter,
+ const double Tarad,
+ const double Tref,
+ const size_t origin[2], /* Tile origin in image plane */
+ const size_t size[2], /* #pixels in X and Y */
+ const size_t spp, /* #samples per pixel */
+ const double pix_sz[2], /* Pixel size in the normalized image plane */
+ struct sdis_mc* estimations)
+{
+ size_t mcode; /* Morton code of the tile pixel */
+ size_t npixels;
+ res_T res = RES_OK;
+ ASSERT(scn && rng && mdm && cam && spp && origin && estimations);
+ ASSERT(size &&size[0] && size[1]);
+ ASSERT(pix_sz && pix_sz[0] > 0 && pix_sz[1] > 0);
+
+ /* Adjust the #pixels to process them wrt a morton order */
+ npixels = round_up_pow2(MMAX(size[0], size[1]));
+ npixels *= npixels;
+
+ FOR_EACH(mcode, 0, npixels) {
+ size_t ipix[2];
+ struct sdis_mc* estimation;
+
+ ipix[0] = morton2D_decode((uint32_t)(mcode>>0));
+ if(ipix[0] >= size[0]) continue;
+ ipix[1] = morton2D_decode((uint32_t)(mcode>>1));
+ if(ipix[1] >= size[1]) continue;
+
+ estimation = estimations + ipix[1]*size[0] + ipix[0];
+ ipix[0] = ipix[0] + origin[0];
+ ipix[1] = ipix[1] + origin[1];
+
+ res = solve_pixel(scn, rng, mdm, cam, time, fp_to_meter, Tarad, Tref, size,
+ spp, pix_sz, estimation);
+ if(res != RES_OK) goto error;
+ }
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
+/*******************************************************************************
+ * Exported functions
+ ******************************************************************************/
+res_T
+sdis_solve_probe
+ (struct sdis_scene* scn,
+ const size_t nrealisations,
+ const double position[3],
+ const double time,
+ const double fp_to_meter,/* Scale factor from floating point unit to meter */
+ const double Tarad, /* Ambient radiative temperature */
+ const double Tref, /* Reference temperature */
+ struct sdis_estimator** out_estimator)
+{
+ const struct sdis_medium* medium = NULL;
+ struct sdis_estimator* estimator = NULL;
+ struct ssp_rng_proxy* rng_proxy = NULL;
+ struct ssp_rng** rngs = NULL;
+ double weight = 0;
+ double sqr_weight = 0;
+ size_t irealisation = 0;
+ size_t N = 0; /* #realisations that do not fail */
+ size_t i;
+ ATOMIC res = RES_OK;
+
+ if(!scn || !nrealisations || !position || time < 0 || fp_to_meter <= 0
+ || Tref < 0 || !out_estimator) {
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ /* Create the proxy RNG */
+ res = ssp_rng_proxy_create(scn->dev->allocator, &ssp_rng_mt19937_64,
+ scn->dev->nthreads, &rng_proxy);
+ if(res != RES_OK) goto error;
+
+ /* Create the per thread RNG */
+ rngs = MEM_CALLOC
+ (scn->dev->allocator, scn->dev->nthreads, sizeof(struct ssp_rng*));
+ if(!rngs) {
+ res = RES_MEM_ERR;
+ goto error;
+ }
+ FOR_EACH(i, 0, scn->dev->nthreads) {
+ res = ssp_rng_proxy_create_rng(rng_proxy, i, rngs+i);
+ if(res != RES_OK) goto error;
+ }
+
+ /* Create the estimator */
+ res = estimator_create(scn->dev, &estimator);
+ if(res != RES_OK) goto error;
+
+ /* Retrieve the medium in which the submitted position lies */
+ res = scene_get_medium(scn, position, &medium);
+ if(res != RES_OK) goto error;
+
+ /* Here we go! Launch the Monte Carlo estimation */
+ #pragma omp parallel for schedule(static) reduction(+:weight,sqr_weight,N)
+ for(irealisation = 0; irealisation < nrealisations; ++irealisation) {
+ res_T res_local;
+ double w;
+ const int ithread = omp_get_thread_num();
+ struct ssp_rng* rng = rngs[ithread];
+
+ if(ATOMIC_GET(&res) != RES_OK) continue; /* An error occured */
+
+ if(scene_is_2d(scn)) {
+ res_local = probe_realisation_2d
+ (scn, rng, medium, position, time, fp_to_meter, Tarad, Tref, &w);
+ } else {
+ res_local = probe_realisation_3d
+ (scn, rng, medium, position, time, fp_to_meter, Tarad, Tref, &w);
+ }
+ if(res_local != RES_OK) {
+ if(res_local != RES_BAD_OP) {
+ ATOMIC_SET(&res, res_local);
+ continue;
+ }
+ } else {
+ weight += w;
+ sqr_weight += w*w;
+ ++N;
+ }
+ }
+
+ estimator->nrealisations = N;
+ estimator->nfailures = nrealisations - N;
+ estimator->temperature.E = weight / (double)N;
+ estimator->temperature.V =
+ sqr_weight / (double)N
+ - estimator->temperature.E * estimator->temperature.E;
+ estimator->temperature.SE = sqrt(estimator->temperature.V / (double)N);
+
+exit:
+ if(rngs) {
+ FOR_EACH(i, 0, scn->dev->nthreads) {
+ if(rngs[i]) SSP(rng_ref_put(rngs[i]));
+ }
+ MEM_RM(scn->dev->allocator, rngs);
+ }
+ if(rng_proxy) SSP(rng_proxy_ref_put(rng_proxy));
+ if(out_estimator) *out_estimator = estimator;
+ return (res_T)res;
+error:
+ if(estimator) {
+ SDIS(estimator_ref_put(estimator));
+ estimator = NULL;
+ }
+ goto exit;
+}
+
+res_T
+sdis_solve_camera
+ (struct sdis_scene* scn,
+ const struct sdis_camera* cam,
+ const double time,
+ const double fp_to_meter, /* Scale from floating point units to meters */
+ const double Tarad, /* In Kelvin */
+ const double Tref, /* In Kelvin */
+ const size_t width, /* #pixels in X */
+ const size_t height, /* #pixels in Y */
+ const size_t spp, /* #samples per pixel */
+ sdis_write_estimations_T writer,
+ void* writer_data)
+{
+ #define TILE_SIZE 32 /* definition in X & Y of a tile */
+ STATIC_ASSERT(IS_POW2(TILE_SIZE), TILE_SIZE_must_be_a_power_of_2);
+
+ const struct sdis_medium* medium = NULL;
+ struct darray_mc* tiles = NULL;
+ struct ssp_rng_proxy* rng_proxy = NULL;
+ struct ssp_rng** rngs = NULL;
+ size_t ntiles_x, ntiles_y, ntiles;
+ double pix_sz[2]; /* Size of a pixel in the normalized image plane */
+ int64_t mcode; /* Morton code of a tile */
+ size_t i;
+ ATOMIC res = RES_OK;
+
+ if(!scn || !cam || time < 0 || fp_to_meter <= 0 || Tref < 0 || !width
+ || !height || !spp || !writer) {
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ if(scene_is_2d(scn)) {
+ log_err(scn->dev, "%s: 2D scene are not supported.\n", FUNC_NAME);
+ goto error;
+ }
+
+ /* Retrieve the medium in which the submitted position lies */
+ res = scene_get_medium(scn, cam->position, &medium);
+ if(res != RES_OK) goto error;
+
+ if(medium != SDIS_MEDIUM_FLUID) {
+ log_err(scn->dev, "%s: the camera position `%g %g %g' is not in a fluid.\n",
+ FUNC_NAME, SPLIT3(cam->position));
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ /* Create the proxy RNG */
+ res = ssp_rng_proxy_create(scn->dev->allocator, &ssp_rng_mt19937_64,
+ scn->dev->nthreads, &rng_proxy);
+ if(res != RES_OK) goto error;
+
+ /* Create the per thread RNG */
+ rngs = MEM_CALLOC
+ (scn->dev->allocator, scn->dev->nthreads, sizeof(struct ssp_rng*));
+ if(!rngs) {
+ res = RES_MEM_ERR;
+ goto error;
+ }
+ FOR_EACH(i, 0, scn->dev->nthreads) {
+ res = ssp_rng_proxy_create_rng(rng_proxy, i, rngs+i);
+ if(res != RES_OK) goto error;
+ }
+
+ /* Allocate per thread buffer of estimations */
+ tiles = darray_tile_data_get(&scn->dev->tiles);
+ ASSERT(darray_tile_size_get(&scn->dev->tiles) == scn->dev->nthreads);
+ FOR_EACH(i, 0, scn->dev->nthreads) {
+ const size_t nestimations = TILE_SIZE * TILE_SIZE;
+ res = darray_mc_resize(tiles+i, nestimations);
+ if(res != RES_OK) goto error;
+ }
+
+ ntiles_x = (width + (TILE_SIZE-1)/*ceil*/)/TILE_SIZE;
+ ntiles_y = (height + (TILE_SIZE-1)/*ceil*/)/TILE_SIZE;
+ ntiles = round_up_pow2(MMAX(ntiles_x, ntiles_y));
+ ntiles *= ntiles;
+
+ pix_sz[0] = 1.0 / (double)width;
+ pix_sz[1] = 1.0 / (double)height;
+
+ omp_set_num_threads((int)scn->dev->nthreads);
+ /*#pragma omp parallel for schedule(static)*/
+ for(mcode = 0; mcode < (int64_t)ntiles; ++mcode) {
+ size_t tile_org[2] = {0, 0};
+ size_t tile_sz[2] = {0, 0};
+ const int ithread = omp_get_thread_num();
+ struct sdis_mc* estimations = NULL;
+ struct ssp_rng* rng = rngs[ithread];
+ res_T res_local = RES_OK;
+
+ if(ATOMIC_GET(&res) != RES_OK) continue;
+
+ tile_org[0] = morton2D_decode((uint32_t)(mcode>>0));
+ if(tile_org[0] >= ntiles_x) continue; /* Discard tile */
+ tile_org[1] = morton2D_decode((uint32_t)(mcode>>1));
+ if(tile_org[1] >= ntiles_y) continue; /* Disaard tile */
+
+ /* Setup the tile coordinates in the image plane */
+ tile_org[0] *= TILE_SIZE;
+ tile_org[1] *= TILE_SIZE;
+ tile_sz[0] = MMIN(TILE_SIZE, width - tile_org[0]);
+ tile_sz[0] = MMIN(TILE_SIZE, width - tile_org[1]);
+
+ /* Fetch the estimation buffer */
+ estimations = darray_mc_data_get(tiles+ithread);
+
+ /* Draw the tile */
+ res_local = solve_tile(scn, rng, medium, cam, time, fp_to_meter, Tarad,
+ Tref, tile_org, tile_sz, spp, pix_sz, estimations);
+ if(res_local != RES_OK) {
+ ATOMIC_SET(&res, res_local);
+ continue;
+ }
+
+ /* Write the estimations */
+ res_local = writer(writer_data, tile_org, tile_sz, estimations);
+ if(res_local != RES_OK) {
+ ATOMIC_SET(&res, res_local);
+ continue;
+ }
+ }
+
+exit:
+ if(rngs) {
+ FOR_EACH(i, 0, scn->dev->nthreads) {
+ if(rngs[i]) SSP(rng_ref_put(rngs[i]));
+ }
+ MEM_RM(scn->dev->allocator, rngs);
+ }
+ if(rng_proxy) SSP(rng_proxy_ref_put(rng_proxy));
+ return (res_T)res;
+error:
+ goto exit;
+}
+
diff --git a/src/sdis_solve_Xd.h b/src/sdis_solve_Xd.h
@@ -0,0 +1,826 @@
+/* Copyright (C) |Meso|Star> 2016-2018 (contact@meso-star.com)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. */
+
+#ifndef SDIS_SOLVE_DIMENSION
+#ifndef SDIS_SOLVE_XD_H
+#define SDIS_SOLVE_XD_H
+
+#include "sdis_device_c.h"
+#include "sdis_interface_c.h"
+#include "sdis_medium_c.h"
+#include "sdis_scene_c.h"
+
+#include <rsys/stretchy_array.h>
+
+#include <star/ssp.h>
+
+/* Emperical scale factor to apply to the upper bound of the ray range in order
+ * to handle numerical imprecisions */
+#define RAY_RANGE_MAX_SCALE 1.0001f
+
+#define BOLTZMANN_CONSTANT 5.6696e-8 /* W/m^2/K^4 */
+
+struct rwalk_context {
+ double Tarad; /* Ambient radiative temperature */
+ double Tref3; /* Reference temperature ^ 3 */
+};
+
+/* Reflect the vector V wrt the normal N. By convention V points outward the
+ * surface. */
+static INLINE float*
+reflect(float res[3], const float V[3], const float N[3])
+{
+ float tmp[3];
+ float cos_V_N;
+ ASSERT(res && V && N);
+ ASSERT(f3_is_normalized(V) && f3_is_normalized(N));
+ cos_V_N = f3_dot(V, N);
+ f3_mulf(tmp, N, 2*cos_V_N);
+ f3_sub(res, tmp, V);
+ return res;
+}
+
+#endif /* SDIS_SOLVE_XD_H */
+#else
+
+#if (SDIS_SOLVE_DIMENSION == 2)
+ #include <rsys/double2.h>
+ #include <rsys/float2.h>
+ #include <star/s2d.h>
+#elif (SDIS_SOLVE_DIMENSION == 3)
+ #include <rsys/double2.h>
+ #include <rsys/double3.h>
+ #include <rsys/float3.h>
+ #include <star/s3d.h>
+#else
+ #error "Invalid SDIS_SOLVE_DIMENSION value."
+#endif
+
+/* Syntactic sugar */
+#define DIM SDIS_SOLVE_DIMENSION
+
+/* Star-XD macros generic to SDIS_SOLVE_DIMENSION */
+#define sXd(Name) CONCAT(CONCAT(CONCAT(s, DIM), d_), Name)
+#define SXD_HIT_NONE CONCAT(CONCAT(S,DIM), D_HIT_NONE)
+#define SXD_HIT_NULL CONCAT(CONCAT(S,DIM), D_HIT_NULL)
+#define SXD_HIT_NULL__ CONCAT(CONCAT(S, DIM), D_HIT_NULL__)
+#define SXD CONCAT(CONCAT(S, DIM), D)
+
+/* Vector macros generic to SDIS_SOLVE_DIMENSION */
+#define dX(Func) CONCAT(CONCAT(CONCAT(d, DIM), _), Func)
+#define fX(Func) CONCAT(CONCAT(CONCAT(f, DIM), _), Func)
+#define fX_set_dX CONCAT(CONCAT(CONCAT(f, DIM), _set_d), DIM)
+#define dX_set_fX CONCAT(CONCAT(CONCAT(d, DIM), _set_f), DIM)
+
+/* Macro making generic its subimitted name to SDIS_SOLVE_DIMENSION */
+#define XD(Name) CONCAT(CONCAT(CONCAT(Name, _), DIM), d)
+
+/* Current state of the random walk */
+struct XD(rwalk) {
+ struct sdis_rwalk_vertex vtx; /* Position and time of the Random walk */
+ const struct sdis_medium* mdm; /* Medium in which the random walk lies */
+ struct sXd(hit) hit; /* Hit of the random walk */
+};
+static const struct XD(rwalk) XD(RWALK_NULL) = {
+ SDIS_RWALK_VERTEX_NULL__, NULL, SXD_HIT_NULL__
+};
+
+struct XD(temperature) {
+ res_T (*func)/* Next function to invoke in order to compute the temperature */
+ (const struct sdis_scene* scn,
+ const double fp_to_meter,
+ const struct rwalk_context* ctx,
+ struct XD(rwalk)* rwalk,
+ struct ssp_rng* rng,
+ struct XD(temperature)* temp);
+ double value; /* Current value of the temperature */
+ int done;
+};
+static const struct XD(temperature) XD(TEMPERATURE_NULL) = { NULL, 0, 0 };
+
+static res_T
+XD(boundary_temperature)
+ (const struct sdis_scene* scn,
+ const double fp_to_meter,
+ const struct rwalk_context* ctx,
+ struct XD(rwalk)* rwalk,
+ struct ssp_rng* rng,
+ struct XD(temperature)* T);
+
+static res_T
+XD(solid_temperature)
+ (const struct sdis_scene* scn,
+ const double fp_to_meter,
+ const struct rwalk_context* ctx,
+ struct XD(rwalk)* rwalk,
+ struct ssp_rng* rng,
+ struct XD(temperature)* T);
+
+static res_T
+XD(fluid_temperature)
+ (const struct sdis_scene* scn,
+ const double fp_to_meter,
+ const struct rwalk_context* ctx,
+ struct XD(rwalk)* rwalk,
+ struct ssp_rng* rng,
+ struct XD(temperature)* T);
+
+static res_T
+XD(radiative_temperature)
+ (const struct sdis_scene* scn,
+ const double fp_to_meter,
+ const struct rwalk_context* ctx,
+ struct XD(rwalk)* rwalk,
+ struct ssp_rng* rng,
+ struct XD(temperature)* T);
+
+/*******************************************************************************
+ * Helper functions
+ ******************************************************************************/
+
+static FINLINE void
+XD(move_pos)(double pos[DIM], const float dir[DIM], const float delta)
+{
+ ASSERT(pos && dir);
+ pos[0] += dir[0] * delta;
+ pos[1] += dir[1] * delta;
+#if(SDIS_SOLVE_DIMENSION == 3)
+ pos[2] += dir[2] * delta;
+#endif
+}
+
+/* Check that the interface fragment is consistent with the current state of
+ * the random walk */
+static INLINE int
+XD(check_rwalk_fragment_consistency)
+ (const struct XD(rwalk)* rwalk,
+ const struct sdis_interface_fragment* frag)
+{
+ double N[DIM];
+ double uv[2] = {0, 0};
+ ASSERT(rwalk && frag);
+ dX(normalize)(N, dX_set_fX(N, rwalk->hit.normal));
+ if( SXD_HIT_NONE(&rwalk->hit)
+ || !dX(eq_eps)(rwalk->vtx.P, frag->P, 1.e-6)
+ || !dX(eq_eps)(N, frag->Ng, 1.e-6)
+ || !( (IS_INF(rwalk->vtx.time) && IS_INF(frag->time))
+ || eq_eps(rwalk->vtx.time, frag->time, 1.e-6))) {
+ return 0;
+ }
+#if (SDIS_SOLVE_DIMENSION == 2)
+ uv[0] = rwalk->hit.u;
+#else
+ d2_set_f2(uv, rwalk->hit.uv);
+#endif
+ return d2_eq_eps(uv, frag->uv, 1.e-6);
+}
+
+res_T
+XD(radiative_temperature)
+ (const struct sdis_scene* scn,
+ const double fp_to_meter,
+ const struct rwalk_context* ctx,
+ struct XD(rwalk)* rwalk,
+ struct ssp_rng* rng,
+ struct XD(temperature)* T)
+{
+ const struct sdis_interface* interf;
+
+ /* The radiative random walk is always perform in 3D. In 2D, the geometry are
+ * assumed to be extruded to the infinty along the Z dimension. */
+ float N[3] = {0, 0, 0};
+ float dir[3] = {0, 0, 0};
+
+ res_T res = RES_OK;
+
+ ASSERT(scn && fp_to_meter > 0 && ctx && rwalk && rng && T);
+ ASSERT(!SXD_HIT_NONE(&rwalk->hit));
+ (void)fp_to_meter;
+
+ /* Fetch the current interface */
+ interf = scene_get_interface(scn, rwalk->hit.prim.prim_id);
+
+ /* Normalize the normal of the interface and ensure that it points toward the
+ * current medium */
+ fX(normalize(N, rwalk->hit.normal));
+ if(interf->medium_back == rwalk->mdm) {
+ fX(minus(N, N));
+ }
+
+ /* Cosine weighted sampling of a direction around the surface normal */
+ ssp_ran_hemisphere_cos_float(rng, N, dir, NULL);
+
+ /* Launch the radiative random walk */
+ for(;;) {
+ struct sdis_interface_fragment frag = SDIS_INTERFACE_FRAGMENT_NULL;
+ const struct sdis_medium* chk_mdm = NULL;
+ double alpha;
+ double epsilon;
+ double r;
+ float pos[DIM];
+ const float range[2] = { 0, FLT_MAX };
+
+ fX_set_dX(pos, rwalk->vtx.P);
+
+ /* Trace the radiative ray */
+#if (SDIS_SOLVE_DIMENSION == 2)
+ SXD(scene_view_trace_ray_3d
+ (scn->sXd(view), pos, dir, range, &rwalk->hit, &rwalk->hit));
+#else
+ SXD(scene_view_trace_ray
+ (scn->sXd(view), pos, dir, range, &rwalk->hit, &rwalk->hit));
+#endif
+ if(SXD_HIT_NONE(&rwalk->hit)) { /* Fetch the ambient radiative temperature */
+ if(ctx->Tarad >= 0) {
+ T->value += ctx->Tarad;
+ T->done = 1;
+ break;
+ } else {
+ log_err(scn->dev,
+"%s: the random walk reaches an invalid ambient radiative temperature of `%gK'\n"
+"at position `%g %g %g'. This may be due to numerical inaccuracies or to\n"
+"inconsistency in the simulated system (eg: unclosed geometry). For systems\n"
+"where the random walks can reach such temperature, one has to setup a valid\n"
+"ambient radiative temperature, i.e. it must be greater or equal to 0.\n",
+ FUNC_NAME,
+ ctx->Tarad,
+ SPLIT3(rwalk->vtx.P));
+ res = RES_BAD_ARG;
+ goto error;
+ }
+ }
+
+ /* Move the random walk to the hit position */
+ XD(move_pos)(rwalk->vtx.P, dir, rwalk->hit.distance);
+
+ /* Fetch the new interface and setup the hit fragment */
+ interf = scene_get_interface(scn, rwalk->hit.prim.prim_id);
+ XD(setup_interface_fragment)(&frag, &rwalk->vtx, &rwalk->hit);
+
+ /* Fetch the interface emissivity */
+ epsilon = interface_get_emissivity(interf, &frag);
+ if(epsilon > 1 && epsilon >= 0) {
+ log_err(scn->dev,
+ "%s: invalid overall emissivity `%g' at position `%g %g %g'.\n",
+ FUNC_NAME, epsilon, SPLIT3(rwalk->vtx.P));
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ /* Switch in boundary temperature ? */
+ r = ssp_rng_canonical(rng);
+ if(r < epsilon) {
+ T->func = XD(boundary_temperature);
+ break;
+ }
+
+ /* Normalize the normal of the interface and ensure that it points toward the
+ * current medium */
+ fX(normalize)(N, rwalk->hit.normal);
+ if(f3_dot(N, dir) > 0) {
+ chk_mdm = interf->medium_back;
+ fX(minus)(N, N);
+ } else {
+ chk_mdm = interf->medium_front;
+ }
+
+ if(chk_mdm != rwalk->mdm) {
+ log_err(scn->dev, "%s: inconsistent medium definition at `%g %g %g'.\n",
+ FUNC_NAME, SPLIT3(rwalk->vtx.P));
+ res = RES_BAD_ARG;
+ goto error;
+ }
+ alpha = interface_get_specular_fraction(interf, &frag);
+ r = ssp_rng_canonical(rng);
+ if(r < alpha) { /* Sample specular part */
+ reflect(dir, f3_minus(dir, dir), N);
+ } else { /* Sample diffuse part */
+ ssp_ran_hemisphere_cos_float(rng, N, dir, NULL);
+ }
+ }
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
+res_T
+XD(fluid_temperature)
+ (const struct sdis_scene* scn,
+ const double fp_to_meter,
+ const struct rwalk_context* ctx,
+ struct XD(rwalk)* rwalk,
+ struct ssp_rng* rng,
+ struct XD(temperature)* T)
+{
+ double tmp;
+ (void)rng, (void)fp_to_meter, (void)ctx;
+ ASSERT(scn && fp_to_meter > 0 && ctx && rwalk && rng && T);
+ ASSERT(rwalk->mdm->type == SDIS_MEDIUM_FLUID);
+
+ tmp = fluid_get_temperature(rwalk->mdm, &rwalk->vtx);
+ if(tmp < 0) {
+ log_err(scn->dev, "%s: unknown fluid temperature is not supported.\n",
+ FUNC_NAME);
+ return RES_BAD_OP;
+ }
+ T->value += tmp;
+ T->done = 1;
+ return RES_OK;
+}
+
+static void
+XD(solid_solid_boundary_temperature)
+ (const struct sdis_scene* scn,
+ const double fp_to_meter,
+ const struct rwalk_context* ctx,
+ const struct sdis_interface_fragment* frag,
+ struct XD(rwalk)* rwalk,
+ struct ssp_rng* rng,
+ struct XD(temperature)* T)
+{
+ const struct sdis_interface* interf = NULL;
+ const struct sdis_medium* solid_front = NULL;
+ const struct sdis_medium* solid_back = NULL;
+ double lambda_front, lambda_back;
+ double delta_front_boundary, delta_back_boundary;
+ double delta_front_boundary_meter, delta_back_boundary_meter;
+ double delta_boundary;
+ double proba;
+ double tmp;
+ double r;
+ float pos[DIM], dir[DIM], range[2];
+ ASSERT(scn && fp_to_meter > 0 && ctx && frag && rwalk && rng && T);
+ ASSERT(XD(check_rwalk_fragment_consistency)(rwalk, frag));
+ (void)frag, (void)ctx;
+
+ /* Retrieve the current boundary media */
+ interf = scene_get_interface(scn, rwalk->hit.prim.prim_id);
+ solid_front = interface_get_medium(interf, SDIS_FRONT);
+ solid_back = interface_get_medium(interf, SDIS_BACK);
+ ASSERT(solid_front->type == SDIS_MEDIUM_SOLID);
+ ASSERT(solid_back->type == SDIS_MEDIUM_SOLID);
+
+ /* Fetch the properties of the media */
+ lambda_front = solid_get_thermal_conductivity(solid_front, &rwalk->vtx);
+ lambda_back = solid_get_thermal_conductivity(solid_back, &rwalk->vtx);
+ delta_front_boundary = solid_get_delta_boundary(solid_front, &rwalk->vtx);
+ delta_back_boundary = solid_get_delta_boundary(solid_back, &rwalk->vtx);
+
+ /* Convert the delta boundary from floating point units to meters */
+ delta_front_boundary_meter = fp_to_meter * delta_front_boundary;
+ delta_back_boundary_meter = fp_to_meter * delta_back_boundary;
+
+ /* Compute the proba to switch in solid0 or in solid1 */
+ tmp = lambda_front / delta_front_boundary_meter;
+ proba = tmp / (tmp + lambda_back / delta_back_boundary_meter);
+
+ r = ssp_rng_canonical(rng);
+ fX(normalize)(dir, rwalk->hit.normal);
+ if(r < proba) { /* Reinject in front */
+ delta_boundary = delta_front_boundary;
+ rwalk->mdm = solid_front;
+ } else { /* Reinject in back */
+ delta_boundary = delta_back_boundary;
+ fX(minus)(dir, dir);
+ rwalk->mdm = solid_back;
+ }
+
+ /* "Reinject" the path into the solid along the surface normal. */
+ fX_set_dX(pos, rwalk->vtx.P);
+ range[0] = 0, range[1] = (float)delta_boundary*RAY_RANGE_MAX_SCALE;
+ SXD(scene_view_trace_ray
+ (scn->sXd(view), pos, dir, range, &rwalk->hit, &rwalk->hit));
+ if(!SXD_HIT_NONE(&rwalk->hit)) delta_boundary = rwalk->hit.distance * 0.5;
+ XD(move_pos)(rwalk->vtx.P, dir, (float)delta_boundary);
+
+ /* Switch in solid random walk */
+ T->func = XD(solid_temperature);
+}
+
+static void
+XD(solid_fluid_boundary_temperature)
+ (const struct sdis_scene* scn,
+ const double fp_to_meter,
+ const struct rwalk_context* ctx,
+ const struct sdis_interface_fragment* frag,
+ struct XD(rwalk)* rwalk,
+ struct ssp_rng* rng,
+ struct XD(temperature)* T)
+{
+ const struct sdis_interface* interf = NULL;
+ const struct sdis_medium* mdm_front = NULL;
+ const struct sdis_medium* mdm_back = NULL;
+ const struct sdis_medium* solid = NULL;
+ const struct sdis_medium* fluid = NULL;
+ double hc;
+ double hr;
+ double epsilon; /* Interface emissivity */
+ double lambda;
+ double fluid_proba;
+ double radia_proba;
+ double delta_boundary;
+ double r;
+ double tmp;
+ float dir[DIM], pos[DIM], range[2];
+
+ ASSERT(scn && fp_to_meter > 0 && rwalk && rng && T && ctx);
+ ASSERT(XD(check_rwalk_fragment_consistency)(rwalk, frag));
+
+ /* Retrieve the solid and the fluid split by the boundary */
+ interf = scene_get_interface(scn, rwalk->hit.prim.prim_id);
+ mdm_front = interface_get_medium(interf, SDIS_FRONT);
+ mdm_back = interface_get_medium(interf, SDIS_BACK);
+ ASSERT(mdm_front->type != mdm_back->type);
+ if(mdm_front->type == SDIS_MEDIUM_SOLID) {
+ solid = mdm_front;
+ fluid = mdm_back;
+ } else {
+ solid = mdm_back;
+ fluid = mdm_front;
+ }
+
+ /* Fetch the solid properties */
+ lambda = solid_get_thermal_conductivity(solid, &rwalk->vtx);
+ delta_boundary = solid_get_delta_boundary(solid, &rwalk->vtx);
+
+ /* Fetch the boundary properties */
+ epsilon = interface_get_emissivity(interf, frag);
+ hc = interface_get_convection_coef(interf, frag);
+
+ /* Compute the radiative coefficient */
+ hr = 4.0 * BOLTZMANN_CONSTANT * ctx->Tref3 * epsilon;
+
+ /* Compute the probas to switch in solid or fluid random walk */
+ tmp = lambda / (delta_boundary*fp_to_meter);
+ fluid_proba = hc / (tmp + hr + hc);
+ radia_proba = hr / (tmp + hr + hc);
+ /*solid_proba = tmp / (tmp + hr + hc);*/
+
+ r = ssp_rng_canonical(rng);
+ if(r < radia_proba) { /* Switch in radiative random walk */
+ rwalk->mdm = fluid;
+ T->func = XD(radiative_temperature);
+ } else if(r < fluid_proba + radia_proba) { /* Switch to fluid random walk */
+ rwalk->mdm = fluid;
+ T->func = XD(fluid_temperature);
+ } else { /* Solid random walk */
+ rwalk->mdm = solid;
+ fX(normalize)(dir, rwalk->hit.normal);
+ if(solid == mdm_back) fX(minus)(dir, dir);
+
+ /* "Reinject" the random walk into the solid */
+ fX_set_dX(pos, rwalk->vtx.P);
+ range[0] = 0, range[1] = (float)delta_boundary*RAY_RANGE_MAX_SCALE;
+ SXD(scene_view_trace_ray
+ (scn->sXd(view), pos, dir, range, &rwalk->hit, &rwalk->hit));
+ if(!SXD_HIT_NONE(&rwalk->hit)) delta_boundary = rwalk->hit.distance * 0.5;
+ XD(move_pos)(rwalk->vtx.P, dir, (float)delta_boundary);
+
+ /* Switch in solid random walk */
+ T->func = XD(solid_temperature);
+ }
+}
+
+res_T
+XD(boundary_temperature)
+ (const struct sdis_scene* scn,
+ const double fp_to_meter,
+ const struct rwalk_context* ctx,
+ struct XD(rwalk)* rwalk,
+ struct ssp_rng* rng,
+ struct XD(temperature)* T)
+{
+ struct sdis_interface_fragment frag = SDIS_INTERFACE_FRAGMENT_NULL;
+ const struct sdis_interface* interf = NULL;
+ const struct sdis_medium* mdm_front = NULL;
+ const struct sdis_medium* mdm_back = NULL;
+ double tmp;
+ ASSERT(scn && fp_to_meter > 0 && ctx && rwalk && rng && T);
+ ASSERT(!SXD_HIT_NONE(&rwalk->hit));
+
+ XD(setup_interface_fragment)(&frag, &rwalk->vtx, &rwalk->hit);
+
+ /* Retrieve the current interface */
+ interf = scene_get_interface(scn, rwalk->hit.prim.prim_id);
+
+ /* Check if the boundary condition is known */
+ tmp = interface_get_temperature(interf, &frag);
+ if(tmp >= 0) {
+ T->value += tmp;
+ T->done = 1;
+ return RES_OK;
+ }
+
+ mdm_front = interface_get_medium(interf, SDIS_FRONT);
+ mdm_back = interface_get_medium(interf, SDIS_BACK);
+
+ if(mdm_front->type == mdm_back->type) {
+ XD(solid_solid_boundary_temperature)
+ (scn, fp_to_meter, ctx, &frag, rwalk, rng, T);
+ } else {
+ XD(solid_fluid_boundary_temperature)
+ (scn, fp_to_meter, ctx, &frag, rwalk, rng, T);
+ }
+ return RES_OK;
+}
+
+res_T
+XD(solid_temperature)
+ (const struct sdis_scene* scn,
+ const double fp_to_meter,
+ const struct rwalk_context* ctx,
+ struct XD(rwalk)* rwalk,
+ struct ssp_rng* rng,
+ struct XD(temperature)* T)
+{
+ double position_start[DIM];
+ const struct sdis_medium* mdm;
+ ASSERT(scn && fp_to_meter > 0 && rwalk && rng && T);
+ ASSERT(rwalk->mdm->type == SDIS_MEDIUM_SOLID);
+ (void)ctx;
+
+ /* Check the random walk consistency */
+ CHK(scene_get_medium(scn, rwalk->vtx.P, &mdm) == RES_OK);
+ if(mdm != rwalk->mdm) {
+ log_err(scn->dev, "%s: invalid solid random walk. "
+ "Unexpected medium at {%g, %g, %g}.\n",
+ FUNC_NAME, SPLIT3(rwalk->vtx.P));
+ return RES_BAD_ARG;
+ }
+ /* Save the submitted position */
+ dX(set)(position_start, rwalk->vtx.P);
+
+ do { /* Solid random walk */
+ struct sXd(hit) hit0, hit1;
+ double lambda; /* Thermal conductivity */
+ double rho; /* Volumic mass */
+ double cp; /* Calorific capacity */
+ double tau, mu;
+ double tmp;
+ float delta, delta_solid; /* Random walk numerical parameter */
+ float range[2];
+ float dir0[DIM], dir1[DIM];
+ float org[DIM];
+
+ /* Check the limit condition */
+ tmp = solid_get_temperature(mdm, &rwalk->vtx);
+ if(tmp >= 0) {
+ T->value += tmp;
+ T->done = 1;
+ return RES_OK;
+ }
+
+ /* Fetch solid properties */
+ delta_solid = (float)solid_get_delta(mdm, &rwalk->vtx);
+ lambda = solid_get_thermal_conductivity(mdm, &rwalk->vtx);
+ rho = solid_get_volumic_mass(mdm, &rwalk->vtx);
+ cp = solid_get_calorific_capacity(mdm, &rwalk->vtx);
+
+#if (SDIS_SOLVE_DIMENSION == 2)
+ /* Sample a direction around 2PI */
+ ssp_ran_circle_uniform_float(rng, dir0, NULL);
+#else
+ /* Sample a direction around 4PI */
+ ssp_ran_sphere_uniform_float(rng, dir0, NULL);
+#endif
+
+ /* Trace a ray along the sampled direction and its opposite to check if a
+ * surface is hit in [0, delta_solid]. */
+ fX_set_dX(org, rwalk->vtx.P);
+ fX(minus)(dir1, dir0);
+ hit0 = hit1 = SXD_HIT_NULL;
+ range[0] = 0.f, range[1] = delta_solid*RAY_RANGE_MAX_SCALE;
+ SXD(scene_view_trace_ray(scn->sXd(view), org, dir0, range, NULL, &hit0));
+ SXD(scene_view_trace_ray(scn->sXd(view), org, dir1, range, NULL, &hit1));
+
+ if(SXD_HIT_NONE(&hit0) && SXD_HIT_NONE(&hit1)) {
+ /* Hit nothing: move along dir0 of the original delta */
+ delta = delta_solid;
+ } else {
+ /* Hit something: move along dir0 of the minimum hit distance */
+ delta = MMIN(hit0.distance, hit1.distance);
+ }
+
+ /* Sample the time */
+ mu = (2*DIM*lambda) / (rho*cp*delta*fp_to_meter*delta*fp_to_meter);
+ tau = ssp_ran_exp(rng, mu);
+ rwalk->vtx.time -= tau;
+
+ /* Check the initial condition */
+ tmp = solid_get_temperature(mdm, &rwalk->vtx);
+ if(tmp >= 0) {
+ T->value += tmp;
+ T->done = 1;
+ return RES_OK;
+ }
+
+ /* Define if the random walk hits something along dir0 */
+ rwalk->hit = hit0.distance > delta ? SXD_HIT_NULL : hit0;
+
+ /* Update the random walk position */
+ XD(move_pos)(rwalk->vtx.P, dir0, delta);
+
+ /* Fetch the current medium */
+ if(SXD_HIT_NONE(&rwalk->hit)) {
+ CHK(scene_get_medium(scn, rwalk->vtx.P, &mdm) == RES_OK);
+ } else {
+ const struct sdis_interface* interf;
+ interf = scene_get_interface(scn, rwalk->hit.prim.prim_id);
+ mdm = interface_get_medium
+ (interf,
+ fX(dot(rwalk->hit.normal, dir0)) < 0 ? SDIS_FRONT : SDIS_BACK);
+ }
+
+ /* Check random walk consistency */
+ if(mdm != rwalk->mdm) {
+ log_err(scn->dev,
+ "%s: inconsistent medium during the solid random walk.\n", FUNC_NAME);
+ if(DIM == 2) {
+ log_err(scn->dev,
+ " start position: %g %g; current position: %g %g\n",
+ SPLIT2(position_start), SPLIT2(rwalk->vtx.P));
+ } else {
+ log_err(scn->dev,
+ " start position: %g %g %g; current position: %g %g %g\n",
+ SPLIT3(position_start), SPLIT3(rwalk->vtx.P));
+ }
+ return RES_BAD_OP;
+ }
+
+ /* Keep going while the solid random walk does not hit an interface */
+ } while(SXD_HIT_NONE(&rwalk->hit));
+
+ T->func = XD(boundary_temperature);
+ return RES_OK;
+}
+
+static res_T
+XD(compute_temperature)
+ (struct sdis_scene* scn,
+ const double fp_to_meter,
+ const struct rwalk_context* ctx,
+ struct XD(rwalk)* rwalk,
+ struct ssp_rng* rng,
+ struct XD(temperature)* T)
+{
+#ifndef NDEBUG
+ struct XD(temperature)* stack = NULL;
+ size_t istack = 0;
+#endif
+ res_T res = RES_OK;
+ ASSERT(scn && fp_to_meter > 0 && ctx && rwalk && rng && T);
+
+ do {
+ res = T->func(scn, fp_to_meter, ctx, rwalk, rng, T);
+ if(res != RES_OK) goto error;
+
+#ifndef NDEBUG
+ sa_push(stack, *T);
+ ++istack;
+#endif
+ } while(!T->done);
+
+exit:
+#ifndef NDEBUG
+ sa_release(stack);
+#endif
+ return res;
+error:
+ goto exit;
+}
+
+static res_T
+XD(probe_realisation)
+ (struct sdis_scene* scn,
+ struct ssp_rng* rng,
+ const struct sdis_medium* medium,
+ const double position[],
+ const double time,
+ const double fp_to_meter,/* Scale factor from floating point unit to meter */
+ const double ambient_radiative_temperature,
+ const double reference_temperature,
+ double* weight)
+{
+ struct rwalk_context ctx;
+ struct XD(rwalk) rwalk = XD(RWALK_NULL);
+ struct XD(temperature) T = XD(TEMPERATURE_NULL);
+ res_T res = RES_OK;
+ ASSERT(medium && position && fp_to_meter > 0 && weight && time >= 0);
+
+ switch(medium->type) {
+ case SDIS_MEDIUM_FLUID: T.func = XD(fluid_temperature); break;
+ case SDIS_MEDIUM_SOLID: T.func = XD(solid_temperature); break;
+ default: FATAL("Unreachable code\n"); break;
+ }
+
+ dX(set)(rwalk.vtx.P, position);
+ rwalk.vtx.time = time;
+ rwalk.hit = SXD_HIT_NULL;
+ rwalk.mdm = medium;
+
+ ctx.Tarad = ambient_radiative_temperature;
+ ctx.Tref3 =
+ reference_temperature
+ * reference_temperature
+ * reference_temperature;
+
+ res = XD(compute_temperature)(scn, fp_to_meter, &ctx, &rwalk, rng, &T);
+ if(res != RES_OK) return res;
+
+ *weight = T.value;
+ return RES_OK;
+}
+
+#if SDIS_SOLVE_DIMENSION == 3
+static res_T
+XD(ray_realisation)
+ (struct sdis_scene* scn,
+ struct ssp_rng* rng,
+ const struct sdis_medium* medium,
+ const double position[],
+ const double direction[],
+ const double time,
+ const double fp_to_meter,
+ const double Tarad,
+ const double Tref,
+ double* weight)
+{
+ struct sXd(hit) hit = SXD_HIT_NULL;
+ const float range[2] = {0, FLT_MAX};
+ float org[3] = {0, 0, 0};
+ float dir[3] = {0, 0, 0};
+ res_T res = RES_OK;
+ ASSERT(scn && position && direction && time>=0 && fp_to_meter>0 && weight);
+ ASSERT(medium && medium->type == SDIS_MEDIUM_FLUID);
+
+ fX_set_dX(org, position);
+ fX_set_dX(dir, direction);
+ SXD(scene_view_trace_ray(scn->sXd(view), org, dir, range, &hit, &hit));
+
+ if(SXD_HIT_NONE(&hit)) {
+ if(Tarad >= 0) {
+ *weight = Tarad;
+ } else {
+ log_err(scn->dev,
+"%s: the ray starting from `%g %g %g' and traced along the `%g %g %g' direction\n"
+"reaches an invalid ambient temperature of `%gK'. One has to provide a valid\n"
+"ambient radiative temperature, i.e. it must be greater or equal to 0.\n",
+ FUNC_NAME, SPLIT3(org), SPLIT3(dir), Tarad);
+ res = RES_BAD_ARG;
+ goto error;
+ }
+ } else {
+ struct rwalk_context ctx;
+ struct XD(rwalk) rwalk = XD(RWALK_NULL);
+ struct XD(temperature) T = XD(TEMPERATURE_NULL);
+
+ dX(set)(rwalk.vtx.P, position);
+ XD(move_pos)(rwalk.vtx.P, dir, hit.distance);
+ rwalk.vtx.time = time;
+ rwalk.hit = hit;
+ rwalk.mdm = medium;
+
+ ctx.Tarad = Tarad;
+ ctx.Tref3 = Tref*Tref*Tref;
+
+ T.func = XD(boundary_temperature);
+ res = XD(compute_temperature)(scn, fp_to_meter, &ctx, &rwalk, rng, &T);
+ if(res != RES_OK) return res;
+
+ *weight = T.value;
+ }
+exit:
+ return res;
+error:
+ goto exit;
+}
+#endif /* SDIS_SOLVE_DIMENSION == 3 */
+
+#undef SDIS_SOLVE_DIMENSION
+#undef DIM
+#undef sXd
+#undef SXD_HIT_NONE
+#undef SXD_HIT_NULL
+#undef SXD_HIT_NULL__
+#undef SXD
+#undef dX
+#undef fX
+#undef fX_set_dX
+#undef XD
+
+#endif /* !SDIS_SOLVE_DIMENSION */
+
diff --git a/src/sdis_solve_probe.c b/src/sdis_solve_probe.c
@@ -1,413 +0,0 @@
-/* Copyright (C) |Meso|Star> 2016-2018 (contact@meso-star.com)
- *
- * This program is free software: you can redistribute it and/or modify
- * it under the terms of the GNU General Public License as published by
- * the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * This program is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU General Public License
- * along with this program. If not, see <http://www.gnu.org/licenses/>. */
-
-#include "sdis.h"
-#include "sdis_camera.h"
-#include "sdis_device_c.h"
-#include "sdis_estimator_c.h"
-#include "sdis_solve_probe_Xd.h"
-
-/* Generate the 2D solver */
-#define SDIS_SOLVE_PROBE_DIMENSION 2
-#include "sdis_solve_probe_Xd.h"
-
-/* Generate the 3D solver */
-#define SDIS_SOLVE_PROBE_DIMENSION 3
-#include "sdis_solve_probe_Xd.h"
-
-#include <star/ssp.h>
-#include <omp.h>
-
-/*******************************************************************************
- * Helper functions
- ******************************************************************************/
-static FINLINE uint16_t
-morton2D_decode(const uint32_t u32)
-{
- uint32_t x = u32 & 0x55555555;
- x = (x | (x >> 1)) & 0x33333333;
- x = (x | (x >> 2)) & 0x0F0F0F0F;
- x = (x | (x >> 4)) & 0x00FF00FF;
- x = (x | (x >> 8)) & 0x0000FFFF;
- return (uint16_t)x;
-}
-
-static res_T
-solve_pixel
- (struct sdis_scene* scn,
- struct ssp_rng* rng,
- const struct sdis_medium* mdm,
- const struct sdis_camera* cam,
- const double time, /* Observation time */
- const double fp_to_meter, /* Scale from floating point units to meters */
- const double Tarad, /* In Kelvin */
- const double Tref, /* In Kelvin */
- const size_t ipix[2], /* Pixel coordinate in the image plane */
- const size_t nrealisations,
- const double pix_sz[2], /* Pixel size in the normalized image plane */
- struct sdis_mc* estimation)
-{
- double weight = 0;
- double sqr_weight = 0;
- size_t N = 0; /* #realisations that do not fail */
- size_t irealisation;
- res_T res = RES_OK;
- ASSERT(scn && mdm && rng && cam && ipix && nrealisations);
- ASSERT(pix_sz && pix_sz[0] > 0 && pix_sz[1] > 0);
-
- FOR_EACH(irealisation, 0, nrealisations) {
- double samp[2]; /* Pixel sample */
- double ray_pos[3];
- double ray_dir[3];
- double w = 0;
-
- /* Generate a sample into the pixel to estimate */
- samp[0] = (double)ipix[0] + ssp_rng_uniform_double(rng, 0, pix_sz[0]);
- samp[1] = (double)ipix[1] + ssp_rng_uniform_double(rng, 0, pix_sz[1]);
-
- /* Generate a ray starting from the camera position and passing through
- * pixel sample */
- camera_ray(cam, samp, ray_pos, ray_dir);
-
- /* Launch the realisation */
- res = ray_realisation_3d(scn, rng, mdm, ray_pos, ray_dir,
- time, fp_to_meter, Tarad, Tref, &w);
- if(res != RES_OK) goto error;
-
- if(res == RES_OK) {
- weight += w;
- sqr_weight += w*w;
- ++N;
- } else if(res != RES_BAD_OP) {
- goto error;
- }
- }
-
- if(N == 0) {
- estimation->E = NaN;
- estimation->SE = NaN;
- estimation->V = NaN;
- } else {
- estimation->E = weight/(double)N;
- estimation->V = sqr_weight/(double)N - estimation->E*estimation->E;
- estimation->SE = sqrt(estimation->V/(double)N);
- }
-
-exit:
- return res;
-error:
- goto exit;
-}
-
-static res_T
-solve_tile
- (struct sdis_scene* scn,
- struct ssp_rng* rng,
- const struct sdis_medium* mdm,
- const struct sdis_camera* cam,
- const double time,
- const double fp_to_meter,
- const double Tarad,
- const double Tref,
- const size_t origin[2], /* Tile origin in image plane */
- const size_t size[2], /* #pixels in X and Y */
- const size_t spp, /* #samples per pixel */
- const double pix_sz[2], /* Pixel size in the normalized image plane */
- struct sdis_mc* estimations)
-{
- size_t mcode; /* Morton code of the tile pixel */
- size_t npixels;
- res_T res = RES_OK;
- ASSERT(scn && rng && mdm && cam && spp && origin && estimations);
- ASSERT(size &&size[0] && size[1]);
- ASSERT(pix_sz && pix_sz[0] > 0 && pix_sz[1] > 0);
-
- /* Adjust the #pixels to process them wrt a morton order */
- npixels = round_up_pow2(MMAX(size[0], size[1]));
- npixels *= npixels;
-
- FOR_EACH(mcode, 0, npixels) {
- size_t ipix[2];
- struct sdis_mc* estimation;
-
- ipix[0] = morton2D_decode((uint32_t)(mcode>>0));
- if(ipix[0] >= size[0]) continue;
- ipix[1] = morton2D_decode((uint32_t)(mcode>>1));
- if(ipix[1] >= size[1]) continue;
-
- estimation = estimations + ipix[1]*size[0] + ipix[0];
- ipix[0] = ipix[0] + origin[0];
- ipix[1] = ipix[1] + origin[1];
-
- res = solve_pixel(scn, rng, mdm, cam, time, fp_to_meter, Tarad, Tref, size,
- spp, pix_sz, estimation);
- if(res != RES_OK) goto error;
- }
-
-exit:
- return res;
-error:
- goto exit;
-}
-
-/*******************************************************************************
- * Exported functions
- ******************************************************************************/
-res_T
-sdis_solve_probe
- (struct sdis_scene* scn,
- const size_t nrealisations,
- const double position[3],
- const double time,
- const double fp_to_meter,/* Scale factor from floating point unit to meter */
- const double Tarad, /* Ambient radiative temperature */
- const double Tref, /* Reference temperature */
- struct sdis_estimator** out_estimator)
-{
- const struct sdis_medium* medium = NULL;
- struct sdis_estimator* estimator = NULL;
- struct ssp_rng_proxy* rng_proxy = NULL;
- struct ssp_rng** rngs = NULL;
- double weight = 0;
- double sqr_weight = 0;
- size_t irealisation = 0;
- size_t N = 0; /* #realisations that do not fail */
- size_t i;
- ATOMIC res = RES_OK;
-
- if(!scn || !nrealisations || !position || time < 0 || fp_to_meter <= 0
- || Tref < 0 || !out_estimator) {
- res = RES_BAD_ARG;
- goto error;
- }
-
- /* Create the proxy RNG */
- res = ssp_rng_proxy_create(scn->dev->allocator, &ssp_rng_mt19937_64,
- scn->dev->nthreads, &rng_proxy);
- if(res != RES_OK) goto error;
-
- /* Create the per thread RNG */
- rngs = MEM_CALLOC
- (scn->dev->allocator, scn->dev->nthreads, sizeof(struct ssp_rng*));
- if(!rngs) {
- res = RES_MEM_ERR;
- goto error;
- }
- FOR_EACH(i, 0, scn->dev->nthreads) {
- res = ssp_rng_proxy_create_rng(rng_proxy, i, rngs+i);
- if(res != RES_OK) goto error;
- }
-
- /* Create the estimator */
- res = estimator_create(scn->dev, &estimator);
- if(res != RES_OK) goto error;
-
- /* Retrieve the medium in which the submitted position lies */
- res = scene_get_medium(scn, position, &medium);
- if(res != RES_OK) goto error;
-
- /* Here we go! Launch the Monte Carlo estimation */
- #pragma omp parallel for schedule(static) reduction(+:weight,sqr_weight,N)
- for(irealisation = 0; irealisation < nrealisations; ++irealisation) {
- res_T res_local;
- double w;
- const int ithread = omp_get_thread_num();
- struct ssp_rng* rng = rngs[ithread];
-
- if(ATOMIC_GET(&res) != RES_OK) continue; /* An error occured */
-
- if(scene_is_2d(scn)) {
- res_local = probe_realisation_2d
- (scn, rng, medium, position, time, fp_to_meter, Tarad, Tref, &w);
- } else {
- res_local = probe_realisation_3d
- (scn, rng, medium, position, time, fp_to_meter, Tarad, Tref, &w);
- }
- if(res_local != RES_OK) {
- if(res_local != RES_BAD_OP) {
- ATOMIC_SET(&res, res_local);
- continue;
- }
- } else {
- weight += w;
- sqr_weight += w*w;
- ++N;
- }
- }
-
- estimator->nrealisations = N;
- estimator->nfailures = nrealisations - N;
- estimator->temperature.E = weight / (double)N;
- estimator->temperature.V =
- sqr_weight / (double)N
- - estimator->temperature.E * estimator->temperature.E;
- estimator->temperature.SE = sqrt(estimator->temperature.V / (double)N);
-
-exit:
- if(rngs) {
- FOR_EACH(i, 0, scn->dev->nthreads) {
- if(rngs[i]) SSP(rng_ref_put(rngs[i]));
- }
- MEM_RM(scn->dev->allocator, rngs);
- }
- if(rng_proxy) SSP(rng_proxy_ref_put(rng_proxy));
- if(out_estimator) *out_estimator = estimator;
- return (res_T)res;
-error:
- if(estimator) {
- SDIS(estimator_ref_put(estimator));
- estimator = NULL;
- }
- goto exit;
-}
-
-res_T
-sdis_solve_camera
- (struct sdis_scene* scn,
- const struct sdis_camera* cam,
- const double time,
- const double fp_to_meter, /* Scale from floating point units to meters */
- const double Tarad, /* In Kelvin */
- const double Tref, /* In Kelvin */
- const size_t width, /* #pixels in X */
- const size_t height, /* #pixels in Y */
- const size_t spp, /* #samples per pixel */
- sdis_write_estimations_T writer,
- void* writer_data)
-{
- #define TILE_SIZE 32 /* definition in X & Y of a tile */
- STATIC_ASSERT(IS_POW2(TILE_SIZE), TILE_SIZE_must_be_a_power_of_2);
-
- const struct sdis_medium* medium = NULL;
- struct darray_mc* tiles = NULL;
- struct ssp_rng_proxy* rng_proxy = NULL;
- struct ssp_rng** rngs = NULL;
- size_t ntiles_x, ntiles_y, ntiles;
- double pix_sz[2]; /* Size of a pixel in the normalized image plane */
- int64_t mcode; /* Morton code of a tile */
- size_t i;
- ATOMIC res = RES_OK;
-
- if(!scn || !cam || time < 0 || fp_to_meter <= 0 || Tref < 0 || !width
- || !height || !spp || !writer) {
- res = RES_BAD_ARG;
- goto error;
- }
-
- if(scene_is_2d(scn)) {
- log_err(scn->dev, "%s: 2D scene are not supported.\n", FUNC_NAME);
- goto error;
- }
-
- /* Retrieve the medium in which the submitted position lies */
- res = scene_get_medium(scn, cam->position, &medium);
- if(res != RES_OK) goto error;
-
- if(medium != SDIS_MEDIUM_FLUID) {
- log_err(scn->dev, "%s: the camera position `%g %g %g' is not in a fluid.\n",
- FUNC_NAME, SPLIT3(cam->position));
- res = RES_BAD_ARG;
- goto error;
- }
-
- /* Create the proxy RNG */
- res = ssp_rng_proxy_create(scn->dev->allocator, &ssp_rng_mt19937_64,
- scn->dev->nthreads, &rng_proxy);
- if(res != RES_OK) goto error;
-
- /* Create the per thread RNG */
- rngs = MEM_CALLOC
- (scn->dev->allocator, scn->dev->nthreads, sizeof(struct ssp_rng*));
- if(!rngs) {
- res = RES_MEM_ERR;
- goto error;
- }
- FOR_EACH(i, 0, scn->dev->nthreads) {
- res = ssp_rng_proxy_create_rng(rng_proxy, i, rngs+i);
- if(res != RES_OK) goto error;
- }
-
- /* Allocate per thread buffer of estimations */
- tiles = darray_tile_data_get(&scn->dev->tiles);
- ASSERT(darray_tile_size_get(&scn->dev->tiles) == scn->dev->nthreads);
- FOR_EACH(i, 0, scn->dev->nthreads) {
- const size_t nestimations = TILE_SIZE * TILE_SIZE;
- res = darray_mc_resize(tiles+i, nestimations);
- if(res != RES_OK) goto error;
- }
-
- ntiles_x = (width + (TILE_SIZE-1)/*ceil*/)/TILE_SIZE;
- ntiles_y = (height + (TILE_SIZE-1)/*ceil*/)/TILE_SIZE;
- ntiles = round_up_pow2(MMAX(ntiles_x, ntiles_y));
- ntiles *= ntiles;
-
- pix_sz[0] = 1.0 / (double)width;
- pix_sz[1] = 1.0 / (double)height;
-
- omp_set_num_threads((int)scn->dev->nthreads);
- /*#pragma omp parallel for schedule(static)*/
- for(mcode = 0; mcode < (int64_t)ntiles; ++mcode) {
- size_t tile_org[2] = {0, 0};
- size_t tile_sz[2] = {0, 0};
- const int ithread = omp_get_thread_num();
- struct sdis_mc* estimations = NULL;
- struct ssp_rng* rng = rngs[ithread];
- res_T res_local = RES_OK;
-
- if(ATOMIC_GET(&res) != RES_OK) continue;
-
- tile_org[0] = morton2D_decode((uint32_t)(mcode>>0));
- if(tile_org[0] >= ntiles_x) continue; /* Discard tile */
- tile_org[1] = morton2D_decode((uint32_t)(mcode>>1));
- if(tile_org[1] >= ntiles_y) continue; /* Disaard tile */
-
- /* Setup the tile coordinates in the image plane */
- tile_org[0] *= TILE_SIZE;
- tile_org[1] *= TILE_SIZE;
- tile_sz[0] = MMIN(TILE_SIZE, width - tile_org[0]);
- tile_sz[0] = MMIN(TILE_SIZE, width - tile_org[1]);
-
- /* Fetch the estimation buffer */
- estimations = darray_mc_data_get(tiles+ithread);
-
- /* Draw the tile */
- res_local = solve_tile(scn, rng, medium, cam, time, fp_to_meter, Tarad,
- Tref, tile_org, tile_sz, spp, pix_sz, estimations);
- if(res_local != RES_OK) {
- ATOMIC_SET(&res, res_local);
- continue;
- }
-
- /* Write the estimations */
- res_local = writer(writer_data, tile_org, tile_sz, estimations);
- if(res_local != RES_OK) {
- ATOMIC_SET(&res, res_local);
- continue;
- }
- }
-
-exit:
- if(rngs) {
- FOR_EACH(i, 0, scn->dev->nthreads) {
- if(rngs[i]) SSP(rng_ref_put(rngs[i]));
- }
- MEM_RM(scn->dev->allocator, rngs);
- }
- if(rng_proxy) SSP(rng_proxy_ref_put(rng_proxy));
- return (res_T)res;
-error:
- goto exit;
-}
-
diff --git a/src/sdis_solve_probe_Xd.h b/src/sdis_solve_probe_Xd.h
@@ -1,826 +0,0 @@
-/* Copyright (C) |Meso|Star> 2016-2018 (contact@meso-star.com)
- *
- * This program is free software: you can redistribute it and/or modify
- * it under the terms of the GNU General Public License as published by
- * the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * This program is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU General Public License
- * along with this program. If not, see <http://www.gnu.org/licenses/>. */
-
-#ifndef SDIS_SOLVE_PROBE_DIMENSION
-#ifndef SDIS_SOLVE_PROBE_XD_H
-#define SDIS_SOLVE_PROBE_XD_H
-
-#include "sdis_device_c.h"
-#include "sdis_interface_c.h"
-#include "sdis_medium_c.h"
-#include "sdis_scene_c.h"
-
-#include <rsys/stretchy_array.h>
-
-#include <star/ssp.h>
-
-/* Emperical scale factor to apply to the upper bound of the ray range in order
- * to handle numerical imprecisions */
-#define RAY_RANGE_MAX_SCALE 1.0001f
-
-#define BOLTZMANN_CONSTANT 5.6696e-8 /* W/m^2/K^4 */
-
-struct rwalk_context {
- double Tarad; /* Ambient radiative temperature */
- double Tref3; /* Reference temperature ^ 3 */
-};
-
-/* Reflect the vector V wrt the normal N. By convention V points outward the
- * surface. */
-static INLINE float*
-reflect(float res[3], const float V[3], const float N[3])
-{
- float tmp[3];
- float cos_V_N;
- ASSERT(res && V && N);
- ASSERT(f3_is_normalized(V) && f3_is_normalized(N));
- cos_V_N = f3_dot(V, N);
- f3_mulf(tmp, N, 2*cos_V_N);
- f3_sub(res, tmp, V);
- return res;
-}
-
-#endif /* SDIS_SOLVE_PROBE_XD_H */
-#else
-
-#if (SDIS_SOLVE_PROBE_DIMENSION == 2)
- #include <rsys/double2.h>
- #include <rsys/float2.h>
- #include <star/s2d.h>
-#elif (SDIS_SOLVE_PROBE_DIMENSION == 3)
- #include <rsys/double2.h>
- #include <rsys/double3.h>
- #include <rsys/float3.h>
- #include <star/s3d.h>
-#else
- #error "Invalid SDIS_SOLVE_PROBE_DIMENSION value."
-#endif
-
-/* Syntactic sugar */
-#define DIM SDIS_SOLVE_PROBE_DIMENSION
-
-/* Star-XD macros generic to SDIS_SOLVE_PROBE_DIMENSION */
-#define sXd(Name) CONCAT(CONCAT(CONCAT(s, DIM), d_), Name)
-#define SXD_HIT_NONE CONCAT(CONCAT(S,DIM), D_HIT_NONE)
-#define SXD_HIT_NULL CONCAT(CONCAT(S,DIM), D_HIT_NULL)
-#define SXD_HIT_NULL__ CONCAT(CONCAT(S, DIM), D_HIT_NULL__)
-#define SXD CONCAT(CONCAT(S, DIM), D)
-
-/* Vector macros generic to SDIS_SOLVE_PROBE_DIMENSION */
-#define dX(Func) CONCAT(CONCAT(CONCAT(d, DIM), _), Func)
-#define fX(Func) CONCAT(CONCAT(CONCAT(f, DIM), _), Func)
-#define fX_set_dX CONCAT(CONCAT(CONCAT(f, DIM), _set_d), DIM)
-#define dX_set_fX CONCAT(CONCAT(CONCAT(d, DIM), _set_f), DIM)
-
-/* Macro making generic its subimitted name to SDIS_SOLVE_PROBE_DIMENSION */
-#define XD(Name) CONCAT(CONCAT(CONCAT(Name, _), DIM), d)
-
-/* Current state of the random walk */
-struct XD(rwalk) {
- struct sdis_rwalk_vertex vtx; /* Position and time of the Random walk */
- const struct sdis_medium* mdm; /* Medium in which the random walk lies */
- struct sXd(hit) hit; /* Hit of the random walk */
-};
-static const struct XD(rwalk) XD(RWALK_NULL) = {
- SDIS_RWALK_VERTEX_NULL__, NULL, SXD_HIT_NULL__
-};
-
-struct XD(temperature) {
- res_T (*func)/* Next function to invoke in order to compute the temperature */
- (const struct sdis_scene* scn,
- const double fp_to_meter,
- const struct rwalk_context* ctx,
- struct XD(rwalk)* rwalk,
- struct ssp_rng* rng,
- struct XD(temperature)* temp);
- double value; /* Current value of the temperature */
- int done;
-};
-static const struct XD(temperature) XD(TEMPERATURE_NULL) = { NULL, 0, 0 };
-
-static res_T
-XD(boundary_temperature)
- (const struct sdis_scene* scn,
- const double fp_to_meter,
- const struct rwalk_context* ctx,
- struct XD(rwalk)* rwalk,
- struct ssp_rng* rng,
- struct XD(temperature)* T);
-
-static res_T
-XD(solid_temperature)
- (const struct sdis_scene* scn,
- const double fp_to_meter,
- const struct rwalk_context* ctx,
- struct XD(rwalk)* rwalk,
- struct ssp_rng* rng,
- struct XD(temperature)* T);
-
-static res_T
-XD(fluid_temperature)
- (const struct sdis_scene* scn,
- const double fp_to_meter,
- const struct rwalk_context* ctx,
- struct XD(rwalk)* rwalk,
- struct ssp_rng* rng,
- struct XD(temperature)* T);
-
-static res_T
-XD(radiative_temperature)
- (const struct sdis_scene* scn,
- const double fp_to_meter,
- const struct rwalk_context* ctx,
- struct XD(rwalk)* rwalk,
- struct ssp_rng* rng,
- struct XD(temperature)* T);
-
-/*******************************************************************************
- * Helper functions
- ******************************************************************************/
-
-static FINLINE void
-XD(move_pos)(double pos[DIM], const float dir[DIM], const float delta)
-{
- ASSERT(pos && dir);
- pos[0] += dir[0] * delta;
- pos[1] += dir[1] * delta;
-#if(SDIS_SOLVE_PROBE_DIMENSION == 3)
- pos[2] += dir[2] * delta;
-#endif
-}
-
-/* Check that the interface fragment is consistent with the current state of
- * the random walk */
-static INLINE int
-XD(check_rwalk_fragment_consistency)
- (const struct XD(rwalk)* rwalk,
- const struct sdis_interface_fragment* frag)
-{
- double N[DIM];
- double uv[2] = {0, 0};
- ASSERT(rwalk && frag);
- dX(normalize)(N, dX_set_fX(N, rwalk->hit.normal));
- if( SXD_HIT_NONE(&rwalk->hit)
- || !dX(eq_eps)(rwalk->vtx.P, frag->P, 1.e-6)
- || !dX(eq_eps)(N, frag->Ng, 1.e-6)
- || !( (IS_INF(rwalk->vtx.time) && IS_INF(frag->time))
- || eq_eps(rwalk->vtx.time, frag->time, 1.e-6))) {
- return 0;
- }
-#if (SDIS_SOLVE_PROBE_DIMENSION == 2)
- uv[0] = rwalk->hit.u;
-#else
- d2_set_f2(uv, rwalk->hit.uv);
-#endif
- return d2_eq_eps(uv, frag->uv, 1.e-6);
-}
-
-res_T
-XD(radiative_temperature)
- (const struct sdis_scene* scn,
- const double fp_to_meter,
- const struct rwalk_context* ctx,
- struct XD(rwalk)* rwalk,
- struct ssp_rng* rng,
- struct XD(temperature)* T)
-{
- const struct sdis_interface* interf;
-
- /* The radiative random walk is always perform in 3D. In 2D, the geometry are
- * assumed to be extruded to the infinty along the Z dimension. */
- float N[3] = {0, 0, 0};
- float dir[3] = {0, 0, 0};
-
- res_T res = RES_OK;
-
- ASSERT(scn && fp_to_meter > 0 && ctx && rwalk && rng && T);
- ASSERT(!SXD_HIT_NONE(&rwalk->hit));
- (void)fp_to_meter;
-
- /* Fetch the current interface */
- interf = scene_get_interface(scn, rwalk->hit.prim.prim_id);
-
- /* Normalize the normal of the interface and ensure that it points toward the
- * current medium */
- fX(normalize(N, rwalk->hit.normal));
- if(interf->medium_back == rwalk->mdm) {
- fX(minus(N, N));
- }
-
- /* Cosine weighted sampling of a direction around the surface normal */
- ssp_ran_hemisphere_cos_float(rng, N, dir, NULL);
-
- /* Launch the radiative random walk */
- for(;;) {
- struct sdis_interface_fragment frag = SDIS_INTERFACE_FRAGMENT_NULL;
- const struct sdis_medium* chk_mdm = NULL;
- double alpha;
- double epsilon;
- double r;
- float pos[DIM];
- const float range[2] = { 0, FLT_MAX };
-
- fX_set_dX(pos, rwalk->vtx.P);
-
- /* Trace the radiative ray */
-#if (SDIS_SOLVE_PROBE_DIMENSION == 2)
- SXD(scene_view_trace_ray_3d
- (scn->sXd(view), pos, dir, range, &rwalk->hit, &rwalk->hit));
-#else
- SXD(scene_view_trace_ray
- (scn->sXd(view), pos, dir, range, &rwalk->hit, &rwalk->hit));
-#endif
- if(SXD_HIT_NONE(&rwalk->hit)) { /* Fetch the ambient radiative temperature */
- if(ctx->Tarad >= 0) {
- T->value += ctx->Tarad;
- T->done = 1;
- break;
- } else {
- log_err(scn->dev,
-"%s: the random walk reaches an invalid ambient radiative temperature of `%gK'\n"
-"at position `%g %g %g'. This may be due to numerical inaccuracies or to\n"
-"inconsistency in the simulated system (eg: unclosed geometry). For systems\n"
-"where the random walks can reach such temperature, one has to setup a valid\n"
-"ambient radiative temperature, i.e. it must be greater or equal to 0.\n",
- FUNC_NAME,
- ctx->Tarad,
- SPLIT3(rwalk->vtx.P));
- res = RES_BAD_ARG;
- goto error;
- }
- }
-
- /* Move the random walk to the hit position */
- XD(move_pos)(rwalk->vtx.P, dir, rwalk->hit.distance);
-
- /* Fetch the new interface and setup the hit fragment */
- interf = scene_get_interface(scn, rwalk->hit.prim.prim_id);
- XD(setup_interface_fragment)(&frag, &rwalk->vtx, &rwalk->hit);
-
- /* Fetch the interface emissivity */
- epsilon = interface_get_emissivity(interf, &frag);
- if(epsilon > 1 && epsilon >= 0) {
- log_err(scn->dev,
- "%s: invalid overall emissivity `%g' at position `%g %g %g'.\n",
- FUNC_NAME, epsilon, SPLIT3(rwalk->vtx.P));
- res = RES_BAD_ARG;
- goto error;
- }
-
- /* Switch in boundary temperature ? */
- r = ssp_rng_canonical(rng);
- if(r < epsilon) {
- T->func = XD(boundary_temperature);
- break;
- }
-
- /* Normalize the normal of the interface and ensure that it points toward the
- * current medium */
- fX(normalize)(N, rwalk->hit.normal);
- if(f3_dot(N, dir) > 0) {
- chk_mdm = interf->medium_back;
- fX(minus)(N, N);
- } else {
- chk_mdm = interf->medium_front;
- }
-
- if(chk_mdm != rwalk->mdm) {
- log_err(scn->dev, "%s: inconsistent medium definition at `%g %g %g'.\n",
- FUNC_NAME, SPLIT3(rwalk->vtx.P));
- res = RES_BAD_ARG;
- goto error;
- }
- alpha = interface_get_specular_fraction(interf, &frag);
- r = ssp_rng_canonical(rng);
- if(r < alpha) { /* Sample specular part */
- reflect(dir, f3_minus(dir, dir), N);
- } else { /* Sample diffuse part */
- ssp_ran_hemisphere_cos_float(rng, N, dir, NULL);
- }
- }
-
-exit:
- return res;
-error:
- goto exit;
-}
-
-res_T
-XD(fluid_temperature)
- (const struct sdis_scene* scn,
- const double fp_to_meter,
- const struct rwalk_context* ctx,
- struct XD(rwalk)* rwalk,
- struct ssp_rng* rng,
- struct XD(temperature)* T)
-{
- double tmp;
- (void)rng, (void)fp_to_meter, (void)ctx;
- ASSERT(scn && fp_to_meter > 0 && ctx && rwalk && rng && T);
- ASSERT(rwalk->mdm->type == SDIS_MEDIUM_FLUID);
-
- tmp = fluid_get_temperature(rwalk->mdm, &rwalk->vtx);
- if(tmp < 0) {
- log_err(scn->dev, "%s: unknown fluid temperature is not supported.\n",
- FUNC_NAME);
- return RES_BAD_OP;
- }
- T->value += tmp;
- T->done = 1;
- return RES_OK;
-}
-
-static void
-XD(solid_solid_boundary_temperature)
- (const struct sdis_scene* scn,
- const double fp_to_meter,
- const struct rwalk_context* ctx,
- const struct sdis_interface_fragment* frag,
- struct XD(rwalk)* rwalk,
- struct ssp_rng* rng,
- struct XD(temperature)* T)
-{
- const struct sdis_interface* interf = NULL;
- const struct sdis_medium* solid_front = NULL;
- const struct sdis_medium* solid_back = NULL;
- double lambda_front, lambda_back;
- double delta_front_boundary, delta_back_boundary;
- double delta_front_boundary_meter, delta_back_boundary_meter;
- double delta_boundary;
- double proba;
- double tmp;
- double r;
- float pos[DIM], dir[DIM], range[2];
- ASSERT(scn && fp_to_meter > 0 && ctx && frag && rwalk && rng && T);
- ASSERT(XD(check_rwalk_fragment_consistency)(rwalk, frag));
- (void)frag, (void)ctx;
-
- /* Retrieve the current boundary media */
- interf = scene_get_interface(scn, rwalk->hit.prim.prim_id);
- solid_front = interface_get_medium(interf, SDIS_FRONT);
- solid_back = interface_get_medium(interf, SDIS_BACK);
- ASSERT(solid_front->type == SDIS_MEDIUM_SOLID);
- ASSERT(solid_back->type == SDIS_MEDIUM_SOLID);
-
- /* Fetch the properties of the media */
- lambda_front = solid_get_thermal_conductivity(solid_front, &rwalk->vtx);
- lambda_back = solid_get_thermal_conductivity(solid_back, &rwalk->vtx);
- delta_front_boundary = solid_get_delta_boundary(solid_front, &rwalk->vtx);
- delta_back_boundary = solid_get_delta_boundary(solid_back, &rwalk->vtx);
-
- /* Convert the delta boundary from floating point units to meters */
- delta_front_boundary_meter = fp_to_meter * delta_front_boundary;
- delta_back_boundary_meter = fp_to_meter * delta_back_boundary;
-
- /* Compute the proba to switch in solid0 or in solid1 */
- tmp = lambda_front / delta_front_boundary_meter;
- proba = tmp / (tmp + lambda_back / delta_back_boundary_meter);
-
- r = ssp_rng_canonical(rng);
- fX(normalize)(dir, rwalk->hit.normal);
- if(r < proba) { /* Reinject in front */
- delta_boundary = delta_front_boundary;
- rwalk->mdm = solid_front;
- } else { /* Reinject in back */
- delta_boundary = delta_back_boundary;
- fX(minus)(dir, dir);
- rwalk->mdm = solid_back;
- }
-
- /* "Reinject" the path into the solid along the surface normal. */
- fX_set_dX(pos, rwalk->vtx.P);
- range[0] = 0, range[1] = (float)delta_boundary*RAY_RANGE_MAX_SCALE;
- SXD(scene_view_trace_ray
- (scn->sXd(view), pos, dir, range, &rwalk->hit, &rwalk->hit));
- if(!SXD_HIT_NONE(&rwalk->hit)) delta_boundary = rwalk->hit.distance * 0.5;
- XD(move_pos)(rwalk->vtx.P, dir, (float)delta_boundary);
-
- /* Switch in solid random walk */
- T->func = XD(solid_temperature);
-}
-
-static void
-XD(solid_fluid_boundary_temperature)
- (const struct sdis_scene* scn,
- const double fp_to_meter,
- const struct rwalk_context* ctx,
- const struct sdis_interface_fragment* frag,
- struct XD(rwalk)* rwalk,
- struct ssp_rng* rng,
- struct XD(temperature)* T)
-{
- const struct sdis_interface* interf = NULL;
- const struct sdis_medium* mdm_front = NULL;
- const struct sdis_medium* mdm_back = NULL;
- const struct sdis_medium* solid = NULL;
- const struct sdis_medium* fluid = NULL;
- double hc;
- double hr;
- double epsilon; /* Interface emissivity */
- double lambda;
- double fluid_proba;
- double radia_proba;
- double delta_boundary;
- double r;
- double tmp;
- float dir[DIM], pos[DIM], range[2];
-
- ASSERT(scn && fp_to_meter > 0 && rwalk && rng && T && ctx);
- ASSERT(XD(check_rwalk_fragment_consistency)(rwalk, frag));
-
- /* Retrieve the solid and the fluid split by the boundary */
- interf = scene_get_interface(scn, rwalk->hit.prim.prim_id);
- mdm_front = interface_get_medium(interf, SDIS_FRONT);
- mdm_back = interface_get_medium(interf, SDIS_BACK);
- ASSERT(mdm_front->type != mdm_back->type);
- if(mdm_front->type == SDIS_MEDIUM_SOLID) {
- solid = mdm_front;
- fluid = mdm_back;
- } else {
- solid = mdm_back;
- fluid = mdm_front;
- }
-
- /* Fetch the solid properties */
- lambda = solid_get_thermal_conductivity(solid, &rwalk->vtx);
- delta_boundary = solid_get_delta_boundary(solid, &rwalk->vtx);
-
- /* Fetch the boundary properties */
- epsilon = interface_get_emissivity(interf, frag);
- hc = interface_get_convection_coef(interf, frag);
-
- /* Compute the radiative coefficient */
- hr = 4.0 * BOLTZMANN_CONSTANT * ctx->Tref3 * epsilon;
-
- /* Compute the probas to switch in solid or fluid random walk */
- tmp = lambda / (delta_boundary*fp_to_meter);
- fluid_proba = hc / (tmp + hr + hc);
- radia_proba = hr / (tmp + hr + hc);
- /*solid_proba = tmp / (tmp + hr + hc);*/
-
- r = ssp_rng_canonical(rng);
- if(r < radia_proba) { /* Switch in radiative random walk */
- rwalk->mdm = fluid;
- T->func = XD(radiative_temperature);
- } else if(r < fluid_proba + radia_proba) { /* Switch to fluid random walk */
- rwalk->mdm = fluid;
- T->func = XD(fluid_temperature);
- } else { /* Solid random walk */
- rwalk->mdm = solid;
- fX(normalize)(dir, rwalk->hit.normal);
- if(solid == mdm_back) fX(minus)(dir, dir);
-
- /* "Reinject" the random walk into the solid */
- fX_set_dX(pos, rwalk->vtx.P);
- range[0] = 0, range[1] = (float)delta_boundary*RAY_RANGE_MAX_SCALE;
- SXD(scene_view_trace_ray
- (scn->sXd(view), pos, dir, range, &rwalk->hit, &rwalk->hit));
- if(!SXD_HIT_NONE(&rwalk->hit)) delta_boundary = rwalk->hit.distance * 0.5;
- XD(move_pos)(rwalk->vtx.P, dir, (float)delta_boundary);
-
- /* Switch in solid random walk */
- T->func = XD(solid_temperature);
- }
-}
-
-res_T
-XD(boundary_temperature)
- (const struct sdis_scene* scn,
- const double fp_to_meter,
- const struct rwalk_context* ctx,
- struct XD(rwalk)* rwalk,
- struct ssp_rng* rng,
- struct XD(temperature)* T)
-{
- struct sdis_interface_fragment frag = SDIS_INTERFACE_FRAGMENT_NULL;
- const struct sdis_interface* interf = NULL;
- const struct sdis_medium* mdm_front = NULL;
- const struct sdis_medium* mdm_back = NULL;
- double tmp;
- ASSERT(scn && fp_to_meter > 0 && ctx && rwalk && rng && T);
- ASSERT(!SXD_HIT_NONE(&rwalk->hit));
-
- XD(setup_interface_fragment)(&frag, &rwalk->vtx, &rwalk->hit);
-
- /* Retrieve the current interface */
- interf = scene_get_interface(scn, rwalk->hit.prim.prim_id);
-
- /* Check if the boundary condition is known */
- tmp = interface_get_temperature(interf, &frag);
- if(tmp >= 0) {
- T->value += tmp;
- T->done = 1;
- return RES_OK;
- }
-
- mdm_front = interface_get_medium(interf, SDIS_FRONT);
- mdm_back = interface_get_medium(interf, SDIS_BACK);
-
- if(mdm_front->type == mdm_back->type) {
- XD(solid_solid_boundary_temperature)
- (scn, fp_to_meter, ctx, &frag, rwalk, rng, T);
- } else {
- XD(solid_fluid_boundary_temperature)
- (scn, fp_to_meter, ctx, &frag, rwalk, rng, T);
- }
- return RES_OK;
-}
-
-res_T
-XD(solid_temperature)
- (const struct sdis_scene* scn,
- const double fp_to_meter,
- const struct rwalk_context* ctx,
- struct XD(rwalk)* rwalk,
- struct ssp_rng* rng,
- struct XD(temperature)* T)
-{
- double position_start[DIM];
- const struct sdis_medium* mdm;
- ASSERT(scn && fp_to_meter > 0 && rwalk && rng && T);
- ASSERT(rwalk->mdm->type == SDIS_MEDIUM_SOLID);
- (void)ctx;
-
- /* Check the random walk consistency */
- CHK(scene_get_medium(scn, rwalk->vtx.P, &mdm) == RES_OK);
- if(mdm != rwalk->mdm) {
- log_err(scn->dev, "%s: invalid solid random walk. "
- "Unexpected medium at {%g, %g, %g}.\n",
- FUNC_NAME, SPLIT3(rwalk->vtx.P));
- return RES_BAD_ARG;
- }
- /* Save the submitted position */
- dX(set)(position_start, rwalk->vtx.P);
-
- do { /* Solid random walk */
- struct sXd(hit) hit0, hit1;
- double lambda; /* Thermal conductivity */
- double rho; /* Volumic mass */
- double cp; /* Calorific capacity */
- double tau, mu;
- double tmp;
- float delta, delta_solid; /* Random walk numerical parameter */
- float range[2];
- float dir0[DIM], dir1[DIM];
- float org[DIM];
-
- /* Check the limit condition */
- tmp = solid_get_temperature(mdm, &rwalk->vtx);
- if(tmp >= 0) {
- T->value += tmp;
- T->done = 1;
- return RES_OK;
- }
-
- /* Fetch solid properties */
- delta_solid = (float)solid_get_delta(mdm, &rwalk->vtx);
- lambda = solid_get_thermal_conductivity(mdm, &rwalk->vtx);
- rho = solid_get_volumic_mass(mdm, &rwalk->vtx);
- cp = solid_get_calorific_capacity(mdm, &rwalk->vtx);
-
-#if (SDIS_SOLVE_PROBE_DIMENSION == 2)
- /* Sample a direction around 2PI */
- ssp_ran_circle_uniform_float(rng, dir0, NULL);
-#else
- /* Sample a direction around 4PI */
- ssp_ran_sphere_uniform_float(rng, dir0, NULL);
-#endif
-
- /* Trace a ray along the sampled direction and its opposite to check if a
- * surface is hit in [0, delta_solid]. */
- fX_set_dX(org, rwalk->vtx.P);
- fX(minus)(dir1, dir0);
- hit0 = hit1 = SXD_HIT_NULL;
- range[0] = 0.f, range[1] = delta_solid*RAY_RANGE_MAX_SCALE;
- SXD(scene_view_trace_ray(scn->sXd(view), org, dir0, range, NULL, &hit0));
- SXD(scene_view_trace_ray(scn->sXd(view), org, dir1, range, NULL, &hit1));
-
- if(SXD_HIT_NONE(&hit0) && SXD_HIT_NONE(&hit1)) {
- /* Hit nothing: move along dir0 of the original delta */
- delta = delta_solid;
- } else {
- /* Hit something: move along dir0 of the minimum hit distance */
- delta = MMIN(hit0.distance, hit1.distance);
- }
-
- /* Sample the time */
- mu = (2*DIM*lambda) / (rho*cp*delta*fp_to_meter*delta*fp_to_meter);
- tau = ssp_ran_exp(rng, mu);
- rwalk->vtx.time -= tau;
-
- /* Check the initial condition */
- tmp = solid_get_temperature(mdm, &rwalk->vtx);
- if(tmp >= 0) {
- T->value += tmp;
- T->done = 1;
- return RES_OK;
- }
-
- /* Define if the random walk hits something along dir0 */
- rwalk->hit = hit0.distance > delta ? SXD_HIT_NULL : hit0;
-
- /* Update the random walk position */
- XD(move_pos)(rwalk->vtx.P, dir0, delta);
-
- /* Fetch the current medium */
- if(SXD_HIT_NONE(&rwalk->hit)) {
- CHK(scene_get_medium(scn, rwalk->vtx.P, &mdm) == RES_OK);
- } else {
- const struct sdis_interface* interf;
- interf = scene_get_interface(scn, rwalk->hit.prim.prim_id);
- mdm = interface_get_medium
- (interf,
- fX(dot(rwalk->hit.normal, dir0)) < 0 ? SDIS_FRONT : SDIS_BACK);
- }
-
- /* Check random walk consistency */
- if(mdm != rwalk->mdm) {
- log_err(scn->dev,
- "%s: inconsistent medium during the solid random walk.\n", FUNC_NAME);
- if(DIM == 2) {
- log_err(scn->dev,
- " start position: %g %g; current position: %g %g\n",
- SPLIT2(position_start), SPLIT2(rwalk->vtx.P));
- } else {
- log_err(scn->dev,
- " start position: %g %g %g; current position: %g %g %g\n",
- SPLIT3(position_start), SPLIT3(rwalk->vtx.P));
- }
- return RES_BAD_OP;
- }
-
- /* Keep going while the solid random walk does not hit an interface */
- } while(SXD_HIT_NONE(&rwalk->hit));
-
- T->func = XD(boundary_temperature);
- return RES_OK;
-}
-
-static res_T
-XD(compute_temperature)
- (struct sdis_scene* scn,
- const double fp_to_meter,
- const struct rwalk_context* ctx,
- struct XD(rwalk)* rwalk,
- struct ssp_rng* rng,
- struct XD(temperature)* T)
-{
-#ifndef NDEBUG
- struct XD(temperature)* stack = NULL;
- size_t istack = 0;
-#endif
- res_T res = RES_OK;
- ASSERT(scn && fp_to_meter > 0 && ctx && rwalk && rng && T);
-
- do {
- res = T->func(scn, fp_to_meter, ctx, rwalk, rng, T);
- if(res != RES_OK) goto error;
-
-#ifndef NDEBUG
- sa_push(stack, *T);
- ++istack;
-#endif
- } while(!T->done);
-
-exit:
-#ifndef NDEBUG
- sa_release(stack);
-#endif
- return res;
-error:
- goto exit;
-}
-
-static res_T
-XD(probe_realisation)
- (struct sdis_scene* scn,
- struct ssp_rng* rng,
- const struct sdis_medium* medium,
- const double position[],
- const double time,
- const double fp_to_meter,/* Scale factor from floating point unit to meter */
- const double ambient_radiative_temperature,
- const double reference_temperature,
- double* weight)
-{
- struct rwalk_context ctx;
- struct XD(rwalk) rwalk = XD(RWALK_NULL);
- struct XD(temperature) T = XD(TEMPERATURE_NULL);
- res_T res = RES_OK;
- ASSERT(medium && position && fp_to_meter > 0 && weight && time >= 0);
-
- switch(medium->type) {
- case SDIS_MEDIUM_FLUID: T.func = XD(fluid_temperature); break;
- case SDIS_MEDIUM_SOLID: T.func = XD(solid_temperature); break;
- default: FATAL("Unreachable code\n"); break;
- }
-
- dX(set)(rwalk.vtx.P, position);
- rwalk.vtx.time = time;
- rwalk.hit = SXD_HIT_NULL;
- rwalk.mdm = medium;
-
- ctx.Tarad = ambient_radiative_temperature;
- ctx.Tref3 =
- reference_temperature
- * reference_temperature
- * reference_temperature;
-
- res = XD(compute_temperature)(scn, fp_to_meter, &ctx, &rwalk, rng, &T);
- if(res != RES_OK) return res;
-
- *weight = T.value;
- return RES_OK;
-}
-
-#if SDIS_SOLVE_PROBE_DIMENSION == 3
-static res_T
-XD(ray_realisation)
- (struct sdis_scene* scn,
- struct ssp_rng* rng,
- const struct sdis_medium* medium,
- const double position[],
- const double direction[],
- const double time,
- const double fp_to_meter,
- const double Tarad,
- const double Tref,
- double* weight)
-{
- struct sXd(hit) hit = SXD_HIT_NULL;
- const float range[2] = {0, FLT_MAX};
- float org[3] = {0, 0, 0};
- float dir[3] = {0, 0, 0};
- res_T res = RES_OK;
- ASSERT(scn && position && direction && time>=0 && fp_to_meter>0 && weight);
- ASSERT(medium && medium->type == SDIS_MEDIUM_FLUID);
-
- fX_set_dX(org, position);
- fX_set_dX(dir, direction);
- SXD(scene_view_trace_ray(scn->sXd(view), org, dir, range, &hit, &hit));
-
- if(SXD_HIT_NONE(&hit)) {
- if(Tarad >= 0) {
- *weight = Tarad;
- } else {
- log_err(scn->dev,
-"%s: the ray starting from `%g %g %g' and traced along the `%g %g %g' direction\n"
-"reaches an invalid ambient temperature of `%gK'. One has to provide a valid\n"
-"ambient radiative temperature, i.e. it must be greater or equal to 0.\n",
- FUNC_NAME, SPLIT3(org), SPLIT3(dir), Tarad);
- res = RES_BAD_ARG;
- goto error;
- }
- } else {
- struct rwalk_context ctx;
- struct XD(rwalk) rwalk = XD(RWALK_NULL);
- struct XD(temperature) T = XD(TEMPERATURE_NULL);
-
- dX(set)(rwalk.vtx.P, position);
- XD(move_pos)(rwalk.vtx.P, dir, hit.distance);
- rwalk.vtx.time = time;
- rwalk.hit = hit;
- rwalk.mdm = medium;
-
- ctx.Tarad = Tarad;
- ctx.Tref3 = Tref*Tref*Tref;
-
- T.func = XD(boundary_temperature);
- res = XD(compute_temperature)(scn, fp_to_meter, &ctx, &rwalk, rng, &T);
- if(res != RES_OK) return res;
-
- *weight = T.value;
- }
-exit:
- return res;
-error:
- goto exit;
-}
-#endif /* SDIS_SOLVE_PROBE_DIMENSION == 3 */
-
-#undef SDIS_SOLVE_PROBE_DIMENSION
-#undef DIM
-#undef sXd
-#undef SXD_HIT_NONE
-#undef SXD_HIT_NULL
-#undef SXD_HIT_NULL__
-#undef SXD
-#undef dX
-#undef fX
-#undef fX_set_dX
-#undef XD
-
-#endif /* !SDIS_SOLVE_PROBE_DIMENSION */
-
