commit 11b6f2b1f333bda4ba7c4e0a53bd5e439b445c01
parent 70b07c7af2c43a8b6f051b9c6028bd29bcd857e7
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Fri, 6 Apr 2018 14:52:58 +0200
Add the sdis_solve_probe_boundary function
Diffstat:
6 files changed, 243 insertions(+), 0 deletions(-)
diff --git a/src/sdis.h b/src/sdis.h
@@ -441,6 +441,18 @@ sdis_solve_probe
struct sdis_estimator** estimator);
SDIS_API res_T
+sdis_solve_probe_boundary
+ (struct sdis_scene* scn,
+ const size_t nrealisations, /* #realisations */
+ const size_t iprim, /* Identifier of the primitive on which the probe lies */
+ const double uv[2], /* Parametric coordinates of the probe onto the primitve */
+ const double time, /* Observation time */
+ const double fp_to_meter, /* Scale from floating point units to meters */
+ const double ambient_radiative_temperature, /* In Kelvin */
+ const double reference_temperature, /* In Kelvin */
+ struct sdis_estimator** estimator);
+
+SDIS_API res_T
sdis_solve_camera
(struct sdis_scene* scn,
const struct sdis_camera* cam, /* Point of view */
diff --git a/src/sdis_medium_c.h b/src/sdis_medium_c.h
@@ -102,6 +102,16 @@ solid_get_delta_boundary
}
static INLINE double
+solid_get_volumic_power
+ (const struct sdis_medium* mdm, const struct sdis_rwalk_vertex* vtx)
+{
+ ASSERT(mdm && mdm->type == SDIS_MEDIUM_SOLID);
+ return mdm->shader.solid.volumic_power
+ ? mdm->shader.solid.volumic_power(vtx, mdm->data)
+ : 0;
+}
+
+static INLINE double
solid_get_temperature
(const struct sdis_medium* mdm, const struct sdis_rwalk_vertex* vtx)
{
diff --git a/src/sdis_scene_c.h b/src/sdis_scene_c.h
@@ -46,6 +46,13 @@ struct sdis_scene {
struct sdis_device* dev;
};
+static FINLINE size_t
+scene_get_primitives_count(const struct sdis_scene* scn)
+{
+ ASSERT(scn);
+ return darray_interf_size_get(&scn->prim_interfaces);
+}
+
extern LOCAL_SYM const struct sdis_interface*
scene_get_interface
(const struct sdis_scene* scene,
diff --git a/src/sdis_solve.c b/src/sdis_solve.c
@@ -267,6 +267,146 @@ error:
}
res_T
+sdis_solve_probe_boundary
+ (struct sdis_scene* scn,
+ const size_t nrealisations, /* #realisations */
+ const size_t iprim, /* Identifier of the primitive on which the probe lies */
+ const double uv[2], /* Parametric coordinates of the probe onto the primitve */
+ 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 */
+ struct sdis_estimator** out_estimator)
+{
+ 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;
+ res_T res = RES_OK;
+
+ if(!scn || !nrealisations || !uv || time < 0 || fp_to_meter <= 0
+ || Tref < 0 || !out_estimator) {
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ /* Check the primitive identifier */
+ if(iprim > scene_get_primitives_count(scn)) {
+ log_err(scn->dev,
+"%s: invalid primitive identifier `%lu'. It must be less than %lu.\n",
+ FUNC_NAME,
+ (unsigned long)iprim,
+ (unsigned long)scene_get_primitives_count(scn));
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ /* Check parametric coordinates */
+ if(scene_is_2d(scn)) {
+ const double v = CLAMP(1.0 - uv[0], 0, 1);
+ if(uv[0] < 0 || uv[0] > 1 || !eq_eps(uv[0] + v, 1, 1.e-6)) {
+ log_err(scn->dev,
+"%s: invalid parametric coordinates %g.\n"
+"u + (1-u) must be equal to 1 with u [0, 1].\n",
+ FUNC_NAME, uv[0]);
+ res = RES_BAD_ARG;
+ goto error;
+ }
+ } else {
+ const double w = CLAMP(1 - uv[0] - uv[1], 0, 1);
+ if(uv[0] < 0 || uv[1] < 0 || uv[0] > 1 || uv[1] > 1
+ || !eq_eps(w + uv[0] + uv[1], 1, 1.e-6)) {
+ log_err(scn->dev,
+"%s: invalid parametric coordinates [%g, %g].\n"
+"u + v + (1-u-v) must be equal to 1 with u and v in [0, 1].\n",
+ FUNC_NAME, uv[0], uv[1]);
+ 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;
+
+ /* Here we go! Launch the Monte Carlo estimation */
+ omp_set_num_threads((int)scn->dev->nthreads);
+ #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 = boundary_realisation_2d
+ (scn, rng, iprim, uv, time, fp_to_meter, Tarad, Tref, &w);
+ } else {
+ res_local = boundary_realisation_3d
+ (scn, rng, iprim, uv, 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;
+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,
diff --git a/src/sdis_solve_Xd.h b/src/sdis_solve_Xd.h
@@ -76,6 +76,8 @@ reflect(float res[3], const float V[3], const float N[3])
#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_POSITION CONCAT(CONCAT(S, DIM), D_POSITION)
+#define SXD_GEOMETRY_NORMAL CONCAT(CONCAT(S, DIM), D_GEOMETRY_NORMAL)
#define SXD CONCAT(CONCAT(S, DIM), D)
/* Vector macros generic to SDIS_SOLVE_DIMENSION */
@@ -270,6 +272,7 @@ XD(trace_radiative_path)
r = ssp_rng_canonical(rng);
if(r < epsilon) {
T->func = XD(boundary_temperature);
+ rwalk->mdm = NULL; /* The random walk is at an interface between 2 media */
break;
}
@@ -298,6 +301,7 @@ XD(trace_radiative_path)
}
}
+
exit:
return res;
error:
@@ -541,6 +545,7 @@ XD(boundary_temperature)
const struct sdis_medium* mdm_back = NULL;
double tmp;
ASSERT(scn && fp_to_meter > 0 && ctx && rwalk && rng && T);
+ ASSERT(rwalk->mdm == NULL);
ASSERT(!SXD_HIT_NONE(&rwalk->hit));
XD(setup_interface_fragment)(&frag, &rwalk->vtx, &rwalk->hit);
@@ -705,6 +710,7 @@ XD(solid_temperature)
} while(SXD_HIT_NONE(&rwalk->hit));
T->func = XD(boundary_temperature);
+ rwalk->mdm = NULL; /* The random walk is at an interface between 2 media */
return RES_OK;
}
@@ -785,6 +791,71 @@ XD(probe_realisation)
return RES_OK;
}
+static res_T
+XD(boundary_realisation)
+ (struct sdis_scene* scn,
+ struct ssp_rng* rng,
+ const size_t iprim,
+ const double uv[DIM],
+ const double time,
+ const double fp_to_meter,
+ const double Tarad,
+ const double Tref,
+ double* weight)
+{
+ struct rwalk_context ctx;
+ struct XD(rwalk) rwalk = XD(RWALK_NULL);
+ struct XD(temperature) T = XD(TEMPERATURE_NULL);
+ struct sXd(attrib) attr;
+#if SDIS_SOLVE_DIMENSION == 2
+ float st;
+#else
+ float st[2];
+#endif
+ res_T res = RES_OK;
+ ASSERT(uv && fp_to_meter > 0 && weight && time >= 0);
+
+ T.func = XD(boundary_temperature);
+
+ rwalk.hit.distance = 0;
+ rwalk.vtx.time = time;
+ rwalk.mdm = NULL; /* The random walk is at an interface between 2 media */
+
+#if SDIS_SOLVE_DIMENSION == 2
+ st = (float)uv[0];
+#else
+ f2_set_d2(st, uv);
+#endif
+
+ /* Fetch the primitive */
+ SXD(scene_view_get_primitive
+ (scn->sXd(view), (unsigned int)iprim, &rwalk.hit.prim));
+
+ /* Retrieve the world space position of the probe onto the primitive */
+ SXD(primitive_get_attrib(&rwalk.hit.prim, SXD_POSITION, st, &attr));
+ dX_set_fX(rwalk.vtx.P, attr.value);
+
+ /* Retrieve the primitive normal */
+ SXD(primitive_get_attrib(&rwalk.hit.prim, SXD_GEOMETRY_NORMAL, st, &attr));
+ fX(set)(rwalk.hit.normal, attr.value);
+
+#if SDIS_SOLVE_DIMENSION==2
+ rwalk.hit.u = st;
+#else
+ f2_set(rwalk.hit.uv, st);
+#endif
+
+ ctx.Tarad = Tarad;
+ ctx.Tref3 = Tref*Tref*Tref;
+
+ 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)
@@ -840,6 +911,8 @@ error:
#undef SXD_HIT_NONE
#undef SXD_HIT_NULL
#undef SXD_HIT_NULL__
+#undef SXD_POSITION
+#undef SXD_GEOMETRY_NORMAL
#undef SXD
#undef dX
#undef fX
diff --git a/src/test_sdis_utils.h b/src/test_sdis_utils.h
@@ -104,6 +104,7 @@ static const struct sdis_solid_shader DUMMY_SOLID_SHADER = {
dummy_medium_getter,
dummy_medium_getter,
dummy_medium_getter,
+ dummy_medium_getter,
dummy_medium_getter
};
