commit b9fed004322bf3939d09eac4231afc6ebefdd8b5
parent 6f004ab27b1a3ae6f469fdadff62337f76579136
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Tue, 14 May 2019 12:26:26 +0200
Merge branch 'feature_reduce_failures' into develop
Diffstat:
15 files changed, 1921 insertions(+), 649 deletions(-)
diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt
@@ -37,7 +37,7 @@ find_package(Star3D 0.6 REQUIRED)
find_package(StarSP 0.8 REQUIRED)
find_package(StarEnc 0.2.2 REQUIRED)
find_package(StarEnc2D 0.2.2 REQUIRED)
-find_package(RSys 0.6.1 REQUIRED)
+find_package(RSys 0.8 REQUIRED)
find_package(OpenMP 2.0 REQUIRED)
set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} ${RCMAKE_SOURCE_DIR})
@@ -181,6 +181,7 @@ if(NOT NO_TEST)
new_test(test_sdis_interface)
new_test(test_sdis_medium)
new_test(test_sdis_scene)
+ new_test(test_sdis_solid_random_walk_robustness)
new_test(test_sdis_solve_camera)
new_test(test_sdis_solve_probe)
new_test(test_sdis_solve_probe2)
@@ -193,7 +194,7 @@ if(NOT NO_TEST)
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)
+ new_test(test_sdis_volumic_power4)
# Additionnal tests
build_test(test_sdis_volumic_power2)
@@ -206,6 +207,7 @@ if(NOT NO_TEST)
add_test(test_sdis_volumic_power3_2d test_sdis_volumic_power3_2d)
endif()
+ target_link_libraries(test_sdis_solid_random_walk_robustness Star3DUT)
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})
diff --git a/src/sdis_green.c b/src/sdis_green.c
@@ -781,7 +781,7 @@ green_function_redux_and_clear
{
size_t i;
res_T res = RES_OK;
- ASSERT(dst && greens && ngreens);
+ ASSERT(dst && greens);
FOR_EACH(i, 0, ngreens) {
res = green_function_merge_and_clear(dst, greens[i]);
diff --git a/src/sdis_heat_path_boundary_Xd.h b/src/sdis_heat_path_boundary_Xd.h
@@ -25,15 +25,15 @@
#include "sdis_Xd_begin.h"
/* Emperical scale factor applied to the challenged reinjection distance. If
- * the distance to reinject is less than this adjusted value, the solver
- * switches from 2D reinjection scheme to the 1D reinjection scheme in order to
- * avoid numerical issues. */
+ * the distance to reinject is less than this adjusted value, the solver will
+ * try to discard the reinjection distance if possible in order to avoid
+ * numerical issues. */
#define REINJECT_DST_MIN_SCALE 0.125f
/*******************************************************************************
* Helper functions
******************************************************************************/
-static FINLINE void
+static void
XD(sample_reinjection_dir)
(const struct XD(rwalk)* rwalk, struct ssp_rng* rng, float dir[DIM])
{
@@ -75,6 +75,355 @@ XD(sample_reinjection_dir)
#endif
}
+#if DIM == 2
+static void
+XD(move_away_primitive_boundaries)
+ (struct XD(rwalk)* rwalk,
+ const double delta)
+{
+ struct sXd(attrib) attr;
+ float pos[DIM];
+ float dir[DIM];
+ float len;
+ const float st = 0.5f;
+ ASSERT(rwalk && !SXD_HIT_NONE(&rwalk->hit) && delta > 0);
+
+ SXD(primitive_get_attrib(&rwalk->hit.prim, SXD_POSITION, st, &attr));
+
+ fX_set_dX(pos, rwalk->vtx.P);
+ fX(sub)(dir, attr.value, pos);
+ len = fX(normalize)(dir, dir);
+ len = MMIN(len, (float)(delta*0.1));
+
+ XD(move_pos)(rwalk->vtx.P, dir, len);
+}
+#else
+/* Move the random walk away from the primitive boundaries to avoid numerical
+ * issues leading to inconsistent random walks. */
+static void
+XD(move_away_primitive_boundaries)
+ (struct XD(rwalk)* rwalk,
+ const double delta)
+{
+ struct s3d_attrib v0, v1, v2; /* Triangle vertices */
+ float E[3][4]; /* 3D edge equations */
+ float dst[3]; /* Distance from current position to edge equation */
+ float N[3]; /* Triangle normal */
+ float P[3]; /* Random walk position */
+ float tmp[3];
+ float min_dst, max_dst;
+ float cos_a1, cos_a2;
+ float len;
+ int imax;
+ int imin;
+ int imid;
+ int i;
+ ASSERT(rwalk && delta > 0 && !S3D_HIT_NONE(&rwalk->hit));
+
+ fX_set_dX(P, rwalk->vtx.P);
+
+ /* Fetch triangle vertices */
+ S3D(triangle_get_vertex_attrib(&rwalk->hit.prim, 0, S3D_POSITION, &v0));
+ S3D(triangle_get_vertex_attrib(&rwalk->hit.prim, 1, S3D_POSITION, &v1));
+ S3D(triangle_get_vertex_attrib(&rwalk->hit.prim, 2, S3D_POSITION, &v2));
+
+ /* Compute the edge vector */
+ f3_sub(E[0], v1.value, v0.value);
+ f3_sub(E[1], v2.value, v1.value);
+ f3_sub(E[2], v0.value, v2.value);
+
+ /* Compute the triangle normal */
+ f3_cross(N, E[1], E[0]);
+
+ /* Compute the 3D edge equation */
+ f3_normalize(E[0], f3_cross(E[0], E[0], N));
+ f3_normalize(E[1], f3_cross(E[1], E[1], N));
+ f3_normalize(E[2], f3_cross(E[2], E[2], N));
+ E[0][3] = -f3_dot(E[0], v0.value);
+ E[1][3] = -f3_dot(E[1], v1.value);
+ E[2][3] = -f3_dot(E[2], v2.value);
+
+ /* Compute the distance from current position to the edges */
+ dst[0] = f3_dot(E[0], P) + E[0][3];
+ dst[1] = f3_dot(E[1], P) + E[1][3];
+ dst[2] = f3_dot(E[2], P) + E[2][3];
+
+ /* Retrieve the min and max distance from random walk position to triangle
+ * edges */
+ min_dst = MMIN(MMIN(dst[0], dst[1]), dst[2]);
+ max_dst = MMAX(MMAX(dst[0], dst[1]), dst[2]);
+
+ /* Sort the edges with respect to their distance to the random walk position */
+ FOR_EACH(i, 0, 3) {
+ if(dst[i] == min_dst) {
+ imin = i;
+ } else if(dst[i] == max_dst) {
+ imax = i;
+ } else {
+ imid = i;
+ }
+ }
+ (void)imax;
+
+ /* TODO if the current position is near a vertex, one should move toward the
+ * farthest edge along its normal to avoid too small displacement */
+
+ /* Compute the distance `dst' from the current position to the edges to move
+ * to, along the normal of the edge from which the random walk is the nearest
+ *
+ * +. cos(a) = d / dst => dst = d / cos_a
+ * / `*.
+ * / | `*.
+ * / dst| a /`*.
+ * / | / `*.
+ * / | / d `*.
+ * / |/ `*.
+ * +---------o----------------+ */
+ cos_a1 = f3_dot(E[imin], f3_minus(tmp, E[imid]));
+ cos_a2 = f3_dot(E[imin], f3_minus(tmp, E[imax]));
+ dst[imid] = cos_a1 > 0 ? dst[imid] / cos_a1 : FLT_MAX;
+ dst[imax] = cos_a2 > 0 ? dst[imax] / cos_a2 : FLT_MAX;
+
+ /* Compute the maximum displacement distance into the triangle along the
+ * normal of the edge from which the random walk is the nearest */
+ len = MMIN(dst[imid], dst[imax]);
+ ASSERT(len != FLT_MAX);
+
+ /* Define the displacement distance as the minimum between 10 percent of
+ * delta and len / 2. */
+ len = MMIN(len*0.5f, (float)(delta*0.1));
+ XD(move_pos)(rwalk->vtx.P, E[imin], len);
+}
+#endif
+
+static res_T
+XD(select_reinjection_dir)
+ (const struct sdis_scene* scn,
+ const struct sdis_medium* mdm, /* Medium into which the reinjection occurs */
+ struct XD(rwalk)* rwalk, /* Current random walk state */
+ const float dir0[DIM], /* Challenged direction */
+ const float dir1[DIM], /* Challanged direction */
+ const double delta, /* Max reinjection distance */
+ float reinject_dir[DIM], /* Selected direction */
+ float* reinject_dst, /* Effective reinjection distance */
+ int can_move, /* Define of the random wal pos can be moved or not */
+ int* move_pos, /* Define if the current random walk was moved. May be NULL */
+ struct sXd(hit)* reinject_hit) /* Hit along the reinjection dir */
+{
+ struct sdis_interface* interf;
+ struct sdis_medium* mdm0;
+ struct sdis_medium* mdm1;
+ struct hit_filter_data filter_data;
+ struct sXd(hit) hit;
+ struct sXd(hit) hit0;
+ struct sXd(hit) hit1;
+ double tmp[DIM];
+ double dst;
+ double dst0;
+ double dst1;
+ const double delta_adjusted = delta * RAY_RANGE_MAX_SCALE;
+ const float* dir;
+ const float reinject_threshold = (float)delta * REINJECT_DST_MIN_SCALE;
+ float org[DIM];
+ float range[2];
+ enum sdis_side side;
+ int iattempt = 0;
+ const int MAX_ATTEMPTS = can_move ? 2 : 1;
+ res_T res = RES_OK;
+ ASSERT(scn && mdm && rwalk && dir0 && dir1 && delta > 0);
+ ASSERT(reinject_dir && reinject_dst && reinject_hit);
+
+ if(move_pos) *move_pos = 0;
+
+ do {
+ f2(range, 0, FLT_MAX);
+ fX_set_dX(org, rwalk->vtx.P);
+ filter_data.XD(hit) = rwalk->hit;
+ filter_data.epsilon = delta * 0.01;
+ SXD(scene_view_trace_ray(scn->sXd(view), org, dir0, range, &filter_data, &hit0));
+ SXD(scene_view_trace_ray(scn->sXd(view), org, dir1, range, &filter_data, &hit1));
+
+ /* Retrieve the medium at the reinjection pos along dir0 */
+ if(SXD_HIT_NONE(&hit0)) {
+ XD(move_pos)(dX(set)(tmp, rwalk->vtx.P), dir0, (float)delta);
+ res = scene_get_medium_in_closed_boundaries(scn, tmp, &mdm0);
+ if(res != RES_OK) goto error;
+ } else {
+ interf = scene_get_interface(scn, hit0.prim.prim_id);
+ side = fX(dot)(dir0, hit0.normal) < 0 ? SDIS_FRONT : SDIS_BACK;
+ mdm0 = interface_get_medium(interf, side);
+ }
+
+ /* Retrieve the medium at the reinjection pos along dir1 */
+ if(SXD_HIT_NONE(&hit1)) {
+ XD(move_pos)(dX(set)(tmp, rwalk->vtx.P), dir1, (float)delta);
+ res = scene_get_medium_in_closed_boundaries(scn, tmp, &mdm1);
+ if(res != RES_OK) goto error;
+ } else {
+ interf = scene_get_interface(scn, hit1.prim.prim_id);
+ side = fX(dot)(dir1, hit1.normal) < 0 ? SDIS_FRONT : SDIS_BACK;
+ mdm1 = interface_get_medium(interf, side);
+ }
+
+ dst0 = dst1 = -1;
+ if(mdm0 == mdm) { /* Check reinjection consistency */
+ if(hit0.distance <= delta_adjusted) {
+ dst0 = hit0.distance;
+ } else {
+ dst0 = delta;
+ hit0 = SXD_HIT_NULL;
+ }
+ }
+ if(mdm1 == mdm) {/* Check reinjection consistency */
+ if(hit1.distance <= delta_adjusted) {
+ dst1 = hit1.distance;
+ } else {
+ dst1 = delta;
+ hit1 = SXD_HIT_NULL;
+ }
+ }
+
+ /* No valid reinjection. Maybe the random walk is near a sharp corner and
+ * thus the ray-tracing misses the enclosure geometry. Another possibility
+ * is that the random walk lies roughly on an edge. In this case, sampled
+ * reinjecton dirs can intersect the primitive on the other side of the
+ * edge. Normally, this primitive should be filtered by the "hit_filter"
+ * function but this may be not the case due to a "threshold effect". In
+ * both situations, try to slightly move away from the primitive boundaries
+ * and retry to find a valid reinjection. */
+ if(dst0 == -1 && dst1 == -1) {
+ XD(move_away_primitive_boundaries)(rwalk, delta);
+ if(move_pos) *move_pos = 1;
+ }
+ } while(dst0 == -1 && dst1 == -1 && ++iattempt < MAX_ATTEMPTS);
+
+ if(dst0 == -1 && dst1 == -1) { /* No valid reinjection */
+ log_err(scn->dev, "%s: no valid reinjection direction at {%g, %g, %g}.\n",
+ FUNC_NAME, SPLIT3(rwalk->vtx.P));
+ res = RES_BAD_OP_IRRECOVERABLE;
+ goto error;
+ }
+
+ if(dst0 == -1) {
+ /* Invalid dir0 -> move along dir1 */
+ dir = dir1;
+ dst = dst1;
+ hit = hit1;
+ } else if(dst1 == -1) {
+ /* Invalid dir1 -> move along dir0 */
+ dir = dir0;
+ dst = dst0;
+ hit = hit0;
+ } else if(dst0 < reinject_threshold && dst1 < reinject_threshold) {
+ /* The displacement along dir0 and dir1 are both below the reinjection
+ * threshold that defines a distance under which the temperature gradients
+ * are ignored. Move along the direction that allows the maximum
+ * displacement. */
+ if(dst0 > dst1) {
+ dir = dir0;
+ dst = dst0;
+ hit = hit0;
+ } else {
+ dir = dir1;
+ dst = dst1;
+ hit = hit1;
+ }
+ } else if(dst0 < reinject_threshold) {
+ /* Ingore dir0 that is bellow the reinject threshold */
+ dir = dir1;
+ dst = dst1;
+ hit = hit1;
+ } else if(dst1 < reinject_threshold) {
+ /* Ingore dir1 that is bellow the reinject threshold */
+ dir = dir0;
+ dst = dst0;
+ hit = hit0;
+ } else {
+ /* All reinjection directions are valid. Choose the first 1 that was
+ * randomly selected by the sample_reinjection_dir procedure and adjust
+ * the displacement distance. */
+ dir = dir0;
+ if(dst0 < dst1) {
+ dst = dst0;
+ hit = hit0;
+ } else {
+ dst = dst1;
+ hit = SXD_HIT_NULL;
+ }
+
+ /* If the displacement distance is too close of a boundary, move to the
+ * boundary in order to avoid numerical uncertainty. */
+ if(!SXD_HIT_NONE(&hit0)
+ && dst0 != dst
+ && eq_eps(dst0, dst, dst0*0.1)) {
+ dst = dst0;
+ hit = hit0;
+ }
+ }
+
+ /* Setup output variable */
+ fX(set)(reinject_dir, dir);
+ *reinject_dst = (float)dst;
+ *reinject_hit = hit;
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
+static res_T
+XD(select_reinjection_dir_and_check_validity)
+ (const struct sdis_scene* scn,
+ const struct sdis_medium* mdm, /* Medium into which the reinjection occurs */
+ struct XD(rwalk)* rwalk, /* Current random walk state */
+ const float dir0[DIM], /* Challenged direction */
+ const float dir1[DIM], /* Challanged direction */
+ const double delta, /* Max reinjection distance */
+ float out_reinject_dir[DIM], /* Selected direction */
+ float* out_reinject_dst, /* Effective reinjection distance */
+ int can_move, /* Define of the random wal pos can be moved or not */
+ int* move_pos, /* Define if the current random walk was moved. May be NULL */
+ int* is_valid, /* Define if the reinjection defines a valid pos */
+ struct sXd(hit)* out_reinject_hit) /* Hit along the reinjection dir */
+{
+ double pos[DIM];
+ struct sdis_medium* reinject_mdm;
+ struct sXd(hit) reinject_hit;
+ float reinject_dir[DIM];
+ float reinject_dst;
+ res_T res = RES_OK;
+ ASSERT(is_valid && out_reinject_dir && out_reinject_dst && out_reinject_hit);
+
+ /* Select a reinjection direction */
+ res = XD(select_reinjection_dir)(scn, mdm, rwalk, dir0, dir1, delta,
+ reinject_dir, &reinject_dst, can_move, move_pos, &reinject_hit);
+ if(res != RES_OK) goto error;
+
+ if(!SXD_HIT_NONE(&reinject_hit)) {
+ *is_valid = 1;
+ } else {
+ /* Check medium consistency at the reinjection position */
+ XD(move_pos)(dX(set)(pos, rwalk->vtx.P), reinject_dir, reinject_dst);
+ res = scene_get_medium_in_closed_boundaries
+ (scn, pos, &reinject_mdm);
+ if(res != RES_OK) goto error;
+
+ *is_valid = reinject_mdm == mdm;
+ }
+
+ if(*is_valid) {
+ fX(set)(out_reinject_dir, reinject_dir);
+ *out_reinject_dst = reinject_dst;
+ *out_reinject_hit = reinject_hit;
+ }
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
/* Check that the interface fragment is consistent with the current state of
* the random walk */
static INLINE int
@@ -111,8 +460,9 @@ XD(solid_solid_boundary_path)
struct ssp_rng* rng,
struct XD(temperature)* T)
{
- struct sXd(hit) hit0, hit1, hit2, hit3;
+ struct sXd(hit) hit0, hit1;
struct sXd(hit)* hit;
+ struct XD(rwalk) rwalk_saved;
struct sdis_interface* interf = NULL;
struct sdis_medium* solid_front = NULL;
struct sdis_medium* solid_back = NULL;
@@ -120,18 +470,21 @@ XD(solid_solid_boundary_path)
double lambda_front, lambda_back;
double delta_front, delta_back;
double delta_boundary_front, delta_boundary_back;
- double delta_boundary;
- double reinject_dst_front, reinject_dst_back;
- double reinject_dst;
double proba;
double tmp;
double r;
double power;
- float range0[2], range1[2];
float dir0[DIM], dir1[DIM], dir2[DIM], dir3[DIM];
+ float dir_front[DIM], dir_back[DIM];
float* dir;
- float pos[DIM];
- int dim = DIM;
+ float reinject_dst_front, reinject_dst_back;
+ float reinject_dst;
+ /* In 2D it is useless to try to resample a reinjection direction since there
+ * is only one possible direction */
+ const int MAX_ATTEMPTS = DIM == 2 ? 1 : 10;
+ int iattempt;
+ int move;
+ int reinjection_is_valid;
res_T res = RES_OK;
ASSERT(scn && fp_to_meter > 0 && ctx && frag && rwalk && rng && T);
ASSERT(XD(check_rwalk_fragment_consistency)(rwalk, frag));
@@ -151,30 +504,57 @@ XD(solid_solid_boundary_path)
/* Note that reinjection distance is *FIXED*. It MUST ensure that the orthogonal
* distance from the boundary to the point to challenge is equal to delta. */
delta_front = solid_get_delta(solid_front, &rwalk->vtx);
- delta_back = solid_get_delta(solid_back, &rwalk->vtx);
+ delta_back = solid_get_delta(solid_back, &rwalk->vtx);
delta_boundary_front = delta_front*sqrt(DIM);
- delta_boundary_back = delta_back *sqrt(DIM);
-
- /* Sample a reinjection direction and reflect it around the normal. Then
- * reflect them on the back side of the interface. */
- XD(sample_reinjection_dir)(rwalk, rng, dir0);
- XD(reflect)(dir2, dir0, rwalk->hit.normal);
- fX(minus)(dir1, dir0);
- fX(minus)(dir3, dir2);
-
- /* Trace the sampled directions on both sides of the interface to adjust the
- * reinjection distance of the random walk . */
- fX_set_dX(pos, rwalk->vtx.P);
- f2(range0, 0, (float)delta_boundary_front*RAY_RANGE_MAX_SCALE);
- f2(range1, 0, (float)delta_boundary_back *RAY_RANGE_MAX_SCALE);
- SXD(scene_view_trace_ray(scn->sXd(view), pos, dir0, range0, &rwalk->hit, &hit0));
- SXD(scene_view_trace_ray(scn->sXd(view), pos, dir1, range1, &rwalk->hit, &hit1));
- SXD(scene_view_trace_ray(scn->sXd(view), pos, dir2, range0, &rwalk->hit, &hit2));
- SXD(scene_view_trace_ray(scn->sXd(view), pos, dir3, range1, &rwalk->hit, &hit3));
+ delta_boundary_back = delta_back *sqrt(DIM);
+
+ rwalk_saved = *rwalk;
+ reinjection_is_valid = 0;
+ iattempt = 0;
+ do {
+ if(iattempt != 0) *rwalk = rwalk_saved;
+
+ /* Sample a reinjection direction and reflect it around the normal. Then
+ * reflect them on the back side of the interface. */
+ XD(sample_reinjection_dir)(rwalk, rng, dir0);
+ XD(reflect)(dir2, dir0, rwalk->hit.normal);
+ fX(minus)(dir1, dir0);
+ fX(minus)(dir3, dir2);
+
+ /* Select the reinjection direction and distance for the front side */
+ res = XD(select_reinjection_dir_and_check_validity)(scn, solid_front, rwalk,
+ dir0, dir2, delta_boundary_front, dir_front, &reinject_dst_front, 1, &move,
+ &reinjection_is_valid, &hit0);
+ if(res != RES_OK) goto error;
+ if(!reinjection_is_valid) continue;
- /* Adjust the reinjection distance */
- reinject_dst_front = MMIN(MMIN(delta_boundary_front, hit0.distance), hit2.distance);
- reinject_dst_back = MMIN(MMIN(delta_boundary_back, hit1.distance), hit3.distance);
+ /* Select the reinjection direction and distance for the back side */
+ res = XD(select_reinjection_dir_and_check_validity)(scn, solid_back, rwalk,
+ dir1, dir3, delta_boundary_back, dir_back, &reinject_dst_back, 1, &move,
+ &reinjection_is_valid, &hit1);
+ if(res != RES_OK) goto error;
+ if(!reinjection_is_valid) continue;
+
+ /* If random walk was moved by the select_reinjection_dir on back side, one
+ * has to rerun the select_reinjection_dir on front side at the new pos */
+ if(move) {
+ res = XD(select_reinjection_dir_and_check_validity)(scn, solid_front,
+ rwalk, dir0, dir2, delta_boundary_front, dir_front, &reinject_dst_front,
+ 0, NULL, &reinjection_is_valid, &hit0);
+ if(res != RES_OK) goto error;
+ if(!reinjection_is_valid) continue;
+ }
+ } while(!reinjection_is_valid && ++iattempt < MAX_ATTEMPTS);
+
+ /* Could not find a valid reinjection */
+ if(iattempt >= MAX_ATTEMPTS) {
+ *rwalk = rwalk_saved;
+ log_err(scn->dev,
+ "%s: could not find a valid soid/solid reinjection at {%g, %g, %g}.\n",
+ FUNC_NAME, SPLIT3(rwalk->vtx.P));
+ res = RES_BAD_OP_IRRECOVERABLE;
+ goto error;
+ }
/* Define the reinjection side. Note that the proba should be :
* Lf/Df' / (Lf/Df' + Lb/Db')
@@ -190,38 +570,15 @@ XD(solid_solid_boundary_path)
proba = (lambda_front/reinject_dst_front)
/ (lambda_front/reinject_dst_front + lambda_back/reinject_dst_back);
if(r < proba) { /* Reinject in front */
- dir = dir0;
+ dir = dir_front;
hit = &hit0;
mdm = solid_front;
reinject_dst = reinject_dst_front;
- delta_boundary = delta_boundary_front;
} else { /* Reinject in back */
- dir = dir1;
+ dir = dir_back;
hit = &hit1;
mdm = solid_back;
reinject_dst = reinject_dst_back;
- delta_boundary = delta_boundary_back;
- }
-
- /* Switch in 1D reinjection scheme */
- if(reinject_dst < delta_boundary * REINJECT_DST_MIN_SCALE) {
- if(dir == dir0) {
- fX(set)(dir, rwalk->hit.normal);
- } else {
- fX(minus)(dir, rwalk->hit.normal);
- }
-
- f2(range0, 0, (float)delta_boundary*RAY_RANGE_MAX_SCALE);
- SXD(scene_view_trace_ray(scn->sXd(view), pos, dir, range0, &rwalk->hit, hit));
- reinject_dst = MMIN(delta_boundary, hit->distance),
- dim = 1;
-
- /* Hit something in 1D. Arbitrarily move the random walk to 0.5 of the hit
- * distance */
- if(!SXD_HIT_NONE(hit)) {
- reinject_dst *= 0.5;
- *hit = SXD_HIT_NULL;
- }
}
/* Handle the volumic power */
@@ -229,7 +586,7 @@ XD(solid_solid_boundary_path)
if(power != SDIS_VOLUMIC_POWER_NONE) {
const double delta_in_meter = reinject_dst * fp_to_meter;
const double lambda = solid_get_thermal_conductivity(mdm, &rwalk->vtx);
- tmp = delta_in_meter * delta_in_meter / (2.0 * dim * lambda);
+ tmp = delta_in_meter * delta_in_meter / (2.0 * DIM * lambda);
T->value += power * tmp;
if(ctx->green_path) {
@@ -238,9 +595,9 @@ XD(solid_solid_boundary_path)
}
}
- /* Reinject */
+ /* Perform reinjection. */
XD(move_pos)(rwalk->vtx.P, dir, (float)reinject_dst);
- if(eq_epsf(hit->distance, (float)reinject_dst, 1.e-4f)) {
+ if(hit->distance == reinject_dst) {
T->func = XD(boundary_path);
rwalk->mdm = NULL;
rwalk->hit = *hit;
@@ -278,8 +635,8 @@ XD(solid_fluid_boundary_path)
struct sdis_medium* mdm_back = NULL;
struct sdis_medium* solid = NULL;
struct sdis_medium* fluid = NULL;
- struct sXd(hit) hit0 = SXD_HIT_NULL;
- struct sXd(hit) hit1 = SXD_HIT_NULL;
+ struct XD(rwalk) rwalk_saved;
+ struct sXd(hit) hit = SXD_HIT_NULL;
struct sdis_interface_fragment frag_fluid;
double hc;
double hr;
@@ -291,10 +648,13 @@ XD(solid_fluid_boundary_path)
double delta_boundary;
double r;
double tmp;
- float pos[DIM];
float dir0[DIM], dir1[DIM];
- float range[2];
- int dim = DIM;
+ float reinject_dst;
+ /* In 2D it is useless to try to resample a reinjection direction since there
+ * is only one possible direction */
+ const int MAX_ATTEMPTS = DIM == 2 ? 1 : 10;
+ int iattempt;
+ int reinjection_is_valid = 0;
res_T res = RES_OK;
ASSERT(scn && fp_to_meter > 0 && rwalk && rng && T && ctx);
ASSERT(XD(check_rwalk_fragment_consistency)(rwalk, frag));
@@ -325,48 +685,44 @@ XD(solid_fluid_boundary_path)
* delta. */
delta_boundary = sqrt(DIM) * delta;
- /* Sample a reinjection direction */
- XD(sample_reinjection_dir)(rwalk, rng, dir0);
+ rwalk_saved = *rwalk;
+ reinjection_is_valid = 0;
+ iattempt = 0;
+ do {
+ if(iattempt != 0) *rwalk = rwalk_saved;
- /* Reflect the sampled direction around the normal */
- XD(reflect)(dir1, dir0, rwalk->hit.normal);
+ /* Sample a reinjection direction */
+ XD(sample_reinjection_dir)(rwalk, rng, dir0);
- if(solid == mdm_back) {
- fX(minus)(dir0, dir0);
- fX(minus)(dir1, dir1);
- }
+ /* Reflect the sampled direction around the normal */
+ XD(reflect)(dir1, dir0, rwalk->hit.normal);
- /* Trace dir0/dir1 to adjust the reinjection distance */
- fX_set_dX(pos, rwalk->vtx.P);
- f2(range, 0, (float)delta_boundary*RAY_RANGE_MAX_SCALE);
- SXD(scene_view_trace_ray(scn->sXd(view), pos, dir0, range, &rwalk->hit, &hit0));
- SXD(scene_view_trace_ray(scn->sXd(view), pos, dir1, range, &rwalk->hit, &hit1));
-
- /* Adjust the delta boundary to the hit distance */
- tmp = MMIN(MMIN(delta_boundary, hit0.distance), hit1.distance);
-
- if(tmp >= delta_boundary * REINJECT_DST_MIN_SCALE) {
- delta_boundary = tmp;
- /* Define the orthogonal dst from the reinjection pos to the interface */
- delta = delta_boundary / sqrt(DIM);
- } else { /* Switch in 1D reinjection scheme. */
- fX(set)(dir0, rwalk->hit.normal);
- if(solid == mdm_back) fX(minus)(dir0, dir0);
- f2(range, 0, (float)delta*RAY_RANGE_MAX_SCALE);
- SXD(scene_view_trace_ray(scn->sXd(view), pos, dir0, range, &rwalk->hit, &hit0));
- delta_boundary = MMIN(hit0.distance, delta);
-
- /* Hit something in 1D. Arbitrarily move the random walk to 0.5 of the hit
- * distance in order to avoid infinite bounces for parallel plane */
- if(!SXD_HIT_NONE(&hit0)) {
- delta_boundary *= 0.5;
- hit0 = SXD_HIT_NULL;
+ if(solid == mdm_back) {
+ fX(minus)(dir0, dir0);
+ fX(minus)(dir1, dir1);
}
- delta = delta_boundary;
- dim = 1;
+ /* Select the solid reinjection direction and distance */
+ res = XD(select_reinjection_dir_and_check_validity)(scn, solid, rwalk,
+ dir0, dir1, delta_boundary, dir0, &reinject_dst, 1, NULL,
+ &reinjection_is_valid, &hit);
+ if(res != RES_OK) goto error;
+
+ } while(!reinjection_is_valid && ++iattempt < MAX_ATTEMPTS);
+
+ /* Could not find a valid reinjecton */
+ if(iattempt >= MAX_ATTEMPTS) {
+ *rwalk = rwalk_saved;
+ log_err(scn->dev,
+ "%s: could not find a valid solid/fluid reinjection at {%g, %g %g}.\n",
+ FUNC_NAME, SPLIT3(rwalk->vtx.P));
+ res = RES_BAD_OP_IRRECOVERABLE;
+ goto error;
}
+ /* Define the orthogonal dst from the reinjection pos to the interface */
+ delta = reinject_dst / sqrt(DIM);
+
/* Fetch the boundary properties */
epsilon = interface_side_get_emissivity(interf, &frag_fluid);
hc = interface_get_convection_coef(interf, frag);
@@ -393,8 +749,8 @@ XD(solid_fluid_boundary_path)
/* Handle the volumic power */
const double power = solid_get_volumic_power(solid, &rwalk->vtx);
if(power != SDIS_VOLUMIC_POWER_NONE) {
- const double delta_in_meter = delta_boundary * fp_to_meter;
- tmp = delta_in_meter * delta_in_meter / (2.0 * dim * lambda);
+ const double delta_in_meter = reinject_dst * fp_to_meter;
+ tmp = delta_in_meter * delta_in_meter / (2.0 * DIM * lambda);
T->value += power * tmp;
if(ctx->green_path) {
@@ -403,13 +759,13 @@ XD(solid_fluid_boundary_path)
}
}
- /* Reinject */
- XD(move_pos)(rwalk->vtx.P, dir0, (float)delta_boundary);
- if(eq_epsf(hit0.distance, (float)delta_boundary, 1.e-4f)) {
+ /* Perform solid reinjection */
+ XD(move_pos)(rwalk->vtx.P, dir0, reinject_dst);
+ if(hit.distance == reinject_dst) {
T->func = XD(boundary_path);
rwalk->mdm = NULL;
- rwalk->hit = hit0;
- rwalk->hit_side = fX(dot)(hit0.normal, dir0) < 0 ? SDIS_FRONT : SDIS_BACK;
+ rwalk->hit = hit;
+ rwalk->hit_side = fX(dot)(hit.normal, dir0) < 0 ? SDIS_FRONT : SDIS_BACK;
} else {
T->func = XD(conductive_path);
rwalk->mdm = solid;
@@ -440,6 +796,7 @@ XD(solid_boundary_with_flux_path)
struct ssp_rng* rng,
struct XD(temperature)* T)
{
+ struct XD(rwalk) rwalk_saved;
struct sdis_interface* interf = NULL;
struct sdis_medium* mdm = NULL;
double lambda;
@@ -448,13 +805,15 @@ XD(solid_boundary_with_flux_path)
double delta_in_meter;
double power;
double tmp;
- struct sXd(hit) hit0;
- struct sXd(hit) hit1;
- float pos[DIM];
+ struct sXd(hit) hit;
float dir0[DIM];
float dir1[DIM];
- float range[2];
- int dim = DIM;
+ float reinject_dst;
+ /* In 2D it is useless to try to resample a reinjection direction since there
+ * is only one possible direction */
+ const int MAX_ATTEMPTS = DIM == 2 ? 1 : 10;
+ int iattempt = 0;
+ int reinjection_is_valid = 0;
res_T res = RES_OK;
ASSERT(frag && phi != SDIS_FLUX_NONE);
ASSERT(XD(check_rwalk_fragment_consistency)(rwalk, frag));
@@ -476,48 +835,43 @@ XD(solid_boundary_with_flux_path)
* distance from the boundary to the point to chalenge is equal to delta. */
delta_boundary = delta * sqrt(DIM);
- /* Sample a reinjection direction */
- XD(sample_reinjection_dir)(rwalk, rng, dir0);
+ rwalk_saved = *rwalk;
+ reinjection_is_valid = 0;
+ iattempt = 0;
+ do {
+ if(iattempt != 0) *rwalk = rwalk_saved;
+ /* Sample a reinjection direction */
+ XD(sample_reinjection_dir)(rwalk, rng, dir0);
- /* Reflect the sampled direction around the normal */
- XD(reflect)(dir1, dir0, rwalk->hit.normal);
+ /* Reflect the sampled direction around the normal */
+ XD(reflect)(dir1, dir0, rwalk->hit.normal);
- if(frag->side == SDIS_BACK) {
- fX(minus)(dir0, dir0);
- fX(minus)(dir1, dir1);
- }
-
- /* Trace dir0/dir1 to adjust the reinjection distance wrt the geometry */
- fX_set_dX(pos, rwalk->vtx.P);
- f2(range, 0, (float)delta_boundary*RAY_RANGE_MAX_SCALE);
- SXD(scene_view_trace_ray(scn->sXd(view), pos, dir0, range, &rwalk->hit, &hit0));
- SXD(scene_view_trace_ray(scn->sXd(view), pos, dir1, range, &rwalk->hit, &hit1));
-
- /* Adjust the delta boundary to the hit distance */
- tmp = MMIN(MMIN(delta_boundary, hit0.distance), hit1.distance);
-
- if(tmp >= delta_boundary * REINJECT_DST_MIN_SCALE) {
- delta_boundary = tmp;
- /* Define the orthogonal dst from the reinjection pos to the interface */
- delta = delta_boundary / sqrt(DIM);
- } else { /* Switch in 1D reinjection scheme. */
- fX(set)(dir0, rwalk->hit.normal);
- if(frag->side == SDIS_BACK) fX(minus)(dir0, dir0);
- f2(range, 0, (float)delta*RAY_RANGE_MAX_SCALE);
- SXD(scene_view_trace_ray(scn->sXd(view), pos, dir0, range, &rwalk->hit, &hit0));
- delta_boundary = MMIN(hit0.distance, delta_boundary);
-
- /* Hit something in 1D. Arbitrarily move the random walk to 0.5 of the hit
- * distance in order to avoid infinite bounces for parallel plane */
- if(!SXD_HIT_NONE(&hit0)) {
- delta_boundary *= 0.5;
- hit0 = SXD_HIT_NULL;
+ if(frag->side == SDIS_BACK) {
+ fX(minus)(dir0, dir0);
+ fX(minus)(dir1, dir1);
}
- delta = delta_boundary;
- dim = 1;
+ /* Select the reinjection direction and distance */
+ res = XD(select_reinjection_dir_and_check_validity)(scn, mdm, rwalk, dir0,
+ dir1, delta_boundary, dir0, &reinject_dst, 1, NULL,
+ &reinjection_is_valid, &hit);
+ if(res != RES_OK) goto error;
+
+ } while(!reinjection_is_valid && ++iattempt < MAX_ATTEMPTS);
+
+ /* Could not find a valid reinjecton */
+ if(iattempt >= MAX_ATTEMPTS) {
+ *rwalk = rwalk_saved;
+ log_err(scn->dev,
+ "%s: could not find a valid solid/fluid with flux reinjection "
+ "at {%g, %g, %g}.\n", FUNC_NAME, SPLIT3(rwalk->vtx.P));
+ res = RES_BAD_OP_IRRECOVERABLE;
+ goto error;
}
+ /* Define the orthogonal dst from the reinjection pos to the interface */
+ delta = reinject_dst / sqrt(DIM);
+
/* Handle the flux */
delta_in_meter = delta*fp_to_meter;
tmp = delta_in_meter / lambda;
@@ -530,8 +884,8 @@ XD(solid_boundary_with_flux_path)
/* Handle the volumic power */
power = solid_get_volumic_power(mdm, &rwalk->vtx);
if(power != SDIS_VOLUMIC_POWER_NONE) {
- delta_in_meter = delta_boundary * fp_to_meter;
- tmp = delta_in_meter * delta_in_meter / (2.0 * dim * lambda);
+ delta_in_meter = reinject_dst * fp_to_meter;
+ tmp = delta_in_meter * delta_in_meter / (2.0 * DIM * lambda);
T->value += power * tmp;
if(ctx->green_path) {
res = green_path_add_power_term(ctx->green_path, mdm, &rwalk->vtx, tmp);
@@ -539,13 +893,14 @@ XD(solid_boundary_with_flux_path)
}
}
- /* Reinject into the solid */
- XD(move_pos)(rwalk->vtx.P, dir0, (float)delta_boundary);
- if(eq_epsf(hit0.distance, (float)delta_boundary, 1.e-4f)) {
+ /* Reinject. If the reinjection move the point too close of a boundary,
+ * assume that the zone is isotherm and move to the boundary. */
+ XD(move_pos)(rwalk->vtx.P, dir0, reinject_dst);
+ if(hit.distance == reinject_dst) {
T->func = XD(boundary_path);
rwalk->mdm = NULL;
- rwalk->hit = hit0;
- rwalk->hit_side = fX(dot)(hit0.normal, dir0) < 0 ? SDIS_FRONT : SDIS_BACK;
+ rwalk->hit = hit;
+ rwalk->hit_side = fX(dot)(hit.normal, dir0) < 0 ? SDIS_FRONT : SDIS_BACK;
} else {
T->func = XD(conductive_path);
rwalk->mdm = mdm;
diff --git a/src/sdis_heat_path_conductive_Xd.h b/src/sdis_heat_path_conductive_Xd.h
@@ -24,6 +24,177 @@
#include "sdis_Xd_begin.h"
+/*******************************************************************************
+ * Helper functions
+ ******************************************************************************/
+/* Sample the next direction to walk toward and compute the distance to travel.
+ * Return the sampled direction `dir0', the distance to travel along this
+ * direction, the hit `hit0' along `dir0' wrt to the returned distance, the
+ * direction `dir1' used to adjust the displacement distance, and the hit
+ * `hit1' along `dir1' used to adjust the displacement distance. */
+static float
+XD(sample_next_step)
+ (struct sdis_scene* scn,
+ struct ssp_rng* rng,
+ const float pos[DIM],
+ const float delta_solid,
+ float dir0[DIM], /* Sampled direction */
+ float dir1[DIM], /* Direction used to adjust delta */
+ struct sXd(hit)* hit0, /* Hit along the sampled direction */
+ struct sXd(hit)* hit1) /* Hit used to adjust delta */
+{
+ struct sXd(hit) hits[2];
+ float dirs[2][DIM];
+ float range[2];
+ float delta;
+ ASSERT(scn && rng && pos && delta_solid>0 && dir0 && dir1 && hit0 && hit1);
+
+ *hit0 = SXD_HIT_NULL;
+ *hit1 = SXD_HIT_NULL;
+
+#if DIM == 2
+ /* Sample a main direction around 2PI */
+ ssp_ran_circle_uniform_float(rng, dirs[0], NULL);
+#else
+ /* Sample a main direction around 4PI */
+ ssp_ran_sphere_uniform_float(rng, dirs[0], NULL);
+#endif
+
+ /* Negate the sampled dir */
+ fX(minus)(dirs[1], dirs[0]);
+
+ /* Use the previously sampled direction to estimate the minimum distance from
+ * `pos' to the scene boundary */
+ f2(range, 0.f, delta_solid*RAY_RANGE_MAX_SCALE);
+ SXD(scene_view_trace_ray(scn->sXd(view), pos, dirs[0], range, NULL, &hits[0]));
+ SXD(scene_view_trace_ray(scn->sXd(view), pos, dirs[1], range, NULL, &hits[1]));
+ if(SXD_HIT_NONE(&hits[0]) && SXD_HIT_NONE(&hits[1])) {
+ delta = delta_solid;
+ } else {
+ delta = MMIN(hits[0].distance, hits[1].distance);
+ }
+
+ if(!SXD_HIT_NONE(&hits[0])
+ && delta != hits[0].distance
+ && eq_eps(hits[0].distance, delta, delta_solid*0.1)) {
+ /* Set delta to the main hit distance if it is roughly equal to it in order
+ * to avoid numerical issues on moving along the main direction. */
+ delta = hits[0].distance;
+ }
+
+ /* Setup outputs */
+ if(delta <= delta_solid*0.1 && hits[1].distance == delta) {
+ /* Snap the random walk to the boundary if delta is too small */
+ fX(set)(dir0, dirs[1]);
+ *hit0 = hits[1];
+ fX(splat)(dir1, (float)INF);
+ *hit1 = SXD_HIT_NULL;
+ } else {
+ fX(set)(dir0, dirs[0]);
+ *hit0 = hits[0];
+ if(delta == hits[0].distance) {
+ fX(set)(dir1, dirs[0]);
+ *hit1 = hits[0];
+ } else if(delta == hits[1].distance) {
+ fX(set)(dir1, dirs[1]);
+ *hit1 = hits[1];
+ } else {
+ fX(splat)(dir1, 0);
+ *hit1 = SXD_HIT_NULL;
+ }
+ }
+
+ return delta;
+}
+
+/* Sample the next direction to walk toward and compute the distance to travel.
+ * If the targeted position does not lie inside the current medium, reject it
+ * and sample a new next step. */
+static res_T
+XD(sample_next_step_robust)
+ (struct sdis_scene* scn,
+ struct sdis_medium* current_mdm,
+ struct ssp_rng* rng,
+ const double pos[DIM],
+ const float delta_solid,
+ float dir0[DIM], /* Sampled direction */
+ float dir1[DIM], /* Direction used to adjust delta */
+ struct sXd(hit)* hit0, /* Hit along the sampled direction */
+ struct sXd(hit)* hit1, /* Hit used to adjust delta */
+ float* out_delta)
+{
+ struct sdis_medium* mdm;
+ float delta;
+ float org[DIM];
+ const size_t MAX_ATTEMPTS = 100;
+ size_t iattempt = 0;
+ res_T res = RES_OK;
+ ASSERT(scn && current_mdm && rng && pos && delta_solid > 0);
+ ASSERT(dir0 && dir1 && hit0 && hit1 && out_delta);
+
+ fX_set_dX(org, pos);
+ do {
+ double pos_next[DIM];
+
+ /* Compute the next step */
+ delta = XD(sample_next_step)
+ (scn, rng, org, delta_solid, dir0, dir1, hit0, hit1);
+
+ /* Retrieve the medium of the next step */
+ if(hit0->distance > delta) {
+ XD(move_pos)(dX(set)(pos_next, pos), dir0, delta);
+ res = scene_get_medium_in_closed_boundaries(scn, pos_next, &mdm);
+ if(res != RES_OK) goto error;
+ } else {
+ struct sdis_interface* interf;
+ enum sdis_side side;
+ interf = scene_get_interface(scn, hit0->prim.prim_id);
+ side = fX(dot)(dir0, hit0->normal) < 0 ? SDIS_FRONT : SDIS_BACK;
+ mdm = interface_get_medium(interf, side);
+ }
+
+ /* Check medium consistency */
+ if(current_mdm != mdm) {
+#if 0
+#if DIM == 2
+ log_err(scn->dev,
+ "%s: inconsistent medium during the solid random walk at {%g, %g}.\n",
+ FUNC_NAME, SPLIT2(pos));
+#else
+ log_err(scn->dev,
+ "%s: inconsistent medium during the solid random walk at {%g, %g, %g}.\n",
+ FUNC_NAME, SPLIT3(pos));
+#endif
+#endif
+ }
+ } while(current_mdm != mdm && ++iattempt < MAX_ATTEMPTS);
+
+ /* Handle error */
+ if(iattempt >= MAX_ATTEMPTS) {
+#if DIM == 2
+ log_err(scn->dev,
+ "%s: could not find a next valid conductive step at {%g, %g}.\n",
+ FUNC_NAME, SPLIT2(pos));
+#else
+ log_err(scn->dev,
+ "%s: could not find a next valid conductive step at {%g, %g, %g}.\n",
+ FUNC_NAME, SPLIT3(pos));
+#endif
+ res = RES_BAD_OP;
+ goto error;
+ }
+
+ *out_delta = delta;
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
+/*******************************************************************************
+ * Local function
+ ******************************************************************************/
res_T
XD(conductive_path)
(struct sdis_scene* scn,
@@ -44,8 +215,8 @@ XD(conductive_path)
(void)ctx, (void)istep;
/* Check the random walk consistency */
- CHK(scene_get_medium(scn, rwalk->vtx.P, NULL, &mdm) == RES_OK);
- if(mdm != rwalk->mdm) {
+ res = scene_get_medium_in_closed_boundaries(scn, rwalk->vtx.P, &mdm);
+ if(res != RES_OK || 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));
res = RES_BAD_OP_IRRECOVERABLE;
@@ -61,7 +232,6 @@ XD(conductive_path)
}
do { /* Solid random walk */
- struct get_medium_info info = GET_MEDIUM_INFO_NULL;
struct sXd(hit) hit0, hit1;
double lambda; /* Thermal conductivity */
double rho; /* Volumic mass */
@@ -70,7 +240,6 @@ XD(conductive_path)
double power_factor = 0;
double power;
float delta, delta_solid; /* Random walk numerical parameter */
- float range[2];
float dir0[DIM], dir1[DIM];
float org[DIM];
@@ -109,39 +278,20 @@ XD(conductive_path)
goto error;
}
-#if DIM == 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;
-
- /* Add the volumic power density to the measured temperature */
- if(power != SDIS_VOLUMIC_POWER_NONE) {
+
+ /* Sample the direction to walk toward and compute the distance to travel */
+ res = XD(sample_next_step_robust)(scn, mdm, rng, rwalk->vtx.P, delta_solid,
+ dir0, dir1, &hit0, &hit1, &delta);
+ if(res != RES_OK) goto error;
+
+ /* Add the volumic power density to the measured temperature */
+ if(power != SDIS_VOLUMIC_POWER_NONE) {
+ if((S3D_HIT_NONE(&hit0) && S3D_HIT_NONE(&hit1))) { /* Hit nothing */
const double delta_in_meter = delta * fp_to_meter;
power_factor = delta_in_meter * delta_in_meter / (2.0 * DIM * lambda);
T->value += power * power_factor;
- }
- } else {
- /* Hit something: move along dir0 of the minimum hit distance */
- delta = MMIN(hit0.distance, hit1.distance);
-
- /* Add the volumic power density to the measured temperature */
- if(power != SDIS_VOLUMIC_POWER_NONE) {
+ } else {
const double delta_s_adjusted = delta_solid * RAY_RANGE_MAX_SCALE;
const double delta_s_in_meter = delta_solid * fp_to_meter;
double h;
@@ -226,10 +376,8 @@ XD(conductive_path)
}
}
- /* Define if the random walk hits something along dir0. Multiply delta by
- * the empirical ray range scale factor to ensure that once moved, the
- * random walk does not lie in the uncertainty zone near the geometry */
- if(hit0.distance > delta * RAY_RANGE_MAX_SCALE) {
+ /* Define if the random walk hits something along dir0 */
+ if(hit0.distance > delta) {
rwalk->hit = SXD_HIT_NULL;
rwalk->hit_side = SDIS_SIDE_NULL__;
} else {
@@ -245,45 +393,6 @@ XD(conductive_path)
(ctx->heat_path, &rwalk->vtx, T->value, SDIS_HEAT_VERTEX_CONDUCTION);
if(res != RES_OK) goto error;
- /* Fetch the current medium */
- if(SXD_HIT_NONE(&rwalk->hit)) {
- CHK(scene_get_medium(scn, rwalk->vtx.P, &info, &mdm) == RES_OK);
- } else {
- const struct sdis_interface* interf;
- interf = scene_get_interface(scn, rwalk->hit.prim.prim_id);
- mdm = interface_get_medium(interf, rwalk->hit_side);
- }
-
- /* 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
- #define VEC_STR "%g %g"
- #define VEC_SPLIT SPLIT2
-#else
- #define VEC_STR "%g %g %g"
- #define VEC_SPLIT SPLIT3
-#endif
- log_err(scn->dev,
- " start position: " VEC_STR "; current position: " VEC_STR "\n",
- VEC_SPLIT(position_start), VEC_SPLIT(rwalk->vtx.P));
- if(SXD_HIT_NONE(&rwalk->hit)) {
- float hit_pos[DIM];
- fX(mulf)(hit_pos, info.ray_dir, info.XD(hit).distance);
- fX(add)(hit_pos, info.ray_org, hit_pos);
- log_err(scn->dev, " ray org: " VEC_STR "; ray dir: " VEC_STR "\n",
- VEC_SPLIT(info.ray_org), VEC_SPLIT(info.ray_dir));
- log_err(scn->dev, " targeted point: " VEC_STR "\n",
- VEC_SPLIT(info.pos_tgt));
- log_err(scn->dev, " hit pos: " VEC_STR "\n", VEC_SPLIT(hit_pos));
- }
-#undef VEC_STR
-#undef VEC_SPLIT
- res = RES_BAD_OP;
- goto error;
- }
-
++istep;
/* Keep going while the solid random walk does not hit an interface */
diff --git a/src/sdis_heat_path_convective_Xd.h b/src/sdis_heat_path_convective_Xd.h
@@ -34,6 +34,7 @@ XD(register_heat_vertex_in_fluid)
const double weight)
{
struct sdis_rwalk_vertex vtx = SDIS_RWALK_VERTEX_NULL;
+ struct hit_filter_data filter_data;
const float empirical_dst = 0.1f;
const float range[2] = {0, FLT_MAX};
float org[DIM];
@@ -50,7 +51,9 @@ XD(register_heat_vertex_in_fluid)
fX(set)(dir, rwalk->hit.normal);
if(rwalk->hit_side == SDIS_BACK) fX(minus)(dir, dir);
- SXD(scene_view_trace_ray(scn->sXd(view), org, dir, range, &rwalk->hit, &hit));
+ filter_data.XD(hit) = rwalk->hit;
+ filter_data.epsilon = 1.e-6;
+ SXD(scene_view_trace_ray(scn->sXd(view), org, dir, range, &filter_data, &hit));
dst = SXD_HIT_NONE(&hit) ? empirical_dst : hit.distance * 0.5f;
vtx = rwalk->vtx;
diff --git a/src/sdis_heat_path_radiative_Xd.h b/src/sdis_heat_path_radiative_Xd.h
@@ -52,6 +52,7 @@ XD(trace_radiative_path)
/* Launch the radiative random walk */
for(;;) {
const struct sdis_interface* interf = NULL;
+ struct hit_filter_data filter_data;
struct sdis_interface_fragment frag = SDIS_INTERFACE_FRAGMENT_NULL;
struct sdis_medium* chk_mdm = NULL;
double alpha;
@@ -63,12 +64,14 @@ XD(trace_radiative_path)
fX_set_dX(pos, rwalk->vtx.P);
/* Trace the radiative ray */
+ filter_data.XD(hit) = rwalk->hit;
+ filter_data.epsilon = 1.e-6;
#if (SDIS_XD_DIMENSION == 2)
SXD(scene_view_trace_ray_3d
- (scn->sXd(view), pos, dir, range, &rwalk->hit, &rwalk->hit));
+ (scn->sXd(view), pos, dir, range, &filter_data, &rwalk->hit));
#else
SXD(scene_view_trace_ray
- (scn->sXd(view), pos, dir, range, &rwalk->hit, &rwalk->hit));
+ (scn->sXd(view), pos, dir, range, &filter_data, &rwalk->hit));
#endif
if(SXD_HIT_NONE(&rwalk->hit)) { /* Fetch the ambient radiative temperature */
rwalk->hit_side = SDIS_SIDE_NULL__;
diff --git a/src/sdis_realisation_Xd.h b/src/sdis_realisation_Xd.h
@@ -45,9 +45,9 @@ XD(compute_temperature)
}* stack = NULL;
size_t istack = 0;
#endif
- /* Maximum accepted #failures before stopping the realisation */
struct sdis_heat_vertex* heat_vtx = NULL;
- const size_t MAX_FAILS = 10;
+ /* Maximum accepted #failures before stopping the realisation */
+ const size_t MAX_FAILS = 1;
res_T res = RES_OK;
ASSERT(scn && fp_to_meter > 0 && ctx && rwalk && rng && T);
diff --git a/src/sdis_scene.c b/src/sdis_scene.c
@@ -373,3 +373,14 @@ scene_get_medium
: scene_get_medium_3d(scn, pos, info, out_medium);
}
+res_T
+scene_get_medium_in_closed_boundaries
+ (const struct sdis_scene* scn,
+ const double pos[],
+ struct sdis_medium** out_medium)
+{
+ return scene_is_2d(scn)
+ ? scene_get_medium_in_closed_boundaries_2d(scn, pos, out_medium)
+ : scene_get_medium_in_closed_boundaries_3d(scn, pos, out_medium);
+}
+
diff --git a/src/sdis_scene_Xd.h b/src/sdis_scene_Xd.h
@@ -22,8 +22,14 @@
#include "sdis_medium_c.h"
#include "sdis_scene_c.h"
+#include <star/ssp.h>
+#include <rsys/float22.h>
+#include <rsys/float33.h>
#include <rsys/rsys.h>
+/* Emperical cos threshold defining if an angle is sharp */
+#define SHARP_ANGLE_COS_THRESOLD -0.70710678 /* ~ cos(3*PI/4) */
+
/*******************************************************************************
* Define the helper functions and the data types used by the scene
* independently of its dimension, i.e. 2D or 3D.
@@ -148,6 +154,7 @@ clear_properties(struct sdis_scene* scn)
/* Vector macros generic to SDIS_SCENE_DIMENSION */
#define fX(Func) CONCAT(CONCAT(CONCAT(f, DIM), _), Func)
#define fX_set_dX CONCAT(CONCAT(CONCAT(f, DIM), _set_d), DIM)
+#define fXX_mulfX CONCAT(CONCAT(CONCAT(CONCAT(f, DIM), DIM), _mulf), DIM)
/* Macro making generic its subimitted name to SDIS_SCENE_DIMENSION */
#define XD(Name) CONCAT(CONCAT(CONCAT(Name, _), DIM), d)
@@ -162,40 +169,261 @@ clear_properties(struct sdis_scene* scn)
* Helper functions
******************************************************************************/
#if DIM == 2
-/* Check that `hit' roughly lies on a vertex. For segments, a simple but
- * approximative way is to test that its position have at least one barycentric
- * coordinate roughly equal to 0 or 1. */
-static FINLINE int
-hit_on_vertex(const struct s2d_hit* hit)
+#define ON_VERTEX_EPSILON 1.e-4f
+/* Check that `hit' roughly lies on a vertex. */
+static INLINE int
+hit_on_vertex
+ (const struct s2d_hit* hit,
+ const float org[3],
+ const float dir[3])
+{
+ struct s2d_attrib v0, v1;
+ float E[2];
+ float hit_pos[2];
+ float segment_len;
+ float hit_len0;
+ float hit_len1;
+ ASSERT(hit && !S2D_HIT_NONE(hit) && org && dir);
+
+ /* Rertieve the segment vertices */
+ S2D(segment_get_vertex_attrib(&hit->prim, 0, S2D_POSITION, &v0));
+ S2D(segment_get_vertex_attrib(&hit->prim, 1, S2D_POSITION, &v1));
+
+ /* Compute the length of the segment */
+ segment_len = f2_len(f2_sub(E, v1.value, v0.value));
+
+ /* Compute the hit position onto the segment */
+ f2_add(hit_pos, org, f2_mulf(hit_pos, dir, hit->distance));
+
+ /* Compute the length from hit position to segment vertices */
+ hit_len0 = f2_len(f2_sub(E, v0.value, hit_pos));
+ hit_len1 = f2_len(f2_sub(E, v1.value, hit_pos));
+
+ if(hit_len0 / segment_len < ON_VERTEX_EPSILON
+ || hit_len1 / segment_len < ON_VERTEX_EPSILON)
+ return 1;
+ return 0;
+}
+
+static int
+hit_shared_vertex
+ (const struct s2d_primitive* seg0,
+ const struct s2d_primitive* seg1,
+ const float pos0[2], /* Tested position onto the segment 0 */
+ const float pos1[2]) /* Tested Position onto the segment 1 */
{
- const float on_vertex_eps = 1.e-4f;
- float v;
- ASSERT(hit && !S2D_HIT_NONE(hit));
- v = 1.f - hit->u;
- return eq_epsf(hit->u, 0.f, on_vertex_eps)
- || eq_epsf(hit->u, 1.f, on_vertex_eps)
- || eq_epsf(v, 0.f, on_vertex_eps)
- || eq_epsf(v, 1.f, on_vertex_eps);
+ struct s2d_attrib seg0_vertices[2]; /* Vertex positions of the segment 0 */
+ struct s2d_attrib seg1_vertices[2]; /* Vertex positions of the segment 1 */
+ float d0[2], d1[2]; /* temporary vector */
+ float seg0_len, seg1_len; /* Length of the segments */
+ float tmp0_len, tmp1_len;
+ float cos_normals;
+ int seg0_vert = -1; /* Id of the shared vertex for the segment 0 */
+ int seg1_vert = -1; /* Id of the shared vertex for the segment 1 */
+ int seg0_ivertex, seg1_ivertex;
+ ASSERT(seg0 && seg1 && pos0 && pos1);
+
+ /* Fetch the vertices of the segment 0 */
+ S2D(segment_get_vertex_attrib(seg0, 0, S2D_POSITION, &seg0_vertices[0]));
+ S2D(segment_get_vertex_attrib(seg0, 1, S2D_POSITION, &seg0_vertices[1]));
+
+ /* Fetch the vertices of the segment 1 */
+ S2D(segment_get_vertex_attrib(seg1, 0, S2D_POSITION, &seg1_vertices[0]));
+ S2D(segment_get_vertex_attrib(seg1, 1, S2D_POSITION, &seg1_vertices[1]));
+
+ /* Look for the vertex shared by the 2 segments */
+ for(seg0_ivertex = 0; seg0_ivertex < 2 && seg0_vert < 0; ++seg0_ivertex) {
+ for(seg1_ivertex = 0; seg1_ivertex < 2 && seg1_vert < 0; ++seg1_ivertex) {
+ const int vertex_eq = f2_eq_eps
+ (seg0_vertices[seg0_ivertex].value,
+ seg1_vertices[seg1_ivertex].value,
+ 1.e-6f);
+ if(vertex_eq) {
+ seg0_vert = seg0_ivertex;
+ seg1_vert = seg1_ivertex;
+ /* We assume that the segments are not degenerated. As a consequence we
+ * can break here since a vertex of the segment 0 can be equal to at most
+ * one vertex of the segment 1 */
+ break;
+ }
+ }}
+
+ /* The segments do not have a common vertex */
+ if(seg0_vert < 0) return 0;
+
+ /* Compute the dirctions from shared vertex to the opposite segment vertex */
+ f2_sub(d0, seg0_vertices[(seg0_vert+1)%2].value, seg0_vertices[seg0_vert].value);
+ f2_sub(d1, seg1_vertices[(seg1_vert+1)%2].value, seg1_vertices[seg1_vert].value);
+
+ /* Compute the cosine between the segments */
+ seg0_len = f2_normalize(d0, d0);
+ seg1_len = f2_normalize(d1, d1);
+ cos_normals = f2_dot(d0, d1);
+
+ /* The angle formed by the 2 segments is sharp. Do not filter the hit */
+ if(cos_normals > SHARP_ANGLE_COS_THRESOLD) return 0;
+
+ /* Compute the length from pos<0|1> to shared vertex */
+ f2_sub(d0, seg0_vertices[seg0_vert].value, pos0);
+ f2_sub(d1, seg1_vertices[seg1_vert].value, pos1);
+ tmp0_len = f2_len(d0);
+ tmp1_len = f2_len(d1);
+
+ return (eq_epsf(seg0_len, 0, 1.e-6f) || tmp0_len/seg0_len < ON_VERTEX_EPSILON)
+ && (eq_epsf(seg1_len, 0, 1.e-6f) || tmp1_len/seg1_len < ON_VERTEX_EPSILON);
}
#else /* DIM == 3 */
-/* Check that `hit' roughly lies on an edge. For triangular primitives, a
- * simple but approximative way is to test that its position have at least one
- * barycentric coordinate roughly equal to 0 or 1. */
-static FINLINE int
-hit_on_edge(const struct s3d_hit* hit)
+#define ON_EDGE_EPSILON 1.e-4f
+/* Check that `hit' roughly lies on an edge. */
+static INLINE int
+hit_on_edge
+ (const struct s3d_hit* hit,
+ const float org[3],
+ const float dir[3])
+{
+ struct s3d_attrib v0, v1, v2;
+ float E0[3], E1[3], N[3];
+ float tri_2area;
+ float hit_2area0;
+ float hit_2area1;
+ float hit_2area2;
+ float hit_pos[3];
+ ASSERT(hit && !S3D_HIT_NONE(hit) && org && dir);
+
+ /* Retrieve the triangle vertices */
+ S3D(triangle_get_vertex_attrib(&hit->prim, 0, S3D_POSITION, &v0));
+ S3D(triangle_get_vertex_attrib(&hit->prim, 1, S3D_POSITION, &v1));
+ S3D(triangle_get_vertex_attrib(&hit->prim, 2, S3D_POSITION, &v2));
+
+ /* Compute the triangle area * 2 */
+ f3_sub(E0, v1.value, v0.value);
+ f3_sub(E1, v2.value, v0.value);
+ tri_2area = f3_len(f3_cross(N, E0, E1));
+
+ /* Compute the hit position */
+ f3_add(hit_pos, org, f3_mulf(hit_pos, dir, hit->distance));
+
+ /* Compute areas */
+ f3_sub(E0, v0.value, hit_pos);
+ f3_sub(E1, v1.value, hit_pos);
+ hit_2area0 = f3_len(f3_cross(N, E0, E1));
+ f3_sub(E0, v1.value, hit_pos);
+ f3_sub(E1, v2.value, hit_pos);
+ hit_2area1 = f3_len(f3_cross(N, E0, E1));
+ f3_sub(E0, v2.value, hit_pos);
+ f3_sub(E1, v0.value, hit_pos);
+ hit_2area2 = f3_len(f3_cross(N, E0, E1));
+
+ if(hit_2area0 / tri_2area < ON_EDGE_EPSILON
+ || hit_2area1 / tri_2area < ON_EDGE_EPSILON
+ || hit_2area2 / tri_2area < ON_EDGE_EPSILON)
+ return 1;
+
+ return 0;
+}
+
+static int
+hit_shared_edge
+ (const struct s3d_primitive* tri0,
+ const struct s3d_primitive* tri1,
+ const float pos0[3], /* Tested position onto the triangle 0 */
+ const float pos1[3]) /* Tested Position onto the triangle 1 */
{
- const float on_edge_eps = 1.e-4f;
- float w;
- ASSERT(hit && !S3D_HIT_NONE(hit));
- w = 1.f - hit->uv[0] - hit->uv[1];
- return eq_epsf(hit->uv[0], 0.f, on_edge_eps)
- || eq_epsf(hit->uv[0], 1.f, on_edge_eps)
- || eq_epsf(hit->uv[1], 0.f, on_edge_eps)
- || eq_epsf(hit->uv[1], 1.f, on_edge_eps)
- || eq_epsf(w, 0.f, on_edge_eps)
- || eq_epsf(w, 1.f, on_edge_eps);
+ struct s3d_attrib tri0_vertices[3]; /* Vertex positions of the triangle 0 */
+ struct s3d_attrib tri1_vertices[3]; /* Vertex positions of the triangle 1 */
+ float E0[3], E1[3]; /* Temporary variables storing triangle edges */
+ float N0[3], N1[3]; /* Temporary Normals */
+ float tri0_2area, tri1_2area; /* 2*area of the submitted triangles */
+ float tmp0_2area, tmp1_2area;
+ float cos_normals;
+ int tri0_edge[2] = {-1, -1}; /* Shared edge vertex ids for the triangle 0 */
+ int tri1_edge[2] = {-1, -1}; /* Shared edge vertex ids for the triangle 1 */
+ int edge_ivertex = 0; /* Temporary variable */
+ int tri0_ivertex, tri1_ivertex;
+ int iv0, iv1, iv2;
+ ASSERT(tri0 && tri1 && pos0 && pos1);
+
+ /* Fetch the vertices of the triangle 0 */
+ S3D(triangle_get_vertex_attrib(tri0, 0, S3D_POSITION, &tri0_vertices[0]));
+ S3D(triangle_get_vertex_attrib(tri0, 1, S3D_POSITION, &tri0_vertices[1]));
+ S3D(triangle_get_vertex_attrib(tri0, 2, S3D_POSITION, &tri0_vertices[2]));
+
+ /* Fetch the vertices of the triangle 1 */
+ S3D(triangle_get_vertex_attrib(tri1, 0, S3D_POSITION, &tri1_vertices[0]));
+ S3D(triangle_get_vertex_attrib(tri1, 1, S3D_POSITION, &tri1_vertices[1]));
+ S3D(triangle_get_vertex_attrib(tri1, 2, S3D_POSITION, &tri1_vertices[2]));
+
+ /* Look for the vertices shared by the 2 triangles */
+ for(tri0_ivertex=0; tri0_ivertex < 3 && edge_ivertex < 2; ++tri0_ivertex) {
+ for(tri1_ivertex=0; tri1_ivertex < 3 && edge_ivertex < 2; ++tri1_ivertex) {
+ const int vertex_eq = f3_eq_eps
+ (tri0_vertices[tri0_ivertex].value,
+ tri1_vertices[tri1_ivertex].value,
+ 1.e-6f);
+ if(vertex_eq) {
+ tri0_edge[edge_ivertex] = tri0_ivertex;
+ tri1_edge[edge_ivertex] = tri1_ivertex;
+ ++edge_ivertex;
+ /* We assume that the triangles are not degenerated. As a consequence we
+ * can break here since a vertex of the triangle 0 can be equal to at
+ * most one vertex of the triangle 1 */
+ break;
+ }
+ }}
+
+ /* The triangles do not have a common edge */
+ if(edge_ivertex < 2) return 0;
+
+ /* Ensure that the vertices of the shared edge are registered in the right
+ * order regarding the triangle vertices, i.e. (0,1), (1,2) or (2,0) */
+ if((tri0_edge[0]+1)%3 != tri0_edge[1]) SWAP(int, tri0_edge[0], tri0_edge[1]);
+ if((tri1_edge[0]+1)%3 != tri1_edge[1]) SWAP(int, tri1_edge[0], tri1_edge[1]);
+
+ /* Compute the shared edge normal lying in the triangle 0 plane */
+ iv0 = tri0_edge[0];
+ iv1 = tri0_edge[1];
+ iv2 = (tri0_edge[1]+1) % 3;
+ f3_sub(E0, tri0_vertices[iv1].value, tri0_vertices[iv0].value);
+ f3_sub(E1, tri0_vertices[iv2].value, tri0_vertices[iv0].value);
+ f3_cross(N0, E0, E1); /* Triangle 0 normal */
+ tri0_2area = f3_len(N0);
+ f3_cross(N0, N0, E0);
+
+ /* Compute the shared edge normal lying in the triangle 1 plane */
+ iv0 = tri1_edge[0];
+ iv1 = tri1_edge[1];
+ iv2 = (tri1_edge[1]+1) % 3;
+ f3_sub(E0, tri1_vertices[iv1].value, tri1_vertices[iv0].value);
+ f3_sub(E1, tri1_vertices[iv2].value, tri1_vertices[iv0].value);
+ f3_cross(N1, E0, E1);
+ tri1_2area = f3_len(N1);
+ f3_cross(N1, N1, E0);
+
+ /* Compute the cosine between the 2 edge normals */
+ f3_normalize(N0, N0);
+ f3_normalize(N1, N1);
+ cos_normals = f3_dot(N0, N1);
+
+ /* The angle formed by the 2 triangles is sharp */
+ if(cos_normals > SHARP_ANGLE_COS_THRESOLD) return 0;
+
+ /* Compute the 2 times the area of the (pos0, shared_edge.vertex0,
+ * shared_edge.vertex1) triangles */
+ f3_sub(E0, tri0_vertices[tri0_edge[0]].value, pos0);
+ f3_sub(E1, tri0_vertices[tri0_edge[1]].value, pos0);
+ tmp0_2area = f3_len(f3_cross(N0, E0, E1));
+
+ /* Compute the 2 times the area of the (pos1, shared_edge.vertex0,
+ * shared_edge.vertex1) triangles */
+ f3_sub(E0, tri1_vertices[tri1_edge[0]].value, pos1);
+ f3_sub(E1, tri1_vertices[tri1_edge[1]].value, pos1);
+ tmp1_2area = f3_len(f3_cross(N1, E0, E1));
+
+ return (eq_epsf(tri0_2area, 0, 1.e-6f) || tmp0_2area/tri0_2area < ON_EDGE_EPSILON)
+ && (eq_epsf(tri1_2area, 0, 1.e-6f) || tmp1_2area/tri1_2area < ON_EDGE_EPSILON);
}
+#undef ON_EDGE_EPSILON
#endif /* DIM == 2 */
/* Avoid self-intersection for a ray starting from a planar primitive, i.e. a
@@ -206,22 +434,25 @@ XD(hit_filter_function)
const float org[DIM],
const float dir[DIM],
void* ray_data,
- void* filter_data)
+ void* global_data)
{
- const struct sXd(hit)* hit_from = ray_data;
- (void)org, (void)dir, (void)filter_data;
+ const struct hit_filter_data* filter_data = ray_data;
+ const struct sXd(hit)* hit_from = &filter_data->XD(hit);
+ (void)org, (void)dir, (void)global_data;
- if(!hit_from || SXD_HIT_NONE(hit_from)) return 0; /* No filtering */
+ if(!ray_data || SXD_HIT_NONE(hit_from)) return 0; /* No filtering */
if(SXD_PRIMITIVE_EQ(&hit_from->prim, &hit->prim)) return 1;
- if(eq_epsf(hit->distance, 0, 1.e-6f)) {
+ if(eq_epsf(hit->distance, 0, (float)filter_data->epsilon)) {
+ float pos[DIM];
+ fX(add)(pos, org, fX(mulf)(pos, dir, hit->distance));
/* If the targeted point is near of the origin, check that it lies on an
- * edge/vertex shared by the 2 primitives. */
+ * edge/vertex shared by the 2 primitives */
#if DIM == 2
- return hit_on_vertex(hit_from) && hit_on_vertex(hit);
+ return hit_shared_vertex(&hit_from->prim, &hit->prim, org, pos);
#else
- return hit_on_edge(hit_from) && hit_on_edge(hit);
+ return hit_shared_edge(&hit_from->prim, &hit->prim, org, pos);
#endif
}
return 0;
@@ -321,7 +552,7 @@ XD(run_analyze)
ASSERT(scn && nprims && indices && interf && nverts && position && out_desc);
res = sencXd(device_create)(scn->dev->logger, scn->dev->allocator,
- scn->dev->nthreads, scn->dev->verbose, &senc);
+ 1/*scn->dev->nthreads*/, scn->dev->verbose, &senc);
if(res != RES_OK) goto error;
res = sencXd(scene_create)(senc,
@@ -801,6 +1032,7 @@ XD(scene_get_medium)
size_t iprim, nprims;
size_t nfailures = 0;
const size_t max_failures = 10;
+ float P[DIM];
/* Range of the parametric coordinate into which positions are challenged */
#if DIM == 2
float st[3];
@@ -821,23 +1053,37 @@ XD(scene_get_medium)
f2(st[2], 5.f/12.f, 5.f/12.f);
#endif
+ fX_set_dX(P, pos);
+
SXD(scene_view_primitives_count(scn->sXd(view), &nprims));
FOR_EACH(iprim, 0, nprims) {
struct sXd(hit) hit;
struct sXd(attrib) attr;
struct sXd(primitive) prim;
+ size_t iprim2;
const float range[2] = {0.f, FLT_MAX};
- float N[DIM], P[DIM], dir[DIM], cos_N_dir;
+ float N[DIM], dir[DIM], cos_N_dir;
size_t istep = 0;
+ /* 1 primitive over 2, take a primitive from the end of the primitive list.
+ * When primitives are sorted in a coherent manner regarding their
+ * position, this strategy avoids to test primitives that are going to be
+ * rejected of the same manner due to possible numerical issues of the
+ * resulting intersection. */
+ if((iprim % 2) == 0) {
+ iprim2 = iprim / 2;
+ } else {
+ iprim2 = nprims - 1 - (iprim / 2);
+ }
+
do {
/* Retrieve a position onto the primitive */
- SXD(scene_view_get_primitive(scn->sXd(view), (unsigned)iprim, &prim));
+ SXD(scene_view_get_primitive(scn->sXd(view), (unsigned)iprim2, &prim));
SXD(primitive_get_attrib(&prim, SXD_POSITION, st[istep], &attr));
/* Trace a ray from the random walk vertex toward the retrieved primitive
* position */
- fX(normalize)(dir, fX(sub)(dir, attr.value, fX_set_dX(P, pos)));
+ fX(normalize)(dir, fX(sub)(dir, attr.value, P));
SXD(scene_view_trace_ray(scn->sXd(view), P, dir, range, NULL, &hit));
/* Unforeseen error. One has to intersect a primitive ! */
@@ -850,9 +1096,9 @@ XD(scene_get_medium)
goto error;
}
}
- /* Discard the hit if it is on a vertex, i.e. between 2 segments, and
- * target a new position onto the current primitive */
- } while((SXD_HIT_NONE(&hit) || HIT_ON_BOUNDARY(&hit))
+ /* Discard the hit if it is on a vertex/edge, and target a new position
+ * onto the current primitive */
+ } while((SXD_HIT_NONE(&hit) || HIT_ON_BOUNDARY(&hit, P, dir))
&& ++istep < nsteps);
/* The hits of all targeted positions on the current primitive are on
@@ -863,7 +1109,7 @@ XD(scene_get_medium)
cos_N_dir = fX(dot)(N, dir);
/* Not too close and not roughly orthognonal */
- if(hit.distance > 1.e-6 || absf(cos_N_dir) > 1.e-1f) {
+ if(hit.distance > 1.e-6 && absf(cos_N_dir) > 1.e-2f) {
const struct sdis_interface* interf;
interf = scene_get_interface(scn, hit.prim.prim_id);
medium = interface_get_medium
@@ -880,6 +1126,18 @@ XD(scene_get_medium)
}
}
+ if(iprim >= nprims) {
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ if(iprim > 10 && iprim > (size_t)((double)nprims * 0.05)) {
+ log_warn(scn->dev,
+ "%s: performance issue. Up to %lu primitives were tested to define the "
+ "current medium at {%g, %g, %g}.\n",
+ FUNC_NAME, (unsigned long)iprim, SPLIT3(P));
+ }
+
exit:
*out_medium = medium;
return res;
@@ -893,6 +1151,83 @@ error:
#endif
goto exit;
}
+
+static INLINE res_T
+XD(scene_get_medium_in_closed_boundaries)
+ (const struct sdis_scene* scn,
+ const double pos[DIM],
+ struct sdis_medium** out_medium)
+{
+ struct sdis_medium* medium = NULL;
+ float P[DIM];
+ float frame[DIM*DIM];
+ float dirs[6][3] = {{1,0,0},{-1,0,0},{0,1,0},{0,-1,0},{0,0,1},{0,0,-1}};
+ int idir;
+ res_T res = RES_OK;
+ ASSERT(scn && pos);
+
+ /* Build a frame that will be used to rotate the main axis by PI/4 around
+ * each axis. This can avoid numerical issues when geometry is discretized
+ * along the main axis */
+#if DIM == 2
+ f22_rotation(frame, (float)PI/4);
+#else
+/* N[0] = N[1] = N[2] = (float)(1.0 / sqrt(3.0));*/
+/* f33_basis(frame, N);*/
+ f33_rotation(frame, (float)PI/4, (float)PI/4, (float)PI/4);
+#endif
+
+ fX_set_dX(P, pos);
+ FOR_EACH(idir, 0, 2*DIM) {
+ struct sXd(hit) hit;
+ float N[DIM];
+ const float range[2] = {0.f, FLT_MAX};
+ float cos_N_dir;
+
+ /* Transform the directions to avoid to be aligned with the axis */
+ fXX_mulfX(dirs[idir], frame, dirs[idir]);
+
+ /* Trace a ray from the random walk vertex toward the retrieved primitive
+ * position */
+ SXD(scene_view_trace_ray(scn->sXd(view), P, dirs[idir], range, NULL, &hit));
+
+ /* Unforeseen error. One has to intersect a primitive ! */
+ if(SXD_HIT_NONE(&hit)) continue;
+
+ /* Discard a hits if it lies on an edge/point */
+ if(HIT_ON_BOUNDARY(&hit, P, dirs[idir])) continue;
+
+ fX(normalize)(N, hit.normal);
+ cos_N_dir = fX(dot)(N, dirs[idir]);
+
+ /* Not too close and not roughly orthogonal */
+ if(hit.distance > 1.e-6 && absf(cos_N_dir) > 1.e-2f) {
+ const struct sdis_interface* interf;
+ interf = scene_get_interface(scn, hit.prim.prim_id);
+ medium = interface_get_medium
+ (interf, cos_N_dir < 0 ? SDIS_FRONT : SDIS_BACK);
+ break;
+ }
+ }
+ if(idir >= 2*DIM) {
+ res = XD(scene_get_medium)(scn, pos, NULL, &medium);
+ if(res != RES_OK) goto error;
+ }
+
+exit:
+ *out_medium = medium;
+ return res;
+error:
+#if DIM == 2
+ log_err(scn->dev, "%s: could not retrieve the medium at {%g, %g}.\n",
+ FUNC_NAME, SPLIT2(pos));
+#else
+ log_err(scn->dev, "%s: could not retrieve the medium at {%g, %g, %g}.\n",
+ FUNC_NAME, SPLIT3(pos));
+#endif
+ goto exit;
+}
+
#undef SDIS_SCENE_DIMENSION
#undef DIM
#undef sencXd
@@ -909,6 +1244,7 @@ error:
#undef SXD_PRIMITIVE_EQ
#undef fX
#undef fX_set_dX
+#undef fXX_mulfX
#undef XD
#undef HIT_ON_BOUNDARY
diff --git a/src/sdis_scene_c.h b/src/sdis_scene_c.h
@@ -31,6 +31,12 @@ struct prim_prop {
unsigned back_enclosure; /* Id of the back facing enclosure */
};
+struct hit_filter_data {
+ struct s2d_hit hit_2d;
+ struct s3d_hit hit_3d;
+ double epsilon; /* Threshold defining roughly equal intersections */
+};
+
struct get_medium_info {
/* Targeted position */
float pos_tgt[3];
@@ -226,6 +232,22 @@ scene_get_medium
struct get_medium_info* info, /* May be NULL */
struct sdis_medium** medium);
+/* This function assumes that the tested position lies into finite enclosure.
+ * The medium into which it lies is thus retrieved by tracing a random ray
+ * around the current position. For possible infinite enclosure, one has to use
+ * the `scene_get_medium' function instead that, in counterpart, can be more
+ * time consuming.
+ *
+ * Note that actually, the function internally calls scene_get_medium if no
+ * valid medium is found with the regular procedure. This may be due to
+ * numerical issues or wrong assumptions on the current medium (its boundaries
+ * are opened to infinity). */
+extern LOCAL_SYM res_T
+scene_get_medium_in_closed_boundaries
+ (const struct sdis_scene* scn,
+ const double position[],
+ struct sdis_medium** medium);
+
static INLINE void
scene_get_enclosure_ids
(const struct sdis_scene* scn,
diff --git a/src/test_sdis_conducto_radiative.c b/src/test_sdis_conducto_radiative.c
@@ -411,7 +411,7 @@ main(int argc, char** argv)
CHK(nfails + nreals == N);
CHK(nfails < N/1000);
- CHK(eq_eps(T.E, ref, 2*T.SE) == 1);
+ CHK(eq_eps(T.E, ref, 3*T.SE) == 1);
/* Check green function */
OK(sdis_solve_probe_green_function(scn, N, pos, 1, -1, Tref, &green));
diff --git a/src/test_sdis_solid_random_walk_robustness.c b/src/test_sdis_solid_random_walk_robustness.c
@@ -0,0 +1,391 @@
+/* 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 <star/s3dut.h>
+#include <rsys/math.h>
+
+#define Tfluid 350.0 /* Temperature of the fluid in Kelvin */
+#define LAMBDA 10.0 /* Thermal conductivity */
+#define Hcoef 1.0 /* Convection coefficient */
+#define Pw 10000.0 /* Volumetric power */
+#define Nreals 10000 /* #realisations */
+#define UNKNOWN -1
+
+/*******************************************************************************
+ * Helper functions
+ ******************************************************************************/
+static void
+compute_aabb
+ (const double* positions,
+ const size_t nvertices,
+ double lower[3],
+ double upper[3])
+{
+ size_t i;
+ CHK(positions && nvertices && lower && upper);
+
+ lower[0] = lower[1] = lower[2] = DBL_MAX;
+ upper[0] = upper[1] = upper[2] =-DBL_MAX;
+
+ FOR_EACH(i, 0, nvertices) {
+ lower[0] = MMIN(lower[0], positions[i*3+0]);
+ lower[1] = MMIN(lower[1], positions[i*3+1]);
+ lower[2] = MMIN(lower[2], positions[i*3+2]);
+ upper[0] = MMAX(upper[0], positions[i*3+0]);
+ upper[1] = MMAX(upper[1], positions[i*3+1]);
+ upper[2] = MMAX(upper[2], positions[i*3+2]);
+ }
+}
+
+static double
+trilinear_temperature(const double pos[3])
+{
+ const double a = 333;
+ const double b = 432;
+ const double c = 579;
+ double x, y, z;
+ CHK(pos[0] >= -10.0 && pos[0] <= 10.0);
+ CHK(pos[1] >= -10.0 && pos[1] <= 10.0);
+ CHK(pos[2] >= -10.0 && pos[2] <= 10.0);
+
+ x = (pos[0] + 10.0) / 20.0;
+ y = (pos[1] + 10.0) / 20.0;
+ z = (pos[2] + 10.0) / 20.0;
+
+ return a*x + b*y + c*z;
+}
+
+static double
+volumetric_temperature(const double pos[3], const double upper[3])
+{
+ const double beta = - 1.0 / 3.0 * Pw / (2*LAMBDA);
+ const double temp =
+ beta * (pos[0]*pos[0] - upper[0]*upper[0])
+ + beta * (pos[1]*pos[1] - upper[1]*upper[1])
+ + beta * (pos[2]*pos[2] - upper[2]*upper[2]);
+ CHK(temp > 0);
+ return temp;
+}
+
+static void
+print_estimation_result
+ (const struct sdis_estimator* estimator, const double Tref)
+{
+ struct sdis_mc T = SDIS_MC_NULL;
+ size_t nreals;
+ size_t nfails;
+
+ OK(sdis_estimator_get_temperature(estimator, &T));
+ OK(sdis_estimator_get_realisation_count(estimator, &nreals));
+ OK(sdis_estimator_get_failure_count(estimator, &nfails));
+ CHK(nfails <= Nreals * 0.0005);
+ printf("T = %g ~ %g +/- %g [%g, %g]; #failures = %lu / %lu\n",
+ Tref, T.E, T.SE, T.E - 3*T.SE, T.E + 3*T.SE,
+ (unsigned long)nfails, (unsigned long)Nreals);
+ CHK(eq_eps(T.E, Tref, T.SE*3));
+}
+
+/*******************************************************************************
+ * Geometry
+ ******************************************************************************/
+struct context {
+ struct s3dut_mesh_data msh;
+ struct sdis_interface* interf;
+};
+
+static void
+get_indices(const size_t itri, size_t ids[3], void* context)
+{
+ const struct context* ctx = context;
+ CHK(ids && ctx && itri < ctx->msh.nprimitives);
+ ids[0] = ctx->msh.indices[itri*3+0];
+ ids[2] = ctx->msh.indices[itri*3+1];
+ ids[1] = ctx->msh.indices[itri*3+2];
+}
+
+static void
+get_position(const size_t ivert, double pos[3], void* context)
+{
+ const struct context* ctx = context;
+ CHK(pos && ctx && ivert < ctx->msh.nvertices);
+ pos[0] = ctx->msh.positions[ivert*3+0];
+ pos[1] = ctx->msh.positions[ivert*3+1];
+ pos[2] = ctx->msh.positions[ivert*3+2];
+}
+
+static void
+get_interface(const size_t itri, struct sdis_interface** bound, void* context)
+{
+ const struct context* ctx = context;
+ CHK(bound && ctx && itri < ctx->msh.nprimitives);
+ *bound = ctx->interf;
+}
+
+/*******************************************************************************
+ * Interface
+ ******************************************************************************/
+enum profile {
+ PROFILE_UNKNOWN,
+ PROFILE_TRILINEAR,
+ PROFILE_VOLUMETRIC_POWER
+};
+struct interf {
+ enum profile profile;
+ double upper[3]; /* Upper bound of the scene */
+ double h;
+};
+
+static double
+interface_get_temperature
+ (const struct sdis_interface_fragment* frag, struct sdis_data* data)
+{
+ const struct interf* interf;
+ double temperature;
+ CHK(data != NULL && frag != NULL);
+ interf = sdis_data_cget(data);
+ switch(interf->profile) {
+ case PROFILE_UNKNOWN:
+ temperature = UNKNOWN;
+ break;
+ case PROFILE_VOLUMETRIC_POWER:
+ temperature = volumetric_temperature(frag->P, interf->upper);
+ break;
+ case PROFILE_TRILINEAR:
+ temperature = trilinear_temperature(frag->P);
+ break;
+ default: FATAL("Unreachable code.\n"); break;
+ }
+ return temperature;
+}
+
+static double
+interface_get_convection_coef
+ (const struct sdis_interface_fragment* frag, struct sdis_data* data)
+{
+ CHK(frag != NULL);
+ return ((const struct interf*)sdis_data_cget(data))->h;;
+}
+
+/*******************************************************************************
+ * Fluid
+ ******************************************************************************/
+static double
+fluid_get_temperature
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ CHK(vtx != NULL);
+ (void)data;
+ return Tfluid;
+}
+
+/*******************************************************************************
+ * Solid API
+ ******************************************************************************/
+struct solid {
+ double delta;
+ double cp;
+ double lambda;
+ double rho;
+ double temperature;
+ double power;
+};
+
+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;
+}
+
+static double
+solid_get_volumetric_power
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ CHK(data != NULL && vtx != NULL);
+ return ((const struct solid*)sdis_data_cget(data))->power;
+}
+
+/*******************************************************************************
+ * Main function
+ ******************************************************************************/
+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* msh = NULL;
+ struct sdis_device* dev = NULL;
+ struct sdis_data* data = NULL;
+ struct sdis_estimator* estimator = NULL;
+ struct sdis_medium* fluid = NULL;
+ struct sdis_medium* solid = NULL;
+ struct sdis_interface* interf = NULL;
+ struct sdis_scene* scn = NULL;
+ struct sdis_fluid_shader fluid_shader = DUMMY_FLUID_SHADER;
+ struct sdis_solid_shader solid_shader = SDIS_SOLID_SHADER_NULL;
+ struct sdis_interface_shader interf_shader = SDIS_INTERFACE_SHADER_NULL;
+ struct interf* interf_param = NULL;
+ struct solid* solid_param = NULL;
+ struct context ctx;
+ double probe[3];
+ double lower[3];
+ double upper[3];
+ double size[3];
+ const double time[2] = {DBL_MAX, DBL_MAX};
+ (void)argc, (void)argv;
+
+ OK(mem_init_proxy_allocator(&allocator, &mem_default_allocator));
+ OK(sdis_device_create(NULL, &allocator, SDIS_NTHREADS_DEFAULT, 1, &dev));
+
+ /* Create the fluid medium */
+ fluid_shader.temperature = fluid_get_temperature;
+ OK(sdis_fluid_create(dev, &fluid_shader, NULL, &fluid));
+
+ /* 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;
+ solid_shader.volumic_power = solid_get_volumetric_power;
+
+ /* Create the solid medium */
+ OK(sdis_data_create
+ (dev, sizeof(struct solid), ALIGNOF(struct solid), NULL, &data));
+ solid_param = sdis_data_get(data);
+ solid_param->delta = 0.01;
+ solid_param->lambda = 10;
+ solid_param->cp = 1.0;
+ solid_param->rho = 1.0;
+ solid_param->temperature = -1;
+ solid_param->power = SDIS_VOLUMIC_POWER_NONE;
+ OK(sdis_solid_create(dev, &solid_shader, data, &solid));
+ OK(sdis_data_ref_put(data));
+
+ /* Create the interface */
+ OK(sdis_data_create
+ (dev, sizeof(struct interf), ALIGNOF(struct interf), NULL, &data));
+ interf_param = sdis_data_get(data);
+ interf_param->h = 1;
+ interf_shader.convection_coef = interface_get_convection_coef;
+ interf_shader.front.temperature = interface_get_temperature;
+ OK(sdis_interface_create(dev, solid, fluid, &interf_shader, data, &interf));
+ OK(sdis_data_ref_put(data));
+
+ /* Create the solid super shape */
+ f0.A = 1; f0.B = 1; f0.M = 20; f0.N0 = 1; f0.N1 = 1; f0.N2 = 5;
+ f1.A = 1; f1.B = 1; f1.M = 7; f1.N0 = 1; f1.N1 = 2; f1.N2 = 5;
+ OK(s3dut_create_super_shape(&allocator, &f0, &f1, 1, 128, 64, &msh));
+ OK(s3dut_mesh_get_data(msh, &ctx.msh));
+
+ compute_aabb(ctx.msh.positions, ctx.msh.nvertices, lower, upper);
+
+ /* Create the scene */
+ ctx.interf = interf;
+ OK(sdis_scene_create(dev, ctx.msh.nprimitives, get_indices, get_interface,
+ ctx.msh.nvertices, get_position, &ctx, &scn));
+ /*dump_mesh(stdout, ctx.msh.positions,
+ ctx.msh.nvertices, ctx.msh.indices, ctx.msh.nprimitives);*/
+
+ /* Compute the delta of the solid random walk */
+ size[0] = upper[0] - lower[0];
+ size[1] = upper[1] - lower[1];
+ size[2] = upper[2] - lower[2];
+ solid_param->delta = MMIN(MMIN(size[0], size[1]), size[2]) / 20.0;
+
+ interf_param->upper[0] = upper[0];
+ interf_param->upper[1] = upper[1];
+ interf_param->upper[2] = upper[2];
+ interf_param->h = 1;
+
+ probe[0] = 0;
+ probe[1] = 0;
+ probe[2] = 0;
+
+ /* Launch probe estimation with trilinear profile set at interfaces */
+ interf_param->profile = PROFILE_TRILINEAR;
+ OK(sdis_solve_probe(scn, Nreals, probe, time, 1.0, -1, 0,
+ SDIS_HEAT_PATH_FAILED, &estimator));
+ print_estimation_result(estimator, trilinear_temperature(probe));
+ OK(sdis_estimator_ref_put(estimator));
+
+ /* Launch probe estimation with volumetric power profile set at interfaces */
+ interf_param->profile = PROFILE_VOLUMETRIC_POWER;
+ solid_param->power = Pw;
+ OK(sdis_solve_probe(scn, Nreals, probe, time, 1.0, -1, 0,
+ SDIS_HEAT_PATH_FAILED, &estimator));
+ print_estimation_result(estimator, volumetric_temperature(probe, upper));
+ solid_param->power = SDIS_VOLUMIC_POWER_NONE;
+ OK(sdis_estimator_ref_put(estimator));
+
+ /* Launch medium integration */
+ interf_param->profile = PROFILE_UNKNOWN;
+ OK(sdis_solve_medium(scn, Nreals, solid, time, 1.0, -1, 0,
+ SDIS_HEAT_PATH_FAILED, &estimator));
+ print_estimation_result(estimator, Tfluid);
+ /*dump_heat_paths(stdout, estimator);*/
+ OK(sdis_estimator_ref_put(estimator));
+
+ /* Release */
+ OK(s3dut_mesh_ref_put(msh));
+ OK(sdis_device_ref_put(dev));
+ OK(sdis_medium_ref_put(fluid));
+ OK(sdis_medium_ref_put(solid));
+ OK(sdis_interface_ref_put(interf));
+ 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_probe3.c b/src/test_sdis_solve_probe3.c
@@ -312,7 +312,7 @@ main(int argc, char** argv)
/* Check the results */
CHK(nfails + nreals == N);
CHK(nfails < N/1000);
- CHK(eq_eps(T.E, ref, 2*T.SE));
+ CHK(eq_eps(T.E, ref, 3*T.SE));
/* Check green function */
OK(sdis_solve_probe_green_function(scn, N, pos, 1.0, -1, 0, &green));
diff --git a/src/test_sdis_volumic_power4.c b/src/test_sdis_volumic_power4.c
@@ -0,0 +1,403 @@
+/* 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/clock_time.h>
+#include <rsys/math.h>
+
+#define Tf 100.0
+#define Power 10000.0
+#define H 50.0
+#define LAMBDA 100.0
+#define DELTA 0.4/*(1.0/2.0)*/
+#define N 100000
+#define LENGTH 10000.0
+
+/*
+ * The 2D scene is a solid slabs stretched along the X dimension to simulate a
+ * 1D case. The slab has a volumic power and has a convective exchange with
+ * surrounding fluid whose temperature is fixed to Tf.
+ *
+ *
+ * _\ Tf
+ * / /
+ * \__/
+ *
+ * ... -----Hboundary----- ...
+ *
+ * Lambda, Power
+ *
+ * ... -----Hboundary----- ...
+ *
+ * _\ Tf
+ * / /
+ * \__/
+ *
+ */
+
+static const double vertices_2d[4/*#vertices*/*2/*#coords per vertex*/] = {
+ LENGTH,-0.5,
+ -LENGTH,-0.5,
+ -LENGTH, 0.5,
+ LENGTH, 0.5
+};
+
+static const double vertices_3d[8/*#vertices*/*3/*#coords per vertex*/] = {
+ -LENGTH,-0.5,-LENGTH,
+ LENGTH,-0.5,-LENGTH,
+ -LENGTH, 0.5,-LENGTH,
+ LENGTH, 0.5,-LENGTH,
+ -LENGTH,-0.5, LENGTH,
+ LENGTH,-0.5, LENGTH,
+ -LENGTH, 0.5, LENGTH,
+ LENGTH, 0.5, LENGTH
+};
+
+/*******************************************************************************
+ * Geometry
+ ******************************************************************************/
+static void
+get_position_2d(const size_t ivert, double pos[2], void* context)
+{
+ (void)context;
+ CHK(pos);
+ pos[0] = vertices_2d[ivert*2+0];
+ pos[1] = vertices_2d[ivert*2+1];
+}
+
+static void
+get_position_3d(const size_t ivert, double pos[3], void* context)
+{
+ (void)context;
+ CHK(pos);
+ pos[0] = vertices_3d[ivert*3+0];
+ pos[1] = vertices_3d[ivert*3+1];
+ pos[2] = vertices_3d[ivert*3+2];
+}
+
+static void
+get_interface(const size_t iprim, struct sdis_interface** bound, void* context)
+{
+ struct sdis_interface** interfaces = context;
+ CHK(context && bound);
+ *bound = interfaces[iprim];
+}
+
+/*******************************************************************************
+ * Solid medium
+ ******************************************************************************/
+struct solid {
+ double cp;
+ double lambda;
+ double rho;
+ double delta;
+ double volumic_power;
+ 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;
+}
+
+static double
+solid_get_volumic_power
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ CHK(data != NULL && vtx != NULL);
+ return ((const struct solid*)sdis_data_cget(data))->volumic_power;
+}
+
+/*******************************************************************************
+ * Fluid medium
+ ******************************************************************************/
+struct fluid {
+ double temperature;
+};
+
+static double
+fluid_get_temperature
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ const struct fluid* fluid;
+ CHK(data != NULL && vtx != NULL);
+ fluid = sdis_data_cget(data);
+ return fluid->temperature;
+}
+
+/*******************************************************************************
+ * Interfaces
+ ******************************************************************************/
+struct interf {
+ double h;
+ double temperature;
+};
+
+static double
+interface_get_convection_coef
+ (const struct sdis_interface_fragment* frag, struct sdis_data* data)
+{
+ CHK(frag && data);
+ return ((const struct interf*)sdis_data_cget(data))->h;
+}
+
+static double
+interface_get_temperature
+ (const struct sdis_interface_fragment* frag, struct sdis_data* data)
+{
+ CHK(frag && data);
+ return ((const struct interf*)sdis_data_cget(data))->temperature;
+}
+
+/*******************************************************************************
+ * Test
+ ******************************************************************************/
+int
+main(int argc, char** argv)
+{
+ char dump[128];
+ struct time t0, t1;
+ struct mem_allocator allocator;
+ struct solid* solid_param = NULL;
+ struct fluid* fluid_param = NULL;
+ struct interf* interf_param = NULL;
+ struct sdis_device* dev = NULL;
+ struct sdis_data* data = NULL;
+ struct sdis_medium* fluid1 = NULL;
+ struct sdis_medium* fluid2 = NULL;
+ struct sdis_medium* solid = NULL;
+ struct sdis_scene* scn_2d = NULL;
+ struct sdis_scene* scn_3d = NULL;
+ struct sdis_estimator* estimator = NULL;
+ struct sdis_fluid_shader fluid_shader = SDIS_FLUID_SHADER_NULL;
+ struct sdis_solid_shader solid_shader = SDIS_SOLID_SHADER_NULL;
+ struct sdis_interface_shader interf_shader = SDIS_INTERFACE_SHADER_NULL;
+ struct sdis_interface* interf_adiabatic = NULL;
+ struct sdis_interface* interf_solid_fluid1 = NULL;
+ struct sdis_interface* interf_solid_fluid2 = NULL;
+ struct sdis_interface* interfaces[12/*#max primitives*/];
+ struct sdis_mc T = SDIS_MC_NULL;
+ size_t nreals, nfails;
+ double pos[3];
+ double time_range[2] = { INF, INF };
+ double Tref;
+ double x;
+ (void)argc, (void)argv;
+
+ OK(mem_init_proxy_allocator(&allocator, &mem_default_allocator));
+ OK(sdis_device_create(NULL, &allocator, SDIS_NTHREADS_DEFAULT, 1, &dev));
+
+ /* Create the fluid medium */
+ fluid_shader.temperature = fluid_get_temperature;
+ fluid_shader.calorific_capacity = dummy_medium_getter;
+ fluid_shader.volumic_mass = dummy_medium_getter;
+
+ OK(sdis_data_create
+ (dev, sizeof(struct fluid), ALIGNOF(struct fluid), NULL, &data));
+ fluid_param = sdis_data_get(data);
+ fluid_param->temperature = Tf;
+ OK(sdis_fluid_create(dev, &fluid_shader, data, &fluid1));
+ OK(sdis_data_ref_put(data));
+
+ OK(sdis_data_create
+ (dev, sizeof(struct fluid), ALIGNOF(struct fluid), NULL, &data));
+ fluid_param = sdis_data_get(data);
+ fluid_param->temperature = Tf;
+ OK(sdis_fluid_create(dev, &fluid_shader, data, &fluid2));
+ 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;
+ solid_shader.volumic_power = solid_get_volumic_power;
+
+ /* Create the solid medium */
+ OK(sdis_data_create
+ (dev, sizeof(struct solid), ALIGNOF(struct solid), NULL, &data));
+ solid_param = sdis_data_get(data);
+ solid_param->cp = 500000;
+ solid_param->rho = 1000;
+ solid_param->lambda = LAMBDA;
+ solid_param->delta = DELTA;
+ solid_param->volumic_power = Power;
+ solid_param->temperature = -1;
+ OK(sdis_solid_create(dev, &solid_shader, data, &solid));
+ OK(sdis_data_ref_put(data));
+
+ /* Setup the interface shader */
+ interf_shader.convection_coef = interface_get_convection_coef;
+ interf_shader.front.temperature = interface_get_temperature;
+
+ /* Create the adiabatic interface */
+ OK(sdis_data_create (dev, sizeof(struct interf), ALIGNOF(struct interf),
+ NULL, &data));
+ interf_param = sdis_data_get(data);
+ interf_param->h = 0;
+ interf_param->temperature = -1;
+ OK(sdis_interface_create(dev, solid, fluid1, &interf_shader, data,
+ &interf_adiabatic));
+ OK(sdis_data_ref_put(data));
+
+ /* Create the solid fluid1 interface */
+ OK(sdis_data_create (dev, sizeof(struct interf), ALIGNOF(struct interf),
+ NULL, &data));
+ interf_param = sdis_data_get(data);
+ interf_param->h = H;
+ interf_param->temperature = -1;
+ OK(sdis_interface_create(dev, solid, fluid1, &interf_shader, data,
+ &interf_solid_fluid1));
+ OK(sdis_data_ref_put(data));
+
+ /* Create the solid fluid2 interface */
+ OK(sdis_data_create (dev, sizeof(struct interf), ALIGNOF(struct interf),
+ NULL, &data));
+ interf_param = sdis_data_get(data);
+ interf_param->h = H;
+ interf_param->temperature = -1;
+ OK(sdis_interface_create(dev, solid, fluid2, &interf_shader, data,
+ &interf_solid_fluid2));
+ OK(sdis_data_ref_put(data));
+
+ /* Release the media */
+ OK(sdis_medium_ref_put(fluid1));
+ OK(sdis_medium_ref_put(fluid2));
+ OK(sdis_medium_ref_put(solid));
+
+#if 0
+ dump_segments(stdout, vertices, nvertices, indices, nsegments);
+ exit(0);
+#endif
+
+ /* Map the interfaces to their square segments */
+ interfaces[0] = interf_solid_fluid2; /* Bottom */
+ interfaces[1] = interf_adiabatic; /* Left */
+ interfaces[2] = interf_solid_fluid1; /* Top */
+ interfaces[3] = interf_adiabatic; /* Right */
+
+ /* Create the 2D scene */
+ OK(sdis_scene_2d_create(dev, square_nsegments, square_get_indices, get_interface,
+ square_nvertices, get_position_2d, interfaces, &scn_2d));
+
+ /* Map the interfaces to their box triangles */
+ interfaces[0] = interfaces[1] = interf_adiabatic; /* Front */
+ interfaces[2] = interfaces[3] = interf_adiabatic; /* Left */
+ interfaces[4] = interfaces[5] = interf_adiabatic; /* Back */
+ interfaces[6] = interfaces[7] = interf_adiabatic; /* Right */
+ interfaces[8] = interfaces[9] = interf_solid_fluid1; /* Top */
+ interfaces[10]= interfaces[11]= interf_solid_fluid2; /* Bottom */
+
+ /* Create the 3D scene */
+ OK(sdis_scene_create(dev, box_ntriangles, box_get_indices, get_interface,
+ box_nvertices, get_position_3d, interfaces, &scn_3d));
+
+ /* Release the interfaces */
+ OK(sdis_interface_ref_put(interf_adiabatic));
+ OK(sdis_interface_ref_put(interf_solid_fluid1));
+ OK(sdis_interface_ref_put(interf_solid_fluid2));
+
+ pos[0] = 0;
+ pos[1] = 0;
+ pos[2] = 0;
+
+ x = pos[1];
+ Tref = -Power / (2*LAMBDA) * x*x + Tf + Power/(2*H) + Power/(8*LAMBDA);
+
+ printf(">>> 2D\n");
+
+ time_current(&t0);
+ OK(sdis_solve_probe(scn_2d, N, pos, time_range, 1.f, -1, 0, 0, &estimator));
+ time_sub(&t0, time_current(&t1), &t0);
+ time_dump(&t0, TIME_ALL, NULL, dump, sizeof(dump));
+ printf("Elapsed time = %s\n", dump);
+
+ OK(sdis_estimator_get_temperature(estimator, &T));
+ OK(sdis_estimator_get_realisation_count(estimator, &nreals));
+ OK(sdis_estimator_get_failure_count(estimator, &nfails));
+ printf("Temperature at (%g %g) = %g ~ %g +/- %g [%g %g]\n",
+ SPLIT2(pos), Tref, T.E, T.SE, T.E-3*T.SE, T.E+3*T.SE);
+ printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N);
+ OK(sdis_estimator_ref_put(estimator));
+ CHK(nfails + nreals == N);
+ CHK(nfails < N/1000);
+ CHK(eq_eps(T.E, Tref, T.SE*3));
+
+ printf("\n>>> 3D\n");
+
+ time_current(&t0);
+ OK(sdis_solve_probe(scn_3d, N, pos, time_range, 1.f, -1, 0, 0, &estimator));
+ time_sub(&t0, time_current(&t1), &t0);
+ time_dump(&t0, TIME_ALL, NULL, dump, sizeof(dump));
+ printf("Elapsed time = %s\n", dump);
+
+ OK(sdis_estimator_get_temperature(estimator, &T));
+ OK(sdis_estimator_get_realisation_count(estimator, &nreals));
+ OK(sdis_estimator_get_failure_count(estimator, &nfails));
+ printf("Temperature at (%g %g) = %g ~ %g +/- %g [%g %g]\n",
+ SPLIT2(pos), Tref, T.E, T.SE, T.E-3*T.SE, T.E+3*T.SE);
+ printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N);
+ OK(sdis_estimator_ref_put(estimator));
+ CHK(nfails + nreals == N);
+ CHK(nfails < N/1000);
+ CHK(eq_eps(T.E, Tref, T.SE*3));
+
+ OK(sdis_scene_ref_put(scn_2d));
+ OK(sdis_scene_ref_put(scn_3d));
+ OK(sdis_device_ref_put(dev));
+
+ check_memory_allocator(&allocator);
+ mem_shutdown_proxy_allocator(&allocator);
+ CHK(mem_allocated_size() == 0);
+ return 0;
+}
+
diff --git a/src/test_sdis_volumic_power4_2d.c b/src/test_sdis_volumic_power4_2d.c
@@ -1,363 +0,0 @@
-/* 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/clock_time.h>
-#include <rsys/math.h>
-
-#define Tf 100.0
-#define Power 10000.0
-#define H 50.0
-#define LAMBDA 100.0
-#define DELTA (1.0/2.0)
-#define N 10000
-
-/*
- * The 2D scene is a solid slabs stretched along the X dimension to simulate a
- * 1D case. The slab has a volumic power and has a convective exchange with
- * surrounding fluid whose temperature is fixed to Tf.
- *
- *
- * _\ Tf
- * / /
- * \__/
- *
- * ... -----Hboundary----- ...
- *
- * Lambda, Power
- *
- * ... -----Hboundary----- ...
- *
- * _\ Tf
- * / /
- * \__/
- *
- */
-
-static const double vertices[4/*#vertices*/*2/*#coords per vertex*/] = {
- -1000000.5,-0.5,
- -1000000.5, 0.5,
- 1000000.5, 0.5,
- 1000000.5,-0.5
-};
-static const size_t nvertices = sizeof(vertices)/sizeof(double[2]);
-
-static const size_t indices[4/*#segments*/*2/*#indices per segment*/]= {
- 0, 1,
- 1, 2,
- 2, 3,
- 3, 0
-};
-static const size_t nsegments = sizeof(indices)/sizeof(size_t[2]);
-
-/*******************************************************************************
- * Geometry
- ******************************************************************************/
-static void
-get_indices(const size_t iseg, size_t ids[2], void* context)
-{
- (void)context;
- CHK(ids);
- ids[0] = indices[iseg*2+0];
- ids[1] = indices[iseg*2+1];
-}
-
-static void
-get_position(const size_t ivert, double pos[2], void* context)
-{
- (void)context;
- CHK(pos);
- pos[0] = vertices[ivert*2+0];
- pos[1] = vertices[ivert*2+1];
-}
-
-static void
-get_interface(const size_t iseg, struct sdis_interface** bound, void* context)
-{
- struct sdis_interface** interfaces = context;
- CHK(context && bound);
- *bound = interfaces[iseg];
-}
-
-/*******************************************************************************
- * Solid medium
- ******************************************************************************/
-struct solid {
- double cp;
- double lambda;
- double rho;
- double delta;
- double volumic_power;
- 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;
-}
-
-static double
-solid_get_volumic_power
- (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
-{
- CHK(data != NULL && vtx != NULL);
- return ((const struct solid*)sdis_data_cget(data))->volumic_power;
-}
-
-/*******************************************************************************
- * Fluid medium
- ******************************************************************************/
-struct fluid {
- double temperature;
-};
-
-static double
-fluid_get_temperature
- (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
-{
- const struct fluid* fluid;
- CHK(data != NULL && vtx != NULL);
- fluid = sdis_data_cget(data);
- return fluid->temperature;
-}
-
-/*******************************************************************************
- * Interfaces
- ******************************************************************************/
-struct interf {
- double h;
- double temperature;
-};
-
-static double
-interface_get_convection_coef
- (const struct sdis_interface_fragment* frag, struct sdis_data* data)
-{
- CHK(frag && data);
- return ((const struct interf*)sdis_data_cget(data))->h;
-}
-
-static double
-interface_get_temperature
- (const struct sdis_interface_fragment* frag, struct sdis_data* data)
-{
- CHK(frag && data);
- return ((const struct interf*)sdis_data_cget(data))->temperature;
-}
-
-/*******************************************************************************
- * Test
- ******************************************************************************/
-int
-main(int argc, char** argv)
-{
- char dump[128];
- struct time t0, t1;
- struct mem_allocator allocator;
- struct solid* solid_param = NULL;
- struct fluid* fluid_param = NULL;
- struct interf* interf_param = NULL;
- struct sdis_device* dev = NULL;
- struct sdis_data* data = NULL;
- struct sdis_medium* fluid1 = NULL;
- struct sdis_medium* fluid2 = NULL;
- struct sdis_medium* solid = NULL;
- struct sdis_scene* scn = NULL;
- struct sdis_estimator* estimator = NULL;
- struct sdis_fluid_shader fluid_shader = SDIS_FLUID_SHADER_NULL;
- struct sdis_solid_shader solid_shader = SDIS_SOLID_SHADER_NULL;
- struct sdis_interface_shader interf_shader = SDIS_INTERFACE_SHADER_NULL;
- struct sdis_interface* interf_adiabatic = NULL;
- struct sdis_interface* interf_solid_fluid1 = NULL;
- struct sdis_interface* interf_solid_fluid2 = NULL;
- struct sdis_interface* interfaces[4/*#segment*/];
- struct sdis_mc T = SDIS_MC_NULL;
- size_t nreals, nfails;
- double pos[2];
- double time_range[2] = { INF, INF };
- double Tref;
- double x;
- (void)argc, (void)argv;
-
- OK(mem_init_proxy_allocator(&allocator, &mem_default_allocator));
- OK(sdis_device_create(NULL, &allocator, SDIS_NTHREADS_DEFAULT, 1, &dev));
-
- /* Create the fluid medium */
- fluid_shader.temperature = fluid_get_temperature;
- fluid_shader.calorific_capacity = dummy_medium_getter;
- fluid_shader.volumic_mass = dummy_medium_getter;
-
- OK(sdis_data_create
- (dev, sizeof(struct fluid), ALIGNOF(struct fluid), NULL, &data));
- fluid_param = sdis_data_get(data);
- fluid_param->temperature = Tf;
- OK(sdis_fluid_create(dev, &fluid_shader, data, &fluid1));
- OK(sdis_data_ref_put(data));
-
- OK(sdis_data_create
- (dev, sizeof(struct fluid), ALIGNOF(struct fluid), NULL, &data));
- fluid_param = sdis_data_get(data);
- fluid_param->temperature = Tf;
- OK(sdis_fluid_create(dev, &fluid_shader, data, &fluid2));
- 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;
- solid_shader.volumic_power = solid_get_volumic_power;
-
- /* Create the solid medium */
- OK(sdis_data_create
- (dev, sizeof(struct solid), ALIGNOF(struct solid), NULL, &data));
- solid_param = sdis_data_get(data);
- solid_param->cp = 500000;
- solid_param->rho = 1000;
- solid_param->lambda = LAMBDA;
- solid_param->delta = DELTA;
- solid_param->volumic_power = Power;
- solid_param->temperature = -1;
- OK(sdis_solid_create(dev, &solid_shader, data, &solid));
- OK(sdis_data_ref_put(data));
-
- /* Setup the interface shader */
- interf_shader.convection_coef = interface_get_convection_coef;
- interf_shader.front.temperature = interface_get_temperature;
-
- /* Create the adiabatic interface */
- OK(sdis_data_create (dev, sizeof(struct interf), ALIGNOF(struct interf),
- NULL, &data));
- interf_param = sdis_data_get(data);
- interf_param->h = 0;
- interf_param->temperature = -1;
- OK(sdis_interface_create(dev, solid, fluid1, &interf_shader, data,
- &interf_adiabatic));
- OK(sdis_data_ref_put(data));
-
- /* Create the solid fluid1 interface */
- OK(sdis_data_create (dev, sizeof(struct interf), ALIGNOF(struct interf),
- NULL, &data));
- interf_param = sdis_data_get(data);
- interf_param->h = H;
- interf_param->temperature = -1;
- OK(sdis_interface_create(dev, solid, fluid1, &interf_shader, data,
- &interf_solid_fluid1));
- OK(sdis_data_ref_put(data));
-
- /* Create the solid fluid2 interface */
- OK(sdis_data_create (dev, sizeof(struct interf), ALIGNOF(struct interf),
- NULL, &data));
- interf_param = sdis_data_get(data);
- interf_param->h = H;
- interf_param->temperature = -1;
- OK(sdis_interface_create(dev, solid, fluid2, &interf_shader, data,
- &interf_solid_fluid2));
- OK(sdis_data_ref_put(data));
-
- /* Release the media */
- OK(sdis_medium_ref_put(fluid1));
- OK(sdis_medium_ref_put(fluid2));
- OK(sdis_medium_ref_put(solid));
-
- /* Map the interfaces to their square segments */
- interfaces[0] = interf_adiabatic;
- interfaces[1] = interf_solid_fluid1;
- interfaces[2] = interf_adiabatic;
- interfaces[3] = interf_solid_fluid2;
-
-#if 0
- dump_segments(stdout, vertices, nvertices, indices, nsegments);
- exit(0);
-#endif
-
- /* Create the scene */
- OK(sdis_scene_2d_create(dev, nsegments, get_indices, get_interface,
- nvertices, get_position, interfaces, &scn));
-
- /* Release the interfaces */
- OK(sdis_interface_ref_put(interf_adiabatic));
- OK(sdis_interface_ref_put(interf_solid_fluid1));
- OK(sdis_interface_ref_put(interf_solid_fluid2));
-
- pos[0] = 0;
- pos[1] = 0.25;
-
- x = pos[1];
- Tref = -Power / (2*LAMBDA) * x*x + Tf + Power/(2*H) + Power/(8*LAMBDA);
-
- time_current(&t0);
- OK(sdis_solve_probe(scn, N, pos, time_range, 1.f, -1, 0, 0, &estimator));
- time_sub(&t0, time_current(&t1), &t0);
- time_dump(&t0, TIME_ALL, NULL, dump, sizeof(dump));
- printf("Elapsed time = %s\n", dump);
-
- OK(sdis_estimator_get_temperature(estimator, &T));
- OK(sdis_estimator_get_realisation_count(estimator, &nreals));
- OK(sdis_estimator_get_failure_count(estimator, &nfails));
- printf("Temperature at (%g %g) = %g ~ %g +/- %g [%g %g]\n",
- SPLIT2(pos), Tref, T.E, T.SE, T.E-3*T.SE, T.E+3*T.SE);
- printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N);
- OK(sdis_estimator_ref_put(estimator));
- CHK(nfails + nreals == N);
- CHK(nfails < N/1000);
- CHK(eq_eps(T.E, Tref, T.SE*3));
-
- OK(sdis_scene_ref_put(scn));
- OK(sdis_device_ref_put(dev));
-
- check_memory_allocator(&allocator);
- mem_shutdown_proxy_allocator(&allocator);
- CHK(mem_allocated_size() == 0);
- return 0;
-}
-
