commit db9607f292dcbc9562a5b95c16b30165f64e1ba2
parent 47d1debedbcdff637fde53217199f06e493a0af7
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Wed, 6 Mar 2019 15:28:33 +0100
Merge branch 'feature_medium_temperature' into develop
Diffstat:
10 files changed, 1277 insertions(+), 12 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)
@@ -184,6 +185,8 @@ if(NOT NO_TEST)
new_test(test_sdis_solve_probe3_2d)
new_test(test_sdis_solve_boundary)
new_test(test_sdis_solve_boundary_flux)
+ new_test(test_sdis_solve_medium)
+ new_test(test_sdis_solve_medium_2d)
new_test(test_sdis_volumic_power)
new_test(test_sdis_volumic_power4_2d)
@@ -198,6 +201,7 @@ if(NOT NO_TEST)
add_test(test_sdis_volumic_power3_2d test_sdis_volumic_power3_2d)
endif()
+ target_link_libraries(test_sdis_solve_medium Star3DUT)
target_link_libraries(test_sdis_solve_probe3 Star3DUT)
target_link_libraries(test_sdis_solve_probe3_2d ${MATH_LIB})
target_link_libraries(test_sdis_solve_camera Star3DUT)
diff --git a/src/sdis.h b/src/sdis.h
@@ -681,6 +681,15 @@ sdis_scene_get_dimension
(const struct sdis_scene* scn,
enum sdis_scene_dimension* dim);
+/* Return the area/volume of occupied by a medium in a 2D/3D scene. Only
+ * enclosed media are handled, i.e. media whose border are explicitly defined
+ * by a geometry. */
+SDIS_API res_T
+sdis_scene_get_medium_spread
+ (struct sdis_scene* scn,
+ const struct sdis_medium* mdm,
+ double* spread);
+
/*******************************************************************************
* An estimator stores the state of a simulation
******************************************************************************/
@@ -913,6 +922,18 @@ sdis_solve_camera
sdis_write_accums_T writer,
void* writer_data);
+SDIS_API 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 ambient_radiative_temperature, /* In Kelvin */
+ const double reference_temperature, /* In Kelvin */
+ const int register_path, /* Combination of enum sdis_heat_path_flag */
+ struct sdis_estimator** estimator);
+
/*******************************************************************************
* Green solvers.
*
@@ -927,7 +948,7 @@ sdis_solve_camera
* media must be constant in time and space too. Furthermore, note that only
* the interfaces/media that had a flux/volumic power during green estimation
* can update their flux/volumic power value for subsequent
- * sdis_green_function_solve invocations : other interfaces/media are
+ * sdis_green_function_solve invocations : others interfaces/media are
* definitely registered against the green function as interfaces/media with no
* flux/volumic power.
*
diff --git a/src/sdis_scene.c b/src/sdis_scene.c
@@ -312,13 +312,45 @@ sdis_scene_release_analysis(struct sdis_scene* scn)
res_T
sdis_scene_get_dimension
- (const struct sdis_scene* scn, enum sdis_scene_dimension* dim)
+ (const struct sdis_scene* scn, enum sdis_scene_dimension* dim)
{
if(!scn || !dim) return RES_BAD_ARG;
*dim = scene_is_2d(scn) ? SDIS_SCENE_2D : SDIS_SCENE_3D;
return RES_OK;
}
+res_T
+sdis_scene_get_medium_spread
+ (struct sdis_scene* scn,
+ const struct sdis_medium* mdm,
+ double* out_spread)
+{
+ struct htable_enclosure_iterator it, end;
+ double spread = 0;
+ res_T res = RES_OK;
+
+ if(!scn || !mdm || !out_spread) {
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ htable_enclosure_begin(&scn->enclosures, &it);
+ htable_enclosure_end(&scn->enclosures, &end);
+ while(!htable_enclosure_iterator_eq(&it, &end)) {
+ const struct enclosure* enc = htable_enclosure_iterator_data_get(&it);
+ htable_enclosure_iterator_next(&it);
+ if(sdis_medium_get_id(mdm) == enc->medium_id) {
+ spread += enc->V;
+ }
+ }
+ *out_spread = spread;
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
/*******************************************************************************
* Local miscellaneous function
******************************************************************************/
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;
@@ -580,7 +580,7 @@ XD(setup_enclosure_geometry)(struct sdis_scene* scn, struct sencXd(enclosure)* e
#endif
CALL(sXd(scene_create)(sXd_dev, &sXd_scn));
CALL(sXd(scene_attach_shape)(sXd_scn, sXd_shape));
- CALL(sXd(scene_view_create)(sXd_scn, SXD_SAMPLE, &enc_data->sXd(view)));
+ CALL(sXd(scene_view_create)(sXd_scn, SXD_SAMPLE|SXD_TRACE, &enc_data->sXd(view)));
/* Compute the S/V ratio */
#if DIM == 2
@@ -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,27 @@ 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 */
+ const int register_paths, /* Combination of enum sdis_heat_path_flag */
+ 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, register_paths, estimator);
+ } else {
+ res = solve_medium_3d(scn, nrealisations, medium, time_range, fp_to_meter, Tarad,
+ Tref, register_paths, 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,353 @@
+/* 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;
+ enc_cumul.enc = enc;
+ enc_cumul.cumul = accum;
+ res = darray_enclosure_cumul_push_back(cumul, &enc_cumul);
+ if(res != RES_OK) goto error;
+ }
+
+ if(darray_enclosure_cumul_size_get(cumul) == 0) {
+ log_err(scn->dev,
+ "%s: there is no enclosure that encompasses the submitted medium.\n",
+ FUNC_NAME);
+ res = RES_BAD_ARG;
+ 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 */
+ const int register_paths, /* Combination of enum sdis_heat_path_flag */
+ 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
+ || !time_range || time_range[0] > time_range[1] || 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];
+ struct sdis_heat_path* pheat_path = NULL;
+ struct sdis_heat_path heat_path;
+ 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);
+
+ /* Prepare path registration if necessary */
+ if(register_paths) {
+ heat_path_init(scn->dev->allocator, &heat_path);
+ pheat_path = &heat_path;
+ }
+
+ /* 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, pheat_path, &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;
+ }
+
+ /* Finalize the registered path */
+ if(pheat_path) {
+ pheat_path->status = res_local == RES_OK
+ ? SDIS_HEAT_PATH_SUCCEED
+ : SDIS_HEAT_PATH_FAILED;
+
+ /* Check if the path must be saved regarding the register_paths mask */
+ if(!(register_paths & (int)pheat_path->status)) {
+ heat_path_release(pheat_path);
+ } else { /* Register the sampled path */
+ res_local = estimator_add_and_release_heat_path(estimator, pheat_path);
+ if(res_local != RES_OK) { ATOMIC_SET(&res, res_local); continue; }
+ }
+ }
+ }
+ 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"
diff --git a/src/test_sdis_solve_medium.c b/src/test_sdis_solve_medium.c
@@ -0,0 +1,423 @@
+/* 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.h"
+#include "test_sdis_utils.h"
+
+#include <rsys/math.h>
+#include <rsys/stretchy_array.h>
+#include <star/s3dut.h>
+
+#include <string.h>
+
+#define Tf0 300.0
+#define Tf1 330.0
+#define N 1000ul /* #realisations */
+
+/*
+ * The scene is composed of 2 super shapes whose temperature is unknown. The
+ * first super shape is surrounded by a fluid whose temperature is Tf0 while
+ * the second one is in fluid whose temperature is Tf1. The temperatures of the
+ * super shape 0 and 1 are thus uniform and equal to Tf0 and Tf1, respectively.
+ *
+ * This program performs 2 tests. In the first one, the super shapes 0 and 1
+ * have different media; the medium solver thus estimates the
+ * temperature of one super shape. In the second test, the scene is updated to
+ * use the same medium for the 2 super shapes. In this case, when invoked on
+ * the right medium, the estimated temperature T is equal to :
+ *
+ * T = Tf0 * V0/(V0 + V1) + Tf1 * V1/(V0 + V1)
+ *
+ * with V0 and V1 the volume of the super shapes 0 and 1, respectively.
+ */
+
+/*******************************************************************************
+ * Geometry
+ ******************************************************************************/
+struct context {
+ struct s3dut_mesh_data msh0;
+ struct s3dut_mesh_data msh1;
+ struct sdis_interface* interf0;
+ struct sdis_interface* interf1;
+};
+
+static void
+get_indices(const size_t itri, size_t ids[3], void* context)
+{
+ const struct context* ctx = context;
+ /* Note that we swap the indices to ensure that the triangle normals point
+ * inward the super shape */
+ if(itri < ctx->msh0.nprimitives) {
+ ids[0] = ctx->msh0.indices[itri*3+0];
+ ids[2] = ctx->msh0.indices[itri*3+1];
+ ids[1] = ctx->msh0.indices[itri*3+2];
+ } else {
+ const size_t itri2 = itri - ctx->msh0.nprimitives;
+ ids[0] = ctx->msh1.indices[itri2*3+0] + ctx->msh0.nvertices;
+ ids[2] = ctx->msh1.indices[itri2*3+1] + ctx->msh0.nvertices;
+ ids[1] = ctx->msh1.indices[itri2*3+2] + ctx->msh0.nvertices;
+ }
+}
+
+static void
+get_position(const size_t ivert, double pos[3], void* context)
+{
+ const struct context* ctx = context;
+ if(ivert < ctx->msh0.nvertices) {
+ pos[0] = ctx->msh0.positions[ivert*3+0] - 2.0;
+ pos[1] = ctx->msh0.positions[ivert*3+1];
+ pos[2] = ctx->msh0.positions[ivert*3+2];
+ } else {
+ const size_t ivert2 = ivert - ctx->msh0.nvertices;
+ pos[0] = ctx->msh1.positions[ivert2*3+0] + 2.0;
+ pos[1] = ctx->msh1.positions[ivert2*3+1];
+ pos[2] = ctx->msh1.positions[ivert2*3+2];
+ }
+}
+
+static void
+get_interface(const size_t itri, struct sdis_interface** bound, void* context)
+{
+ const struct context* ctx = context;
+ *bound = itri < ctx->msh0.nprimitives ? ctx->interf0 : ctx->interf1;
+}
+
+/*******************************************************************************
+ * Fluid medium
+ ******************************************************************************/
+struct fluid {
+ double temperature;
+};
+
+static double
+fluid_get_temperature
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ CHK(data != NULL && vtx != NULL);
+ return ((const struct fluid*)sdis_data_cget(data))->temperature;
+}
+
+/*******************************************************************************
+ * Solid medium
+ ******************************************************************************/
+struct solid {
+ double cp;
+ double lambda;
+ double rho;
+ double delta;
+ double temperature;
+};
+
+static double
+solid_get_calorific_capacity
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ CHK(data != NULL && vtx != NULL);
+ return ((const struct solid*)sdis_data_cget(data))->cp;
+}
+
+static double
+solid_get_thermal_conductivity
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ CHK(data != NULL && vtx != NULL);
+ return ((const struct solid*)sdis_data_cget(data))->lambda;
+}
+
+static double
+solid_get_volumic_mass
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ CHK(data != NULL && vtx != NULL);
+ return ((const struct solid*)sdis_data_cget(data))->rho;
+}
+
+static double
+solid_get_delta
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ CHK(data != NULL && vtx != NULL);
+ return ((const struct solid*)sdis_data_cget(data))->delta;
+}
+
+static double
+solid_get_temperature
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ CHK(data != NULL && vtx != NULL);
+ return ((const struct solid*)sdis_data_cget(data))->temperature;
+}
+
+struct interf {
+ double hc;
+ double epsilon;
+ double specular_fraction;
+};
+
+static double
+interface_get_convection_coef
+ (const struct sdis_interface_fragment* frag, struct sdis_data* data)
+{
+ CHK(data != NULL && frag != NULL);
+ return ((const struct interf*)sdis_data_cget(data))->hc;
+}
+
+static double
+interface_get_emissivity
+ (const struct sdis_interface_fragment* frag, struct sdis_data* data)
+{
+ CHK(data != NULL && frag != NULL);
+ return ((const struct interf*)sdis_data_cget(data))->epsilon;
+}
+
+static double
+interface_get_specular_fraction
+ (const struct sdis_interface_fragment* frag, struct sdis_data* data)
+{
+ CHK(data != NULL && frag != NULL);
+ return ((const struct interf*)sdis_data_cget(data))->specular_fraction;
+}
+
+/*******************************************************************************
+ * Test
+ ******************************************************************************/
+int
+main(int argc, char** argv)
+{
+ struct mem_allocator allocator;
+ struct s3dut_super_formula f0 = S3DUT_SUPER_FORMULA_NULL;
+ struct s3dut_super_formula f1 = S3DUT_SUPER_FORMULA_NULL;
+ struct s3dut_mesh* msh0 = NULL;
+ struct s3dut_mesh* msh1 = NULL;
+ struct sdis_mc T = SDIS_MC_NULL;
+ struct sdis_device* dev = NULL;
+ struct sdis_medium* solid0 = NULL;
+ struct sdis_medium* solid1 = NULL;
+ struct sdis_medium* fluid0 = NULL;
+ struct sdis_medium* fluid1 = NULL;
+ struct sdis_interface* solid0_fluid0 = NULL;
+ struct sdis_interface* solid0_fluid1 = NULL;
+ struct sdis_interface* solid1_fluid1 = NULL;
+ struct sdis_scene* scn = NULL;
+ struct sdis_data* data = NULL;
+ struct sdis_estimator* estimator = NULL;
+ struct fluid* fluid_param = NULL;
+ struct solid* solid_param = NULL;
+ struct interf* interface_param = NULL;
+ struct sdis_fluid_shader fluid_shader = DUMMY_FLUID_SHADER;
+ struct sdis_solid_shader solid_shader = DUMMY_SOLID_SHADER;
+ struct sdis_interface_shader interface_shader = SDIS_INTERFACE_SHADER_NULL;
+ struct context ctx;
+ const double trange[2] = {0, INF};
+ double ref;
+ double v, v0, v1;
+ size_t nreals;
+ size_t nfails;
+ size_t ntris;
+ size_t nverts;
+ (void)argc, (void)argv;
+
+ OK(mem_init_proxy_allocator(&allocator, &mem_default_allocator));
+ OK(sdis_device_create(NULL, &allocator, SDIS_NTHREADS_DEFAULT, 1, &dev));
+
+ fluid_shader.temperature = fluid_get_temperature;
+
+ /* Create the fluid0 medium */
+ OK(sdis_data_create
+ (dev, sizeof(struct fluid), ALIGNOF(struct fluid), NULL, &data));
+ fluid_param = sdis_data_get(data);
+ fluid_param->temperature = Tf0;
+ OK(sdis_fluid_create(dev, &fluid_shader, data, &fluid0));
+ OK(sdis_data_ref_put(data));
+
+ /* Create the fluid1 medium */
+ OK(sdis_data_create
+ (dev, sizeof(struct fluid), ALIGNOF(struct fluid), NULL, &data));
+ fluid_param = sdis_data_get(data);
+ fluid_param->temperature = Tf1;
+ OK(sdis_fluid_create(dev, &fluid_shader, data, &fluid1));
+ OK(sdis_data_ref_put(data));
+
+ /* Setup the solid shader */
+ solid_shader.calorific_capacity = solid_get_calorific_capacity;
+ solid_shader.thermal_conductivity = solid_get_thermal_conductivity;
+ solid_shader.volumic_mass = solid_get_volumic_mass;
+ solid_shader.delta_solid = solid_get_delta;
+ solid_shader.temperature = solid_get_temperature;
+
+ /* Create the solid0 medium */
+ OK(sdis_data_create
+ (dev, sizeof(struct solid), ALIGNOF(struct solid), NULL, &data));
+ solid_param = sdis_data_get(data);
+ solid_param->cp = 1.0;
+ solid_param->lambda = 0.1;
+ solid_param->rho = 1.0;
+ solid_param->delta = 1.0/20.0;
+ solid_param->temperature = -1; /* Unknown temperature */
+ OK(sdis_solid_create(dev, &solid_shader, data, &solid0));
+ OK(sdis_data_ref_put(data));
+
+ /* Create the solid1 medium */
+ OK(sdis_data_create
+ (dev, sizeof(struct solid), ALIGNOF(struct solid), NULL, &data));
+ solid_param = sdis_data_get(data);
+ solid_param->cp = 1.0;
+ solid_param->lambda = 1.0;
+ solid_param->rho = 1.0;
+ solid_param->delta = 1.0/20.0;
+ solid_param->temperature = -1; /* Unknown temperature */
+ OK(sdis_solid_create(dev, &solid_shader, data, &solid1));
+ OK(sdis_data_ref_put(data));
+
+ /* Create the interfaces */
+ OK(sdis_data_create(dev, sizeof(struct interf),
+ ALIGNOF(struct interf), NULL, &data));
+ interface_param = sdis_data_get(data);
+ interface_param->hc = 0.5;
+ interface_param->epsilon = 0;
+ interface_param->specular_fraction = 0;
+ interface_shader.convection_coef = interface_get_convection_coef;
+ interface_shader.front = SDIS_INTERFACE_SIDE_SHADER_NULL;
+ interface_shader.back.temperature = NULL;
+ interface_shader.back.emissivity = interface_get_emissivity;
+ interface_shader.back.specular_fraction = interface_get_specular_fraction;
+ OK(sdis_interface_create
+ (dev, solid0, fluid0, &interface_shader, data, &solid0_fluid0));
+ OK(sdis_interface_create
+ (dev, solid0, fluid1, &interface_shader, data, &solid0_fluid1));
+ OK(sdis_interface_create
+ (dev, solid1, fluid1, &interface_shader, data, &solid1_fluid1));
+ OK(sdis_data_ref_put(data));
+
+ /* Create the mesh0 */
+ f0.A = 1; f0.B = 1; f0.M = 3; f0.N0 = 1; f0.N1 = 1; f0.N2 = 2;
+ f1.A = 1; f1.B = 1; f1.M = 10; f1.N0 = 1; f1.N1 = 1; f1.N2 = 3;
+ OK(s3dut_create_super_shape(&allocator, &f0, &f1, 1, 64, 32, &msh0));
+ OK(s3dut_mesh_get_data(msh0, &ctx.msh0));
+
+ /* Create the mesh1 */
+ f0.A = 1; f0.B = 1; f0.M = 10; f0.N0 = 1; f0.N1 = 1; f0.N2 = 5;
+ f1.A = 1; f1.B = 1; f1.M = 1; f1.N0 = 1; f1.N1 = 1; f1.N2 = 1;
+ OK(s3dut_create_super_shape(&allocator, &f0, &f1, 1, 64, 32, &msh1));
+ OK(s3dut_mesh_get_data(msh1, &ctx.msh1));
+
+ /* Create the scene */
+ ctx.interf0 = solid0_fluid0;
+ ctx.interf1 = solid1_fluid1;
+ ntris = ctx.msh0.nprimitives + ctx.msh1.nprimitives;
+ nverts = ctx.msh0.nvertices + ctx.msh1.nvertices;
+#if 0
+ {
+ double* vertices = NULL;
+ size_t* indices = NULL;
+ size_t i;
+ CHK(vertices = MEM_CALLOC(&allocator, nverts*3, sizeof(*vertices)));
+ CHK(indices = MEM_CALLOC(&allocator, ntris*3, sizeof(*indices)));
+ FOR_EACH(i, 0, ntris) get_indices(i, indices + i*3, &ctx);
+ FOR_EACH(i, 0, nverts) get_position(i, vertices + i*3, &ctx);
+ dump_mesh(stdout, vertices, nverts, indices, ntris);
+ MEM_RM(&allocator, vertices);
+ MEM_RM(&allocator, indices);
+ }
+#endif
+
+ OK(sdis_scene_create(dev, ntris, get_indices, get_interface, nverts,
+ get_position, &ctx, &scn));
+
+ BA(sdis_scene_get_medium_spread(NULL, solid0, &v0));
+ BA(sdis_scene_get_medium_spread(scn, NULL, &v0));
+ BA(sdis_scene_get_medium_spread(scn, solid0, NULL));
+ OK(sdis_scene_get_medium_spread(scn, solid0, &v0));
+ CHK(v0 > 0);
+ OK(sdis_scene_get_medium_spread(scn, solid1, &v1));
+ CHK(v1 > 0);
+ OK(sdis_scene_get_medium_spread(scn, fluid0, &v));
+ CHK(v == 0);
+ OK(sdis_scene_get_medium_spread(scn, fluid1, &v));
+ CHK(v == 0);
+
+ BA(sdis_solve_medium(NULL, N, solid0, trange, 1.f, -1, 0, 0, &estimator));
+ BA(sdis_solve_medium(scn, 0, solid0, trange, 1.f, -1, 0, 0, &estimator));
+ BA(sdis_solve_medium(scn, N, NULL, trange, 1.f, -1, 0, 0, &estimator));
+ BA(sdis_solve_medium(scn, N, solid0, NULL, 1.f, -1, 0, 0, &estimator));
+ BA(sdis_solve_medium(scn, N, solid0, trange, 0.f, -1, 0, 0, &estimator));
+ BA(sdis_solve_medium(scn, N, solid0, trange, 1.f, -1, 0, 0, NULL));
+ OK(sdis_solve_medium(scn, N, solid0, trange, 1.f, -1, 0, 0, &estimator));
+
+ OK(sdis_estimator_get_realisation_count(estimator, &nreals));
+ OK(sdis_estimator_get_failure_count(estimator, &nfails));
+ OK(sdis_estimator_get_temperature(estimator, &T));
+ printf("Shape0 temperature = "STR(Tf0)" ~ %g +/- %g\n", T.E, T.SE);
+ printf("#failures = %lu/%lu\n", nfails, N);
+ CHK(eq_eps(T.E, Tf0, T.SE));
+ CHK(nreals + nfails == N);
+ OK(sdis_estimator_ref_put(estimator));
+
+ OK(sdis_solve_medium(scn, N, solid1, trange, 1.f, -1, 0, 0, &estimator));
+ OK(sdis_estimator_get_realisation_count(estimator, &nreals));
+ OK(sdis_estimator_get_failure_count(estimator, &nfails));
+ OK(sdis_estimator_get_temperature(estimator, &T));
+ printf("Shape1 temperature = "STR(Tf1)" ~ %g +/- %g\n", T.E, T.SE);
+ printf("#failures = %lu/%lu\n", nfails, N);
+ CHK(eq_eps(T.E, Tf1, T.SE));
+ CHK(nreals + nfails == N);
+ OK(sdis_estimator_ref_put(estimator));
+
+#if 0
+ OK(sdis_solve_medium(scn, 1, solid1, trange, 1.f, -1, 0,
+ SDIS_HEAT_PATH_ALL, &estimator));
+ dump_heat_paths(stderr, estimator);
+ OK(sdis_estimator_ref_put(estimator));
+#endif
+
+ /* Create a new scene with the same medium in the 2 super shapes */
+ OK(sdis_scene_ref_put(scn));
+ ctx.interf0 = solid0_fluid0;
+ ctx.interf1 = solid0_fluid1;
+ OK(sdis_scene_create(dev, ntris, get_indices, get_interface, nverts,
+ get_position, &ctx, &scn));
+
+ OK(sdis_scene_get_medium_spread(scn, solid0, &v));
+ CHK(eq_eps(v, v0+v1, 1.e-6));
+
+ BA(sdis_solve_medium(scn, N, solid1, trange, 1.f, -1, 0, 0, &estimator));
+ OK(sdis_solve_medium(scn, 10000, solid0, trange, 1.f, -1, 0, 0, &estimator));
+ OK(sdis_estimator_get_temperature(estimator, &T));
+ OK(sdis_estimator_get_realisation_count(estimator, &nreals));
+ OK(sdis_estimator_get_failure_count(estimator, &nfails));
+ ref = Tf0 * v0/v + Tf1 * v1/v;
+ printf("Shape0 + Shape1 temperature = %g ~ %g +/- %g\n", ref, T.E, T.SE);
+ printf("#failures = %lu/10000\n", nfails);
+ CHK(eq_eps(T.E, ref, T.SE*3));
+ OK(sdis_estimator_ref_put(estimator));
+
+ /* Release */
+ OK(s3dut_mesh_ref_put(msh0));
+ OK(s3dut_mesh_ref_put(msh1));
+ OK(sdis_device_ref_put(dev));
+ OK(sdis_medium_ref_put(fluid0));
+ OK(sdis_medium_ref_put(fluid1));
+ OK(sdis_medium_ref_put(solid0));
+ OK(sdis_medium_ref_put(solid1));
+ OK(sdis_interface_ref_put(solid0_fluid0));
+ OK(sdis_interface_ref_put(solid0_fluid1));
+ OK(sdis_interface_ref_put(solid1_fluid1));
+ OK(sdis_scene_ref_put(scn));
+
+ check_memory_allocator(&allocator);
+ mem_shutdown_proxy_allocator(&allocator);
+ CHK(mem_allocated_size() == 0);
+ return 0;
+}
diff --git a/src/test_sdis_solve_medium_2d.c b/src/test_sdis_solve_medium_2d.c
@@ -0,0 +1,397 @@
+/* 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.h"
+#include "test_sdis_utils.h"
+
+#include <rsys/stretchy_array.h>
+#include <rsys/math.h>
+
+#include <string.h>
+
+#define Tf0 300.0
+#define Tf1 330.0
+#define N 1000ul /* #realisations */
+
+/*
+ * The scene is composed of a square and a disk whose temperature is unknown.
+ * The square is surrounded by a fluid whose temperature is Tf0 while the disk
+ * is in a fluid whose temperature is Tf1. The temperature of the square
+ * and the disk are thus uniform and equal to Tf0 and Tf1, respectively.
+ *
+ * # # Tf1 +---------+
+ * # # _\ | | _\ Tf0
+ * # # / / | | / /
+ * # # \__/ | | \__/
+ * # # | |
+ * # # +---------+
+ *
+ * This program performs 2 tests. In the first one, the square and the disk
+ * have different media; the medium solver estimates the temperature of
+ * the square or the one of the disk. In the second test, the scene is updated
+ * to use the same medium for the 2 shapes. When invoked on
+ * the right medium, the estimated temperature T is equal to :
+ *
+ * T = Tf0 * A0/(A0 + A1) + Tf1 * A1/(A0 + A1)
+ *
+ * with A0 and A1 the area of the shape and the area of the disk, respectively.
+ */
+
+/*******************************************************************************
+ * Geometry
+ ******************************************************************************/
+struct context {
+ const double* positions;
+ const size_t* indices;
+ size_t nsegments_interf0; /* #segments of the interface 0 */
+ struct sdis_interface* interf0;
+ struct sdis_interface* interf1;
+};
+
+static void
+get_indices(const size_t iseg, size_t ids[2], void* context)
+{
+ const struct context* ctx = context;
+ ids[0] = ctx->indices[iseg*2+0];
+ ids[1] = ctx->indices[iseg*2+1];
+}
+
+static void
+get_position(const size_t ivert, double pos[2], void* context)
+{
+ const struct context* ctx = context;
+ pos[0] = ctx->positions[ivert*2+0];
+ pos[1] = ctx->positions[ivert*2+1];
+}
+
+static void
+get_interface(const size_t iseg, struct sdis_interface** bound, void* context)
+{
+ const struct context* ctx = context;
+ *bound = iseg < ctx->nsegments_interf0 ? ctx->interf0 : ctx->interf1;
+}
+
+/*******************************************************************************
+ * Fluid medium
+ ******************************************************************************/
+struct fluid {
+ double temperature;
+};
+
+static double
+fluid_get_temperature
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ CHK(data != NULL && vtx != NULL);
+ return ((const struct fluid*)sdis_data_cget(data))->temperature;
+}
+
+/*******************************************************************************
+ * Solid medium
+ ******************************************************************************/
+struct solid {
+ double cp;
+ double lambda;
+ double rho;
+ double delta;
+ double temperature;
+};
+
+static double
+solid_get_calorific_capacity
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ CHK(data != NULL && vtx != NULL);
+ return ((const struct solid*)sdis_data_cget(data))->cp;
+}
+
+static double
+solid_get_thermal_conductivity
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ CHK(data != NULL && vtx != NULL);
+ return ((const struct solid*)sdis_data_cget(data))->lambda;
+}
+
+static double
+solid_get_volumic_mass
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ CHK(data != NULL && vtx != NULL);
+ return ((const struct solid*)sdis_data_cget(data))->rho;
+}
+
+static double
+solid_get_delta
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ CHK(data != NULL && vtx != NULL);
+ return ((const struct solid*)sdis_data_cget(data))->delta;
+}
+
+static double
+solid_get_temperature
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ CHK(data != NULL && vtx != NULL);
+ return ((const struct solid*)sdis_data_cget(data))->temperature;
+}
+
+struct interf {
+ double hc;
+ double epsilon;
+ double specular_fraction;
+};
+
+static double
+interface_get_convection_coef
+ (const struct sdis_interface_fragment* frag, struct sdis_data* data)
+{
+ CHK(data != NULL && frag != NULL);
+ return ((const struct interf*)sdis_data_cget(data))->hc;
+}
+
+static double
+interface_get_emissivity
+ (const struct sdis_interface_fragment* frag, struct sdis_data* data)
+{
+ CHK(data != NULL && frag != NULL);
+ return ((const struct interf*)sdis_data_cget(data))->epsilon;
+}
+
+static double
+interface_get_specular_fraction
+ (const struct sdis_interface_fragment* frag, struct sdis_data* data)
+{
+ CHK(data != NULL && frag != NULL);
+ return ((const struct interf*)sdis_data_cget(data))->specular_fraction;
+}
+
+/*******************************************************************************
+ * Test
+ ******************************************************************************/
+int
+main(int argc, char** argv)
+{
+ struct mem_allocator allocator;
+ struct sdis_mc T = SDIS_MC_NULL;
+ struct sdis_device* dev = NULL;
+ struct sdis_medium* solid0 = NULL;
+ struct sdis_medium* solid1 = NULL;
+ struct sdis_medium* fluid0 = NULL;
+ struct sdis_medium* fluid1 = NULL;
+ struct sdis_interface* solid0_fluid0 = NULL;
+ struct sdis_interface* solid0_fluid1 = NULL;
+ struct sdis_interface* solid1_fluid1 = NULL;
+ struct sdis_scene* scn = NULL;
+ struct sdis_data* data = NULL;
+ struct sdis_estimator* estimator = NULL;
+ struct fluid* fluid_param = NULL;
+ struct solid* solid_param = NULL;
+ struct interf* interface_param = NULL;
+ struct sdis_fluid_shader fluid_shader = DUMMY_FLUID_SHADER;
+ struct sdis_solid_shader solid_shader = DUMMY_SOLID_SHADER;
+ struct sdis_interface_shader interface_shader = SDIS_INTERFACE_SHADER_NULL;
+ struct context ctx;
+ const double trange[2] = {0, INF};
+ double a, a0, a1;
+ double ref;
+ double* positions = NULL;
+ size_t* indices = NULL;
+ size_t nverts;
+ size_t nreals;
+ size_t nfails;
+ size_t i;
+ (void)argc, (void)argv;
+
+ OK(mem_init_proxy_allocator(&allocator, &mem_default_allocator));
+ OK(sdis_device_create(NULL, &allocator, SDIS_NTHREADS_DEFAULT, 1, &dev));
+
+ fluid_shader.temperature = fluid_get_temperature;
+
+ /* Create the fluid0 medium */
+ OK(sdis_data_create
+ (dev, sizeof(struct fluid), ALIGNOF(struct fluid), NULL, &data));
+ fluid_param = sdis_data_get(data);
+ fluid_param->temperature = Tf0;
+ OK(sdis_fluid_create(dev, &fluid_shader, data, &fluid0));
+ OK(sdis_data_ref_put(data));
+
+ /* Create the fluid1 medium */
+ OK(sdis_data_create
+ (dev, sizeof(struct fluid), ALIGNOF(struct fluid), NULL, &data));
+ fluid_param = sdis_data_get(data);
+ fluid_param->temperature = Tf1;
+ OK(sdis_fluid_create(dev, &fluid_shader, data, &fluid1));
+ OK(sdis_data_ref_put(data));
+
+ /* Setup the solid shader */
+ solid_shader.calorific_capacity = solid_get_calorific_capacity;
+ solid_shader.thermal_conductivity = solid_get_thermal_conductivity;
+ solid_shader.volumic_mass = solid_get_volumic_mass;
+ solid_shader.delta_solid = solid_get_delta;
+ solid_shader.temperature = solid_get_temperature;
+
+ /* Create the solid0 medium */
+ OK(sdis_data_create
+ (dev, sizeof(struct solid), ALIGNOF(struct solid), NULL, &data));
+ solid_param = sdis_data_get(data);
+ solid_param->cp = 1.0;
+ solid_param->lambda = 0.1;
+ solid_param->rho = 1.0;
+ solid_param->delta = 1.0/20.0;
+ solid_param->temperature = -1; /* Unknown temperature */
+ OK(sdis_solid_create(dev, &solid_shader, data, &solid0));
+ OK(sdis_data_ref_put(data));
+
+ /* Create the solid1 medium */
+ OK(sdis_data_create
+ (dev, sizeof(struct solid), ALIGNOF(struct solid), NULL, &data));
+ solid_param = sdis_data_get(data);
+ solid_param->cp = 1.0;
+ solid_param->lambda = 1.0;
+ solid_param->rho = 1.0;
+ solid_param->delta = 1.0/20.0;
+ solid_param->temperature = -1; /* Unknown temperature */
+ OK(sdis_solid_create(dev, &solid_shader, data, &solid1));
+ OK(sdis_data_ref_put(data));
+
+ /* Create the interfaces */
+ OK(sdis_data_create(dev, sizeof(struct interf),
+ ALIGNOF(struct interf), NULL, &data));
+ interface_param = sdis_data_get(data);
+ interface_param->hc = 0.5;
+ interface_param->epsilon = 0;
+ interface_param->specular_fraction = 0;
+ interface_shader.convection_coef = interface_get_convection_coef;
+ interface_shader.front = SDIS_INTERFACE_SIDE_SHADER_NULL;
+ interface_shader.back.temperature = NULL;
+ interface_shader.back.emissivity = interface_get_emissivity;
+ interface_shader.back.specular_fraction = interface_get_specular_fraction;
+ OK(sdis_interface_create
+ (dev, solid0, fluid0, &interface_shader, data, &solid0_fluid0));
+ OK(sdis_interface_create
+ (dev, solid0, fluid1, &interface_shader, data, &solid0_fluid1));
+ OK(sdis_interface_create
+ (dev, solid1, fluid1, &interface_shader, data, &solid1_fluid1));
+ OK(sdis_data_ref_put(data));
+
+ /* Setup the square geometry */
+ sa_add(positions, square_nvertices*2);
+ sa_add(indices, square_nsegments*2);
+ memcpy(positions, square_vertices, square_nvertices*sizeof(double[2]));
+ memcpy(indices, square_indices, square_nsegments*sizeof(size_t[2]));
+
+ /* Transate the square in X */
+ FOR_EACH(i, 0, square_nvertices) positions[i*2] += 2;
+
+ /* Setup a disk */
+ nverts = 64;
+ FOR_EACH(i, 0, nverts) {
+ const double theta = (double)i * (2*PI)/(double)nverts;
+ const double r = 1; /* Radius */
+ const double x = cos(theta) * r - 2/* X translation */;
+ const double y = sin(theta) * r + 0.5/* Y translation */;
+ sa_push(positions, x);
+ sa_push(positions, y);
+ }
+ FOR_EACH(i, 0, nverts) {
+ const size_t i0 = i + square_nvertices;
+ const size_t i1 = (i+1) % nverts + square_nvertices;
+ /* Flip the ids to ensure that the normals point inward the disk */
+ sa_push(indices, i1);
+ sa_push(indices, i0);
+ }
+
+ /* Create the scene */
+ ctx.positions = positions;
+ ctx.indices = indices;
+ ctx.nsegments_interf0 = square_nsegments;
+ ctx.interf0 = solid0_fluid0;
+ ctx.interf1 = solid1_fluid1;
+ OK(sdis_scene_2d_create(dev, sa_size(indices)/2, get_indices, get_interface,
+ sa_size(positions)/2, get_position, &ctx, &scn));
+
+ OK(sdis_scene_get_medium_spread(scn, solid0, &a0));
+ CHK(eq_eps(a0, 1.0, 1.e-6));
+ OK(sdis_scene_get_medium_spread(scn, solid1, &a1));
+ /* Rough estimation since the disk is coarsely discretized */
+ CHK(eq_eps(a1, PI, 1.e-1));
+
+ /* Estimate the temperature of the square */
+ OK(sdis_solve_medium(scn, N, solid0, trange, 1.f, -1, 0, 0, &estimator));
+ OK(sdis_estimator_get_temperature(estimator, &T));
+ OK(sdis_estimator_get_realisation_count(estimator, &nreals));
+ OK(sdis_estimator_get_failure_count(estimator, &nfails));
+ printf("Square temperature = "STR(Tf0)" ~ %g +/- %g\n", T.E, T.SE);
+ printf("#failures = %lu / %lu\n", nfails, N);
+ CHK(eq_eps(T.E, Tf0, T.SE));
+ CHK(nreals + nfails == N);
+ OK(sdis_estimator_ref_put(estimator));
+
+ /* Estimate the temperature of the disk */
+ OK(sdis_solve_medium(scn, N, solid1, trange, 1.f, -1, 0, 0, &estimator));
+ OK(sdis_estimator_get_temperature(estimator, &T));
+ OK(sdis_estimator_get_realisation_count(estimator, &nreals));
+ OK(sdis_estimator_get_failure_count(estimator, &nfails));
+ printf("Disk temperature = "STR(Tf1)" ~ %g +/- %g\n", T.E, T.SE);
+ printf("#failures = %lu / %lu\n", nfails, N);
+ CHK(eq_eps(T.E, Tf1, T.SE));
+ CHK(nreals + nfails == N);
+ OK(sdis_estimator_ref_put(estimator));
+
+ /* Create a new scene with the same medium for the disk and the square */
+ OK(sdis_scene_ref_put(scn));
+ ctx.interf0 = solid0_fluid0;
+ ctx.interf1 = solid0_fluid1;
+ OK(sdis_scene_2d_create(dev, sa_size(indices)/2, get_indices, get_interface,
+ sa_size(positions)/2, get_position, &ctx, &scn));
+
+ OK(sdis_scene_get_medium_spread(scn, solid0, &a));
+ CHK(eq_eps(a, a0+a1, 1.e-6));
+
+ /* Estimate the temperature of the square and disk shapes */
+ BA(sdis_solve_medium(scn, N, solid1, trange, 1.f, -1, 0, 0, &estimator));
+ OK(sdis_solve_medium(scn, 10000, solid0, trange, 1.f, -1, 0, 0, &estimator));
+ OK(sdis_estimator_get_temperature(estimator, &T));
+ OK(sdis_estimator_get_realisation_count(estimator, &nreals));
+ OK(sdis_estimator_get_failure_count(estimator, &nfails));
+ ref = Tf0 * a0/a + Tf1 * a1/a;
+ printf("Square + Disk temperature = %g ~ %g +/- %g\n", ref, T.E, T.SE);
+ printf("#failures = %lu / 10000\n", nfails);
+ CHK(eq_eps(T.E, ref, 3*T.SE));
+ CHK(nreals + nfails == 10000);
+ OK(sdis_estimator_ref_put(estimator));
+
+ /* Release */
+ OK(sdis_device_ref_put(dev));
+ OK(sdis_medium_ref_put(solid0));
+ OK(sdis_medium_ref_put(solid1));
+ OK(sdis_medium_ref_put(fluid0));
+ OK(sdis_medium_ref_put(fluid1));
+ OK(sdis_interface_ref_put(solid0_fluid0));
+ OK(sdis_interface_ref_put(solid0_fluid1));
+ OK(sdis_interface_ref_put(solid1_fluid1));
+ OK(sdis_scene_ref_put(scn));
+
+ sa_release(positions);
+ sa_release(indices);
+
+ check_memory_allocator(&allocator);
+ mem_shutdown_proxy_allocator(&allocator);
+ CHK(mem_allocated_size() == 0);
+ return 0;
+}
+
