commit 2296ae2a67e8103733d6438c2f9db0c854fc8b94
parent 47d1debedbcdff637fde53217199f06e493a0af7
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Tue, 5 Mar 2019 12:37:14 +0100
Implement the medium solver (untested version)
Diffstat:
7 files changed, 374 insertions(+), 9 deletions(-)
diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt
@@ -93,6 +93,7 @@ set(SDIS_FILES_INC
sdis_scene_c.h
sdis_scene_Xd.h
sdis_solve_Xd.h
+ sdis_solve_medium_Xd.h
sdis_Xd_begin.h
sdis_Xd_end.h)
diff --git a/src/sdis.h b/src/sdis.h
@@ -913,6 +913,17 @@ sdis_solve_camera
sdis_write_accums_T writer,
void* writer_data);
+SDIS_API res_T
+sdis_solve_medium
+ (struct sdis_scene* sanc,
+ const size_t nrealisations, /* #realisations */
+ struct sdis_medium* medium, /* Medium to solve */
+ const double time_range[2], /* 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);
+
/*******************************************************************************
* Green solvers.
*
diff --git a/src/sdis_scene_Xd.h b/src/sdis_scene_Xd.h
@@ -558,7 +558,7 @@ XD(setup_enclosure_geometry)(struct sdis_scene* scn, struct sencXd(enclosure)* e
enc_data = htable_enclosure_find(&scn->enclosures, &header.enclosure_id);
ASSERT(enc_data != NULL);
- /* Setup the vertex data */
+ /* Setup the vertex data */
vdata.usage = SXD_POSITION;
#if DIM == 2
vdata.type = S2D_FLOAT2;
@@ -590,11 +590,12 @@ XD(setup_enclosure_geometry)(struct sdis_scene* scn, struct sencXd(enclosure)* e
CALL(s3d_scene_view_compute_area(enc_data->s3d_view, &S));
CALL(s3d_scene_view_compute_volume(enc_data->s3d_view, &V));
#endif
+ enc_data->V = V;
enc_data->S_over_V = S / V;
ASSERT(enc_data->S_over_V >= 0);
#undef CALL
- /* Set enclosure hc upper bound regardless of its media being a fluid */
+ /* Set enclosure hc upper bound regardless of its media being a fluid */
p_ub = htable_d_find(&scn->tmp_hc_ub, &header.enclosure_id);
ASSERT(p_ub);
enc_data->hc_upper_bound = *p_ub;
@@ -612,6 +613,9 @@ XD(setup_enclosure_geometry)(struct sdis_scene* scn, struct sencXd(enclosure)* e
#endif
}
+ /* Setup the medium id of the enclosure */
+ SENCXD(enclosure_get_medium(enc, 0, &enc_data->medium_id));
+
exit:
enclosure_release(&enc_dummy);
if(sXd_shape) SXD(shape_ref_put(sXd_shape));
@@ -642,7 +646,6 @@ XD(setup_enclosures)(struct sdis_scene* scn, struct sencXd(descriptor)* desc)
#else
struct senc_enclosure_header header;
#endif
- const struct sdis_medium* mdm;
SENCXD(descriptor_get_enclosure(desc, ienc, &enc));
SENCXD(enclosure_get_header(enc, &header));
@@ -687,13 +690,10 @@ XD(setup_enclosures)(struct sdis_scene* scn, struct sencXd(descriptor)* desc)
}
SENCXD(enclosure_get_medium(enc, 0, &enclosed_medium));
- if(res != RES_OK) goto error;
ASSERT(enclosed_medium < darray_medium_size_get(&scn->media));
- mdm = darray_medium_cdata_get(&scn->media)[enclosed_medium];
- ASSERT(mdm);
- /* Silently discard the solid and infinite enclosures */
- if(mdm->type == SDIS_FLUID && !header.is_infinite) {
+ /* Silently discard infinite enclosures */
+ if(!header.is_infinite) {
res = XD(setup_enclosure_geometry)(scn, enc);
if(res != RES_OK) goto error;
}
diff --git a/src/sdis_scene_c.h b/src/sdis_scene_c.h
@@ -78,6 +78,9 @@ struct enclosure {
double hc_upper_bound;
double S_over_V; /* in 3D = surface/volume; in 2D = perimeter/area */
+ double V; /* 3D = volume; 2D = area; */
+
+ unsigned medium_id;
};
static INLINE void
@@ -88,6 +91,7 @@ enclosure_init(struct mem_allocator* allocator, struct enclosure* enc)
enc->s3d_view = NULL;
darray_uint_init(allocator, &enc->local2global);
enc->S_over_V = 0;
+ enc->V = 0;
enc->hc_upper_bound = 0;
}
@@ -111,6 +115,7 @@ enclosure_copy(struct enclosure* dst, const struct enclosure* src)
dst->s2d_view = src->s2d_view;
}
dst->S_over_V = src->S_over_V;
+ dst->V = src->V;
dst->hc_upper_bound = src->hc_upper_bound;
return darray_uint_copy(&dst->local2global, &src->local2global);
}
@@ -132,6 +137,7 @@ enclosure_copy_and_release(struct enclosure* dst, struct enclosure* src)
src->s2d_view = NULL;
}
dst->S_over_V = src->S_over_V;
+ dst->V = src->V;
dst->hc_upper_bound = src->hc_upper_bound;
return RES_OK;
}
diff --git a/src/sdis_solve.c b/src/sdis_solve.c
@@ -25,6 +25,12 @@
#define SDIS_XD_DIMENSION 3
#include "sdis_solve_Xd.h"
+/* Generate the medium solvers */
+#define SDIS_XD_DIMENSION 2
+#include "sdis_solve_medium_Xd.h"
+#define SDIS_XD_DIMENSION 3
+#include "sdis_solve_medium_Xd.h"
+
#include <star/ssp.h>
#include <omp.h>
@@ -430,3 +436,26 @@ error:
goto exit;
}
+res_T
+sdis_solve_medium
+ (struct sdis_scene* scn,
+ const size_t nrealisations, /* #realisations */
+ struct sdis_medium* medium, /* Medium to solve */
+ const double time_range[2], /* 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** estimator)
+{
+ res_T res = RES_OK;
+ if(!scn) return RES_BAD_ARG;
+ if(scene_is_2d(scn)) {
+ res = solve_medium_2d(scn, nrealisations, medium, time_range, fp_to_meter, Tarad,
+ Tref, estimator);
+ } else {
+ res = solve_medium_3d(scn, nrealisations, medium, time_range, fp_to_meter, Tarad,
+ Tref, estimator);
+ }
+ return res;
+}
+
diff --git a/src/sdis_solve_Xd.h b/src/sdis_solve_Xd.h
@@ -33,7 +33,6 @@ struct XD(boundary_context) {
static const struct XD(boundary_context) XD(BOUNDARY_CONTEXT_NULL) = {
NULL, NULL
};
-
/*******************************************************************************
* Helper functions
******************************************************************************/
diff --git a/src/sdis_solve_medium_Xd.h b/src/sdis_solve_medium_Xd.h
@@ -0,0 +1,319 @@
+/* Copyright (C) 2016-2019 |Meso|Star> (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_device_c.h"
+#include "sdis_estimator_c.h"
+#include "sdis_realisation.h"
+#include "sdis_scene_c.h"
+
+#include <rsys/algorithm.h>
+#include <rsys/dynamic_array.h>
+
+#include "sdis_Xd_begin.h"
+
+#ifndef SDIS_SOLVE_MEDIUM_XD_H
+#define SDIS_SOLVE_MEDIUM_XD_H
+
+/*
+ * Define the data structures and functions that are not generic to the
+ * SDIS_XD_DIMENSION parameter
+ */
+
+struct enclosure_cumul {
+ const struct enclosure* enc;
+ double cumul;
+};
+
+/* Define the darray_enclosure_cumul dynamic array */
+#define DARRAY_NAME enclosure_cumul
+#define DARRAY_DATA struct enclosure_cumul
+#include <rsys/dynamic_array.h>
+
+/*******************************************************************************
+ * Helper functions
+ ******************************************************************************/
+static INLINE int
+cmp_double_to_enc_cumuls(const void* a, const void* b)
+{
+ const double key = *(const double*)a;
+ const struct enclosure_cumul* enc_cumul = (const struct enclosure_cumul*)b;
+ if(key < enc_cumul->cumul) return -1;
+ if(key > enc_cumul->cumul) return +1;
+ return 0;
+}
+
+static res_T
+compute_medium_enclosure_cumulative
+ (struct sdis_scene* scn,
+ const struct sdis_medium* mdm,
+ struct darray_enclosure_cumul* cumul)
+{
+ struct htable_enclosure_iterator it, end;
+ double accum = 0;
+ res_T res = RES_OK;
+ ASSERT(scn && mdm && cumul);
+
+ darray_enclosure_cumul_clear(cumul);
+
+ htable_enclosure_begin(&scn->enclosures, &it);
+ htable_enclosure_end(&scn->enclosures, &end);
+ while(!htable_enclosure_iterator_eq(&it, &end)) {
+ struct enclosure_cumul enc_cumul;
+ const struct enclosure* enc = htable_enclosure_iterator_data_get(&it);
+ htable_enclosure_iterator_next(&it);
+
+ if(sdis_medium_get_id(mdm) != enc->medium_id) continue;
+
+ accum += enc->V + accum;
+ enc_cumul.enc = enc;
+ enc_cumul.cumul = accum;
+ res = darray_enclosure_cumul_push_back(cumul, &enc_cumul);
+ if(res != RES_OK) goto error;
+ }
+exit:
+ return res;
+error:
+ darray_enclosure_cumul_clear(cumul);
+ goto exit;
+}
+
+static const struct enclosure*
+sample_medium_enclosure
+ (const struct darray_enclosure_cumul* cumul, struct ssp_rng* rng)
+{
+ const struct enclosure_cumul* enc_cumuls = NULL;
+ const struct enclosure_cumul* enc_cumul_found = NULL;
+ double r;
+ size_t i;
+ size_t sz;
+ ASSERT(cumul && rng && darray_enclosure_cumul_size_get(cumul));
+
+ sz = darray_enclosure_cumul_size_get(cumul);
+ enc_cumuls = darray_enclosure_cumul_cdata_get(cumul);
+ if(sz == 1) {
+ enc_cumul_found = enc_cumuls;
+ } else {
+ /* Generate an uniform random number in [0, cumul[ */
+ r = ssp_rng_canonical(rng);
+ r = r * enc_cumuls[sz-1].cumul;
+
+ enc_cumul_found = search_lower_bound
+ (&r, enc_cumuls, sz, sizeof(*enc_cumuls), cmp_double_to_enc_cumuls);
+ ASSERT(enc_cumul_found);
+
+ /* search_lower_bound returns the first entry that is not less than `r'.
+ * The following code discards entries that are also equal to `r'. */
+ i = (size_t)(enc_cumul_found - enc_cumuls);
+ while(enc_cumuls[i].cumul == r && i < sz) ++i;
+ ASSERT(i < sz);
+
+ enc_cumul_found = enc_cumuls + i;
+ }
+ return enc_cumul_found->enc;
+}
+
+#endif /* !SDIS_SOLVE_MEDIUM_XD_H */
+
+/*******************************************************************************
+ * Helper functions generic to the SDIS_XD_DIMENSION parameter
+ ******************************************************************************/
+static res_T
+XD(sample_enclosure_position)
+ (const struct enclosure* enc,
+ struct ssp_rng* rng,
+ double pos[DIM])
+{
+ const size_t MAX_NCHALLENGES = 1000;
+ float lower[DIM], upper[DIM];
+ size_t ichallenge;
+ size_t i;
+ res_T res = RES_OK;
+ ASSERT(enc && rng && pos);
+
+ SXD(scene_view_get_aabb(enc->sXd(view), lower, upper));
+
+ FOR_EACH(i, 0, DIM) {
+ if(lower[i] > upper[i] || eq_epsf(lower[i], upper[i], 1.e-6f)) {
+ res = RES_BAD_ARG; /* Degenerated enclosure */
+ goto error;
+ }
+ }
+
+ FOR_EACH(ichallenge, 0, MAX_NCHALLENGES) {
+ struct sXd(hit) hit = SXD_HIT_NULL;
+ const float dir[3] = {1,0,0};
+ const float range[2] = {0, FLT_MAX};
+ float org[DIM];
+
+ /* Generate an uniform position into the enclosure AABB */
+ FOR_EACH(i, 0, DIM) {
+ org[i] = ssp_rng_uniform_float(rng, lower[i], upper[i]);
+ }
+
+ /* Check that pos lies into the enclosure; trace a ray and check that it
+ * hits something and that the normal points towards the traced ray
+ * direction (enclosure normals point inword the enclosure) */
+ SXD(scene_view_trace_ray(enc->sXd(view), org, dir, range, NULL, &hit));
+ if(!SXD_HIT_NONE(&hit) && fX(dot)(dir, hit.normal) < 0) {
+ dX_set_fX(pos, org);
+ break;
+ }
+ }
+
+ if(ichallenge >= MAX_NCHALLENGES) {
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
+/*******************************************************************************
+ * Local function
+ ******************************************************************************/
+static res_T
+XD(solve_medium)
+ (struct sdis_scene* scn,
+ const size_t nrealisations,
+ struct sdis_medium* mdm,
+ const double time_range[2],
+ 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)
+{
+ struct darray_enclosure_cumul cumul;
+ struct ssp_rng_proxy* rng_proxy = NULL;
+ struct ssp_rng** rngs = NULL;
+ struct sdis_estimator* estimator = NULL;
+ struct accum { double sum, sum2; size_t naccums; }* accums = NULL;
+ int64_t irealisation;
+ int cumul_is_init = 0;
+ size_t i;
+ ATOMIC res = RES_OK;
+
+ if(!scn || !mdm || !nrealisations || nrealisations > INT64_MAX
+ || fp_to_meter <= 0 || Tref < 0 || !out_estimator) {
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+#if SDIS_XD_DIMENSION == 2
+ if(scene_is_2d(scn) == 0) { res = RES_BAD_ARG; goto error; }
+#else
+ if(scene_is_2d(scn) != 0) { res = RES_BAD_ARG; goto error; }
+#endif
+
+ /* 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(*rngs));
+ 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 per thread accumulator */
+ accums = MEM_CALLOC(scn->dev->allocator, scn->dev->nthreads, sizeof(*accums));
+ if(!accums) { res = RES_MEM_ERR; goto error; }
+
+ /* Create the estimator */
+ res = estimator_create(scn->dev, SDIS_ESTIMATOR_TEMPERATURE, &estimator);
+ if(res != RES_OK) goto error;
+
+ /* Compute the enclosure cumulative */
+ darray_enclosure_cumul_init(scn->dev->allocator, &cumul);
+ cumul_is_init = 1;
+ res = compute_medium_enclosure_cumulative(scn, mdm, &cumul);
+ if(res != RES_OK) goto error;
+
+ omp_set_num_threads((int)scn->dev->nthreads);
+ #pragma omp parallel for schedule(static)
+ for(irealisation = 0; irealisation < (int64_t)nrealisations; ++irealisation) {
+ const int ithread = omp_get_thread_num();
+ const struct enclosure* enc = NULL;
+ struct accum* accum = accums + ithread;
+ struct ssp_rng* rng = rngs[ithread];
+ double weight;
+ double time;
+ double pos[DIM];
+ res_T res_local = RES_OK;
+
+ if(ATOMIC_GET(&res) != RES_OK) continue; /* An error occurred */
+
+ /* Sample the time */
+ time = sample_time(rng, time_range);
+
+ /* Uniformly Sample an enclosure that surround the submitted medium and
+ * uniformly sample a position into it */
+ enc = sample_medium_enclosure(&cumul, rng);
+ res_local = XD(sample_enclosure_position)(enc, rng, pos);
+ if(res_local != RES_OK) {
+ log_err(scn->dev, "%s: could not sample a medium position.\n", FUNC_NAME);
+ ATOMIC_SET(&res, res_local);
+ continue;
+ }
+
+ /* Run a probe realisation */
+ res_local = XD(probe_realisation)((size_t)irealisation, scn, rng, mdm, pos,
+ time, fp_to_meter, Tarad, Tref, NULL, NULL, &weight);
+ if(res_local != RES_OK) {
+ if(res_local != RES_BAD_OP) { ATOMIC_SET(&res, res_local); continue; }
+ } else {
+ accum->sum += weight;
+ accum->sum2 += weight*weight;
+ ++accum->naccums;
+ }
+ }
+ if(res != RES_OK) goto error;
+
+ /* Merge the per thread accumulators into the accumulator of the thread 0 */
+ FOR_EACH(i, 1, scn->dev->nthreads) {
+ accums[0].sum += accums[i].sum;
+ accums[0].sum2 += accums[i].sum2;
+ accums[0].naccums += accums[i].naccums;
+ }
+ ASSERT(accums[0].naccums <= nrealisations);
+
+ /* Setup the estimated temperature */
+ estimator_setup_realisations_count(estimator, nrealisations, accums[0].naccums);
+ estimator_setup_temperature(estimator, accums[0].sum, accums[0].sum2);
+
+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(accums) MEM_RM(scn->dev->allocator, accums);
+ if(cumul_is_init) darray_enclosure_cumul_release(&cumul);
+ 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;
+}
+
+#include "sdis_Xd_end.h"
