commit c79bea63847a5caf927c3e7251b0183a0a444316
parent 34a3152c88013b25186945e52a8678cb41bc5b7c
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Tue, 4 Jun 2024 15:43:41 +0200
Merge branch 'feature_track_enclosures' into develop
Diffstat:
28 files changed, 1731 insertions(+), 1540 deletions(-)
diff --git a/src/sdis_Xd_begin.h b/src/sdis_Xd_begin.h
@@ -76,44 +76,3 @@
/* Macro making generic its submitted name to SDIS_XD_DIMENSION */
#define XD(Name) CONCAT(CONCAT(CONCAT(Name, _), DIM), d)
-
-/* Generate the generic data structures and constants */
-#if (SDIS_XD_DIMENSION == 2 && !defined(SDIS_2D_H)) \
-|| (SDIS_XD_DIMENSION == 3 && !defined(SDIS_3D_H))
- #if SDIS_XD_DIMENSION == 2
- #define SDIS_2D_H
- #else
- #define SDIS_3D_H
- #endif
-
-struct rwalk_context;
-
-/* Current state of the random walk */
-struct XD(rwalk) {
- struct sdis_rwalk_vertex vtx; /* Position and time of the Random walk */
- struct sdis_medium* mdm; /* Medium in which the random walk lies */
- struct sXd(hit) hit; /* Hit of the random walk */
-
- /* Direction along which the random walk reached the radiative environment */
- double dir[3];
-
- double elapsed_time;
- enum sdis_side hit_side;
-};
-static const struct XD(rwalk) XD(RWALK_NULL) = {
- SDIS_RWALK_VERTEX_NULL__, NULL, SXD_HIT_NULL__, {0,0,0}, 0, SDIS_SIDE_NULL__
-};
-
-struct XD(temperature) {
- res_T (*func)/* Next function to invoke in order to compute the temperature */
- (struct sdis_scene* scn,
- struct rwalk_context* ctx,
- struct XD(rwalk)* rwalk,
- struct ssp_rng* rng,
- struct XD(temperature)* temp);
- double value; /* Current value of the temperature */
- int done;
-};
-static const struct XD(temperature) XD(TEMPERATURE_NULL) = { NULL, 0, 0 };
-
-#endif /* SDIX_<2|3>D_H */
diff --git a/src/sdis_heat_path.h b/src/sdis_heat_path.h
@@ -17,6 +17,7 @@
#define SDIS_HEAT_PATH_H
#include "sdis.h"
+#include "sdis_scene_c.h"
#include <rsys/dynamic_array.h>
#include <rsys/dynamic_array_size_t.h>
@@ -24,12 +25,8 @@
/* Forward declarations */
struct green_path_handle;
-struct rwalk_2d;
-struct rwalk_3d;
struct sdis_scene;
struct ssp_rng;
-struct temperature_2d;
-struct temperature_3d;
/*******************************************************************************
* Context of a random walk, i.e. its data concerning the current system and the
@@ -87,6 +84,46 @@ get_picard_order(const struct rwalk_context* ctx)
}
/*******************************************************************************
+ * 2D/3D random walk and associated temperature, i.e. current state of the
+ * sampled path
+ ******************************************************************************/
+struct rwalk {
+ struct sdis_rwalk_vertex vtx; /* Position and time of the Random walk */
+ unsigned enc_id; /* Id of the enclosure in which the random walk lies */
+ struct s2d_hit hit_2d;
+ struct s3d_hit hit_3d;
+
+ /* Direction along which the random walk reached the radiative environment */
+ double dir[3];
+
+ double elapsed_time;
+ enum sdis_side hit_side;
+};
+#define RWALK_NULL__ { \
+ SDIS_RWALK_VERTEX_NULL__, \
+ ENCLOSURE_ID_NULL, \
+ S2D_HIT_NULL__, \
+ S3D_HIT_NULL__, \
+ {0,0,0}, \
+ 0, \
+ SDIS_SIDE_NULL__ \
+}
+static const struct rwalk RWALK_NULL = RWALK_NULL__;
+
+struct temperature {
+ res_T (*func)/* Next function to invoke in order to compute the temperature */
+ (struct sdis_scene* scn,
+ struct rwalk_context* ctx,
+ struct rwalk* rwalk,
+ struct ssp_rng* rng,
+ struct temperature* temp);
+ double value; /* Current value of the temperature */
+ int done;
+};
+#define TEMPERATURE_NULL__ {NULL,0,0}
+static const struct temperature TEMPERATURE_NULL = TEMPERATURE_NULL__;
+
+/*******************************************************************************
* Heat path data structure used to record the geometry of sampled paths
******************************************************************************/
/* Generate the dynamic array of heat vertices */
@@ -261,34 +298,34 @@ trace_radiative_path_2d
(struct sdis_scene* scn,
const float ray_dir[3],
struct rwalk_context* ctx,
- struct rwalk_2d* rwalk,
+ struct rwalk* rwalk,
struct ssp_rng* rng,
- struct temperature_2d* temperature);
+ struct temperature* temperature);
extern LOCAL_SYM res_T
trace_radiative_path_3d
(struct sdis_scene* scn,
const float ray_dir[3],
struct rwalk_context* ctx,
- struct rwalk_3d* rwalk,
+ struct rwalk* rwalk,
struct ssp_rng* rng,
- struct temperature_3d* temperature);
+ struct temperature* temperature);
extern LOCAL_SYM res_T
radiative_path_2d
(struct sdis_scene* scn,
struct rwalk_context* ctx,
- struct rwalk_2d* rwalk,
+ struct rwalk* rwalk,
struct ssp_rng* rng,
- struct temperature_2d* temperature);
+ struct temperature* temperature);
extern LOCAL_SYM res_T
radiative_path_3d
(struct sdis_scene* scn,
struct rwalk_context* ctx,
- struct rwalk_3d* rwalk,
+ struct rwalk* rwalk,
struct ssp_rng* rng,
- struct temperature_3d* temperature);
+ struct temperature* temperature);
/*******************************************************************************
* Convective path
@@ -297,17 +334,17 @@ extern LOCAL_SYM res_T
convective_path_2d
(struct sdis_scene* scn,
struct rwalk_context* ctx,
- struct rwalk_2d* rwalk,
+ struct rwalk* rwalk,
struct ssp_rng* rng,
- struct temperature_2d* temperature);
+ struct temperature* temperature);
extern LOCAL_SYM res_T
convective_path_3d
(struct sdis_scene* scn,
struct rwalk_context* ctx,
- struct rwalk_3d* rwalk,
+ struct rwalk* rwalk,
struct ssp_rng* rng,
- struct temperature_3d* temperature);
+ struct temperature* temperature);
/*******************************************************************************
* Conductive path
@@ -316,17 +353,17 @@ extern LOCAL_SYM res_T
conductive_path_2d
(struct sdis_scene* scn,
struct rwalk_context* ctx,
- struct rwalk_2d* rwalk,
+ struct rwalk* rwalk,
struct ssp_rng* rng,
- struct temperature_2d* temperature);
+ struct temperature* temperature);
extern LOCAL_SYM res_T
conductive_path_3d
(struct sdis_scene* scn,
struct rwalk_context* ctx,
- struct rwalk_3d* rwalk,
+ struct rwalk* rwalk,
struct ssp_rng* rng,
- struct temperature_3d* temperature);
+ struct temperature* temperature);
/*******************************************************************************
* Boundary sub-path
@@ -335,16 +372,16 @@ extern LOCAL_SYM res_T
boundary_path_2d
(struct sdis_scene* scn,
struct rwalk_context* ctx,
- struct rwalk_2d* rwalk,
+ struct rwalk* rwalk,
struct ssp_rng* rng,
- struct temperature_2d* temperature);
+ struct temperature* temperature);
extern LOCAL_SYM res_T
boundary_path_3d
(struct sdis_scene* scn,
struct rwalk_context* ctx,
- struct rwalk_3d* rwalk,
+ struct rwalk* rwalk,
struct ssp_rng* rng,
- struct temperature_3d* temperature);
+ struct temperature* temperature);
#endif /* SDIS_HEAT_PATH_H */
diff --git a/src/sdis_heat_path_boundary_Xd.h b/src/sdis_heat_path_boundary_Xd.h
@@ -25,15 +25,67 @@
#include "sdis_Xd_begin.h"
/*******************************************************************************
+ * Helper functions
+ ******************************************************************************/
+/* This function checks whether the random walk is on a boundary and, if so,
+ * verifies that the temperature of the medium attached to the interface is
+ * known. This medium can be different from the medium of the enclosure. Indeed,
+ * the enclosure can contain several media used to set the temperatures of
+ * several boundary conditions. Hence this function, which queries the medium on
+ * the trajectory coming from a boundary */
+static res_T
+XD(handle_known_medium_temperature)
+ (struct sdis_scene* scn,
+ struct rwalk_context* ctx,
+ struct rwalk* rwalk,
+ struct temperature* T)
+{
+ struct sdis_interface* interf = NULL;
+ struct sdis_medium* mdm = NULL;
+ double temperature = SDIS_TEMPERATURE_NONE;
+ res_T res = RES_OK;
+ ASSERT(scn && ctx && rwalk && T);
+
+ /* Not at an interface */
+ if(SXD_HIT_NONE(&rwalk->XD(hit))) return RES_OK; /* Nothing to do */
+
+ interf = scene_get_interface(scn, rwalk->XD(hit).prim.prim_id);
+ mdm = rwalk->hit_side==SDIS_FRONT ? interf->medium_front: interf->medium_back;
+
+ temperature = medium_get_temperature(mdm, &rwalk->vtx);
+
+ /* Check if the temperature is known */
+ if(SDIS_TEMPERATURE_IS_UNKNOWN(temperature)) goto exit;
+
+ T->value += temperature;
+ T->done = 1;
+
+ if(ctx->green_path) {
+ res = green_path_set_limit_vertex
+ (ctx->green_path, mdm, &rwalk->vtx, rwalk->elapsed_time);
+ if(res != RES_OK) goto error;
+ }
+
+ if(ctx->heat_path) {
+ heat_path_get_last_vertex(ctx->heat_path)->weight = T->value;
+ }
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
+/*******************************************************************************
* Local functions
******************************************************************************/
res_T
XD(boundary_path)
(struct sdis_scene* scn,
struct rwalk_context* ctx,
- struct XD(rwalk)* rwalk,
+ struct rwalk* rwalk,
struct ssp_rng* rng,
- struct XD(temperature)* T)
+ struct temperature* T)
{
struct sdis_interface_fragment frag = SDIS_INTERFACE_FRAGMENT_NULL;
struct sdis_interface* interf = NULL;
@@ -42,15 +94,16 @@ XD(boundary_path)
double tmp;
res_T res = RES_OK;
ASSERT(scn && ctx && rwalk && rng && T);
- ASSERT(rwalk->mdm == NULL);
- ASSERT(!SXD_HIT_NONE(&rwalk->hit));
+ ASSERT(rwalk->enc_id == ENCLOSURE_ID_NULL);
+ ASSERT(!SXD_HIT_NONE(&rwalk->XD(hit)));
- XD(setup_interface_fragment)(&frag, &rwalk->vtx, &rwalk->hit, rwalk->hit_side);
+ XD(setup_interface_fragment)
+ (&frag, &rwalk->vtx, &rwalk->XD(hit), rwalk->hit_side);
- fX(normalize)(rwalk->hit.normal, rwalk->hit.normal);
+ fX(normalize)(rwalk->XD(hit).normal, rwalk->XD(hit).normal);
/* Retrieve the current interface */
- interf = scene_get_interface(scn, rwalk->hit.prim.prim_id);
+ interf = scene_get_interface(scn, rwalk->XD(hit).prim.prim_id);
/* Check if the boundary temperature is known */
tmp = interface_side_get_temperature(interf, &frag);
@@ -82,6 +135,32 @@ XD(boundary_path)
}
if(res != RES_OK) goto error;
+#if 1
+ if(T->done) goto exit;
+
+ /* Handling limit boundary condition, i.e. the trajectory originates from a
+ * boundary and the medium temperature is known (e.g. Robin's condition). To
+ * simplify data description, we allow in such a situation, to define several
+ * medium on the same enclosure, each with a fixed temperature, i.e. different
+ * conditions are defined for the different interfaces that detour the
+ * enclosure. As a result, no path can be sampled in this enclosure, which is
+ * beyond the system boundary. The boundary medium must therefore be
+ * interrogated from the interface.
+ *
+ * Note that we check this boundary condition with convective paths to handle
+ * Robin's boundary conditions. But we also make this check when passing
+ * through conduction when there's no reason why a solid should have a fixed
+ * temperature and not its boundary: it should be a Dirichlet condition.
+ * Although it's not physical, such systems can still be defined
+ * computationally, and in fact it's also a handy way of testing
+ * border cases */
+ if(T->func == XD(convective_path) || T->func == XD(conductive_path)) {
+ res = XD(handle_known_medium_temperature)(scn, ctx, rwalk, T);
+ if(res != RES_OK) goto error;
+ if(T->done) goto exit; /* That's all folks */
+ }
+#endif
+
exit:
return res;
error:
@@ -89,4 +168,3 @@ error:
}
#include "sdis_Xd_end.h"
-
diff --git a/src/sdis_heat_path_boundary_Xd_c.h b/src/sdis_heat_path_boundary_Xd_c.h
@@ -26,18 +26,18 @@
#include "sdis_Xd_begin.h"
struct XD(find_reinjection_ray_args) {
- const struct sdis_medium* solid; /* Medium into which the reinjection occurs */
- const struct XD(rwalk)* rwalk; /* Current random walk state */
+ const struct rwalk* rwalk; /* Current random walk state */
float dir0[DIM]; /* Challenged ray direction */
float dir1[DIM]; /* Challenged ray direction */
double distance; /* Maximum reinjection distance */
+ unsigned solid_enc_id; /* Enclosure id into which the reinjection occurs */
/* Define if the random walk position can be moved or not to find a valid
* reinjection direction */
int can_move;
};
static const struct XD(find_reinjection_ray_args)
-XD(FIND_REINJECTION_RAY_ARGS_NULL) = { NULL, NULL, {0}, {0}, 0, 0 };
+XD(FIND_REINJECTION_RAY_ARGS_NULL) = { NULL, {0}, {0}, 0, ENCLOSURE_ID_NULL, 0 };
struct XD(reinjection_ray) {
double org[DIM]; /* Origin of the reinjection */
@@ -55,80 +55,150 @@ XD(REINJECTION_RAY_NULL) = { {0}, {0}, 0, SXD_HIT_NULL__, 0 };
/*******************************************************************************
* Helper functions
******************************************************************************/
-static INLINE int
+static INLINE res_T
XD(check_find_reinjection_ray_args)
- (const struct XD(find_reinjection_ray_args)* args)
+ (struct sdis_scene* scn,
+ const struct XD(find_reinjection_ray_args)* args)
{
- return args
- && args->solid
- && args->rwalk
- && args->distance > 0
- && fX(is_normalized)(args->dir0)
- && fX(is_normalized)(args->dir1);
+ const struct enclosure* enc = NULL;
+ struct sdis_medium* mdm = NULL;
+ res_T res = RES_OK;
+ ASSERT(scn);
+
+ /* Check pointers */
+ if(!args || !args->rwalk) return RES_BAD_ARG;
+
+ /* Check distance */
+ if(args->distance <= 0) return RES_BAD_ARG;
+
+ /* Check directions */
+ if(!fX(is_normalized)(args->dir0) || !fX(is_normalized)(args->dir1)) {
+ return RES_BAD_ARG;
+ }
+
+ /* Check enclosure id */
+ if(args->solid_enc_id == ENCLOSURE_ID_NULL) {
+ return RES_BAD_ARG;
+ }
+ enc = scene_get_enclosure(scn, args->solid_enc_id);
+
+ /* Check the enclosure */
+ enc = scene_get_enclosure(scn, args->solid_enc_id);
+ if(enc->medium_id != MEDIUM_ID_MULTI) {
+ if((res = scene_get_enclosure_medium(scn, enc, &mdm)) != RES_OK) return res;
+ if(sdis_medium_get_type(mdm) != SDIS_SOLID) {
+ res = RES_BAD_ARG;
+ }
+ }
+
+ return RES_OK;
}
-static INLINE int
+static INLINE res_T
XD(check_sample_reinjection_step_args)
- (const struct XD(sample_reinjection_step_args)* args)
+ (struct sdis_scene* scn,
+ const struct sample_reinjection_step_args* args)
{
- return args
- && args->rng
- && args->solid
- && args->solid->type == SDIS_SOLID
- && args->rwalk
- && args->distance > 0
- && (unsigned)args->side < SDIS_SIDE_NULL__;
+ const struct enclosure* enc = NULL;
+ struct sdis_medium* mdm = NULL;
+ res_T res = RES_OK;
+ ASSERT(scn);
+
+ /* Check pointers */
+ if(!args || !args->rng || !args->rwalk) return RES_BAD_ARG;
+
+ /* Check distance */
+ if(args->distance <= 0) return RES_BAD_ARG;
+
+ /* Check side */
+ if((unsigned)args->side >= SDIS_SIDE_NULL__) return RES_BAD_ARG;
+
+ /* Check enclosure id */
+ if(args->solid_enc_id == ENCLOSURE_ID_NULL) {
+ return RES_BAD_ARG;
+ }
+
+ /* Check the enclosure */
+ enc = scene_get_enclosure(scn, args->solid_enc_id);
+ if(enc->medium_id != MEDIUM_ID_MULTI) {
+ if((res = scene_get_enclosure_medium(scn, enc, &mdm)) != RES_OK) return res;
+ if(sdis_medium_get_type(mdm) != SDIS_SOLID) {
+ return RES_BAD_ARG;
+ }
+ }
+
+ return RES_OK;
}
-static INLINE int
-XD(check_reinjection_step)(const struct XD(reinjection_step)* step)
+static INLINE res_T
+XD(check_reinjection_step)(const struct reinjection_step* step)
{
- return step
- && fX(is_normalized)(step->direction)
- && step->distance > 0;
+ /* Check pointer */
+ if(!step) return RES_BAD_ARG;
+
+ /* Check direction */
+ if(!fX(is_normalized)(step->direction)) return RES_BAD_ARG;
+
+ /* Check distance */
+ if(step->distance <= 0) return RES_BAD_ARG;
+
+ return RES_OK;
}
-static INLINE int
-XD(check_solid_reinjection_args)(const struct XD(solid_reinjection_args)* args)
+static INLINE res_T
+XD(check_solid_reinjection_args)(const struct solid_reinjection_args* args)
{
- return args
- && XD(check_reinjection_step)(args->reinjection)
- && args->rng
- && args->rwalk
- && args->rwalk_ctx
- && args->T
- && args->fp_to_meter > 0;
+ /* Check pointers */
+ if(!args || !args->rng || !args->rwalk || !args->rwalk_ctx || !args->T)
+ return RES_BAD_ARG;
+
+ /* Check unit */
+ if(args->fp_to_meter <= 0) return RES_BAD_ARG;
+
+ return XD(check_reinjection_step)(args->reinjection);
}
/* Check that the interface fragment is consistent with the current state of
* the random walk */
-static INLINE int
+static INLINE res_T
XD(check_rwalk_fragment_consistency)
- (const struct XD(rwalk)* rwalk,
+ (const struct rwalk* rwalk,
const struct sdis_interface_fragment* frag)
{
double N[DIM];
double uv[2] = {0, 0};
ASSERT(rwalk && frag);
- dX(normalize)(N, dX_set_fX(N, rwalk->hit.normal));
- if( SXD_HIT_NONE(&rwalk->hit)
- || !dX(eq_eps)(rwalk->vtx.P, frag->P, 1.e-6)
- || !dX(eq_eps)(N, frag->Ng, 1.e-6)
- || !( (IS_INF(rwalk->vtx.time) && IS_INF(frag->time))
- || eq_eps(rwalk->vtx.time, frag->time, 1.e-6))) {
- return 0;
+
+ /* Check intersection */
+ if(SXD_HIT_NONE(&rwalk->XD(hit))) return RES_BAD_ARG;
+
+ /* Check positions */
+ if(!dX(eq_eps)(rwalk->vtx.P, frag->P, 1.e-6)) return RES_BAD_ARG;
+
+ /* Check normals */
+ dX(normalize)(N, dX_set_fX(N, rwalk->XD(hit).normal));
+ if(!dX(eq_eps)(N, frag->Ng, 1.e-6)) return RES_BAD_ARG;
+
+ /* Check time */
+ if(!eq_eps(rwalk->vtx.time, frag->time, 1.e-6)
+ && !(IS_INF(rwalk->vtx.time) && IS_INF(frag->time))) {
+ return RES_BAD_ARG;
}
+
+ /* Check parametric coordinates */
#if (SDIS_XD_DIMENSION == 2)
- uv[0] = rwalk->hit.u;
+ uv[0] = rwalk->XD(hit).u;
#else
- d2_set_f2(uv, rwalk->hit.uv);
+ d2_set_f2(uv, rwalk->XD(hit).uv);
#endif
- return d2_eq_eps(uv, frag->uv, 1.e-6);
+ if(!d2_eq_eps(uv, frag->uv, 1.e-6)) return RES_BAD_ARG;
+
+ return RES_OK;
}
static void
XD(sample_reinjection_dir)
- (const struct XD(rwalk)* rwalk,
+ (const struct rwalk* rwalk,
struct ssp_rng* rng,
float dir[DIM])
{
@@ -145,15 +215,15 @@ XD(sample_reinjection_dir)
* discard the sqrt(2)/2 constant. */
const uint64_t r = ssp_rng_uniform_uint64(rng, 0, 1);
ASSERT(rwalk && rng && dir);
- ASSERT(!SXD_HIT_NONE(&rwalk->hit));
- ASSERT(!rwalk->mdm);
+ ASSERT(!SXD_HIT_NONE(&rwalk->XD(hit)));
+ ASSERT(rwalk->enc_id == ENCLOSURE_ID_NULL);
if(r) {
- dir[0] = rwalk->hit.normal[0] - rwalk->hit.normal[1];
- dir[1] = rwalk->hit.normal[0] + rwalk->hit.normal[1];
+ dir[0] = rwalk->XD(hit).normal[0] - rwalk->XD(hit).normal[1];
+ dir[1] = rwalk->XD(hit).normal[0] + rwalk->XD(hit).normal[1];
} else {
- dir[0] = rwalk->hit.normal[0] + rwalk->hit.normal[1];
- dir[1] =-rwalk->hit.normal[0] + rwalk->hit.normal[1];
+ dir[0] = rwalk->XD(hit).normal[0] + rwalk->XD(hit).normal[1];
+ dir[1] =-rwalk->XD(hit).normal[0] + rwalk->XD(hit).normal[1];
}
f2_normalize(dir, dir);
#else
@@ -162,17 +232,17 @@ XD(sample_reinjection_dir)
* radius of its base is 1. */
float frame[9];
ASSERT(rwalk && rng && dir);
- ASSERT(!SXD_HIT_NONE(&rwalk->hit));
- ASSERT(!rwalk->mdm);
- ASSERT(fX(is_normalized)(rwalk->hit.normal));
+ ASSERT(!SXD_HIT_NONE(&rwalk->XD(hit)));
+ ASSERT(rwalk->enc_id == ENCLOSURE_ID_NULL);
+ ASSERT(fX(is_normalized)(rwalk->XD(hit).normal));
ssp_ran_circle_uniform_float(rng, dir, NULL);
dir[2] = (float)(1.0/sqrt(2));
- f33_basis(frame, rwalk->hit.normal);
+ f33_basis(frame, rwalk->XD(hit).normal);
f33_mulf3(dir, frame, dir);
f3_normalize(dir, dir);
- ASSERT(eq_epsf(f3_dot(dir, rwalk->hit.normal), (float)(1.0/sqrt(3)), 1.e-4f));
+ ASSERT(eq_epsf(f3_dot(dir, rwalk->XD(hit).normal), (float)(1.0/sqrt(3)), 1.e-4f));
#endif
}
@@ -180,7 +250,7 @@ XD(sample_reinjection_dir)
#if DIM == 2
static void
XD(move_away_primitive_boundaries)
- (const struct XD(rwalk)* rwalk,
+ (const struct rwalk* rwalk,
const double delta,
double position[DIM]) /* Position to move */
{
@@ -189,9 +259,9 @@ XD(move_away_primitive_boundaries)
float dir[DIM];
float len;
const float st = 0.5f;
- ASSERT(rwalk && !SXD_HIT_NONE(&rwalk->hit) && delta > 0);
+ ASSERT(rwalk && !SXD_HIT_NONE(&rwalk->XD(hit)) && delta > 0);
- SXD(primitive_get_attrib(&rwalk->hit.prim, SXD_POSITION, st, &attr));
+ SXD(primitive_get_attrib(&rwalk->XD(hit).prim, SXD_POSITION, st, &attr));
fX_set_dX(pos, position);
fX(sub)(dir, attr.value, pos);
@@ -205,7 +275,7 @@ XD(move_away_primitive_boundaries)
* numerical issues leading to inconsistent random walks. */
static void
XD(move_away_primitive_boundaries)
- (const struct XD(rwalk)* rwalk,
+ (const struct rwalk* rwalk,
const double delta,
double position[DIM])
{
@@ -222,14 +292,14 @@ XD(move_away_primitive_boundaries)
int imin = 0;
int imid = 0;
int i;
- ASSERT(rwalk && delta > 0 && !S3D_HIT_NONE(&rwalk->hit));
+ ASSERT(rwalk && delta > 0 && !S3D_HIT_NONE(&rwalk->XD(hit)));
fX_set_dX(P, position);
/* 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));
+ S3D(triangle_get_vertex_attrib(&rwalk->XD(hit).prim, 0, S3D_POSITION, &v0));
+ S3D(triangle_get_vertex_attrib(&rwalk->XD(hit).prim, 1, S3D_POSITION, &v1));
+ S3D(triangle_get_vertex_attrib(&rwalk->XD(hit).prim, 2, S3D_POSITION, &v2));
/* Compute the edge vector */
f3_sub(E[0], v1.value, v0.value);
@@ -315,10 +385,10 @@ XD(find_reinjection_ray)
/* # attempts to find a ray direction */
int MAX_ATTEMPTS = 1;
- /* Physical properties */
- struct sdis_interface* interf;
- struct sdis_medium* mdm0;
- struct sdis_medium* mdm1;
+ /* Enclosures */
+ unsigned enc_ids[2] = {ENCLOSURE_ID_NULL, ENCLOSURE_ID_NULL};
+ unsigned enc0_id = ENCLOSURE_ID_NULL;
+ unsigned enc1_id = ENCLOSURE_ID_NULL;
struct hit_filter_data filter_data = HIT_FILTER_DATA_NULL;
struct sXd(hit) hit;
@@ -338,7 +408,7 @@ XD(find_reinjection_ray)
res_T res = RES_OK;
ASSERT(scn && args && ray);
- ASSERT(XD(check_find_reinjection_ray_args)(args));
+ ASSERT(XD(check_find_reinjection_ray_args)(scn, args) == RES_OK);
*ray = XD(REINJECTION_RAY_NULL);
MAX_ATTEMPTS = args->can_move ? 2 : 1;
@@ -350,39 +420,39 @@ XD(find_reinjection_ray)
do {
fX_set_dX(org, ray->org);
- filter_data.XD(hit) = args->rwalk->hit;
+ filter_data.XD(hit) = args->rwalk->XD(hit);
filter_data.epsilon = args->distance * 0.01;
SXD(scene_view_trace_ray
(scn->sXd(view), org, args->dir0, range, &filter_data, &hit0));
SXD(scene_view_trace_ray
(scn->sXd(view), org, args->dir1, range, &filter_data, &hit1));
- /* Retrieve the medium at the reinjection pos along dir0 */
- if(SXD_HIT_NONE(&hit0)) {
+ /* Retrieve the enclosure at the reinjection pos along dir0 */
+ if(!SXD_HIT_NONE(&hit0)) {
+ scene_get_enclosure_ids(scn, hit0.prim.prim_id, enc_ids);
+ side = fX(dot)(args->dir0, hit0.normal) < 0 ? SDIS_FRONT : SDIS_BACK;
+ enc0_id = enc_ids[side];
+ } else {
XD(move_pos)(dX(set)(tmp, ray->org), args->dir0, (float)args->distance);
- res = scene_get_medium_in_closed_boundaries(scn, tmp, &mdm0);
- if(res == RES_BAD_OP) { mdm0 = NULL; res = RES_OK; }
+ res = scene_get_enclosure_id_in_closed_boundaries(scn, tmp, &enc0_id);
+ if(res == RES_BAD_OP) { enc0_id = ENCLOSURE_ID_NULL; res = RES_OK; }
if(res != RES_OK) goto error;
- } else {
- interf = scene_get_interface(scn, hit0.prim.prim_id);
- side = fX(dot)(args->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)) {
+ /* Retrieve the enclosure at the reinjection pos along dir1 */
+ if(!SXD_HIT_NONE(&hit1)) {
+ scene_get_enclosure_ids(scn, hit1.prim.prim_id, enc_ids);
+ side = fX(dot)(args->dir1, hit1.normal) < 0 ? SDIS_FRONT : SDIS_BACK;
+ enc1_id = enc_ids[side];
+ } else {
XD(move_pos)(dX(set)(tmp, ray->org), args->dir1, (float)args->distance);
- res = scene_get_medium_in_closed_boundaries(scn, tmp, &mdm1);
- if(res == RES_BAD_OP) { mdm1 = NULL; res = RES_OK; }
+ res = scene_get_enclosure_id_in_closed_boundaries(scn, tmp, &enc1_id);
+ if(res == RES_BAD_OP) { enc1_id = ENCLOSURE_ID_NULL; res = RES_OK; }
if(res != RES_OK) goto error;
- } else {
- interf = scene_get_interface(scn, hit1.prim.prim_id);
- side = fX(dot)(args->dir1, hit1.normal) < 0 ? SDIS_FRONT : SDIS_BACK;
- mdm1 = interface_get_medium(interf, side);
}
dst0 = dst1 = -1;
- if(mdm0 == args->solid) { /* Check reinjection consistency */
+ if(enc0_id == args->solid_enc_id) { /* Check reinjection consistency */
if(hit0.distance <= dst_adjusted) {
dst0 = hit0.distance;
} else {
@@ -390,7 +460,7 @@ XD(find_reinjection_ray)
hit0 = SXD_HIT_NULL;
}
}
- if(mdm1 == args->solid) { /* Check reinjection consistency */
+ if(enc1_id == args->solid_enc_id) { /* Check reinjection consistency */
if(hit1.distance <= dst_adjusted) {
dst1 = hit1.distance;
} else {
@@ -503,23 +573,25 @@ XD(find_reinjection_ray_and_check_validity)
struct XD(reinjection_ray)* ray)
{
double pos[DIM];
- struct sdis_medium* reinject_mdm;
res_T res = RES_OK;
ASSERT(scn && args && ray);
- ASSERT(XD(check_find_reinjection_ray_args)(args));
+ ASSERT(XD(check_find_reinjection_ray_args)(scn, args) == RES_OK);
/* Select a reinjection direction */
res = XD(find_reinjection_ray)(scn, args, ray);
if(res != RES_OK) goto error;
if(SXD_HIT_NONE(&ray->hit)) {
- /* Check medium consistency at the reinjection position */
+ unsigned enc_id = ENCLOSURE_ID_NULL;
+
+ /* Obtain the enclosure in which the reinjection position lies */
XD(move_pos)(dX(set)(pos, ray->org), ray->dir, (float)ray->dst);
- res = scene_get_medium_in_closed_boundaries(scn, pos, &reinject_mdm);
+ res = scene_get_enclosure_id_in_closed_boundaries(scn, pos, &enc_id);
if(res != RES_OK) goto error;
- if(reinject_mdm != args->solid) {
+ /* Check enclosure consistency at the reinjection position */
+ if(enc_id != args->solid_enc_id) {
res = RES_BAD_OP;
goto error;
}
@@ -535,9 +607,9 @@ static res_T
XD(handle_volumic_power)
(struct sdis_medium* solid,
struct rwalk_context* rwalk_ctx,
- struct XD(rwalk)* rwalk,
+ struct rwalk* rwalk,
const double reinject_dst_m,
- struct XD(temperature)* T)
+ struct temperature* T)
{
double power;
double lambda;
@@ -592,22 +664,22 @@ error:
res_T
XD(sample_reinjection_step_solid_fluid)
(struct sdis_scene* scn,
- const struct XD(sample_reinjection_step_args)* args,
- struct XD(reinjection_step)* step)
+ const struct sample_reinjection_step_args* args,
+ struct reinjection_step* step)
{
/* Input/output data of the function finding a valid reinjection ray */
struct XD(find_reinjection_ray_args) find_reinject_ray_args =
XD(FIND_REINJECTION_RAY_ARGS_NULL);
struct XD(reinjection_ray) ray = XD(REINJECTION_RAY_NULL);
+ /* Enclosures */
+ const struct enclosure* solid_enc = NULL;
+ unsigned enc_ids[2] = {ENCLOSURE_ID_NULL, ENCLOSURE_ID_NULL};
+
/* 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;
- /* Enclosure */
- const struct enclosure* solid_enclosure = NULL;
- unsigned enc_ids[2];
-
/* Miscellaneous variables */
float dir0[DIM]; /* Sampled direction */
float dir1[DIM]; /* Sampled direction reflected */
@@ -616,23 +688,23 @@ XD(sample_reinjection_step_solid_fluid)
/* Pre-conditions */
ASSERT(scn && args && step);
- ASSERT(XD(check_sample_reinjection_step_args)(args));
+ ASSERT(XD(check_sample_reinjection_step_args)(scn, args) == RES_OK);
/* Initialise the reinjection step */
- *step = XD(REINJECTION_STEP_NULL);
+ *step = REINJECTION_STEP_NULL;
/* Check whether the solid enclosure is part of the system, or whether it is
* only there to describe boundary conditions. In the latter case, the
* enclosure has no geometrical existence and it is sufficient to return a
* valid reinjection step which will be used to select the next step. Note
* that if the trajectory passes through the solid enclosure, it will stop,
- * i.e. the temperature of the solid should be fixed. If it doesn't, it's an
+ * i.e. the temperature of the solid should be fixed. If it doesn't, it's a
* error */
- scene_get_enclosure_ids(scn, args->rwalk->hit.prim.prim_id, enc_ids);
- solid_enclosure = scene_get_enclosure(scn, enc_ids[args->side]);
- if(solid_enclosure->medium_id == ENCLOSURE_MULTI_MEDIA) {
- step->hit = SXD_HIT_NULL;
- fX(normalize)(step->direction, args->rwalk->hit.normal);
+ scene_get_enclosure_ids(scn, args->rwalk->XD(hit).prim.prim_id, enc_ids);
+ solid_enc = scene_get_enclosure(scn, enc_ids[args->side]);
+ if(solid_enc->medium_id == MEDIUM_ID_MULTI) {
+ step->XD(hit) = SXD_HIT_NULL;
+ fX(normalize)(step->direction, args->rwalk->XD(hit).normal);
if(args->side == SDIS_BACK) fX(minus)(step->direction, step->direction);
step->distance = (float)args->distance;
goto exit; /* That's all folks! */
@@ -644,7 +716,7 @@ XD(sample_reinjection_step_solid_fluid)
XD(sample_reinjection_dir)(args->rwalk, args->rng, dir0);
/* Reflect the sampled direction around the normal */
- XD(reflect)(dir1, dir0, args->rwalk->hit.normal);
+ XD(reflect)(dir1, dir0, args->rwalk->XD(hit).normal);
/* Flip the sampled directions if one wants to reinject to back side */
if(args->side == SDIS_BACK) {
@@ -653,7 +725,7 @@ XD(sample_reinjection_step_solid_fluid)
}
/* Find the reinjection step */
- find_reinject_ray_args.solid = args->solid;
+ find_reinject_ray_args.solid_enc_id = args->solid_enc_id;
find_reinject_ray_args.rwalk = args->rwalk;
find_reinject_ray_args.distance = args->distance;
find_reinject_ray_args.can_move = 1;
@@ -676,7 +748,7 @@ XD(sample_reinjection_step_solid_fluid)
}
/* Setup the reinjection step */
- step->hit = ray.hit;
+ step->XD(hit) = ray.hit;
step->distance = ray.dst;
fX(set)(step->direction, ray.dir);
@@ -687,7 +759,7 @@ XD(sample_reinjection_step_solid_fluid)
/* Post-conditions */
ASSERT(dX(eq)(args->rwalk->vtx.P, ray.org));
- ASSERT(XD(check_reinjection_step)(step));
+ ASSERT(XD(check_reinjection_step)(step) == RES_OK);
exit:
return res;
@@ -698,10 +770,10 @@ error:
res_T
XD(sample_reinjection_step_solid_solid)
(struct sdis_scene* scn,
- const struct XD(sample_reinjection_step_args)* args_frt,
- const struct XD(sample_reinjection_step_args)* args_bck,
- struct XD(reinjection_step)* step_frt,
- struct XD(reinjection_step)* step_bck)
+ const struct sample_reinjection_step_args* args_frt,
+ const struct sample_reinjection_step_args* args_bck,
+ struct reinjection_step* step_frt,
+ struct reinjection_step* step_bck)
{
/* Input/output data of the function finding a valid reinjection ray */
struct XD(find_reinjection_ray_args) find_reinject_ray_frt_args =
@@ -715,7 +787,7 @@ XD(sample_reinjection_step_solid_solid)
double rwalk_pos_backup[DIM];
/* Variables shared by the two side */
- struct XD(rwalk)* rwalk = NULL;
+ struct rwalk* rwalk = NULL;
struct ssp_rng* rng = NULL;
/* In 2D it is useless to try to resample a reinjection direction since there
@@ -723,26 +795,26 @@ XD(sample_reinjection_step_solid_solid)
const int MAX_ATTEMPTS = DIM == 2 ? 1 : 10;
/* Enclosure */
- const struct enclosure* enclosure_bck = NULL;
- const struct enclosure* enclosure_frt = NULL;
- int multiple_media_bck = 0;
- int multiple_media_frt = 0;
unsigned enc_ids[2];
+ const struct enclosure* enc_frt = NULL;
+ const struct enclosure* enc_bck = NULL;
float dir_frt_samp[DIM]; /* Sampled direction */
float dir_frt_refl[DIM]; /* Sampled direction reflected */
float dir_bck_samp[DIM]; /* Negated sampled direction */
float dir_bck_refl[DIM]; /* Negated sampled direction reflected */
+ int multi_frt = 0;
+ int multi_bck = 0;
int iattempt = 0; /* #attempts to find a reinjection dir */
res_T res = RES_OK;
/* Pre-conditions */
ASSERT(scn && args_frt && args_bck && step_frt && step_bck);
- ASSERT(XD(check_sample_reinjection_step_args)(args_frt));
- ASSERT(XD(check_sample_reinjection_step_args)(args_bck));
+ ASSERT(XD(check_sample_reinjection_step_args)(scn, args_frt) == RES_OK);
+ ASSERT(XD(check_sample_reinjection_step_args)(scn, args_bck) == RES_OK);
ASSERT(args_frt->side == SDIS_FRONT);
ASSERT(args_bck->side == SDIS_BACK);
- ASSERT(SXD_PRIMITIVE_EQ(&args_frt->rwalk->hit.prim, &args_bck->rwalk->hit.prim));
+ ASSERT(SXD_PRIMITIVE_EQ(&args_frt->rwalk->XD(hit).prim, &args_bck->rwalk->XD(hit).prim));
rng = args_frt->rng;
rwalk = args_frt->rwalk;
@@ -750,8 +822,8 @@ XD(sample_reinjection_step_solid_solid)
ASSERT(args_bck->rwalk == rwalk);
/* Initialise the reinjection steps */
- *step_frt = XD(REINJECTION_STEP_NULL);
- *step_bck = XD(REINJECTION_STEP_NULL);
+ *step_frt = REINJECTION_STEP_NULL;
+ *step_bck = REINJECTION_STEP_NULL;
/* Check whether the solid enclosure is part of the system, or whether it is
* only there to describe boundary conditions. In the latter case, the
@@ -760,25 +832,23 @@ XD(sample_reinjection_step_solid_solid)
* that if the trajectory passes through the solid enclosure, it will stop,
* i.e. the temperature of the solid should be fixed. If it doesn't, it's an
* error */
- scene_get_enclosure_ids(scn, args_frt->rwalk->hit.prim.prim_id, enc_ids);
- enclosure_frt = scene_get_enclosure(scn, enc_ids[SDIS_FRONT]);
- enclosure_bck = scene_get_enclosure(scn, enc_ids[SDIS_BACK]);
- if(enclosure_frt->medium_id == ENCLOSURE_MULTI_MEDIA) {
- step_frt->hit = SXD_HIT_NULL;
- fX(normalize)(step_frt->direction, args_frt->rwalk->hit.normal);
+ scene_get_enclosure_ids(scn, args_frt->rwalk->XD(hit).prim.prim_id, enc_ids);
+ enc_frt = scene_get_enclosure(scn, enc_ids[SDIS_FRONT]);
+ enc_bck = scene_get_enclosure(scn, enc_ids[SDIS_BACK]);
+ if(enc_frt->medium_id == MEDIUM_ID_MULTI) {
+ step_frt->XD(hit) = SXD_HIT_NULL;
+ fX(normalize)(step_frt->direction, args_frt->rwalk->XD(hit).normal);
step_frt->distance = (float)args_frt->distance;
- multiple_media_frt = 1;
+ multi_frt = 1;
}
- if(enclosure_bck->medium_id == ENCLOSURE_MULTI_MEDIA) {
- step_bck->hit = SXD_HIT_NULL;
- fX(normalize)(step_bck->direction, args_bck->rwalk->hit.normal);
- fX(minus)(step_bck->direction, step_bck->direction);
+ if(enc_bck->medium_id == MEDIUM_ID_MULTI) {
+ step_bck->XD(hit) = SXD_HIT_NULL;
+ fX(normalize)(step_bck->direction, args_bck->rwalk->XD(hit).normal);
step_bck->distance = (float)args_bck->distance;
- multiple_media_bck = 1;
+ multi_bck = 1;
}
- if(multiple_media_frt && multiple_media_bck)
- goto exit; /* That's all folks! */
+ if(multi_frt && multi_bck) goto exit; /* That's all folks */
dX(set)(rwalk_pos_backup, rwalk->vtx.P);
iattempt = 0;
@@ -789,16 +859,15 @@ XD(sample_reinjection_step_solid_solid)
/* 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, dir_frt_samp);
- XD(reflect)(dir_frt_refl, dir_frt_samp, rwalk->hit.normal);
+ XD(reflect)(dir_frt_refl, dir_frt_samp, rwalk->XD(hit).normal);
fX(minus)(dir_bck_samp, dir_frt_samp);
fX(minus)(dir_bck_refl, dir_frt_refl);
/* Reject the sampling of the re-injection step if it has already been
* defined, i.e. if the enclosure is a limit condition */
- if(!multiple_media_frt) {
-
+ if(!multi_frt) {
/* Find the reinjection ray for the front side */
- find_reinject_ray_frt_args.solid = args_frt->solid;
+ find_reinject_ray_frt_args.solid_enc_id = args_frt->solid_enc_id;
find_reinject_ray_frt_args.rwalk = args_frt->rwalk;
find_reinject_ray_frt_args.distance = args_frt->distance;
find_reinject_ray_frt_args.can_move = 1;
@@ -815,10 +884,9 @@ XD(sample_reinjection_step_solid_solid)
/* Reject the sampling of the re-injection step if it has already been
* defined, i.e. if the enclosure is a limit condition */
- if(!multiple_media_bck) {
-
+ if(!multi_bck) {
/* Select the reinjection direction and distance for the back side */
- find_reinject_ray_bck_args.solid = args_bck->solid;
+ find_reinject_ray_bck_args.solid_enc_id = args_bck->solid_enc_id;
find_reinject_ray_bck_args.rwalk = args_bck->rwalk;
find_reinject_ray_bck_args.distance = args_bck->distance;
find_reinject_ray_bck_args.can_move = 1;
@@ -834,7 +902,7 @@ XD(sample_reinjection_step_solid_solid)
/* If random walk was moved to find a valid rinjection ray on back side,
* one has to find a valid reinjection on front side from the new pos */
- if(ray_bck.position_was_moved && !multiple_media_frt) {
+ if(ray_bck.position_was_moved) {
find_reinject_ray_frt_args.can_move = 0;
res = XD(find_reinjection_ray_and_check_validity)
(scn, &find_reinject_ray_frt_args, &ray_frt);
@@ -858,17 +926,17 @@ XD(sample_reinjection_step_solid_solid)
}
/* Setup the front and back reinjection steps */
- if(!multiple_media_frt) {
- step_frt->hit = ray_frt.hit;
+ if(!multi_frt) {
+ step_frt->XD(hit) = ray_frt.hit;
step_frt->distance = ray_frt.dst;
fX(set)(step_frt->direction, ray_frt.dir);
- ASSERT(XD(check_reinjection_step)(step_frt)); /* Post-condition */
+ ASSERT(XD(check_reinjection_step)(step_frt) == RES_OK); /* Post-condition */
}
- if(!multiple_media_bck) {
- step_bck->hit = ray_bck.hit;
+ if(!multi_bck) {
+ step_bck->XD(hit) = ray_bck.hit;
step_bck->distance = ray_bck.dst;
fX(set)(step_bck->direction, ray_bck.dir);
- ASSERT(XD(check_reinjection_step)(step_bck)); /* Post-condition */
+ ASSERT(XD(check_reinjection_step)(step_bck) == RES_OK); /* Post-condition */
}
exit:
@@ -879,17 +947,28 @@ error:
res_T
XD(solid_reinjection)
- (struct sdis_medium* solid,
- struct XD(solid_reinjection_args)* args)
+ (struct sdis_scene* scn,
+ const unsigned solid_enc_id,
+ struct solid_reinjection_args* args)
{
+ /* Properties */
struct solid_props props = SOLID_PROPS_NULL;
+ struct sdis_medium* solid = NULL;
+ const struct enclosure* enc = NULL;
+
double reinject_dst_m; /* Reinjection distance in meters */
double mu;
res_T res = RES_OK;
- ASSERT(solid && XD(check_solid_reinjection_args)(args));
+ ASSERT(XD(check_solid_reinjection_args)(args) == RES_OK);
+ ASSERT(solid_enc_id != ENCLOSURE_ID_NULL);
reinject_dst_m = args->reinjection->distance * args->fp_to_meter;
+ /* Get the enclosure medium properties */
+ enc = scene_get_enclosure(scn, solid_enc_id);
+ res = scene_get_enclosure_medium(scn, enc, &solid);
+ if(res != RES_OK) goto error;
+ ASSERT(sdis_medium_get_type(solid) == SDIS_SOLID);
res = solid_get_properties(solid, &args->rwalk->vtx, &props);
if(res != RES_OK) goto error;
@@ -899,10 +978,10 @@ XD(solid_reinjection)
if(res != RES_OK) goto error;
/* Time rewind */
- args->rwalk->mdm = solid; /* Medium into which the time is rewind */
+ args->rwalk->enc_id = solid_enc_id; /* Enclosure into which the time is rewind */
mu = (2*DIM*props.lambda)/(props.rho*props.cp*reinject_dst_m*reinject_dst_m);
- res = XD(time_rewind)
- (mu, props.t0, args->rng, args->rwalk, args->rwalk_ctx, args->T);
+ res = time_rewind
+ (scn, mu, props.t0, args->rng, args->rwalk, args->rwalk_ctx, args->T);
if(res != RES_OK) goto error;
/* Test if a limit condition was reached */
@@ -915,18 +994,18 @@ XD(solid_reinjection)
args->reinjection->distance);
/* The random walk is in the solid */
- if(args->reinjection->hit.distance != args->reinjection->distance) {
+ if(args->reinjection->XD(hit).distance != args->reinjection->distance) {
args->T->func = XD(conductive_path);
- args->rwalk->mdm = solid;
- args->rwalk->hit = SXD_HIT_NULL;
+ args->rwalk->enc_id = solid_enc_id;
+ args->rwalk->XD(hit) = SXD_HIT_NULL;
args->rwalk->hit_side = SDIS_SIDE_NULL__;
/* The random walk is at a boundary */
} else {
args->T->func = XD(boundary_path);
- args->rwalk->mdm = NULL;
- args->rwalk->hit = args->reinjection->hit;
- if(fX(dot)(args->reinjection->hit.normal, args->reinjection->direction) < 0) {
+ args->rwalk->enc_id = ENCLOSURE_ID_NULL;
+ args->rwalk->XD(hit) = args->reinjection->XD(hit);
+ if(fX(dot)(args->reinjection->XD(hit).normal, args->reinjection->direction) < 0) {
args->rwalk->hit_side = SDIS_FRONT;
} else {
args->rwalk->hit_side = SDIS_BACK;
@@ -952,7 +1031,7 @@ res_T
XD(handle_net_flux)
(const struct sdis_scene* scn,
const struct handle_net_flux_args* args,
- struct XD(temperature)* T)
+ struct temperature* T)
{
double flux_term;
double phi;
@@ -1013,7 +1092,7 @@ XD(check_Tref)
if((Bound) < 0) { \
log_err(scn->dev, \
"%s: the "Name" temperature cannot be negative " \
- "to sample a radiative path (T"Name" = %g K).\n", \
+ "to sample a radiative path -- T"Name" = %g K\n", \
func_name, (Bound)); \
return RES_BAD_OP_IRRECOVERABLE; \
} \
@@ -1032,7 +1111,7 @@ XD(check_Tref)
if(SDIS_TEMPERATURE_IS_UNKNOWN(Tref)) {
log_err(scn->dev,
- "%s: the reference temperature is unknown at `"FORMAT_VECX". "
+ "%s: the reference temperature is unknown at ("FORMAT_VECX"). "
"Sampling a radiative path requires a valid reference temperature field.\n",
func_name, SPLITX(pos));
return RES_BAD_OP_IRRECOVERABLE;
@@ -1040,7 +1119,7 @@ XD(check_Tref)
if(Tref < 0) {
log_err(scn->dev,
- "%s: the reference temperature is negative at `"FORMAT_VECX" (Tref = %g K). "
+ "%s: the reference temperature is negative at ("FORMAT_VECX") and Tref = %g K. "
"Sampling a radiative path requires a known, positive reference "
"temperature field.\n",
func_name, SPLITX(pos), Tref);
@@ -1049,7 +1128,7 @@ XD(check_Tref)
if(Tref < scn->tmin || scn->tmax < Tref) {
log_err(scn->dev,
- "%s: invalid reference temperature at `"FORMAT_VECX"' (Tref = %g K). "
+ "%s: invalid reference temperature at ("FORMAT_VECX") and Tref=%g K. "
"It must be included in the provided temperature range "
"(Tmin = %g K; Tmax = %g K)\n",
func_name, SPLITX(pos), Tref, scn->tmin, scn->tmax);
diff --git a/src/sdis_heat_path_boundary_Xd_handle_external_net_flux.h b/src/sdis_heat_path_boundary_Xd_handle_external_net_flux.h
@@ -173,7 +173,7 @@ static INLINE res_T
XD(check_handle_external_net_flux_args)
(const struct sdis_scene* scn,
const char* func_name,
- const struct XD(handle_external_net_flux_args)* args)
+ const struct handle_external_net_flux_args* args)
{
int net_flux = 0;
res_T res = RES_OK;
@@ -181,7 +181,7 @@ XD(check_handle_external_net_flux_args)
/* Handle bugs */
ASSERT(scn && func_name && args);
ASSERT(args->interf && args->frag);
- ASSERT(!SXD_HIT_NONE(args->hit));
+ ASSERT(!SXD_HIT_NONE(args->XD(hit)));
ASSERT(args->h_cond >= 0 && args->h_conv >= 0 && args->h_radi >= 0);
ASSERT(args->h_cond + args->h_conv + args->h_radi > 0);
@@ -453,8 +453,8 @@ res_T
XD(handle_external_net_flux)
(const struct sdis_scene* scn,
struct ssp_rng* rng,
- const struct XD(handle_external_net_flux_args)* args,
- struct XD(temperature)* T)
+ const struct handle_external_net_flux_args* args,
+ struct temperature* T)
{
/* Terms to be registered in the green function */
struct sdis_green_external_flux_terms green =
@@ -516,13 +516,13 @@ XD(handle_external_net_flux)
* interface side */
cos_theta = d3_dot(N, src_sample.dir);
if(cos_theta > 0) {
- Ld = XD(direct_contribution)(scn, &src_sample, frag.P, args->hit);
+ Ld = XD(direct_contribution)(scn, &src_sample, frag.P, args->XD(hit));
incident_flux_direct = cos_theta * Ld / src_sample.pdf; /* [W/m^2] */
}
/* Calculate the incident diffuse flux [W/m^2] */
res = XD(compute_incident_diffuse_flux)
- (scn, rng, frag.P, N, frag.time, args->hit, &incident_flux_diffuse);
+ (scn, rng, frag.P, N, frag.time, args->XD(hit), &incident_flux_diffuse);
if(res != RES_OK) goto error;
/* Calculate the incident flux without the part scattered by the environment.
diff --git a/src/sdis_heat_path_boundary_Xd_solid_fluid_picard1.h b/src/sdis_heat_path_boundary_Xd_solid_fluid_picard1.h
@@ -31,8 +31,8 @@
static INLINE res_T
XD(rwalk_get_Tref)
(const struct sdis_scene* scn,
- const struct XD(rwalk)* rwalk,
- const struct XD(temperature)* T,
+ const struct rwalk* rwalk,
+ const struct temperature* T,
double* out_Tref)
{
double Tref = SDIS_TEMPERATURE_NONE;
@@ -52,16 +52,16 @@ XD(rwalk_get_Tref)
} else {
struct sdis_interface_fragment frag;
struct sdis_interface* interf = NULL;
- ASSERT(!SXD_HIT_NONE(&rwalk->hit));
+ ASSERT(!SXD_HIT_NONE(&rwalk->XD(hit)));
/* Fetch the interface where the random walk ends */
- interf = scene_get_interface(scn, rwalk->hit.prim.prim_id);
+ interf = scene_get_interface(scn, rwalk->XD(hit).prim.prim_id);
ASSERT(rwalk->hit_side!=SDIS_FRONT || interf->medium_front->type==SDIS_FLUID);
ASSERT(rwalk->hit_side!=SDIS_BACK || interf->medium_back->type==SDIS_FLUID);
/* Fragment on the fluid side of the boundary onto which the rwalk ends */
XD(setup_interface_fragment)
- (&frag, &rwalk->vtx, &rwalk->hit, rwalk->hit_side);
+ (&frag, &rwalk->vtx, &rwalk->XD(hit), rwalk->hit_side);
Tref = interface_side_get_reference_temperature(interf, &frag);
}
@@ -85,30 +85,32 @@ XD(solid_fluid_boundary_picard1_path)
(struct sdis_scene* scn,
struct rwalk_context* ctx,
const struct sdis_interface_fragment* frag,
- struct XD(rwalk)* rwalk,
+ struct rwalk* rwalk,
struct ssp_rng* rng,
- struct XD(temperature)* T)
+ struct temperature* T)
{
/* Input argument used to handle the net flux */
struct handle_net_flux_args handle_net_flux_args = HANDLE_NET_FLUX_ARGS_NULL;
/* Input argument used to handle the external net flux */
- struct XD(handle_external_net_flux_args) handle_external_net_flux_args =
- XD(HANDLE_EXTERNAL_NET_FLUX_ARGS_NULL);
+ struct handle_external_net_flux_args handle_external_net_flux_args =
+ HANDLE_EXTERNAL_NET_FLUX_ARGS_NULL;
/* Input/output arguments of the function used to sample a reinjection */
- struct XD(sample_reinjection_step_args) samp_reinject_step_args =
- XD(SAMPLE_REINJECTION_STEP_ARGS_NULL);
- struct XD(reinjection_step) reinject_step =
- XD(REINJECTION_STEP_NULL);
+ struct sample_reinjection_step_args samp_reinject_step_args =
+ SAMPLE_REINJECTION_STEP_ARGS_NULL;
+ struct reinjection_step reinject_step = REINJECTION_STEP_NULL;
/* Temperature and random walk state of the sampled radiative path */
- struct XD(temperature) T_s;
- struct XD(rwalk) rwalk_s;
+ struct temperature T_s;
+ struct rwalk rwalk_s;
/* Fragment on the fluid side of the boundary */
struct sdis_interface_fragment frag_fluid;
+ /* The enclosures split by the boundary */
+ unsigned enc_ids[2] = {ENCLOSURE_ID_NULL, ENCLOSURE_ID_NULL};
+
/* Data attached to the boundary */
struct sdis_interface* interf = NULL;
struct sdis_medium* solid = NULL;
@@ -136,10 +138,10 @@ XD(solid_fluid_boundary_picard1_path)
res_T res = RES_OK;
ASSERT(scn && rwalk && rng && T && ctx);
- ASSERT(XD(check_rwalk_fragment_consistency)(rwalk, frag));
+ ASSERT(XD(check_rwalk_fragment_consistency)(rwalk, frag) == RES_OK);
/* Retrieve the solid and the fluid split by the boundary */
- interf = scene_get_interface(scn, rwalk->hit.prim.prim_id);
+ interf = scene_get_interface(scn, rwalk->XD(hit).prim.prim_id);
solid = interface_get_medium(interf, SDIS_FRONT);
fluid = interface_get_medium(interf, SDIS_BACK);
solid_side = SDIS_FRONT;
@@ -150,6 +152,9 @@ XD(solid_fluid_boundary_picard1_path)
ASSERT(fluid->type == SDIS_FLUID);
}
+ /* Get the enclosures split by the boundary */
+ scene_get_enclosure_ids(scn, rwalk->XD(hit).prim.prim_id, enc_ids);
+
/* Setup a fragment for the fluid side */
frag_fluid = *frag;
frag_fluid.side = fluid_side;
@@ -178,7 +183,7 @@ XD(solid_fluid_boundary_picard1_path)
/* Sample a reinjection step */
samp_reinject_step_args.rng = rng;
- samp_reinject_step_args.solid = solid;
+ samp_reinject_step_args.solid_enc_id = enc_ids[solid_side];
samp_reinject_step_args.rwalk = rwalk;
samp_reinject_step_args.distance = delta_boundary;
samp_reinject_step_args.side = solid_side;
@@ -222,7 +227,7 @@ XD(solid_fluid_boundary_picard1_path)
/* Handle the external net flux if any */
handle_external_net_flux_args.interf = interf;
handle_external_net_flux_args.frag = frag;
- handle_external_net_flux_args.hit = &rwalk->hit;
+ handle_external_net_flux_args.XD(hit) = &rwalk->XD(hit);
handle_external_net_flux_args.green_path = ctx->green_path;
handle_external_net_flux_args.picard_order = get_picard_order(ctx);
handle_external_net_flux_args.h_cond = h_cond;
@@ -248,15 +253,15 @@ XD(solid_fluid_boundary_picard1_path)
/* Switch in convective path */
if(r < p_conv) {
T->func = XD(convective_path);
- rwalk->mdm = fluid;
+ rwalk->enc_id = enc_ids[fluid_side];
rwalk->hit_side = fluid_side;
break;
}
/* Switch in conductive path */
if(r < p_conv + p_cond) {
- struct XD(solid_reinjection_args) solid_reinject_args =
- XD(SOLID_REINJECTION_ARGS_NULL);
+ struct solid_reinjection_args solid_reinject_args =
+ SOLID_REINJECTION_ARGS_NULL;
/* Perform the reinjection into the solid */
solid_reinject_args.reinjection = &reinject_step;
@@ -265,7 +270,7 @@ XD(solid_fluid_boundary_picard1_path)
solid_reinject_args.rng = rng;
solid_reinject_args.T = T;
solid_reinject_args.fp_to_meter = scn->fp_to_meter;
- res = XD(solid_reinjection)(solid, &solid_reinject_args);
+ res = XD(solid_reinjection)(scn, enc_ids[solid_side], &solid_reinject_args);
if(res != RES_OK) goto error;
break;
}
@@ -282,7 +287,7 @@ XD(solid_fluid_boundary_picard1_path)
/* Sample a radiative path and get the Tref at its end. */
T_s = *T;
rwalk_s = *rwalk;
- rwalk_s.mdm = fluid;
+ rwalk_s.enc_id = enc_ids[fluid_side];
rwalk_s.hit_side = fluid_side;
res = XD(radiative_path)(scn, ctx, &rwalk_s, rng, &T_s);
if(res != RES_OK) goto error;
diff --git a/src/sdis_heat_path_boundary_Xd_solid_fluid_picardN.h b/src/sdis_heat_path_boundary_Xd_solid_fluid_picardN.h
@@ -66,23 +66,23 @@ error:
static INLINE res_T
XD(sample_path)
(struct sdis_scene* scn,
- const struct XD(rwalk)* rwalk_from,
+ const struct rwalk* rwalk_from,
struct rwalk_context* ctx,
struct ssp_rng* rng,
- struct XD(temperature)* T)
+ struct temperature* T)
{
- struct XD(rwalk) rwalk = XD(RWALK_NULL);
+ struct rwalk rwalk = RWALK_NULL;
res_T res = RES_OK;
ASSERT(rwalk_from && rng && T);
/* Clean-up the output variable */
- *T = XD(TEMPERATURE_NULL);
+ *T = TEMPERATURE_NULL;
T->func = XD(boundary_path);
/* Init the random walk */
rwalk.vtx = rwalk_from->vtx;
- rwalk.mdm = rwalk_from->mdm;
- rwalk.hit = rwalk_from->hit;
+ rwalk.enc_id = rwalk_from->enc_id;
+ rwalk.XD(hit) = rwalk_from->XD(hit);
rwalk.hit_side = rwalk_from->hit_side;
/* Start the registration of a new heat path */
@@ -128,15 +128,14 @@ XD(solid_fluid_boundary_picardN_path)
(struct sdis_scene* scn,
struct rwalk_context* ctx,
const struct sdis_interface_fragment* frag,
- struct XD(rwalk)* rwalk,
+ struct rwalk* rwalk,
struct ssp_rng* rng,
- struct XD(temperature)* T)
+ struct temperature* T)
{
/* Input/output arguments of the function used to sample a reinjection */
- struct XD(sample_reinjection_step_args) samp_reinject_step_args =
- XD(SAMPLE_REINJECTION_STEP_ARGS_NULL);
- struct XD(reinjection_step) reinject_step =
- XD(REINJECTION_STEP_NULL);
+ struct sample_reinjection_step_args samp_reinject_step_args =
+ SAMPLE_REINJECTION_STEP_ARGS_NULL;
+ struct reinjection_step reinject_step = REINJECTION_STEP_NULL;
/* Fragment on the fluid side of the boundary */
struct sdis_interface_fragment frag_fluid;
@@ -145,6 +144,9 @@ XD(solid_fluid_boundary_picardN_path)
struct sdis_heat_vertex hvtx = SDIS_HEAT_VERTEX_NULL;
struct sdis_heat_vertex hvtx_s = SDIS_HEAT_VERTEX_NULL;
+ /* The enclosures split by the boundary */
+ unsigned enc_ids[2] = {ENCLOSURE_ID_NULL, ENCLOSURE_ID_NULL};
+
/* Data attached to the boundary */
struct sdis_interface* interf = NULL;
struct sdis_medium* solid = NULL;
@@ -172,8 +174,7 @@ XD(solid_fluid_boundary_picardN_path)
res_T res = RES_OK;
ASSERT(scn && rwalk && rng && T && ctx);
- ASSERT(XD(check_rwalk_fragment_consistency)(rwalk, frag));
-
+ ASSERT(XD(check_rwalk_fragment_consistency)(rwalk, frag) == RES_OK);
/* Fetch the Min/max temperature */
Tmin = ctx->Tmin;
@@ -184,7 +185,7 @@ XD(solid_fluid_boundary_picardN_path)
That3 = ctx->That3;
/* Retrieve the solid and the fluid split by the boundary */
- interf = scene_get_interface(scn, rwalk->hit.prim.prim_id);
+ interf = scene_get_interface(scn, rwalk->XD(hit).prim.prim_id);
solid = interface_get_medium(interf, SDIS_FRONT);
fluid = interface_get_medium(interf, SDIS_BACK);
solid_side = SDIS_FRONT;
@@ -195,6 +196,9 @@ XD(solid_fluid_boundary_picardN_path)
ASSERT(fluid->type == SDIS_FLUID);
}
+ /* Get the enclosures split by the boundary */
+ scene_get_enclosure_ids(scn, rwalk->XD(hit).prim.prim_id, enc_ids);
+
/* Check that no net flux is set for this interface since the provided
* picardN algorithm does not handle it */
res = check_net_flux(scn, interf, frag, get_picard_order(ctx));
@@ -219,7 +223,7 @@ XD(solid_fluid_boundary_picardN_path)
/* Sample a reinjection step */
samp_reinject_step_args.rng = rng;
- samp_reinject_step_args.solid = solid;
+ samp_reinject_step_args.solid_enc_id = enc_ids[solid_side];
samp_reinject_step_args.rwalk = rwalk;
samp_reinject_step_args.distance = delta_boundary;
samp_reinject_step_args.side = solid_side;
@@ -256,8 +260,8 @@ XD(solid_fluid_boundary_picardN_path)
/* Null collision main loop */
for(;;) {
/* Temperature and random walk state of the sampled radiative path */
- struct XD(temperature) T_s;
- struct XD(rwalk) rwalk_s;
+ struct temperature T_s;
+ struct rwalk rwalk_s;
double h_radi, h_radi_min, h_radi_max; /* Radiative coefficients */
double p_radi, p_radi_min, p_radi_max; /* Radiative probas */
@@ -272,15 +276,15 @@ XD(solid_fluid_boundary_picardN_path)
/* Switch in convective path */
if(r < p_conv) {
T->func = XD(convective_path);
- rwalk->mdm = fluid;
+ rwalk->enc_id = enc_ids[fluid_side];
rwalk->hit_side = fluid_side;
break;
}
/* Switch in conductive path */
if(r < p_conv + p_cond) {
- struct XD(solid_reinjection_args) solid_reinject_args =
- XD(SOLID_REINJECTION_ARGS_NULL);
+ struct solid_reinjection_args solid_reinject_args =
+ SOLID_REINJECTION_ARGS_NULL;
/* Perform the reinjection into the solid */
solid_reinject_args.reinjection = &reinject_step;
@@ -289,7 +293,7 @@ XD(solid_fluid_boundary_picardN_path)
solid_reinject_args.rng = rng;
solid_reinject_args.T = T;
solid_reinject_args.fp_to_meter = scn->fp_to_meter;
- res = XD(solid_reinjection)(solid, &solid_reinject_args);
+ res = XD(solid_reinjection)(scn, enc_ids[solid_side], &solid_reinject_args);
if(res != RES_OK) goto error;
break;
}
@@ -303,7 +307,7 @@ XD(solid_fluid_boundary_picardN_path)
/* Sample a radiative path */
T_s = *T;
rwalk_s = *rwalk;
- rwalk_s.mdm = fluid;
+ rwalk_s.enc_id = enc_ids[fluid_side];
rwalk_s.hit_side = fluid_side;
res = XD(radiative_path)(scn, ctx, &rwalk_s, rng, &T_s);
if(res != RES_OK) goto error;
@@ -340,9 +344,9 @@ XD(solid_fluid_boundary_picardN_path)
} (void)0
#define COMPUTE_TEMPERATURE(Result, RWalk, Temp) { \
- struct XD(temperature) T_p; \
+ struct temperature T_p; \
if((Temp)->done) { /* Ambient radiative temperature */ \
- ASSERT(SXD_HIT_NONE(&(RWalk)->hit)); \
+ ASSERT(SXD_HIT_NONE(&(RWalk)->XD(hit))); \
T_p = *(Temp); \
} else { \
res = XD(sample_path)(scn, RWalk, ctx, rng, &T_p); \
diff --git a/src/sdis_heat_path_boundary_Xd_solid_solid.h b/src/sdis_heat_path_boundary_Xd_solid_solid.h
@@ -33,28 +33,29 @@ XD(solid_solid_boundary_path)
(struct sdis_scene* scn,
struct rwalk_context* ctx,
const struct sdis_interface_fragment* frag,
- struct XD(rwalk)* rwalk,
+ struct rwalk* rwalk,
struct ssp_rng* rng,
- struct XD(temperature)* T)
+ struct temperature* T)
{
/* Input/output arguments of the function used to sample a reinjection */
- struct XD(sample_reinjection_step_args) samp_reinject_step_frt_args =
- XD(SAMPLE_REINJECTION_STEP_ARGS_NULL);
- struct XD(sample_reinjection_step_args) samp_reinject_step_bck_args =
- XD(SAMPLE_REINJECTION_STEP_ARGS_NULL);
- struct XD(reinjection_step) reinject_step_frt = XD(REINJECTION_STEP_NULL);
- struct XD(reinjection_step) reinject_step_bck = XD(REINJECTION_STEP_NULL);
- struct XD(reinjection_step)* reinject_step = NULL;
+ struct sample_reinjection_step_args samp_reinject_step_frt_args =
+ SAMPLE_REINJECTION_STEP_ARGS_NULL;
+ struct sample_reinjection_step_args samp_reinject_step_bck_args =
+ SAMPLE_REINJECTION_STEP_ARGS_NULL;
+ struct reinjection_step reinject_step_frt = REINJECTION_STEP_NULL;
+ struct reinjection_step reinject_step_bck = REINJECTION_STEP_NULL;
+ struct reinjection_step* reinject_step = NULL;
/* Reinjection arguments */
- struct XD(solid_reinjection_args) solid_reinject_args =
- XD(SOLID_REINJECTION_ARGS_NULL);
+ struct solid_reinjection_args solid_reinject_args =
+ SOLID_REINJECTION_ARGS_NULL;
/* Data attached to the boundary */
+ unsigned enc_ids[2] = {ENCLOSURE_ID_NULL, ENCLOSURE_ID_NULL};
struct sdis_interface* interf = NULL;
struct sdis_medium* solid_frt = NULL;
struct sdis_medium* solid_bck = NULL;
- struct sdis_medium* solid = NULL;
+ unsigned solid_enc_id = ENCLOSURE_ID_NULL; /* Solid to re-inject */
double lambda_frt;
double lambda_bck;
@@ -67,11 +68,12 @@ XD(solid_solid_boundary_path)
res_T res = RES_OK;
ASSERT(scn && ctx && frag && rwalk && rng && T);
- ASSERT(XD(check_rwalk_fragment_consistency)(rwalk, frag));
+ ASSERT(XD(check_rwalk_fragment_consistency)(rwalk, frag) == RES_OK);
(void)frag, (void)ctx;
- /* Retrieve the two solids split by the boundary */
- interf = scene_get_interface(scn, rwalk->hit.prim.prim_id);
+ /* Retrieve the two enclosures and associated media split by the boundary */
+ scene_get_enclosure_ids(scn, rwalk->XD(hit).prim.prim_id, enc_ids);
+ interf = scene_get_interface(scn, rwalk->XD(hit).prim.prim_id);
solid_frt = interface_get_medium(interf, SDIS_FRONT);
solid_bck = interface_get_medium(interf, SDIS_BACK);
ASSERT(solid_frt->type == SDIS_SOLID);
@@ -93,8 +95,8 @@ XD(solid_solid_boundary_path)
/* Sample a front/back reinjection steps */
samp_reinject_step_frt_args.rng = rng;
samp_reinject_step_bck_args.rng = rng;
- samp_reinject_step_frt_args.solid = solid_frt;
- samp_reinject_step_bck_args.solid = solid_bck;
+ samp_reinject_step_frt_args.solid_enc_id = enc_ids[SDIS_FRONT];
+ samp_reinject_step_bck_args.solid_enc_id = enc_ids[SDIS_BACK];
samp_reinject_step_frt_args.rwalk = rwalk;
samp_reinject_step_bck_args.rwalk = rwalk;
samp_reinject_step_frt_args.distance = delta_boundary_frt;
@@ -102,7 +104,7 @@ XD(solid_solid_boundary_path)
samp_reinject_step_frt_args.side = SDIS_FRONT;
samp_reinject_step_bck_args.side = SDIS_BACK;
res = XD(sample_reinjection_step_solid_solid)
- (scn,
+ (scn,
&samp_reinject_step_frt_args,
&samp_reinject_step_bck_args,
&reinject_step_frt,
@@ -150,10 +152,10 @@ XD(solid_solid_boundary_path)
if(r < proba) { /* Reinject in front */
reinject_step = &reinject_step_frt;
- solid = solid_frt;
+ solid_enc_id = enc_ids[SDIS_FRONT];
} else { /* Reinject in back */
reinject_step = &reinject_step_bck;
- solid = solid_bck;
+ solid_enc_id = enc_ids[SDIS_BACK];
}
/* Perform the reinjection into the solid */
@@ -163,7 +165,7 @@ XD(solid_solid_boundary_path)
solid_reinject_args.rwalk_ctx = ctx;
solid_reinject_args.T = T;
solid_reinject_args.fp_to_meter = scn->fp_to_meter;
- res = XD(solid_reinjection)(solid, &solid_reinject_args);
+ res = XD(solid_reinjection)(scn, solid_enc_id, &solid_reinject_args);
if(res != RES_OK) goto error;
exit:
diff --git a/src/sdis_heat_path_boundary_c.h b/src/sdis_heat_path_boundary_c.h
@@ -21,123 +21,92 @@
#include <rsys/rsys.h>
/* Forward declarations */
-struct rwalk_2d;
-struct rwalk_3d;
+struct rwalk;
struct sdis_scene;
struct sdis_medium;
/*******************************************************************************
* Sample a reinjection step
******************************************************************************/
-struct sample_reinjection_step_args_2d {
+struct sample_reinjection_step_args {
struct ssp_rng* rng; /* Random number generator to use */
- const struct sdis_medium* solid; /* Solid in which to reinject */
- struct rwalk_2d* rwalk; /* Current state of the random walk */
+ struct rwalk* rwalk; /* Current state of the random walk */
double distance; /* Maximum Reinjection distance */
+ unsigned solid_enc_id; /* Enclosured Id of the solid in which to reinject */
enum sdis_side side; /* Side of the boundary to re-inject */
};
-struct sample_reinjection_step_args_3d {
- struct ssp_rng* rng; /* Random number generator to use */
- const struct sdis_medium* solid; /* Medium in which to reinject */
- struct rwalk_3d* rwalk; /* Current random walk state */
- double distance; /* Maximum Reinjection distance */
- enum sdis_side side; /* Side of the boundary to re-inject */
-};
-
-struct reinjection_step_2d {
- struct s2d_hit hit; /* Intersection along the reinjection direction */
- float direction[2]; /* Reinjection direction */
- float distance; /* Reinjection distance */
-};
+#define SAMPLE_REINJECTION_STEP_ARGS_NULL__ \
+ {NULL, NULL, -1, ENCLOSURE_ID_NULL, SDIS_SIDE_NULL__}
+static const struct sample_reinjection_step_args
+SAMPLE_REINJECTION_STEP_ARGS_NULL = SAMPLE_REINJECTION_STEP_ARGS_NULL__;
-struct reinjection_step_3d {
- struct s3d_hit hit; /* Intersection along the reinjection direction */
+struct reinjection_step {
+ struct s2d_hit hit_2d; /* 2D Intersection along the reinjection direction */
+ struct s3d_hit hit_3d; /* 3D Intersection along the reinjection direction */
float direction[3]; /* Reinjection direction */
float distance; /* Reinjection distance */
};
-#define SAMPLE_REINJECTION_STEP_ARGS_NULL___2d \
- {NULL, NULL, NULL, -1, SDIS_SIDE_NULL__}
-#define SAMPLE_REINJECTION_STEP_ARGS_NULL___3d \
- {NULL, NULL, NULL, -1, SDIS_SIDE_NULL__}
-static const struct sample_reinjection_step_args_2d
-SAMPLE_REINJECTION_STEP_ARGS_NULL_2d = SAMPLE_REINJECTION_STEP_ARGS_NULL___2d;
-static const struct sample_reinjection_step_args_3d
-SAMPLE_REINJECTION_STEP_ARGS_NULL_3d = SAMPLE_REINJECTION_STEP_ARGS_NULL___3d;
-
-#define REINJECTION_STEP_NULL___2d {S2D_HIT_NULL__, {0,0}, 0}
-#define REINJECTION_STEP_NULL___3d {S3D_HIT_NULL__, {0,0,0}, 0}
-static const struct reinjection_step_2d
-REINJECTION_STEP_NULL_2d = REINJECTION_STEP_NULL___2d;
-static const struct reinjection_step_3d
-REINJECTION_STEP_NULL_3d = REINJECTION_STEP_NULL___3d;
+#define REINJECTION_STEP_NULL__ {S2D_HIT_NULL__, S3D_HIT_NULL__, {0,0,0}, 0}
+static const struct reinjection_step REINJECTION_STEP_NULL =
+ REINJECTION_STEP_NULL__;
extern LOCAL_SYM res_T
sample_reinjection_step_solid_fluid_2d
(struct sdis_scene* scn,
- const struct sample_reinjection_step_args_2d* args,
- struct reinjection_step_2d* step);
+ const struct sample_reinjection_step_args* args,
+ struct reinjection_step* step);
extern LOCAL_SYM res_T
sample_reinjection_step_solid_fluid_3d
(struct sdis_scene* scn,
- const struct sample_reinjection_step_args_3d* args,
- struct reinjection_step_3d *step);
+ const struct sample_reinjection_step_args* args,
+ struct reinjection_step *step);
extern LOCAL_SYM res_T
sample_reinjection_step_solid_solid_2d
(struct sdis_scene* scn,
- const struct sample_reinjection_step_args_2d* args_front,
- const struct sample_reinjection_step_args_2d* args_back,
- struct reinjection_step_2d* step_front,
- struct reinjection_step_2d* step_back);
+ const struct sample_reinjection_step_args* args_front,
+ const struct sample_reinjection_step_args* args_back,
+ struct reinjection_step* step_front,
+ struct reinjection_step* step_back);
extern LOCAL_SYM res_T
sample_reinjection_step_solid_solid_3d
(struct sdis_scene* scn,
- const struct sample_reinjection_step_args_3d* args_front,
- const struct sample_reinjection_step_args_3d* args_back,
- struct reinjection_step_3d* step_front,
- struct reinjection_step_3d* step_back);
+ const struct sample_reinjection_step_args* args_front,
+ const struct sample_reinjection_step_args* args_back,
+ struct reinjection_step* step_front,
+ struct reinjection_step* step_back);
/*******************************************************************************
* Reinject the random walk into a solid
******************************************************************************/
-struct solid_reinjection_args_2d {
- const struct reinjection_step_2d* reinjection; /* Reinjection to do */
+struct solid_reinjection_args {
+ const struct reinjection_step* reinjection; /* Reinjection to do */
struct rwalk_context* rwalk_ctx;
- struct rwalk_2d* rwalk; /* Current state of the random walk */
+ struct rwalk* rwalk; /* Current state of the random walk */
struct ssp_rng* rng; /* Random number generator */
- struct temperature_2d* T;
+ struct temperature* T;
double fp_to_meter;
};
-struct solid_reinjection_args_3d {
- const struct reinjection_step_3d* reinjection; /* Reinjection to do */
- struct rwalk_context* rwalk_ctx;
- struct rwalk_3d* rwalk; /* Current state of the random walk */
- struct ssp_rng* rng; /* Random number generator */
- struct temperature_3d* T;
- double fp_to_meter;
-};
-
-#define SOLID_REINJECTION_ARGS_NULL___2d {NULL,NULL,NULL,NULL,NULL,0}
-#define SOLID_REINJECTION_ARGS_NULL___3d {NULL,NULL,NULL,NULL,NULL,0}
-static const struct solid_reinjection_args_2d SOLID_REINJECTION_ARGS_NULL_2d =
- SOLID_REINJECTION_ARGS_NULL___2d;
-static const struct solid_reinjection_args_3d SOLID_REINJECTION_ARGS_NULL_3d =
- SOLID_REINJECTION_ARGS_NULL___3d;
+#define SOLID_REINJECTION_ARGS_NULL__ {NULL,NULL,NULL,NULL,NULL,0}
+static const struct solid_reinjection_args SOLID_REINJECTION_ARGS_NULL =
+ SOLID_REINJECTION_ARGS_NULL__;
extern LOCAL_SYM res_T
solid_reinjection_2d
- (struct sdis_medium* solid,
- struct solid_reinjection_args_2d* args);
+ (struct sdis_scene* scn,
+ const unsigned solid_enc_id,
+ struct solid_reinjection_args* args);
extern LOCAL_SYM res_T
solid_reinjection_3d
- (struct sdis_medium* solid,
- struct solid_reinjection_args_3d* args);
+ (struct sdis_scene* scn,
+ const unsigned solid_enc_id,
+ struct solid_reinjection_args* args);
/*******************************************************************************
* Handle net flux
@@ -160,21 +129,22 @@ extern LOCAL_SYM res_T
handle_net_flux_2d
(const struct sdis_scene* scn,
const struct handle_net_flux_args* args,
- struct temperature_2d* T);
+ struct temperature* T);
extern LOCAL_SYM res_T
handle_net_flux_3d
(const struct sdis_scene* scn,
const struct handle_net_flux_args* args,
- struct temperature_3d* T);
+ struct temperature* T);
/*******************************************************************************
* Handle external flux
******************************************************************************/
-struct handle_external_net_flux_args_2d {
+struct handle_external_net_flux_args {
struct sdis_interface* interf;
const struct sdis_interface_fragment* frag;
- const struct s2d_hit* hit;
+ const struct s2d_hit* hit_2d;
+ const struct s3d_hit* hit_3d;
struct green_path_handle* green_path; /* Store the propagator */
struct sdis_heat_path* heat_path; /* Save paths */
@@ -199,26 +169,23 @@ struct handle_external_net_flux_args_3d {
double h_radi; /* Radiative coefficient */
};
-#define HANDLE_EXTERNAL_NET_FLUX_ARGS_NULL___2d {NULL,NULL,NULL,NULL,NULL,0,0,0,0}
-#define HANDLE_EXTERNAL_NET_FLUX_ARGS_NULL___3d {NULL,NULL,NULL,NULL,NULL,0,0,0,0}
-static const struct handle_external_net_flux_args_2d
-HANDLE_EXTERNAL_NET_FLUX_ARGS_NULL_2d = HANDLE_EXTERNAL_NET_FLUX_ARGS_NULL___2d;
-static const struct handle_external_net_flux_args_3d
-HANDLE_EXTERNAL_NET_FLUX_ARGS_NULL_3d = HANDLE_EXTERNAL_NET_FLUX_ARGS_NULL___3d;
+#define HANDLE_EXTERNAL_NET_FLUX_ARGS_NULL__ {NULL,NULL,NULL,NULL,NULL,NULL,0,0,0,0}
+static const struct handle_external_net_flux_args
+HANDLE_EXTERNAL_NET_FLUX_ARGS_NULL = HANDLE_EXTERNAL_NET_FLUX_ARGS_NULL__;
extern LOCAL_SYM res_T
handle_external_net_flux_2d
(const struct sdis_scene* scn,
struct ssp_rng* rng,
- const struct handle_external_net_flux_args_2d* args,
- struct temperature_2d* T);
+ const struct handle_external_net_flux_args* args,
+ struct temperature* T);
extern LOCAL_SYM res_T
handle_external_net_flux_3d
(const struct sdis_scene* scn,
struct ssp_rng* rng,
- const struct handle_external_net_flux_args_3d* args,
- struct temperature_3d* T);
+ const struct handle_external_net_flux_args* args,
+ struct temperature* T);
/*******************************************************************************
* Miscellaneous functions
@@ -246,9 +213,9 @@ solid_boundary_with_flux_path_2d
struct rwalk_context* ctx,
const struct sdis_interface_fragment* frag,
const double phi,
- struct rwalk_2d* rwalk,
+ struct rwalk* rwalk,
struct ssp_rng* rng,
- struct temperature_2d* T);
+ struct temperature* T);
extern LOCAL_SYM res_T
solid_boundary_with_flux_path_3d
@@ -256,62 +223,62 @@ solid_boundary_with_flux_path_3d
struct rwalk_context* ctx,
const struct sdis_interface_fragment* frag,
const double phi,
- struct rwalk_3d* rwalk,
+ struct rwalk* rwalk,
struct ssp_rng* rng,
- struct temperature_3d* T);
+ struct temperature* T);
extern LOCAL_SYM res_T
solid_fluid_boundary_picard1_path_2d
(struct sdis_scene* scn,
struct rwalk_context* ctx,
const struct sdis_interface_fragment* frag,
- struct rwalk_2d* rwalk,
+ struct rwalk* rwalk,
struct ssp_rng* rng,
- struct temperature_2d* T);
+ struct temperature* T);
extern LOCAL_SYM res_T
solid_fluid_boundary_picard1_path_3d
(struct sdis_scene* scn,
struct rwalk_context* ctx,
const struct sdis_interface_fragment* frag,
- struct rwalk_3d* rwalk,
+ struct rwalk* rwalk,
struct ssp_rng* rng,
- struct temperature_3d* T);
+ struct temperature* T);
extern LOCAL_SYM res_T
solid_fluid_boundary_picardN_path_2d
(struct sdis_scene* scn,
struct rwalk_context* ctx,
const struct sdis_interface_fragment* frag,
- struct rwalk_2d* rwalk,
+ struct rwalk* rwalk,
struct ssp_rng* rng,
- struct temperature_2d* T);
+ struct temperature* T);
extern LOCAL_SYM res_T
solid_fluid_boundary_picardN_path_3d
(struct sdis_scene* scn,
struct rwalk_context* ctx,
const struct sdis_interface_fragment* frag,
- struct rwalk_3d* rwalk,
+ struct rwalk* rwalk,
struct ssp_rng* rng,
- struct temperature_3d* T);
+ struct temperature* T);
extern LOCAL_SYM res_T
solid_solid_boundary_path_2d
(struct sdis_scene* scn,
struct rwalk_context* ctx,
const struct sdis_interface_fragment* frag,
- struct rwalk_2d* rwalk,
+ struct rwalk* rwalk,
struct ssp_rng* rng,
- struct temperature_2d* T);
+ struct temperature* T);
extern LOCAL_SYM res_T
solid_solid_boundary_path_3d
(struct sdis_scene* scn,
struct rwalk_context* ctx,
const struct sdis_interface_fragment* frag,
- struct rwalk_3d* rwalk,
+ struct rwalk* rwalk,
struct ssp_rng* rng,
- struct temperature_3d* T);
+ struct temperature* T);
#endif /* SDIS_HEAT_PATH_BOUNDARY_C_H */
diff --git a/src/sdis_heat_path_conductive.c b/src/sdis_heat_path_conductive.c
@@ -75,9 +75,9 @@ res_T
conductive_path_2d
(struct sdis_scene* scn,
struct rwalk_context* ctx,
- struct rwalk_2d* rwalk,
+ struct rwalk* rwalk,
struct ssp_rng* rng,
- struct temperature_2d* T)
+ struct temperature* T)
{
res_T res = RES_OK;
ASSERT(ctx);
@@ -98,9 +98,9 @@ res_T
conductive_path_3d
(struct sdis_scene* scn,
struct rwalk_context* ctx,
- struct rwalk_3d* rwalk,
+ struct rwalk* rwalk,
struct ssp_rng* rng,
- struct temperature_3d* T)
+ struct temperature* T)
{
res_T res = RES_OK;
ASSERT(ctx);
diff --git a/src/sdis_heat_path_conductive_c.h b/src/sdis_heat_path_conductive_c.h
@@ -19,15 +19,13 @@
#include <rsys/rsys.h>
/* Forward declarations */
-struct rwalk_2d;
-struct rwalk_3d;
+struct rwalk;
struct rwalk_context;
struct sdis_device;
struct sdis_scene;
struct solid_props;
struct ssp_rng;
-struct temperature_2d;
-struct temperature_3d;
+struct temperature;
extern LOCAL_SYM res_T
check_solid_constant_properties
@@ -44,17 +42,17 @@ extern LOCAL_SYM res_T
conductive_path_delta_sphere_2d
(struct sdis_scene* scn,
struct rwalk_context* ctx,
- struct rwalk_2d* rwalk,
+ struct rwalk* rwalk,
struct ssp_rng* rng,
- struct temperature_2d* T);
+ struct temperature* T);
extern LOCAL_SYM res_T
conductive_path_delta_sphere_3d
(struct sdis_scene* scn,
struct rwalk_context* ctx,
- struct rwalk_3d* rwalk,
+ struct rwalk* rwalk,
struct ssp_rng* rng,
- struct temperature_3d* T);
+ struct temperature* T);
/*******************************************************************************
* Conductive paths using the walk on sphere algorithm
@@ -63,16 +61,16 @@ extern LOCAL_SYM res_T
conductive_path_wos_2d
(struct sdis_scene* scn,
struct rwalk_context* ctx,
- struct rwalk_2d* rwalk,
+ struct rwalk* rwalk,
struct ssp_rng* rng,
- struct temperature_2d* T);
+ struct temperature* T);
extern LOCAL_SYM res_T
conductive_path_wos_3d
(struct sdis_scene* scn,
struct rwalk_context* ctx,
- struct rwalk_3d* rwalk,
+ struct rwalk* rwalk,
struct ssp_rng* rng,
- struct temperature_3d* T);
+ struct temperature* T);
#endif /* SDIS_HEAT_PATH_CONDUCTIVE_C_H */
diff --git a/src/sdis_heat_path_conductive_delta_sphere_Xd.h b/src/sdis_heat_path_conductive_delta_sphere_Xd.h
@@ -111,7 +111,7 @@ XD(sample_next_step)
static res_T
XD(sample_next_step_robust)
(struct sdis_scene* scn,
- struct sdis_medium* current_mdm,
+ const unsigned current_enc_id,
struct ssp_rng* rng,
const double pos[DIM],
const float delta_solid,
@@ -121,13 +121,14 @@ XD(sample_next_step_robust)
struct sXd(hit)* hit1, /* Hit used to adjust delta */
float* out_delta)
{
- struct sdis_medium* mdm;
+ unsigned enc_id = ENCLOSURE_ID_NULL;
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(scn && rng && pos && delta_solid > 0);
+ ASSERT(current_enc_id != ENCLOSURE_ID_NULL);
ASSERT(dir0 && dir1 && hit0 && hit1 && out_delta);
fX_set_dX(org, pos);
@@ -141,44 +142,30 @@ XD(sample_next_step_robust)
/* 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_BAD_OP) { mdm = NULL; res = RES_OK; }
+ res = scene_get_enclosure_id_in_closed_boundaries(scn, pos_next, &enc_id);
+ if(res == RES_BAD_OP) { enc_id = ENCLOSURE_ID_NULL; res = RES_OK; }
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);
+ unsigned enc_ids[2] = {ENCLOSURE_ID_NULL, ENCLOSURE_ID_NULL};
+ scene_get_enclosure_ids(scn, hit0->prim.prim_id, enc_ids);
+ enc_id = fX(dot)(dir0, hit0->normal) < 0 ? enc_ids[0] : enc_ids[1];
}
/* Check medium consistency */
- if(current_mdm != mdm) {
+ if(current_enc_id != enc_id) {
#if 0
-#if DIM == 2
- log_warn(scn->dev,
- "%s: inconsistent medium during the solid random walk at {%g, %g}.\n",
- FUNC_NAME, SPLIT2(pos));
-#else
log_warn(scn->dev,
- "%s: inconsistent medium during the solid random walk at {%g, %g, %g}.\n",
- FUNC_NAME, SPLIT3(pos));
-#endif
+ "%s: inconsistent medium during the solid random walk -- "
+ "pos=("FORMAT_VECX")\n", FUNC_NAME, SPLITX(pos));
#endif
}
- } while(current_mdm != mdm && ++iattempt < MAX_ATTEMPTS);
+ } while(current_enc_id != enc_id && ++iattempt < MAX_ATTEMPTS);
/* Handle error */
if(iattempt >= MAX_ATTEMPTS) {
-#if DIM == 2
- log_warn(scn->dev,
- "%s: could not find a next valid conductive step at {%g, %g}.\n",
- FUNC_NAME, SPLIT2(pos));
-#else
log_warn(scn->dev,
- "%s: could not find a next valid conductive step at {%g, %g, %g}.\n",
- FUNC_NAME, SPLIT3(pos));
-#endif
+ "%s: could not find a next valid conductive -- pos=("FORMAT_VECX")\n",
+ FUNC_NAME, SPLITX(pos));
res = RES_BAD_OP;
goto error;
}
@@ -239,7 +226,7 @@ XD(handle_volumic_power)
(const struct sdis_scene* scn,
const struct XD(handle_volumic_power_args)* args,
double* out_power_term,
- struct XD(temperature)* T)
+ struct temperature* T)
{
double power_term = 0;
res_T res = RES_OK;
@@ -332,28 +319,43 @@ res_T
XD(conductive_path_delta_sphere)
(struct sdis_scene* scn,
struct rwalk_context* ctx,
- struct XD(rwalk)* rwalk,
+ struct rwalk* rwalk,
struct ssp_rng* rng,
- struct XD(temperature)* T)
+ struct temperature* T)
{
- double position_start[DIM];
+ /* Enclosure/medium in which the conductive path starts */
+ struct sdis_medium* mdm = NULL;
+ unsigned enc_id = ENCLOSURE_ID_NULL;
+
+ /* Physical properties */
struct solid_props props_ref = SOLID_PROPS_NULL;
double green_power_term = 0;
- struct sdis_medium* mdm;
+
+ /* Miscellaneous */
+ double position_start[DIM] = {0};
size_t istep = 0; /* Help for debug */
res_T res = RES_OK;
+
+ /* Check pre-conditions */
ASSERT(scn && rwalk && rng && T);
- ASSERT(rwalk->mdm->type == SDIS_SOLID);
- (void)ctx, (void)istep;
+
+ (void)ctx, (void)istep; /* Avoid "unsued variable" warnings */
+
+ res = scene_get_enclosure_id_in_closed_boundaries(scn, rwalk->vtx.P, &enc_id);
+ if(res != RES_OK) goto error;
+ res = scene_get_enclosure_medium(scn, scene_get_enclosure(scn, enc_id), &mdm);
+ if(res != RES_OK) goto error;
+ ASSERT(sdis_medium_get_type(mdm) == SDIS_SOLID);
/* Check the random walk consistency */
- res = scene_get_medium_in_closed_boundaries(scn, rwalk->vtx.P, &mdm);
- if(res != RES_OK || mdm != rwalk->mdm) {
+ if(enc_id != rwalk->enc_id) {
log_err(scn->dev, "%s: invalid solid random walk. "
- "Unexpected medium at {%g, %g, %g}.\n", FUNC_NAME, SPLIT3(rwalk->vtx.P));
+ "Unexpected enclosure -- pos=("FORMAT_VECX")\n",
+ FUNC_NAME, SPLITX(rwalk->vtx.P));
res = RES_BAD_OP_IRRECOVERABLE;
goto error;
}
+
/* Save the submitted position */
dX(set)(position_start, rwalk->vtx.P);
@@ -396,7 +398,7 @@ XD(conductive_path_delta_sphere)
if(ctx->green_path) {
res = green_path_set_limit_vertex
- (ctx->green_path, rwalk->mdm, &rwalk->vtx, rwalk->elapsed_time);
+ (ctx->green_path, mdm, &rwalk->vtx, rwalk->elapsed_time);
if(res != RES_OK) goto error;
}
@@ -410,7 +412,7 @@ XD(conductive_path_delta_sphere)
fX_set_dX(org, rwalk->vtx.P);
/* Sample the direction to walk toward and compute the distance to travel */
- res = XD(sample_next_step_robust)(scn, mdm, rng, rwalk->vtx.P,
+ res = XD(sample_next_step_robust)(scn, enc_id, rng, rwalk->vtx.P,
(float)props.delta, dir0, dir1, &hit0, &hit1, &delta);
if(res != RES_OK) goto error;
@@ -436,16 +438,16 @@ XD(conductive_path_delta_sphere)
/* Rewind the time */
delta_m = delta * scn->fp_to_meter;
mu = (2*DIM*props.lambda)/(props.rho*props.cp*delta_m*delta_m);
- res = XD(time_rewind)(mu, props.t0, rng, rwalk, ctx, T);
+ res = time_rewind(scn, mu, props.t0, rng, rwalk, ctx, T);
if(res != RES_OK) goto error;
if(T->done) break; /* Limit condition was reached */
/* Define if the random walk hits something along dir0 */
if(hit0.distance > delta) {
- rwalk->hit = SXD_HIT_NULL;
+ rwalk->XD(hit) = SXD_HIT_NULL;
rwalk->hit_side = SDIS_SIDE_NULL__;
} else {
- rwalk->hit = hit0;
+ rwalk->XD(hit) = hit0;
rwalk->hit_side = fX(dot)(hit0.normal, dir0) < 0 ? SDIS_FRONT : SDIS_BACK;
}
@@ -460,17 +462,17 @@ XD(conductive_path_delta_sphere)
++istep;
/* Keep going while the solid random walk does not hit an interface */
- } while(SXD_HIT_NONE(&rwalk->hit));
+ } while(SXD_HIT_NONE(&rwalk->XD(hit)));
/* Register the power term for the green function */
if(ctx->green_path && props_ref.power != SDIS_VOLUMIC_POWER_NONE) {
res = green_path_add_power_term
- (ctx->green_path, rwalk->mdm, &rwalk->vtx, green_power_term);
+ (ctx->green_path, mdm, &rwalk->vtx, green_power_term);
if(res != RES_OK) goto error;
}
T->func = XD(boundary_path);
- rwalk->mdm = NULL; /* The random walk is at an interface between 2 media */
+ rwalk->enc_id = ENCLOSURE_ID_NULL; /* At the interface between 2 media */
exit:
return res;
diff --git a/src/sdis_heat_path_conductive_wos_Xd.h b/src/sdis_heat_path_conductive_wos_Xd.h
@@ -23,24 +23,67 @@
#include "sdis_Xd_begin.h"
/*******************************************************************************
+ * Non generic helper functions
+ ******************************************************************************/
+#ifndef SDIS_HEAT_PATH_CONDUCTIVE_WOS_XD_H
+#define SDIS_HEAT_PATH_CONDUCTIVE_WOS_XD_H
+
+static res_T
+update_green_path
+ (struct green_path_handle* green_path,
+ struct rwalk* rwalk,
+ struct sdis_medium* mdm,
+ const struct solid_props* props,
+ const double power_term,
+ const struct temperature* T)
+{
+ res_T res = RES_OK;
+ ASSERT(mdm && props && T);
+
+ /* Is the green function estimated? */
+ if(!green_path) goto exit;
+
+ /* Save power term for green function if any */
+ if(props->power != SDIS_VOLUMIC_POWER_NONE) {
+ res = green_path_add_power_term(green_path, mdm, &rwalk->vtx, power_term);
+ if(res != RES_OK) goto error;
+ }
+
+ /* Set the green path limit to the current position if the initial condition
+ * has been reached */
+ if(T->done) {
+ res = green_path_set_limit_vertex
+ (green_path, mdm, &rwalk->vtx, rwalk->elapsed_time);
+ if(res != RES_OK) goto error;
+ }
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
+#endif /* SDIS_HEAT_PATH_CONDUCTIVE_WOS_XD_H */
+
+/*******************************************************************************
* Helper function
******************************************************************************/
static res_T
-XD(check_medium_consistency)
+XD(check_enclosure_consistency)
(struct sdis_scene* scn,
- const struct XD(rwalk)* rwalk)
+ const struct rwalk* rwalk)
{
- struct sdis_medium* mdm = NULL;
+ unsigned enc_id = ENCLOSURE_ID_NULL;
res_T res = RES_OK;
ASSERT(rwalk);
- res = scene_get_medium_in_closed_boundaries(scn, rwalk->vtx.P, &mdm);
+ res = scene_get_enclosure_id_in_closed_boundaries(scn, rwalk->vtx.P, &enc_id);
if(res != RES_OK) goto error;
- /* Check medium consistency */
- if(mdm != rwalk->mdm) {
+ /* Check enclosure consistency */
+ if(enc_id != rwalk->enc_id) {
log_err(scn->dev,
- "%s:%s: invalid solid walk. Unexpected medium (position: "FORMAT_VECX").\n",
+ "%s:%s: invalid solid walk. Unexpected enclosure -- pos=("FORMAT_VECX")\n",
__FILE__, FUNC_NAME, SPLITX(rwalk->vtx.P));
res = RES_BAD_OP_IRRECOVERABLE;
goto error;
@@ -55,12 +98,13 @@ error:
static res_T
XD(time_travel)
(struct sdis_scene* scn,
- struct XD(rwalk)* rwalk,
+ struct rwalk* rwalk,
struct ssp_rng* rng,
+ struct sdis_medium* mdm,
const double alpha, /* Diffusivity, i.e. lambda/(rho*cp) */
const double t0, /* Initial time [s] */
double* distance, /* Displacement [m/fp_to_meter] */
- struct XD(temperature)* T)
+ struct temperature* T)
{
double dir[DIM] = {0};
double dst = 0; /* Distance [m] */
@@ -120,14 +164,14 @@ XD(time_travel)
dX(add)(rwalk->vtx.P, rwalk->vtx.P, dir);
/* Fetch the initial temperature */
- temperature = medium_get_temperature(rwalk->mdm, &rwalk->vtx);
+ temperature = medium_get_temperature(mdm, &rwalk->vtx);
if(SDIS_TEMPERATURE_IS_UNKNOWN(temperature)) {
log_err(scn->dev,
"%s:%s: the path reaches the initial condition but the "
- "%s temperature remains unknown -- position=%g, %g, %g\n",
+ "%s temperature remains unknown -- pos=("FORMAT_VECX")\n",
__FILE__, FUNC_NAME,
- medium_type_to_string(sdis_medium_get_type(rwalk->mdm)),
- SPLIT3(rwalk->vtx.P));
+ medium_type_to_string(sdis_medium_get_type(mdm)),
+ SPLITX(rwalk->vtx.P));
res = RES_BAD_ARG;
goto error;
}
@@ -142,6 +186,32 @@ error:
goto exit;
}
+static res_T
+XD(handle_volumic_power_wos)
+ (struct sdis_scene* scn,
+ const struct solid_props* props,
+ const double distance, /* [m/fp_to_meter] */
+ double* power_term,
+ struct temperature* T)
+{
+ double dst = distance * scn->fp_to_meter; /* [m] */
+ double term = 0;
+ res_T res = RES_OK;
+ ASSERT(scn && props && distance >= 0 && power_term && T);
+
+ if(props->power == SDIS_VOLUMIC_POWER_NONE) goto exit;
+
+ /* No displacement => no power density */
+ if(distance == 0) goto exit;
+
+ term = dst*dst / (2*DIM*props->lambda);
+ T->value += props->power * term;
+
+exit:
+ *power_term = term;
+ return res;
+}
+
#if DIM == 2
static INLINE enum sdis_side
compute_hit_side_2d
@@ -203,16 +273,15 @@ compute_hit_side_3d
}
#endif
-/* Verify that the submitted position is in the expected medium */
+/* Verify that the submitted position is in the expected enclosure */
static res_T
XD(check_diffusion_position)
(struct sdis_scene* scn,
- const struct sdis_medium* expected_medium,
+ const unsigned expected_enc_id,
const double pos[DIM])
{
- struct sdis_interface* interf = NULL;
- struct sdis_medium* mdm = NULL;
- enum sdis_side side;
+ unsigned enc_ids[2] = {ENCLOSURE_ID_NULL, ENCLOSURE_ID_NULL};
+ enum sdis_side side = SDIS_SIDE_NULL__;
struct sXd(hit) hit = SXD_HIT_NULL;
float wos_pos[DIM] = {0};
@@ -220,7 +289,8 @@ XD(check_diffusion_position)
res_T res = RES_OK;
/* Check pre-conditions */
- ASSERT(scn && expected_medium && pos);
+ ASSERT(scn && pos);
+ ASSERT(expected_enc_id != ENCLOSURE_ID_NULL);
/* Look for the nearest surface within 1 mm of the position to be checked. By
* limiting the search radius we speed up the closest point query. If no
@@ -239,11 +309,10 @@ XD(check_diffusion_position)
SXD(scene_view_closest_point(scn->sXd(view), wos_pos, wos_radius, NULL, &hit));
if(SXD_HIT_NONE(&hit)) goto exit;
- /* Fetch interface properties and check path consistency */
- interf = scene_get_interface(scn, hit.prim.prim_id);
+ /* Check path consistency */
+ scene_get_enclosure_ids(scn, hit.prim.prim_id, enc_ids);
side = XD(compute_hit_side)(&hit, pos);
- mdm = side == SDIS_FRONT ? interf->medium_front : interf->medium_back;
- if(mdm != expected_medium) {
+ if(enc_ids[side] != expected_enc_id) {
res = RES_BAD_ARG;
goto error;
}
@@ -258,15 +327,14 @@ static res_T
XD(setup_hit_wos)
(struct sdis_scene* scn,
const struct sXd(hit)* hit,
- struct XD(rwalk)* rwalk)
+ struct rwalk* rwalk)
{
/* Geometry */
struct sXd(primitive) prim;
struct sXd(attrib) attr;
/* Properties */
- struct sdis_interface* interf = NULL;
- struct sdis_medium* mdm = NULL;
+ unsigned enc_ids[2] = {ENCLOSURE_ID_NULL, ENCLOSURE_ID_NULL};
enum sdis_side side = SDIS_SIDE_NULL__;
/* Miscellaneous */
@@ -288,13 +356,12 @@ XD(setup_hit_wos)
dX_set_fX(tgt, attr.value);
side = XD(compute_hit_side)(hit, rwalk->vtx.P);
- /* Fetch interface properties and check path consistency */
- interf = scene_get_interface(scn, hit->prim.prim_id);
- mdm = side == SDIS_FRONT ? interf->medium_front : interf->medium_back;
- if(mdm != rwalk->mdm) {
+ /* Check path consistency */
+ scene_get_enclosure_ids(scn, hit->prim.prim_id, enc_ids);
+ if(enc_ids[side] != rwalk->enc_id) {
log_err(scn->dev,
- "%s:%s: the conductive path has reached an invalid interface; "
- "unexpected medium (position: "FORMAT_VECX"; side: %s).\n",
+ "%s:%s: the conductive path has reached an invalid interface. "
+ "Unexpected enclosure -- pos=("FORMAT_VECX"), side=%s\n",
__FILE__, FUNC_NAME, SPLITX(tgt), side == SDIS_FRONT ? "front" : "back");
res = RES_BAD_OP_IRRECOVERABLE;
goto error;
@@ -305,7 +372,7 @@ XD(setup_hit_wos)
* the interface. So we can't yet assume that the random walk has left the
* current medium */
dX(set)(rwalk->vtx.P, tgt);
- rwalk->hit = *hit;
+ rwalk->XD(hit) = *hit;
rwalk->hit_side = side;
exit:
@@ -320,11 +387,10 @@ XD(setup_hit_rt)
const double pos[DIM],
const double dir[DIM],
const struct sXd(hit)* hit,
- struct XD(rwalk)* rwalk)
+ struct rwalk* rwalk)
{
/* Properties */
- struct sdis_interface* interf = NULL;
- struct sdis_medium* mdm = NULL;
+ unsigned enc_ids[2] = {ENCLOSURE_ID_NULL, ENCLOSURE_ID_NULL};
enum sdis_side side = SDIS_SIDE_NULL__;
/* Miscellaneous */
@@ -344,12 +410,11 @@ XD(setup_hit_rt)
side = dX(dot)(N, dir) > 0 ? SDIS_BACK : SDIS_FRONT;
/* Fetch interface properties and check path consistency */
- interf = scene_get_interface(scn, hit->prim.prim_id);
- mdm = side == SDIS_FRONT ? interf->medium_front : interf->medium_back;
- if(mdm != rwalk->mdm) {
+ scene_get_enclosure_ids(scn, hit->prim.prim_id, enc_ids);
+ if(enc_ids[side] != rwalk->enc_id) {
log_err(scn->dev,
- "%s:%s: the conductive path has reached an invalid interface; "
- "unexpected medium (position: "FORMAT_VECX"; side: %s).\n",
+ "%s:%s: the conductive path has reached an invalid interface. "
+ "Unexpected enclosure -- pos=("FORMAT_VECX"), side=%s\n",
__FILE__, FUNC_NAME, SPLITX(tgt), side == SDIS_FRONT ? "front" : "back");
res = RES_BAD_OP_IRRECOVERABLE;
goto error;
@@ -360,7 +425,7 @@ XD(setup_hit_rt)
* the interface. So we can't yet assume that the random walk has left the
* current medium */
dX(set)(rwalk->vtx.P, tgt);
- rwalk->hit = *hit;
+ rwalk->XD(hit) = *hit;
rwalk->hit_side = side;
exit:
@@ -372,7 +437,7 @@ error:
static res_T
XD(sample_next_position)
(struct sdis_scene* scn,
- struct XD(rwalk)* rwalk,
+ struct rwalk* rwalk,
struct ssp_rng* rng,
double* distance) /* Displacement distance */
{
@@ -425,7 +490,7 @@ XD(sample_next_position)
* The next diffusion step would then detect an error. This is why we use a
* new function based on the same geometric operator used in the present
* algorithm. */
- res = XD(check_diffusion_position)(scn, rwalk->mdm, pos);
+ res = XD(check_diffusion_position)(scn, rwalk->enc_id, pos);
/* Diffusion position is valid => move the path to the new position */
if(res == RES_OK) {
@@ -459,8 +524,8 @@ XD(sample_next_position)
* we don't care to save it. */
if(SXD_HIT_NONE(&hit)) {
log_err(scn->dev,
- "%s:%s: unable to find the next diffusion position "
- "(position: "FORMAT_VECX"; direction: "FORMAT_VECX"; distance: %g\n",
+ "%s:%s: unable to find the next diffusion position -- "
+ "position=("FORMAT_VECX"), direction=("FORMAT_VECX"), distance=%g\n",
__FILE__, FUNC_NAME, SPLITX(pos), SPLITX(dir), wos_distance);
res = RES_BAD_OP_IRRECOVERABLE;
goto error;
@@ -478,67 +543,6 @@ error:
goto exit;
}
-static res_T
-XD(handle_volumic_power_wos)
- (struct sdis_scene* scn,
- const struct solid_props* props,
- const double distance, /* [m/fp_to_meter] */
- double* power_term,
- struct XD(temperature)* T)
-{
- double dst = distance * scn->fp_to_meter; /* [m] */
- double term = 0;
- res_T res = RES_OK;
- ASSERT(scn && props && distance >= 0 && power_term && T);
-
- if(props->power == SDIS_VOLUMIC_POWER_NONE) goto exit;
-
- /* No displacement => no power density */
- if(distance == 0) goto exit;
-
- term = dst*dst / (2*DIM*props->lambda);
- T->value += props->power * term;
-
-exit:
- *power_term = term;
- return res;
-}
-
-static res_T
-XD(update_green_path)
- (struct green_path_handle* green_path,
- struct XD(rwalk)* rwalk,
- struct sdis_medium* mdm,
- const struct solid_props* props,
- const double power_term,
- const struct XD(temperature)* T)
-{
- res_T res = RES_OK;
- ASSERT(mdm && props && T);
-
- /* Is the green function estimated? */
- if(!green_path) goto exit;
-
- /* Save power term for green function if any */
- if(props->power != SDIS_VOLUMIC_POWER_NONE) {
- res = green_path_add_power_term(green_path, mdm, &rwalk->vtx, power_term);
- if(res != RES_OK) goto error;
- }
-
- /* Set the green path limit to the current position if the initial condition
- * has been reached */
- if(T->done) {
- res = green_path_set_limit_vertex
- (green_path, mdm, &rwalk->vtx, rwalk->elapsed_time);
- if(res != RES_OK) goto error;
- }
-
-exit:
- return res;
-error:
- goto exit;
-}
-
/*******************************************************************************
* Local function
******************************************************************************/
@@ -546,11 +550,12 @@ res_T
XD(conductive_path_wos)
(struct sdis_scene* scn,
struct rwalk_context* ctx,
- struct XD(rwalk)* rwalk,
+ struct rwalk* rwalk,
struct ssp_rng* rng,
- struct XD(temperature)* T)
+ struct temperature* T)
{
/* Properties */
+ const struct enclosure* enc = NULL;
struct sdis_medium* mdm = NULL;
struct solid_props props_ref = SOLID_PROPS_NULL;
struct solid_props props = SOLID_PROPS_NULL;
@@ -566,18 +571,18 @@ XD(conductive_path_wos)
/* Check pre-conditions */
ASSERT(scn && ctx && rwalk && rng && T);
- ASSERT(sdis_medium_get_type(rwalk->mdm) == SDIS_SOLID);
/* Is green evaluated evaluated */
green = ctx->green_path != NULL;
- /* Keep track of the solid. After conduction, a boundary may have been
- * reached, so the random walk medium is NULL. However, this medium is still
- * needed to update the green path. Hence this backup */
- mdm = rwalk->mdm;
+ res = XD(check_enclosure_consistency)(scn, rwalk);
+ if(res != RES_OK) goto error;
- res = XD(check_medium_consistency)(scn, rwalk);
+ /* Get the enclosure medium */
+ enc = scene_get_enclosure(scn, rwalk->enc_id);
+ res = scene_get_enclosure_medium(scn, enc, &mdm);
if(res != RES_OK) goto error;
+ ASSERT(sdis_medium_get_type(mdm) == SDIS_SOLID);
/* Retrieve the solid properties at the current position. Use them to verify
* that those that are supposed to be constant by the conductive random walk
@@ -588,7 +593,7 @@ XD(conductive_path_wos)
* position. By comparing them to the properties along the random walk, we
* thus verify that the properties are constant throughout the random walk
* with respect to the properties of the reinjected position. */
- solid_get_properties(rwalk->mdm, &rwalk->vtx, &props_ref);
+ solid_get_properties(mdm, &rwalk->vtx, &props_ref);
props = props_ref;
/* The algorithm assumes that lambda, rho and cp are constants. The
@@ -620,7 +625,7 @@ XD(conductive_path_wos)
if(res != RES_OK) goto error;
/* Going back in time */
- res = XD(time_travel)(scn, rwalk, rng, alpha, props.t0, &dst, T);
+ res = XD(time_travel)(scn, rwalk, rng, mdm, alpha, props.t0, &dst, T);
if(res != RES_OK) goto error;
/* Add the volumic power density */
@@ -639,16 +644,16 @@ XD(conductive_path_wos)
}
/* The path reaches a boundary */
- if(!SXD_HIT_NONE(&rwalk->hit)) {
+ if(!SXD_HIT_NONE(&rwalk->XD(hit))) {
T->func = XD(boundary_path);
- rwalk->mdm = NULL;
+ rwalk->enc_id = ENCLOSURE_ID_NULL;
break;
}
#undef REGISTER_VERTEX
/* Retreive and check solid properties at the new position */
- res = solid_get_properties(rwalk->mdm, &rwalk->vtx, &props);
+ res = solid_get_properties(mdm, &rwalk->vtx, &props);
if(res != RES_OK) goto error;
res = check_solid_constant_properties(scn->dev, green, wos, &props_ref, &props);
if(res != RES_OK) goto error;
@@ -657,7 +662,7 @@ XD(conductive_path_wos)
}
/* Save green function data */
- res = XD(update_green_path)
+ res = update_green_path
(ctx->green_path, rwalk, mdm, &props_ref, green_power_term, T);
exit:
diff --git a/src/sdis_heat_path_convective_Xd.h b/src/sdis_heat_path_convective_Xd.h
@@ -66,57 +66,19 @@ error:
* Helper functions
******************************************************************************/
static res_T
-XD(register_heat_vertex_in_fluid)
- (struct sdis_scene* scn,
- struct rwalk_context* ctx,
- struct XD(rwalk)* rwalk,
- const double weight)
-{
- struct sdis_rwalk_vertex vtx = SDIS_RWALK_VERTEX_NULL;
- struct hit_filter_data filter_data = HIT_FILTER_DATA_NULL;
- const float empirical_dst = 0.1f;
- const float range[2] = {0, FLT_MAX};
- float org[DIM];
- float dir[DIM];
- float pos[DIM];
- float dst;
- struct sXd(hit) hit;
-
- if(!ctx->heat_path) return RES_OK;
-
- ASSERT(!SXD_HIT_NONE(&rwalk->hit));
-
- fX_set_dX(org, rwalk->vtx.P);
- fX(set)(dir, rwalk->hit.normal);
- if(rwalk->hit_side == SDIS_BACK) fX(minus)(dir, dir);
-
- 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;
- fX(add)(pos, org, fX(mulf)(dir, dir, dst));
- dX_set_fX(vtx.P, pos);
-
- return register_heat_vertex(ctx->heat_path, &vtx, weight,
- SDIS_HEAT_VERTEX_CONVECTION, (int)ctx->nbranchings);
-}
-
-static res_T
XD(handle_known_fluid_temperature)
- (struct sdis_scene* scn,
- struct rwalk_context* ctx,
- struct XD(rwalk)* rwalk,
- struct XD(temperature)* T)
+ (struct rwalk_context* ctx,
+ struct rwalk* rwalk,
+ struct sdis_medium* mdm,
+ struct temperature* T)
{
double temperature;
int known_temperature;
res_T res = RES_OK;
- ASSERT(scn && ctx && rwalk && T);
- ASSERT(sdis_medium_get_type(rwalk->mdm) == SDIS_FLUID);
+ ASSERT(ctx && rwalk && T);
+ ASSERT(sdis_medium_get_type(mdm) == SDIS_FLUID);
- temperature = fluid_get_temperature(rwalk->mdm, &rwalk->vtx);
+ temperature = fluid_get_temperature(mdm, &rwalk->vtx);
/* Check if the temperature is known */
known_temperature = SDIS_TEMPERATURE_IS_KNOWN(temperature);
@@ -127,12 +89,13 @@ XD(handle_known_fluid_temperature)
if(ctx->green_path) {
res = green_path_set_limit_vertex
- (ctx->green_path, rwalk->mdm, &rwalk->vtx, rwalk->elapsed_time);
+ (ctx->green_path, mdm, &rwalk->vtx, rwalk->elapsed_time);
if(res != RES_OK) goto error;
}
- res = XD(register_heat_vertex_in_fluid)(scn, ctx, rwalk, T->value);
- if(res != RES_OK) goto error;
+ if(ctx->heat_path) {
+ heat_path_get_last_vertex(ctx->heat_path)->weight = T->value;
+ }
exit:
return res;
@@ -143,7 +106,7 @@ error:
static res_T
XD(handle_convective_path_startup)
(struct sdis_scene* scn,
- struct XD(rwalk)* rwalk,
+ struct rwalk* rwalk,
int* path_starts_in_fluid)
{
const float range[2] = {FLT_MIN, FLT_MAX};
@@ -151,9 +114,8 @@ XD(handle_convective_path_startup)
float org[DIM] = {0};
res_T res = RES_OK;
ASSERT(scn && rwalk && path_starts_in_fluid);
- ASSERT(sdis_medium_get_type(rwalk->mdm) == SDIS_FLUID);
- *path_starts_in_fluid = SXD_HIT_NONE(&rwalk->hit);
+ *path_starts_in_fluid = SXD_HIT_NONE(&rwalk->XD(hit));
if(*path_starts_in_fluid == 0) goto exit; /* Nothing to do */
dir[DIM-1] = 1;
@@ -161,8 +123,8 @@ XD(handle_convective_path_startup)
/* Init the path hit field required to define the current enclosure and
* fetch the interface data */
- SXD(scene_view_trace_ray(scn->sXd(view), org, dir, range, NULL, &rwalk->hit));
- if(SXD_HIT_NONE(&rwalk->hit)) {
+ SXD(scene_view_trace_ray(scn->sXd(view), org, dir, range, NULL, &rwalk->XD(hit)));
+ if(SXD_HIT_NONE(&rwalk->XD(hit))) {
log_err(scn->dev,
"%s: the position %g %g %g lies in the surrounding fluid whose "
"temperature must be known.\n", FUNC_NAME, SPLIT3(rwalk->vtx.P));
@@ -170,7 +132,8 @@ XD(handle_convective_path_startup)
goto error;
}
- rwalk->hit_side = fX(dot)(rwalk->hit.normal, dir) < 0 ? SDIS_FRONT : SDIS_BACK;
+ rwalk->hit_side = fX(dot)(rwalk->XD(hit).normal, dir) < 0
+ ? SDIS_FRONT : SDIS_BACK;
exit:
return res;
@@ -179,54 +142,28 @@ error:
}
static res_T
-XD(fetch_fluid_enclosure)
+XD(check_enclosure)
(struct sdis_scene* scn,
- struct XD(rwalk)* rwalk,
- const struct enclosure** out_enclosure)
+ const struct rwalk* rwalk,
+ const struct enclosure* enc)
{
- const struct sdis_interface* interf;
- const struct enclosure* enc;
- unsigned enc_ids[2];
- unsigned enc_id;
res_T res = RES_OK;
- ASSERT(scn && rwalk && out_enclosure);
- ASSERT(sdis_medium_get_type(rwalk->mdm) == SDIS_FLUID);
- ASSERT(!SXD_HIT_NONE(&rwalk->hit));
-
- /* Fetch the current interface and its associated enclosures */
- interf = scene_get_interface(scn, rwalk->hit.prim.prim_id);
- scene_get_enclosure_ids(scn, rwalk->hit.prim.prim_id, enc_ids);
-
- /* Find the enclosure identifier of the current medium */
- ASSERT(interf->medium_front != interf->medium_back);
- if(rwalk->mdm == interf->medium_front) {
- enc_id = enc_ids[0];
- ASSERT(rwalk->hit_side == SDIS_FRONT);
- } else {
- ASSERT(rwalk->mdm == interf->medium_back);
- enc_id = enc_ids[1];
- ASSERT(rwalk->hit_side == SDIS_BACK);
- }
+ ASSERT(scn && rwalk && enc);
- /* Fetch the enclosure data */
- enc = scene_get_enclosure(scn, enc_id);
- ASSERT(enc != NULL);
- if(enc->medium_id == ENCLOSURE_MULTI_MEDIA) {
+ if(enc->medium_id == MEDIUM_ID_MULTI) {
/* The enclosures with multiple media are used to describe limit
* conditions and therefore they cannot be fetched */
log_err(scn->dev,
- "%s: enclosure with multiple media at {%g, %g, %g}. "
- "Path should be reached a limit condition before.\n",
- FUNC_NAME, rwalk->vtx.P[0], rwalk->vtx.P[1], DIM==3 ? rwalk->vtx.P[2]:0);
+ "%s: enclosure with multiple media at ("FORMAT_VECX"). "
+ "The path should be reached a limit condition before.\n",
+ FUNC_NAME, SPLITX(rwalk->vtx.P));
res = RES_BAD_ARG;
goto error;
}
exit:
- *out_enclosure = enc;
return res;
error:
- enc = NULL;
goto exit;
}
@@ -237,14 +174,22 @@ res_T
XD(convective_path)
(struct sdis_scene* scn,
struct rwalk_context* ctx,
- struct XD(rwalk)* rwalk,
+ struct rwalk* rwalk,
struct ssp_rng* rng,
- struct XD(temperature)* T)
+ struct temperature* T)
{
- struct sXd(attrib) attr_P, attr_N;
+ /* Properties */
struct fluid_props props_ref = FLUID_PROPS_NULL;
- const struct sdis_interface* interf;
- const struct enclosure* enc;
+ const struct sdis_interface* interf = NULL;
+ struct sdis_medium* mdm = NULL;
+
+ /* Enclosure */
+ unsigned enc_ids[2] = {ENCLOSURE_ID_NULL, ENCLOSURE_ID_NULL};
+ const struct enclosure* enc = NULL;
+
+ /* Miscellaneous */
+ struct sXd(attrib) attr_P, attr_N;
+ struct sXd(hit)* rwalk_hit = NULL;
double r;
#if SDIS_XD_DIMENSION == 2
float st;
@@ -253,11 +198,19 @@ XD(convective_path)
#endif
int path_starts_in_fluid;
res_T res = RES_OK;
- (void)rng, (void)ctx;
+
ASSERT(scn && ctx && rwalk && rng && T);
- ASSERT(rwalk->mdm->type == SDIS_FLUID);
+ (void)rng, (void)ctx; /* Avoid "unsued variable" warnings */
- res = XD(handle_known_fluid_temperature)(scn, ctx, rwalk, T);
+ rwalk_hit = &rwalk->XD(hit);
+
+ /* Get the enclosure medium */
+ enc = scene_get_enclosure(scn, rwalk->enc_id);
+ if((res = XD(check_enclosure)(scn, rwalk, enc)) != RES_OK) goto error;
+ if((res = scene_get_enclosure_medium(scn, enc, &mdm)) != RES_OK) goto error;
+ ASSERT(sdis_medium_get_type(mdm) == SDIS_FLUID);
+
+ res = XD(handle_known_fluid_temperature)(ctx, rwalk, mdm, T);
if(res != RES_OK) goto error;
if(T->done) goto exit; /* The fluid temperature is known */
@@ -266,13 +219,10 @@ XD(convective_path)
res = XD(handle_convective_path_startup)(scn, rwalk, &path_starts_in_fluid);
if(res != RES_OK) goto error;
- res = XD(fetch_fluid_enclosure)(scn, rwalk, &enc);
- if(res != RES_OK) goto error;
-
/* Retrieve the fluid properties at the current position. Use them to verify
* that those that are supposed to be constant by the convective random walk
* remain the same. */
- res = fluid_get_properties(rwalk->mdm, &rwalk->vtx, &props_ref);
+ res = fluid_get_properties(mdm, &rwalk->vtx, &props_ref);
if(res != RES_OK) goto error;
/* The hc upper bound can be 0 if h is uniformly 0. In that case the result
@@ -280,7 +230,7 @@ XD(convective_path)
if(enc->hc_upper_bound == 0) {
ASSERT(path_starts_in_fluid); /* Cannot be in the fluid without starting there. */
rwalk->vtx.time = props_ref.t0;
- res = XD(handle_known_fluid_temperature)(scn, ctx, rwalk, T);
+ res = XD(handle_known_fluid_temperature)(ctx, rwalk, mdm, T);
if(res != RES_OK) goto error;
if(T->done) {
goto exit; /* Stop the random walk */
@@ -300,7 +250,7 @@ XD(convective_path)
double mu;
/* Fetch fluid properties */
- res = fluid_get_properties(rwalk->mdm, &rwalk->vtx, &props);
+ res = fluid_get_properties(mdm, &rwalk->vtx, &props);
if(res != RES_OK) goto error;
res = check_fluid_constant_properties(scn->dev, &props_ref, &props);
@@ -308,7 +258,7 @@ XD(convective_path)
/* Sample the time using the upper bound. */
mu = enc->hc_upper_bound / (props.rho * props.cp) * enc->S_over_V;
- res = XD(time_rewind)(mu, props.t0, rng, rwalk, ctx, T);
+ res = time_rewind(scn, mu, props.t0, rng, rwalk, ctx, T);
if(res != RES_OK) goto error;
if(T->done) break; /* Limit condition was reached */
@@ -318,44 +268,47 @@ XD(convective_path)
(enc->sXd(view),
ssp_rng_canonical_float(rng),
ssp_rng_canonical_float(rng),
- &prim, &rwalk->hit.u));
- st = rwalk->hit.u;
+ &prim, &rwalk_hit->u));
+ st = rwalk_hit->u;
#else
SXD(scene_view_sample
(enc->sXd(view),
ssp_rng_canonical_float(rng),
ssp_rng_canonical_float(rng),
ssp_rng_canonical_float(rng),
- &prim, rwalk->hit.uv));
- f2_set(st, rwalk->hit.uv);
+ &prim, rwalk_hit->uv));
+ f2_set(st, rwalk_hit->uv);
#endif
/* Map the sampled primitive id from the enclosure space to the scene
* space. Note that the overall scene has only one shape. As a consequence
* neither the geom_id nor the inst_id needs to be updated */
- rwalk->hit.prim.prim_id = enclosure_local2global_prim_id(enc, prim.prim_id);
+ rwalk_hit->prim.prim_id = enclosure_local2global_prim_id(enc, prim.prim_id);
- SXD(primitive_get_attrib(&rwalk->hit.prim, SXD_POSITION, st, &attr_P));
- SXD(primitive_get_attrib(&rwalk->hit.prim, SXD_GEOMETRY_NORMAL, st, &attr_N));
+ SXD(primitive_get_attrib(&rwalk_hit->prim, SXD_POSITION, st, &attr_P));
+ SXD(primitive_get_attrib(&rwalk_hit->prim, SXD_GEOMETRY_NORMAL, st, &attr_N));
dX_set_fX(rwalk->vtx.P, attr_P.value);
- fX(set)(rwalk->hit.normal, attr_N.value);
+ fX(set)(rwalk_hit->normal, attr_N.value);
- /* Fetch the interface of the sampled point. */
- interf = scene_get_interface(scn, rwalk->hit.prim.prim_id);
- if(rwalk->mdm == interf->medium_front) {
- rwalk->hit_side = SDIS_FRONT;
- } else if(rwalk->mdm == interf->medium_back) {
+ /* Define the interface side */
+ scene_get_enclosure_ids(scn, rwalk_hit->prim.prim_id, enc_ids);
+ if(rwalk->enc_id == enc_ids[SDIS_BACK]) {
rwalk->hit_side = SDIS_BACK;
+ } else if(rwalk->enc_id == enc_ids[SDIS_FRONT]) {
+ rwalk->hit_side = SDIS_FRONT;
} else {
- FATAL("Unexpected fluid interface.\n");
+ FATAL("Unexpected fluid interface\n");
}
+ /* Get the interface properties */
+ interf = scene_get_interface(scn, rwalk_hit->prim.prim_id);
+
/* Register the new vertex against the heat path */
res = register_heat_vertex(ctx->heat_path, &rwalk->vtx, T->value,
SDIS_HEAT_VERTEX_CONVECTION, (int)ctx->nbranchings);
if(res != RES_OK) goto error;
/* Setup the fragment of the sampled position into the enclosure. */
- XD(setup_interface_fragment)(&frag, &rwalk->vtx, &rwalk->hit, rwalk->hit_side);
+ XD(setup_interface_fragment)(&frag, &rwalk->vtx, rwalk_hit, rwalk->hit_side);
/* Fetch the convection coefficient of the sampled position */
hc = interface_get_convection_coef(interf, &frag);
@@ -374,9 +327,9 @@ XD(convective_path)
}
}
- rwalk->hit.distance = 0;
+ rwalk_hit->distance = 0;
T->func = XD(boundary_path);
- rwalk->mdm = NULL; /* The random walk is at an interface between 2 media */
+ rwalk->enc_id = ENCLOSURE_ID_NULL; /* Interface between 2 enclosures */
exit:
return res;
diff --git a/src/sdis_heat_path_radiative_Xd.h b/src/sdis_heat_path_radiative_Xd.h
@@ -35,9 +35,9 @@ XD(trace_radiative_path)
(struct sdis_scene* scn,
const float ray_dir[3],
struct rwalk_context* ctx,
- struct XD(rwalk)* rwalk,
+ struct rwalk* rwalk,
struct ssp_rng* rng,
- struct XD(temperature)* T)
+ struct temperature* T)
{
/* The radiative random walk is always performed in 3D. In 2D, the geometry
* are assumed to be extruded to the infinity along the Z dimension. */
@@ -56,9 +56,11 @@ XD(trace_radiative_path)
/* Launch the radiative random walk */
for(;;) {
const struct sdis_interface* interf = NULL;
+ struct sdis_medium* chk_mdm = NULL;
struct hit_filter_data filter_data = HIT_FILTER_DATA_NULL;
struct sdis_interface_fragment frag = SDIS_INTERFACE_FRAGMENT_NULL;
- struct sdis_medium* chk_mdm = NULL;
+ unsigned enc_ids[2] = {ENCLOSURE_ID_NULL, ENCLOSURE_ID_NULL};
+ unsigned chk_enc_id = ENCLOSURE_ID_NULL;
double alpha;
double epsilon;
double r;
@@ -68,16 +70,16 @@ XD(trace_radiative_path)
fX_set_dX(pos, rwalk->vtx.P);
/* Trace the radiative ray */
- filter_data.XD(hit) = rwalk->hit;
+ filter_data.XD(hit) = rwalk->XD(hit);
filter_data.epsilon = 1.e-6;
#if (SDIS_XD_DIMENSION == 2)
SXD(scene_view_trace_ray_3d
- (scn->sXd(view), pos, dir, range, &filter_data, &rwalk->hit));
+ (scn->sXd(view), pos, dir, range, &filter_data, &rwalk->XD(hit)));
#else
SXD(scene_view_trace_ray
- (scn->sXd(view), pos, dir, range, &filter_data, &rwalk->hit));
+ (scn->sXd(view), pos, dir, range, &filter_data, &rwalk->XD(hit)));
#endif
- if(SXD_HIT_NONE(&rwalk->hit)) { /* Fetch the ambient radiative temperature */
+ if(SXD_HIT_NONE(&rwalk->XD(hit))) { /* Fetch the ambient radiative temperature */
struct sdis_radiative_ray ray = SDIS_RADIATIVE_RAY_NULL;
double trad = 0; /* [K] */
@@ -127,11 +129,11 @@ XD(trace_radiative_path)
}
/* Define the hit side */
- rwalk->hit_side = fX(dot)(dir, rwalk->hit.normal) < 0
+ rwalk->hit_side = fX(dot)(dir, rwalk->XD(hit).normal) < 0
? SDIS_FRONT : SDIS_BACK;
/* Move the random walk to the hit position */
- XD(move_pos)(rwalk->vtx.P, dir, rwalk->hit.distance);
+ XD(move_pos)(rwalk->vtx.P, dir, rwalk->XD(hit).distance);
/* Register the random walk vertex against the heat path */
res = register_heat_vertex(ctx->heat_path, &rwalk->vtx, T->value,
@@ -139,8 +141,8 @@ XD(trace_radiative_path)
if(res != RES_OK) goto error;
/* Fetch the new interface and setup the hit fragment */
- interf = scene_get_interface(scn, rwalk->hit.prim.prim_id);
- XD(setup_interface_fragment)(&frag, &rwalk->vtx, &rwalk->hit, rwalk->hit_side);
+ interf = scene_get_interface(scn, rwalk->XD(hit).prim.prim_id);
+ XD(setup_interface_fragment)(&frag, &rwalk->vtx, &rwalk->XD(hit), rwalk->hit_side);
/* Fetch the interface emissivity */
epsilon = interface_side_get_emissivity(interf, SDIS_INTERN_SOURCE_ID, &frag);
@@ -156,35 +158,46 @@ XD(trace_radiative_path)
r = ssp_rng_canonical(rng);
if(r < epsilon) {
T->func = XD(boundary_path);
- rwalk->mdm = NULL; /* The random walk is at an interface between 2 media */
+ rwalk->enc_id = ENCLOSURE_ID_NULL; /* Interface between 2 enclosures */
break;
}
/* Normalize the normal of the interface and ensure that it points toward the
* current medium */
- fX(normalize)(N, rwalk->hit.normal);
- if(rwalk->hit_side == SDIS_BACK){
- chk_mdm = interf->medium_back;
- fX(minus)(N, N);
- } else {
- chk_mdm = interf->medium_front;
+ fX(normalize)(N, rwalk->XD(hit).normal);
+ if(rwalk->hit_side == SDIS_BACK) fX(minus)(N, N);
+
+ /* Check that the radiative path still lies in the same enclosure */
+ scene_get_enclosure_ids(scn, rwalk->XD(hit).prim.prim_id, enc_ids);
+ chk_enc_id = rwalk->hit_side == SDIS_FRONT ? enc_ids[0] : enc_ids[1];
+ if(chk_enc_id != rwalk->enc_id) {
+ log_warn(scn->dev,
+ "%s: the radiative path has escaped from its cavity -- pos=(%g, %g, %g)\n",
+ FUNC_NAME, SPLIT3(rwalk->vtx.P));
+ res = RES_BAD_OP;
+ goto error;
}
- if(chk_mdm != rwalk->mdm) {
- /* To ease the setting of models, the external enclosure is allowed to be
- * incoherent regarding media. Here a radiative path is allowed to join
- * 2 different fluids. */
- const int outside = scene_is_outside
- (scn, rwalk->hit_side, rwalk->hit.prim.prim_id);
- if(outside && chk_mdm->type == SDIS_FLUID) {
- rwalk->mdm = chk_mdm;
- } else {
- log_warn(scn->dev, "%s: inconsistent medium definition at `%g %g %g'.\n",
- FUNC_NAME, SPLIT3(rwalk->vtx.P));
- res = RES_BAD_OP;
- goto error;
- }
+ /* Verify that the intersection, although in the same enclosure, touches the
+ * interface of a fluid. We verify this by interface, since a radiative path
+ * can be traced in an enclosure containing several media used to describe a
+ * set of boundary conditions.
+ *
+ * If the enclosure is good but the media type is not, this means that the
+ * radiative path is sampled in the wrong media. This is not a numerical
+ * problem, but a user problem: trying to sample a radiative path in a solid
+ * when semi-transparent solids are not yet supported by Stardis. This error
+ * is therefore fatal for the calculation */
+ chk_mdm = rwalk->hit_side == SDIS_FRONT
+ ? interf->medium_front : interf->medium_back;
+ if(sdis_medium_get_type(chk_mdm)) {
+ log_err(scn->dev,
+ "%s: a radiative path cannot evolve in a solid -- pos=(%g, %g, %g)\n",
+ FUNC_NAME, SPLIT3(rwalk->vtx.P));
+ res = RES_BAD_OP_IRRECOVERABLE;
+ goto error;
}
+
alpha = interface_side_get_specular_fraction(interf, SDIS_INTERN_SOURCE_ID, &frag);
r = ssp_rng_canonical(rng);
if(r < alpha) { /* Sample specular part */
@@ -204,9 +217,9 @@ res_T
XD(radiative_path)
(struct sdis_scene* scn,
struct rwalk_context* ctx,
- struct XD(rwalk)* rwalk,
+ struct rwalk* rwalk,
struct ssp_rng* rng,
- struct XD(temperature)* T)
+ struct temperature* T)
{
/* The radiative random walk is always performed in 3D. In 2D, the geometry
* are assumed to be extruded to the infinity along the Z dimension. */
@@ -215,11 +228,11 @@ XD(radiative_path)
res_T res = RES_OK;
ASSERT(scn && ctx && rwalk && rng && T);
- ASSERT(!SXD_HIT_NONE(&rwalk->hit));
+ ASSERT(!SXD_HIT_NONE(&rwalk->XD(hit)));
/* Normalize the normal of the interface and ensure that it points toward the
* current medium */
- fX(normalize(N, rwalk->hit.normal));
+ fX(normalize(N, rwalk->XD(hit).normal));
if(rwalk->hit_side == SDIS_BACK) {
fX(minus(N, N));
}
diff --git a/src/sdis_misc.c b/src/sdis_misc.c
@@ -14,12 +14,85 @@
* along with this program. If not, see <http://www.gnu.org/licenses/>. */
#include "sdis.h"
+#include "sdis_heat_path.h"
+#include "sdis_log.h"
+#include "sdis_medium_c.h"
+#include "sdis_misc.h"
+#include "sdis_green.h"
-/* Generate the generic functions */
-#define SDIS_XD_DIMENSION 2
-#include "sdis_misc_Xd.h"
-#define SDIS_XD_DIMENSION 3
-#include "sdis_misc_Xd.h"
+#include <star/ssp.h>
+
+res_T
+time_rewind
+ (struct sdis_scene* scn,
+ const double mu,
+ const double t0,
+ struct ssp_rng* rng,
+ struct rwalk* rwalk,
+ const struct rwalk_context* ctx,
+ struct temperature* T)
+{
+ const struct enclosure* enc = NULL;
+ struct sdis_medium* mdm = NULL;
+ double temperature;
+ double tau;
+ res_T res = RES_OK;
+ ASSERT(scn && rwalk && ctx && rng && T);
+
+ /* Get the current medium */
+ enc = scene_get_enclosure(scn, rwalk->enc_id);
+ res = scene_get_enclosure_medium(scn, enc, &mdm);
+ if(res != RES_OK) goto error;
+
+ /* Sample the time using the upper bound. */
+ tau = ssp_ran_exp(rng, mu);
+
+ /* Increment the elapsed time */
+ ASSERT(rwalk->vtx.time >= t0);
+ rwalk->elapsed_time += MMIN(tau, rwalk->vtx.time - t0);
+
+ if(IS_INF(rwalk->vtx.time)) goto exit; /* Steady computation */
+
+ /* Time rewind */
+ rwalk->vtx.time = MMAX(rwalk->vtx.time - tau, t0); /* Time rewind */
+
+ /* The path does not reach the limit condition */
+ if(rwalk->vtx.time > t0) goto exit;
+
+ /* Fetch the initial temperature */
+ temperature = medium_get_temperature(mdm, &rwalk->vtx);
+ if(SDIS_TEMPERATURE_IS_UNKNOWN(temperature)) {
+ log_err(mdm->dev, "the path reaches the limit condition but the "
+ "%s temperature remains unknown -- position=%g, %g, %g\n",
+ medium_type_to_string(sdis_medium_get_type(mdm)),
+ SPLIT3(rwalk->vtx.P));
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ /* Update temperature */
+ T->value += temperature;
+ T->done = 1;
+
+ if(ctx->heat_path) {
+ /* Update the registered vertex data */
+ struct sdis_heat_vertex* vtx;
+ vtx = heat_path_get_last_vertex(ctx->heat_path);
+ vtx->time = rwalk->vtx.time;
+ vtx->weight = T->value;
+ }
+
+ if(ctx->green_path) {
+ res = green_path_set_limit_vertex
+ (ctx->green_path, mdm, &rwalk->vtx, rwalk->elapsed_time);
+ if(res != RES_OK) goto error;
+ }
+
+exit:
+ return res;
+error:
+ goto exit;
+}
res_T
check_primitive_uv_2d(struct sdis_device* dev, const double param_coord[])
diff --git a/src/sdis_misc.h b/src/sdis_misc.h
@@ -148,22 +148,14 @@ register_heat_vertex
}
extern LOCAL_SYM res_T
-time_rewind_2d
- (const double mu,
+time_rewind
+ (struct sdis_scene* scn,
+ const double mu,
const double t0, /* Initial time */
struct ssp_rng* rng,
- struct rwalk_2d* rwalk,
+ struct rwalk* rwalk,
const struct rwalk_context* ctx,
- struct temperature_2d* T);
-
-extern LOCAL_SYM res_T
-time_rewind_3d
- (const double mu,
- const double t0, /* Initial time */
- struct ssp_rng* rng,
- struct rwalk_3d* rwalk,
- const struct rwalk_context* ctx,
- struct temperature_3d* T);
+ struct temperature* T);
/* Check the validity of the parametric coordinate onto a 2D primitive. If it
* is invalid, the function prints an error message and return RES_BAD_ARG. */
diff --git a/src/sdis_misc_Xd.h b/src/sdis_misc_Xd.h
@@ -1,90 +0,0 @@
-/* Copyright (C) 2016-2024 |Méso|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_heat_path.h"
-#include "sdis_log.h"
-#include "sdis_medium_c.h"
-#include "sdis_misc.h"
-#include "sdis_green.h"
-
-#include <star/ssp.h>
-
-#include "sdis_Xd_begin.h"
-
-res_T
-XD(time_rewind)
- (const double mu,
- const double t0,
- struct ssp_rng* rng,
- struct XD(rwalk)* rwalk,
- const struct rwalk_context* ctx,
- struct XD(temperature)* T)
-{
- double temperature;
- double tau;
- res_T res = RES_OK;
- ASSERT(rwalk && rng && T);
-
- /* Sample the time using the upper bound. */
- tau = ssp_ran_exp(rng, mu);
-
- /* Increment the elapsed time */
- ASSERT(rwalk->vtx.time >= t0);
- rwalk->elapsed_time += MMIN(tau, rwalk->vtx.time - t0);
-
- if(IS_INF(rwalk->vtx.time)) goto exit; /* Steady computation */
-
- /* Time rewind */
- rwalk->vtx.time = MMAX(rwalk->vtx.time - tau, t0); /* Time rewind */
-
- /* The path does not reach the limit condition */
- if(rwalk->vtx.time > t0) goto exit;
-
- /* Fetch the initial temperature */
- temperature = medium_get_temperature(rwalk->mdm, &rwalk->vtx);
- if(SDIS_TEMPERATURE_IS_UNKNOWN(temperature)) {
- log_err(rwalk->mdm->dev, "the path reaches the limit condition but the "
- "%s temperature remains unknown -- position=%g, %g, %g\n",
- medium_type_to_string(sdis_medium_get_type(rwalk->mdm)),
- SPLIT3(rwalk->vtx.P));
- res = RES_BAD_ARG;
- goto error;
- }
-
- /* Update temperature */
- T->value += temperature;
- T->done = 1;
-
- if(ctx->heat_path) {
- /* Update the registered vertex data */
- struct sdis_heat_vertex* vtx;
- vtx = heat_path_get_last_vertex(ctx->heat_path);
- vtx->time = rwalk->vtx.time;
- vtx->weight = T->value;
- }
-
- if(ctx->green_path) {
- res = green_path_set_limit_vertex(ctx->green_path, rwalk->mdm,
- &rwalk->vtx, rwalk->elapsed_time);
- if(res != RES_OK) goto error;
- }
-
-exit:
- return res;
-error:
- goto exit;
-}
-
-#include "sdis_Xd_end.h"
diff --git a/src/sdis_realisation.c b/src/sdis_realisation.c
@@ -19,16 +19,30 @@
/*******************************************************************************
* Helper functions
******************************************************************************/
-static INLINE int
+static INLINE res_T
check_ray_realisation_args(const struct ray_realisation_args* args)
{
- return args
- && args->rng
- && args->medium
- && args->medium->type == SDIS_FLUID
- && args->time >= 0
- && args->picard_order > 0
- && (unsigned)args->diff_algo < SDIS_DIFFUSION_ALGORITHMS_COUNT__;
+ /* Check pointers */
+ if(!args || !args->rng) return RES_BAD_ARG;
+
+ if(args->time < 0) return RES_BAD_ARG;
+ if(args->picard_order <= 0) return RES_BAD_ARG;
+
+ if((unsigned)args->diff_algo >= SDIS_DIFFUSION_ALGORITHMS_COUNT__) {
+ return RES_BAD_ARG;
+ }
+
+ /* Check the enclosure identifier. Only its validity is checked, not the fact
+ * that the enclosure is a fluid. Even though Stardis doesn't allow you to
+ * sample a radiative path in a solid, we don't query the medium of the
+ * enclosure since it may contain several: querying the medium will therefore
+ * return an error. The type of medium is checked later, when sampling the
+ * radiative path, when it reaches an interface whose medium must be a fluid*/
+ if(args->enc_id == ENCLOSURE_ID_NULL) {
+ return RES_BAD_ARG;
+ }
+
+ return RES_OK;
}
/*******************************************************************************
@@ -47,17 +61,17 @@ ray_realisation_3d
double* weight)
{
struct rwalk_context ctx = RWALK_CONTEXT_NULL;
- struct rwalk_3d rwalk = RWALK_NULL_3d;
- struct temperature_3d T = TEMPERATURE_NULL_3d;
+ struct rwalk rwalk = RWALK_NULL;
+ struct temperature T = TEMPERATURE_NULL;
float dir[3];
res_T res = RES_OK;
- ASSERT(scn && weight && check_ray_realisation_args(args));
+ ASSERT(scn && weight && check_ray_realisation_args(args) == RES_OK);
d3_set(rwalk.vtx.P, args->position);
rwalk.vtx.time = args->time;
- rwalk.hit = S3D_HIT_NULL;
+ rwalk.hit_3d = S3D_HIT_NULL;
rwalk.hit_side = SDIS_SIDE_NULL__;
- rwalk.mdm = args->medium;
+ rwalk.enc_id = args->enc_id;
ctx.heat_path = args->heat_path;
ctx.Tmin = scn->tmin;
@@ -69,7 +83,7 @@ ray_realisation_3d
ctx.max_branchings = args->picard_order - 1;
ctx.irealisation = args->irealisation;
ctx.diff_algo = args->diff_algo;
-
+
f3_set_d3(dir, args->direction);
/* Register the starting position against the heat path */
diff --git a/src/sdis_realisation.h b/src/sdis_realisation.h
@@ -22,11 +22,13 @@
#include <rsys/rsys.h>
/* Forward declarations */
+struct bound_flux_result;
struct green_path_handle;
+struct rwalk;
struct sdis_heat_path;
struct sdis_scene;
struct ssp_rng;
-struct bound_flux_result;
+struct temperature;
enum flux_flag {
FLUX_FLAG_CONVECTIVE = BIT(FLUX_CONVECTIVE),
@@ -41,24 +43,24 @@ extern LOCAL_SYM res_T
sample_coupled_path_2d
(struct sdis_scene* scn,
struct rwalk_context* ctx,
- struct rwalk_2d* rwalk,
+ struct rwalk* rwalk,
struct ssp_rng* rng,
- struct temperature_2d* T);
+ struct temperature* T);
extern LOCAL_SYM res_T
sample_coupled_path_3d
(struct sdis_scene* scn,
struct rwalk_context* ctx,
- struct rwalk_3d* rwalk,
+ struct rwalk* rwalk,
struct ssp_rng* rng,
- struct temperature_3d* T);
+ struct temperature* T);
/*******************************************************************************
* Realisation at a given position and time IN a medium
******************************************************************************/
struct probe_realisation_args {
struct ssp_rng* rng;
- struct sdis_medium* medium; /* Medium into which the realisation starts */
+ unsigned enc_id; /* Enclosure into which the realisation starts */
double position[3]; /* Probe position */
double time; /* Observation time */
size_t picard_order; /* Picard order to estimate radiative temperature */
@@ -69,7 +71,7 @@ struct probe_realisation_args {
};
#define PROBE_REALISATION_ARGS_NULL__ { \
NULL, /* RNG */ \
- NULL, /* Medium */ \
+ ENCLOSURE_ID_NULL, /* Enclosure */ \
{0,0,0}, /* Position */ \
-1, /* Observation time */ \
0, /* Picard order */ \
@@ -180,7 +182,7 @@ boundary_flux_realisation_3d
******************************************************************************/
struct ray_realisation_args {
struct ssp_rng* rng;
- struct sdis_medium* medium; /* Medium into which the realisation starts */
+ unsigned enc_id; /* Enclosure into which the realisation starts */
double position[3]; /* Ray position */
double direction[3]; /* Ray direction */
double time; /* Observation time */
@@ -191,7 +193,7 @@ struct ray_realisation_args {
};
#define RAY_REALISATION_ARGS_NULL__ { \
NULL, /* RNG */ \
- NULL, /* Medium */ \
+ ENCLOSURE_ID_NULL, /* Enclosure */ \
{0,0,0}, /* Position */ \
{0,0,0}, /* Direction */ \
-1, /* Observation time */ \
diff --git a/src/sdis_realisation_Xd.h b/src/sdis_realisation_Xd.h
@@ -32,18 +32,19 @@
#ifndef SDIS_REALISATION_XD_H
#define SDIS_REALISATION_XD_H
-static INLINE int
+static INLINE res_T
check_probe_realisation_args(const struct probe_realisation_args* args)
{
return args
&& args->rng
- && args->medium
+ && args->enc_id != ENCLOSURE_ID_NULL
&& args->time >= 0
&& args->picard_order > 0
- && (unsigned)args->diff_algo < SDIS_DIFFUSION_ALGORITHMS_COUNT__;
+ && (unsigned)args->diff_algo < SDIS_DIFFUSION_ALGORITHMS_COUNT__
+ ? RES_OK : RES_BAD_ARG;
}
-static INLINE int
+static INLINE res_T
check_boundary_realisation_args(const struct boundary_realisation_args* args)
{
return args
@@ -55,10 +56,11 @@ check_boundary_realisation_args(const struct boundary_realisation_args* args)
&& args->time >= 0
&& args->picard_order > 0
&& (args->side == SDIS_FRONT || args->side == SDIS_BACK)
- && (unsigned)args->diff_algo < SDIS_DIFFUSION_ALGORITHMS_COUNT__;
+ && (unsigned)args->diff_algo < SDIS_DIFFUSION_ALGORITHMS_COUNT__
+ ? RES_OK : RES_BAD_ARG;
}
-static INLINE int
+static INLINE res_T
check_boundary_flux_realisation_args
(const struct boundary_flux_realisation_args* args)
{
@@ -71,7 +73,8 @@ check_boundary_flux_realisation_args
&& args->time >= 0
&& args->picard_order > 0
&& (args->solid_side == SDIS_FRONT || args->solid_side == SDIS_BACK)
- && (unsigned)args->diff_algo < SDIS_DIFFUSION_ALGORITHMS_COUNT__;
+ && (unsigned)args->diff_algo < SDIS_DIFFUSION_ALGORITHMS_COUNT__
+ ? RES_OK : RES_BAD_ARG;
}
#endif /* SDIS_REALISATION_XD_H */
@@ -82,15 +85,15 @@ res_T
XD(sample_coupled_path)
(struct sdis_scene* scn,
struct rwalk_context* ctx,
- struct XD(rwalk)* rwalk,
+ struct rwalk* rwalk,
struct ssp_rng* rng,
- struct XD(temperature)* T)
+ struct temperature* T)
{
#ifndef NDEBUG
/* Stack that saves the state of each recursion steps. */
struct entry {
- struct XD(temperature) temperature;
- struct XD(rwalk) rwalk;
+ struct temperature temperature;
+ struct rwalk rwalk;
}* stack = NULL;
size_t istack = 0;
#endif
@@ -109,8 +112,8 @@ XD(sample_coupled_path)
while(!T->done) {
/* Save the current random walk state */
- const struct XD(rwalk) rwalk_bkp = *rwalk;
- const struct XD(temperature) T_bkp = *T;
+ const struct rwalk rwalk_bkp = *rwalk;
+ const struct temperature T_bkp = *T;
size_t nfails = 0; /* #failures */
#ifndef NDEBUG
@@ -171,27 +174,39 @@ XD(probe_realisation)
struct probe_realisation_args* args,
double* weight)
{
+ /* Starting enclosure/medium */
+ const struct enclosure* enc = NULL;
+ struct sdis_medium* mdm = NULL;
+
+ /* Random walk */
struct rwalk_context ctx = RWALK_CONTEXT_NULL;
- struct XD(rwalk) rwalk = XD(RWALK_NULL);
- struct XD(temperature) T = XD(TEMPERATURE_NULL);
+ struct rwalk rwalk = RWALK_NULL;
+ struct temperature T = TEMPERATURE_NULL;
+
+ /* Miscellaneous */
enum sdis_heat_vertex_type type;
double t0;
double (*get_initial_temperature)
(const struct sdis_medium* mdm,
const struct sdis_rwalk_vertex* vtx);
res_T res = RES_OK;
- ASSERT(scn && weight && check_probe_realisation_args(args));
+ ASSERT(scn && weight && check_probe_realisation_args(args) == RES_OK);
- switch(args->medium->type) {
+ /* Get the enclosure medium */
+ enc = scene_get_enclosure(scn, args->enc_id);
+ res = scene_get_enclosure_medium(scn, enc, &mdm);
+ if(res != RES_OK) goto error;
+
+ switch(sdis_medium_get_type(mdm)) {
case SDIS_FLUID:
T.func = XD(convective_path);
get_initial_temperature = fluid_get_temperature;
- t0 = fluid_get_t0(args->medium);
+ t0 = fluid_get_t0(mdm);
break;
case SDIS_SOLID:
T.func = XD(conductive_path);
get_initial_temperature = solid_get_temperature;
- t0 = solid_get_t0(args->medium);
+ t0 = solid_get_t0(mdm);
break;
default: FATAL("Unreachable code\n"); break;
}
@@ -200,7 +215,7 @@ XD(probe_realisation)
rwalk.vtx.time = args->time;
/* Register the starting position against the heat path */
- type = args->medium->type == SDIS_SOLID
+ type = sdis_medium_get_type(mdm) == SDIS_SOLID
? SDIS_HEAT_VERTEX_CONDUCTION
: SDIS_HEAT_VERTEX_CONVECTION;
res = register_heat_vertex(args->heat_path, &rwalk.vtx, 0, type, 0);
@@ -210,7 +225,7 @@ XD(probe_realisation)
double tmp;
/* Check the initial condition. */
rwalk.vtx.time = t0;
- tmp = get_initial_temperature(args->medium, &rwalk.vtx);
+ tmp = get_initial_temperature(mdm, &rwalk.vtx);
if(SDIS_TEMPERATURE_IS_KNOWN(tmp)) {
*weight = tmp;
goto exit;
@@ -224,8 +239,8 @@ XD(probe_realisation)
goto error;
}
- rwalk.hit = SXD_HIT_NULL;
- rwalk.mdm = args->medium;
+ rwalk.XD(hit) = SXD_HIT_NULL;
+ rwalk.enc_id = args->enc_id;
ctx.green_path = args->green_path;
ctx.heat_path = args->heat_path;
@@ -258,8 +273,8 @@ XD(boundary_realisation)
double* weight)
{
struct rwalk_context ctx = RWALK_CONTEXT_NULL;
- struct XD(rwalk) rwalk = XD(RWALK_NULL);
- struct XD(temperature) T = XD(TEMPERATURE_NULL);
+ struct rwalk rwalk = RWALK_NULL;
+ struct temperature T = TEMPERATURE_NULL;
struct sXd(attrib) attr;
#if SDIS_XD_DIMENSION == 2
float st;
@@ -267,13 +282,13 @@ XD(boundary_realisation)
float st[2];
#endif
res_T res = RES_OK;
- ASSERT(scn && weight && check_boundary_realisation_args(args));
+ ASSERT(scn && weight && check_boundary_realisation_args(args) == RES_OK);
T.func = XD(boundary_path);
rwalk.hit_side = args->side;
- rwalk.hit.distance = 0;
+ rwalk.XD(hit).distance = 0;
rwalk.vtx.time = args->time;
- rwalk.mdm = NULL; /* The random walk is at an interface between 2 media */
+ rwalk.enc_id = ENCLOSURE_ID_NULL; /* At an interface between 2 enclosures */
#if SDIS_XD_DIMENSION == 2
st = (float)args->uv[0];
@@ -283,20 +298,20 @@ XD(boundary_realisation)
/* Fetch the primitive */
SXD(scene_view_get_primitive
- (scn->sXd(view), (unsigned int)args->iprim, &rwalk.hit.prim));
+ (scn->sXd(view), (unsigned int)args->iprim, &rwalk.XD(hit).prim));
/* Retrieve the world space position of the probe onto the primitive */
- SXD(primitive_get_attrib(&rwalk.hit.prim, SXD_POSITION, st, &attr));
+ SXD(primitive_get_attrib(&rwalk.XD(hit).prim, SXD_POSITION, st, &attr));
dX_set_fX(rwalk.vtx.P, attr.value);
/* Retrieve the primitive normal */
- SXD(primitive_get_attrib(&rwalk.hit.prim, SXD_GEOMETRY_NORMAL, st, &attr));
- fX(set)(rwalk.hit.normal, attr.value);
+ SXD(primitive_get_attrib(&rwalk.XD(hit).prim, SXD_GEOMETRY_NORMAL, st, &attr));
+ fX(set)(rwalk.XD(hit).normal, attr.value);
#if SDIS_XD_DIMENSION==2
- rwalk.hit.u = st;
+ rwalk.XD(hit).u = st;
#else
- f2_set(rwalk.hit.uv, st);
+ f2_set(rwalk.XD(hit).uv, st);
#endif
res = register_heat_vertex(args->heat_path, &rwalk.vtx, 0/*weight*/,
@@ -332,14 +347,14 @@ XD(boundary_flux_realisation)
struct boundary_flux_realisation_args* args,
struct bound_flux_result* result)
{
+ /* Random walk */
struct rwalk_context ctx = RWALK_CONTEXT_NULL;
- struct XD(rwalk) rwalk;
- struct XD(temperature) T;
+ struct rwalk rwalk = RWALK_NULL;
+ struct temperature T = TEMPERATURE_NULL;
+
+ /* Boundary */
struct sXd(attrib) attr;
struct sXd(primitive) prim;
- struct sdis_interface* interf = NULL;
- struct sdis_medium* fluid_mdm = NULL;
-
#if SDIS_XD_DIMENSION == 2
float st;
#else
@@ -347,13 +362,17 @@ XD(boundary_flux_realisation)
#endif
double P[SDIS_XD_DIMENSION];
float N[SDIS_XD_DIMENSION];
+ unsigned enc_ids[2] = {ENCLOSURE_ID_NULL, ENCLOSURE_ID_NULL};
+ enum sdis_side fluid_side;
+
+ /* Miscellaneous */
double Tmin, Tmin2, Tmin3;
double That, That2, That3;
- enum sdis_side fluid_side;
res_T res = RES_OK;
char compute_radiative;
char compute_convective;
- ASSERT(scn && result && check_boundary_flux_realisation_args(args));
+
+ ASSERT(scn && result && check_boundary_flux_realisation_args(args) == RES_OK);
#if SDIS_XD_DIMENSION == 2
#define SET_PARAM(Dest, Src) (Dest).u = (Src);
@@ -386,14 +405,14 @@ XD(boundary_flux_realisation)
SXD(primitive_get_attrib(&prim, SXD_GEOMETRY_NORMAL, st, &attr));
fX(set)(N, attr.value);
- #define RESET_WALK(Side, Mdm) { \
- rwalk = XD(RWALK_NULL); \
+ #define RESET_WALK(Side, EncId) { \
+ rwalk = RWALK_NULL; \
rwalk.hit_side = (Side); \
- rwalk.hit.distance = 0; \
+ rwalk.XD(hit).distance = 0; \
rwalk.vtx.time = args->time; \
- rwalk.mdm = (Mdm); \
- rwalk.hit.prim = prim; \
- SET_PARAM(rwalk.hit, st); \
+ rwalk.enc_id = (EncId); \
+ rwalk.XD(hit).prim = prim; \
+ SET_PARAM(rwalk.XD(hit), st); \
ctx.Tmin = Tmin; \
ctx.Tmin3 = Tmin3; \
ctx.That = That; \
@@ -403,24 +422,23 @@ XD(boundary_flux_realisation)
ctx.irealisation = args->irealisation; \
ctx.diff_algo = args->diff_algo; \
dX(set)(rwalk.vtx.P, P); \
- fX(set)(rwalk.hit.normal, N); \
- T = XD(TEMPERATURE_NULL); \
+ fX(set)(rwalk.XD(hit).normal, N); \
+ T = TEMPERATURE_NULL; \
} (void)0
/* Compute boundary temperature */
- RESET_WALK(args->solid_side, NULL);
+ RESET_WALK(args->solid_side, ENCLOSURE_ID_NULL);
T.func = XD(boundary_path);
res = XD(sample_coupled_path)(scn, &ctx, &rwalk, args->rng, &T);
if(res != RES_OK) return res;
result->Tboundary = T.value;
- /* Fetch the fluid medium */
- interf = scene_get_interface(scn, (unsigned)args->iprim);
- fluid_mdm = interface_get_medium(interf, fluid_side);
+ /* Get the enclosures */
+ scene_get_enclosure_ids(scn, (unsigned)args->iprim, enc_ids);
/* Compute radiative temperature */
if(compute_radiative) {
- RESET_WALK(fluid_side, fluid_mdm);
+ RESET_WALK(fluid_side, enc_ids[fluid_side]);
T.func = XD(radiative_path);
res = XD(sample_coupled_path)(scn, &ctx, &rwalk, args->rng, &T);
if(res != RES_OK) return res;
@@ -430,7 +448,7 @@ XD(boundary_flux_realisation)
/* Compute fluid temperature */
if(compute_convective) {
- RESET_WALK(fluid_side, fluid_mdm);
+ RESET_WALK(fluid_side, enc_ids[fluid_side]);
T.func = XD(convective_path);
res = XD(sample_coupled_path)(scn, &ctx, &rwalk, args->rng, &T);
if(res != RES_OK) return res;
diff --git a/src/sdis_scene.c b/src/sdis_scene.c
@@ -427,25 +427,57 @@ scene_get_interface(const struct sdis_scene* scn, const unsigned iprim)
}
res_T
-scene_get_medium
+scene_get_enclosure_id
(struct sdis_scene* scn,
const double pos[],
- struct sdis_medium** out_medium)
+ unsigned* enc_id)
{
return scene_is_2d(scn)
- ? scene_get_medium_2d(scn, pos, out_medium)
- : scene_get_medium_3d(scn, pos, out_medium);
+ ? scene_get_enclosure_id_2d(scn, pos, enc_id)
+ : scene_get_enclosure_id_3d(scn, pos, enc_id);
}
res_T
-scene_get_medium_in_closed_boundaries
+scene_get_enclosure_id_in_closed_boundaries
(struct sdis_scene* scn,
const double pos[],
- struct sdis_medium** out_medium)
+ unsigned* enc_id)
{
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);
+ ? scene_get_enclosure_id_in_closed_boundaries_2d(scn, pos, enc_id)
+ : scene_get_enclosure_id_in_closed_boundaries_3d(scn, pos, enc_id);
+}
+
+res_T
+scene_get_enclosure_medium
+ (struct sdis_scene* scn,
+ const struct enclosure* enc,
+ struct sdis_medium** out_mdm)
+{
+ struct sdis_medium* mdm = NULL;
+ res_T res = RES_OK;
+
+ ASSERT(scn && enc && out_mdm);
+
+ /* Check that the enclosure doesn't surround multiple media */
+ if(enc->medium_id == MEDIUM_ID_MULTI) {
+ log_warn(scn->dev,
+ "%s: invalid medium request. The enclosure includes several media.\n",
+ FUNC_NAME);
+ res = RES_BAD_OP;
+ goto error;
+ }
+
+ /* Obtain enclosure medium */
+ ASSERT(enc->medium_id < darray_medium_size_get(&scn->media));
+ mdm = darray_medium_data_get(&scn->media)[enc->medium_id];
+
+error:
+ *out_mdm = mdm;
+ goto exit;
+exit:
+ mdm = NULL;
+ return res;
}
res_T
diff --git a/src/sdis_scene_Xd.h b/src/sdis_scene_Xd.h
@@ -823,7 +823,7 @@ XD(register_enclosure)(struct sdis_scene* scn, struct sencXd(enclosure)* enc)
/* Setup the medium id of the enclosure */
if(header.enclosed_media_count > 1) {
- enc_data->medium_id = ENCLOSURE_MULTI_MEDIA;
+ enc_data->medium_id = MEDIUM_ID_MULTI;
} else {
SENCXD(enclosure_get_medium(enc, 0, &enc_data->medium_id));
}
@@ -1094,13 +1094,12 @@ error:
/*******************************************************************************
* Local functions
******************************************************************************/
-static INLINE res_T
-XD(scene_get_medium)
+static res_T
+XD(scene_get_enclosure_id)
(struct sdis_scene* scn,
const double pos[DIM],
- struct sdis_medium** out_medium)
+ unsigned* out_enc_id)
{
- struct sdis_medium* medium = NULL;
size_t iprim, nprims;
float P[DIM];
/* Range of the parametric coordinate into which positions are challenged */
@@ -1110,6 +1109,7 @@ XD(scene_get_medium)
float st[3][2];
#endif
size_t nsteps = 3;
+ unsigned enc_id = ENCLOSURE_ID_NULL;
res_T res = RES_OK;
ASSERT(scn && pos);
@@ -1162,8 +1162,8 @@ XD(scene_get_medium)
} while((SXD_HIT_NONE(&hit) || HIT_ON_BOUNDARY(&hit, P, dir))
&& ++istep < nsteps);
- /* No valid intersection is found on the current primitive. Challenge
- * another. */
+ /* No valid intersection is found on the current primitive.
+ * Challenge another. */
if(istep >= nsteps) continue;
fX(normalize)(N, hit.normal);
@@ -1171,18 +1171,17 @@ XD(scene_get_medium)
/* Not too close and not roughly orthognonal */
if(hit.distance > 1.e-6 && absf(cos_N_dir) > 1.e-2f) {
- const struct enclosure* enclosure = NULL;
unsigned enc_ids[2];
- const struct sdis_interface* interf;
-
- interf = scene_get_interface(scn, hit.prim.prim_id);
scene_get_enclosure_ids(scn, hit.prim.prim_id, enc_ids);
- break;
+ enc_id = cos_N_dir < 0 ? enc_ids[0] : enc_ids[1];
+
+ break; /* That's all folks */
}
}
if(iprim >= nprims) {
- log_warn(scn->dev, "%s: could not retrieve the medium at {%g, %g, %g}.\n",
+ log_warn(scn->dev,
+ "%s: cannot retrieve current enclosure at {%g, %g, %g}.\n",
FUNC_NAME, P[0], P[1], DIM == 3 ? P[2] : 0);
res = RES_BAD_OP;
goto error;
@@ -1190,31 +1189,32 @@ XD(scene_get_medium)
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",
+ "%s: performance issue. Up to %lu primitives were tested to find "
+ "current enclosure at {%g, %g, %g}.\n",
FUNC_NAME, (unsigned long)iprim, P[0], P[1], DIM == 3 ? P[2] : 0);
}
exit:
- *out_medium = medium;
+ *out_enc_id = enc_id;
return res;
error:
+ enc_id = ENCLOSURE_ID_NULL;
goto exit;
}
-static INLINE res_T
-XD(scene_get_medium_in_closed_boundaries)
+static res_T
+XD(scene_get_enclosure_id_in_closed_boundaries)
(struct sdis_scene* scn,
const double pos[DIM],
- struct sdis_medium** out_medium)
+ unsigned* out_enc_id)
{
- 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}};
+ float frame[DIM*DIM];
+ float P[DIM];
+ unsigned enc_id = ENCLOSURE_ID_NULL;
int idir;
res_T res = RES_OK;
- ASSERT(scn && pos);
+ ASSERT(scn && pos && out_enc_id);
/* 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
@@ -1222,8 +1222,6 @@ XD(scene_get_medium_in_closed_boundaries)
#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
@@ -1252,22 +1250,23 @@ XD(scene_get_medium_in_closed_boundaries)
/* 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;
+ unsigned enc_ids[2];
+ scene_get_enclosure_ids(scn, hit.prim.prim_id, enc_ids);
+ enc_id = cos_N_dir < 0 ? enc_ids[0] : enc_ids[1];
+
+ break; /* That's all folks */
}
}
- if(idir >= 2*DIM) {
- res = XD(scene_get_medium)(scn, pos, &medium);
+ if(idir >= 2*DIM) { /* Fallback to scene_get_enclosure_id function */
+ res = XD(scene_get_enclosure_id)(scn, pos, &enc_id);
if(res != RES_OK) goto error;
}
exit:
- *out_medium = medium;
+ *out_enc_id = enc_id;
return res;
error:
+ enc_id = ENCLOSURE_ID_NULL;
goto exit;
}
diff --git a/src/sdis_scene_c.h b/src/sdis_scene_c.h
@@ -26,6 +26,9 @@
#include <limits.h>
+#define MEDIUM_ID_MULTI UINT_MAX
+#define ENCLOSURE_ID_NULL UINT_MAX
+
struct prim_prop {
struct sdis_interface* interf;
unsigned front_enclosure; /* Id of the front facing enclosure */
@@ -71,8 +74,6 @@ medium_init(struct mem_allocator* allocator, struct sdis_medium** medium)
*medium = NULL;
}
-#define ENCLOSURE_MULTI_MEDIA UINT_MAX
-
struct enclosure {
struct s2d_scene_view* s2d_view;
struct s3d_scene_view* s3d_view;
@@ -97,7 +98,7 @@ enclosure_init(struct mem_allocator* allocator, struct enclosure* enc)
enc->S_over_V = 0;
enc->V = 0;
enc->hc_upper_bound = 0;
- enc->medium_id = ENCLOSURE_MULTI_MEDIA;
+ enc->medium_id = MEDIUM_ID_MULTI;
}
static INLINE void
@@ -229,25 +230,31 @@ scene_get_interface
const unsigned iprim);
extern LOCAL_SYM res_T
-scene_get_medium
+scene_get_enclosure_id
(struct sdis_scene* scene,
const double position[],
- struct sdis_medium** medium);
+ unsigned* enclosure_id);
-/* 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.
+/* This function assumes that the position under test lies within a finite
+ * enclosure. The enclosure in which it is located is therefore retrieved by
+ * tracing a random ray around the current position. For infinite enclosures,
+ * you need to use the `scene_get_enclosure_id' function, which in turn may take
+ * longer.
*
- * 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). */
+ * Note that the function actually calls scene_get_enclosure internally if no
+ * valid enclosure is found with the normal procedure. This may be due to
+ * numerical problems or incorrect assumptions about the current enclosure (its
+ * limits are open to infinity). */
extern LOCAL_SYM res_T
-scene_get_medium_in_closed_boundaries
- (struct sdis_scene* scn,
+scene_get_enclosure_id_in_closed_boundaries
+ (struct sdis_scene* scene,
const double position[],
+ unsigned* enclosure_id);
+
+extern LOCAL_SYM res_T
+scene_get_enclosure_medium
+ (struct sdis_scene* scene,
+ const struct enclosure* enclosure,
struct sdis_medium** medium);
extern LOCAL_SYM res_T
diff --git a/src/sdis_solve_camera.c b/src/sdis_solve_camera.c
@@ -84,7 +84,7 @@ static res_T
solve_pixel
(struct sdis_scene* scn,
struct ssp_rng* rng,
- struct sdis_medium* mdm,
+ const unsigned enc_id,
const struct sdis_camera* cam,
const double time_range[2], /* Observation time */
const size_t ipix[2], /* Pixel coordinate in the image space */
@@ -99,7 +99,7 @@ solve_pixel
struct sdis_heat_path* pheat_path = NULL;
size_t irealisation;
res_T res = RES_OK;
- ASSERT(scn && mdm && rng && cam && ipix && nrealisations);
+ ASSERT(scn && rng && cam && ipix && nrealisations);
ASSERT(pix_sz && pix_sz[0] > 0 && pix_sz[1] > 0);
ASSERT(pixel && time_range);
@@ -136,7 +136,7 @@ solve_pixel
/* Launch the realisation */
realis_args.rng = rng;
- realis_args.medium = mdm;
+ realis_args.enc_id = enc_id;
realis_args.time = time;
realis_args.picard_order = picard_order;
realis_args.heat_path = pheat_path;
@@ -195,7 +195,7 @@ static res_T
solve_tile
(struct sdis_scene* scn,
struct ssp_rng* rng,
- struct sdis_medium* mdm,
+ const unsigned enc_id,
const struct sdis_camera* cam,
const double time_range[2],
const size_t tile_org[2], /* Origin of the tile in pixel space */
@@ -211,7 +211,7 @@ solve_tile
size_t mcode; /* Morton code of the tile pixel */
size_t npixels;
res_T res = RES_OK;
- ASSERT(scn && rng && mdm && cam && spp);
+ ASSERT(scn && rng && cam && spp);
ASSERT(tile_size && tile_size[0] && tile_size[1]);
ASSERT(pix_sz && pix_sz[0] > 0 && pix_sz[1] > 0 && time_range);
@@ -241,8 +241,8 @@ solve_tile
estimator = estimator_buffer_grab(buf, ipix_image[0], ipix_image[1]);
}
res = solve_pixel
- (scn, rng, mdm, cam, time_range, ipix_image, spp, register_paths, pix_sz,
- picard_order, diff_algo, estimator, pixel);
+ (scn, rng, enc_id, cam, time_range, ipix_image, spp, register_paths,
+ pix_sz, picard_order, diff_algo, estimator, pixel);
if(res != RES_OK) goto error;
}
@@ -504,12 +504,14 @@ sdis_solve_camera
/* Stardis variables */
struct sdis_estimator_buffer* buf = NULL;
- struct sdis_medium* medium = NULL;
/* Random number generators */
struct ssp_rng_proxy* rng_proxy = NULL;
struct ssp_rng** per_thread_rng = NULL;
+ /* Enclosure & medium in which the probe lies */
+ unsigned enc_id = ENCLOSURE_ID_NULL;
+
/* Miscellaneous */
size_t ntiles_x, ntiles_y, ntiles, ntiles_adjusted;
size_t ntiles_proc; /* #tiles for the current proc */
@@ -537,17 +539,9 @@ sdis_solve_camera
}
/* Retrieve the medium in which the submitted position lies */
- res = scene_get_medium(scn, args->cam->position, &medium);
+ res = scene_get_enclosure_id(scn, args->cam->position, &enc_id);
if(res != RES_OK) goto error;
- if(medium->type != SDIS_FLUID) {
- log_err(scn->dev,
- "%s: the camera position `%g %g %g' must be in a fluid medium.\n",
- FUNC_NAME, SPLIT3(args->cam->position));
- res = RES_BAD_ARG;
- goto error;
- }
-
/* Create the per thread RNGs */
res = create_per_thread_rng
(scn->dev, args->rng_state, args->rng_type, &rng_proxy, &per_thread_rng);
@@ -658,7 +652,7 @@ sdis_solve_camera
/* Draw the tile */
res_local = solve_tile
- (scn, rng, medium, args->cam, args->time_range, tile_org, tile_sz,
+ (scn, rng, enc_id, args->cam, args->time_range, tile_org, tile_sz,
args->spp, register_paths, pix_sz, args->picard_order, args->diff_algo,
buf, tile);
if(res_local != RES_OK) {
diff --git a/src/sdis_solve_medium_Xd.h b/src/sdis_solve_medium_Xd.h
@@ -39,7 +39,7 @@
*/
struct enclosure_cumul {
- const struct enclosure* enc;
+ unsigned enc_id;
double cumul;
};
@@ -79,12 +79,13 @@ compute_medium_enclosure_cumulative
while(!htable_enclosure_iterator_eq(&it, &end)) {
struct enclosure_cumul enc_cumul;
const struct enclosure* enc = htable_enclosure_iterator_data_get(&it);
+ const unsigned* enc_id = htable_enclosure_iterator_key_get(&it);
htable_enclosure_iterator_next(&it);
if(sdis_medium_get_id(mdm) != enc->medium_id) continue;
accum += enc->V;
- enc_cumul.enc = enc;
+ enc_cumul.enc_id = *enc_id;
enc_cumul.cumul = accum;
res = darray_enclosure_cumul_push_back(cumul, &enc_cumul);
if(res != RES_OK) goto error;
@@ -105,7 +106,7 @@ error:
goto exit;
}
-static const struct enclosure*
+static unsigned
sample_medium_enclosure
(const struct darray_enclosure_cumul* cumul, struct ssp_rng* rng)
{
@@ -137,7 +138,7 @@ sample_medium_enclosure
enc_cumul_found = enc_cumuls + i;
}
- return enc_cumul_found->enc;
+ return enc_cumul_found->enc_id;
}
static INLINE res_T
@@ -217,17 +218,22 @@ check_compute_power_args(const struct sdis_compute_power_args* args)
******************************************************************************/
static res_T
XD(sample_enclosure_position)
- (const struct enclosure* enc,
+ (struct sdis_scene* scn,
+ const unsigned enc_id,
struct ssp_rng* rng,
double pos[DIM])
{
const size_t MAX_NCHALLENGES = 1000;
+
+ const struct enclosure* enc = NULL;
float lower[DIM], upper[DIM];
size_t ichallenge;
size_t i;
res_T res = RES_OK;
- ASSERT(enc && rng && pos);
+ ASSERT(scn && rng && pos);
+ ASSERT(enc_id != ENCLOSURE_ID_NULL);
+ enc = scene_get_enclosure(scn, enc_id);
SXD(scene_view_get_aabb(enc->sXd(view), lower, upper));
FOR_EACH(i, 0, DIM) {
@@ -396,13 +402,13 @@ XD(solve_medium)
struct accum* acc_time = &per_thread_acc_time[ithread];
struct green_path_handle* pgreen_path = NULL;
struct green_path_handle green_path = GREEN_PATH_HANDLE_NULL;
- const struct enclosure* enc = NULL;
struct sdis_heat_path* pheat_path = NULL;
struct sdis_heat_path heat_path;
double weight;
double time;
double pos[DIM];
size_t n;
+ unsigned enc_id = ENCLOSURE_ID_NULL;
int pcent;
res_T res_local = RES_OK;
res_T res_simul = RES_OK;
@@ -425,8 +431,8 @@ XD(solve_medium)
/* Uniformly Sample an enclosure that surround the submitted medium and
* uniformly sample a position into it */
- enc = sample_medium_enclosure(&cumul, rng);
- res_local = XD(sample_enclosure_position)(enc, rng, pos);
+ enc_id = sample_medium_enclosure(&cumul, rng);
+ res_local = XD(sample_enclosure_position)(scn, enc_id, rng, pos);
if(res_local != RES_OK) {
log_err(scn->dev, "%s: could not sample a medium position.\n", FUNC_NAME);
ATOMIC_SET(&res, res_local);
@@ -435,7 +441,7 @@ XD(solve_medium)
/* Run a probe realisation */
realis_args.rng = rng;
- realis_args.medium = args->medium;
+ realis_args.enc_id = enc_id;
realis_args.time = time;
realis_args.picard_order = args->picard_order;
realis_args.green_path = pgreen_path;
@@ -685,10 +691,10 @@ XD(compute_power)
struct ssp_rng* rng = per_thread_rng[ithread];
struct accum* acc_mpow = &per_thread_acc_mpow[ithread];
struct accum* acc_time = &per_thread_acc_time[ithread];
- const struct enclosure* enc = NULL;
double power = 0;
double usec = 0;
size_t n = 0;
+ unsigned enc_id = ENCLOSURE_ID_NULL;
int pcent = 0;
res_T res_local = RES_OK;
@@ -702,8 +708,8 @@ XD(compute_power)
/* Uniformly Sample an enclosure that surround the submitted medium and
* uniformly sample a position into it */
- enc = sample_medium_enclosure(&cumul, rng);
- res_local = XD(sample_enclosure_position)(enc, rng, vtx.P);
+ enc_id = sample_medium_enclosure(&cumul, rng);
+ res_local = XD(sample_enclosure_position)(scn, enc_id, rng, vtx.P);
if(res_local != RES_OK) {
log_err(scn->dev, "%s: could not sample a medium position.\n", FUNC_NAME);
ATOMIC_SET(&res, res_local);
diff --git a/src/sdis_solve_probe_Xd.h b/src/sdis_solve_probe_Xd.h
@@ -116,7 +116,9 @@ XD(solve_one_probe)
struct accum* acc_temp,
struct accum* acc_time)
{
- struct sdis_medium* medium = NULL; /* Medium in which the probe lies */
+ /* Enclosure in which the probe lies */
+ unsigned enc_id = ENCLOSURE_ID_NULL;
+
size_t irealisation = 0;
res_T res = RES_OK;
ASSERT(scn && rng && check_solve_probe_args(args) == RES_OK);
@@ -126,7 +128,7 @@ XD(solve_one_probe)
*acc_time = ACCUM_NULL;
/* Retrieve the medium in which the submitted position lies */
- res = scene_get_medium(scn, args->position, &medium);
+ res = scene_get_enclosure_id(scn, args->position, &enc_id);
if(res != RES_OK) goto error;
FOR_EACH(irealisation, 0, args->nrealisations) {
@@ -144,7 +146,7 @@ XD(solve_one_probe)
/* Run a realisation */
realis_args.rng = rng;
- realis_args.medium = medium;
+ realis_args.enc_id = enc_id;
realis_args.time = time;
realis_args.picard_order = args->picard_order;
realis_args.irealisation = irealisation;
@@ -191,7 +193,6 @@ XD(solve_probe)
size_t nthreads = 0;
/* Stardis variables */
- struct sdis_medium* medium = NULL;
struct sdis_estimator* estimator = NULL;
struct sdis_green_function* green = NULL;
struct sdis_green_function** per_thread_green = NULL;
@@ -200,6 +201,10 @@ XD(solve_probe)
struct ssp_rng_proxy* rng_proxy = NULL;
struct ssp_rng** per_thread_rng = NULL;
+ /* Enclosure in which the probe lies */
+ const struct enclosure* enc = NULL;
+ unsigned enc_id = ENCLOSURE_ID_NULL;
+
/* Miscellaneous */
struct accum* per_thread_acc_temp = NULL;
struct accum* per_thread_acc_time = NULL;
@@ -255,10 +260,19 @@ XD(solve_probe)
if(!per_thread_acc_temp) { res = RES_MEM_ERR; goto error; }
if(!per_thread_acc_time) { res = RES_MEM_ERR; goto error; }
- /* Retrieve the medium in which the submitted position lies */
- res = scene_get_medium(scn, args->position, &medium);
+ /* Retrieve the enclosure in which the submitted position lies */
+ res = scene_get_enclosure_id(scn, args->position, &enc_id);
if(res != RES_OK) goto error;
+ /* Check that the enclosure does not contain multiple materials */
+ enc = scene_get_enclosure(scn, enc_id);
+ if(enc->medium_id == MEDIUM_ID_MULTI) {
+ log_err(scn->dev, "%s: probe is in an enclosure with several media "
+ "-- pos=("FORMAT_VECX")\n", FUNC_NAME, SPLITX(args->position));
+ res = RES_BAD_OP;
+ goto error;
+ }
+
/* Create the per thread green function */
if(out_green) {
res = create_per_thread_green_function
@@ -327,7 +341,7 @@ XD(solve_probe)
/* Invoke the probe realisation */
realis_args.rng = rng;
- realis_args.medium = medium;
+ realis_args.enc_id = enc_id;
realis_args.time = time;
realis_args.picard_order = args->picard_order;
realis_args.green_path = pgreen_path;
diff --git a/src/test_sdis_solve_camera.c b/src/test_sdis_solve_camera.c
@@ -33,27 +33,38 @@
#define SPP 32 /* #Samples per pixel, i.e. #realisations per pixel */
/*
- * The scene is composed of a solid cube whose temperature is unknown. The
- * emissivity of the cube is 1 and its convection coefficient with the
- * surrounding fluid at 300K is 0.1. At the center of the cube there is a spherical
- * fluid cavity whose temperature is 350K. The convection coefficient between
- * the solid and the cavity is 1 and the emissivity of this interface is null.
- * The ambient radiative temperature of the system is 300K.
+ * The scene consists of a solid cube whose temperature is unknown. The
+ * emissivity of the cube is 1 and the convection coefficient with the
+ * surrounding fluid at 300 K is 0.1. At the center of the cube is a spherical
+ * cavity of fluid with a temperature of 350 K. The convection coefficient
+ * between the solid and the cavity is 1, and the emissivity of this interface
+ * is zero. The ambient radiative temperature of the system is 300 K.
*
- * In this test, we compute the radiative temperature that reaches a camera
- * that looks the cube and we `dump' a normalized image of the result.
+ * Finally, a parallelepiped below the cube symbolizes the ground. The
+ * temperature of its Robin condition is 280 K. This geometry verifies that a
+ * camera can draw a scene in an enclosure containing several media, such as
+ * those used to define several boundary conditions.
*
- * (1,1,1)
- * +----------------+
- * /' # # /|
- * +----*--------*--+ |
- * | ' # # | |
- * | ' # 350K # | |
- * | ' # # | |
- * | +.....#..#.....|.+
- * |/ |/
- * +----------------+
- * (-1,-1,-1)
+ * In this test, we calculate the radiative temperature that reaches a camera
+ * looking at the cube and produce an image of the result written to the
+ * standard output in htrdr-image(5) format.
+ *
+ *
+ * +----------------+
+ * /' # # /|
+ * +----*--------*--+ | __\ 300 K
+ * | ' # # | | / /
+ * | ' # 350K # | | \__/
+ * | ' # # | |
+ * | +.....#..#.....|.+
+ * |/ |/
+ * +----------------+ 280K
+ * +---------------------------------__\------+
+ * / / / /|
+ * / \__/ / +
+ * +------------------------------------------+ /
+ * | |/
+ * +------------------------------------------+
*/
/*******************************************************************************
@@ -82,6 +93,52 @@ struct geometry {
static const struct geometry GEOMETRY_NULL = {NULL, NULL, NULL};
static void
+geometry_add_shape
+ (struct geometry* geom,
+ const double* positions,
+ const size_t nverts,
+ const size_t* indices,
+ const size_t nprims,
+ const double transform[12], /* May be NULL <=> no transformation */
+ struct sdis_interface* interf)
+{
+ struct map_interf* geom_interf = NULL;
+ size_t nverts_prev = 0;
+ size_t i;
+
+ CHK(geom != NULL);
+ CHK(positions != NULL);
+ CHK(indices != NULL);
+ CHK(nverts != 0);
+ CHK(nprims != 0);
+ CHK(interf != NULL);
+
+ /* Save the previous number of vertices/primitives of the geometry */
+ nverts_prev = sa_size(geom->positions) / 3;
+
+ /* Add the vertices */
+ FOR_EACH(i, 0, nverts) {
+ double* pos = sa_add(geom->positions, 3);
+ d3_set(pos, positions + i*3);
+ if(transform) {
+ d33_muld3(pos, transform, pos);
+ d3_add(pos, transform+9, pos);
+ }
+ }
+
+ /* Add the indices */
+ FOR_EACH(i, 0, nprims) {
+ sa_push(geom->indices, indices[i*3+0] + nverts_prev);
+ sa_push(geom->indices, indices[i*3+1] + nverts_prev);
+ sa_push(geom->indices, indices[i*3+2] + nverts_prev);
+ }
+
+ geom_interf = sa_add(geom->interfaces, 1);
+ geom_interf->key = sa_size(geom->indices) / 3 - 1;
+ geom_interf->interf = interf;
+}
+
+static void
geometry_release(struct geometry* geom)
{
CHK(geom != NULL);
@@ -143,43 +200,59 @@ geometry_get_interface
*bound = interf->interf;
}
-/*******************************************************************************
- * Fluid medium
- ******************************************************************************/
-struct fluid {
- double cp;
- double rho;
- double temperature;
-};
-static const struct fluid FLUID_NULL = {0, 0, SDIS_TEMPERATURE_NONE};
-
-static double
-fluid_get_calorific_capacity
- (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+static void
+add_cube(struct geometry* geom, struct sdis_interface* interf)
{
- CHK(data != NULL && vtx != NULL);
- return ((const struct fluid*)sdis_data_cget(data))->cp;
+ struct s3dut_mesh_data msh_data;
+ struct s3dut_mesh* msh = NULL;
+ CHK(geom);
+
+ OK(s3dut_create_cuboid(NULL, 2, 2, 2, &msh));
+ OK(s3dut_mesh_get_data(msh, &msh_data));
+ geometry_add_shape(geom, msh_data.positions, msh_data.nvertices,
+ msh_data.indices, msh_data.nprimitives, NULL, interf);
+ OK(s3dut_mesh_ref_put(msh));
}
-static double
-fluid_get_volumic_mass
- (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+static void
+add_sphere(struct geometry* geom, struct sdis_interface* interf)
{
- CHK(data != NULL && vtx != NULL);
- return ((const struct fluid*)sdis_data_cget(data))->rho;
+ struct s3dut_mesh_data msh_data;
+ struct s3dut_mesh* msh = NULL;
+ CHK(geom);
+
+ OK(s3dut_create_sphere(NULL, 0.5, 32, 16, &msh));
+ OK(s3dut_mesh_get_data(msh, &msh_data));
+ geometry_add_shape(geom, msh_data.positions, msh_data.nvertices,
+ msh_data.indices, msh_data.nprimitives, NULL, interf);
+ OK(s3dut_mesh_ref_put(msh));
}
-static double
-fluid_get_temperature
- (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+static void
+add_ground(struct geometry* geom, struct sdis_interface* interf)
{
- CHK(data != NULL && vtx != NULL);
- return ((const struct fluid*)sdis_data_cget(data))->temperature;
+ struct s3dut_mesh_data msh_data;
+ struct s3dut_mesh* msh = NULL;
+ const double transform[12] = {1,0,0, 0,1,0, 0,0,1, 0,0,-4};
+ CHK(geom);
+
+ OK(s3dut_create_cuboid(NULL, 10, 10, 2, &msh));
+ OK(s3dut_mesh_get_data(msh, &msh_data));
+ geometry_add_shape(geom, msh_data.positions, msh_data.nvertices,
+ msh_data.indices, msh_data.nprimitives, transform, interf);
+ OK(s3dut_mesh_ref_put(msh));
}
/*******************************************************************************
- * Solid medium
+ * Media
******************************************************************************/
+struct fluid {
+ double cp;
+ double rho;
+ double temperature;
+};
+static const struct fluid FLUID_NULL = {0, 0, SDIS_TEMPERATURE_NONE};
+
struct solid {
double cp;
double lambda;
@@ -189,48 +262,87 @@ struct solid {
};
static const struct solid SOLID_NULL = {0, 0, 0, 0, SDIS_TEMPERATURE_NONE};
-static double
-solid_get_calorific_capacity
- (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+#define MEDIUM_GETTER(Type, Prop) \
+ static double \
+ Type##_get_##Prop \
+ (const struct sdis_rwalk_vertex* vtx, \
+ struct sdis_data* data) \
+ { \
+ CHK(data && vtx); \
+ return ((const struct Type*)sdis_data_cget(data))->Prop; \
+ }
+/* Fluid getters */
+MEDIUM_GETTER(fluid, cp)
+MEDIUM_GETTER(fluid, rho)
+MEDIUM_GETTER(fluid, temperature)
+/* Solid getters */
+MEDIUM_GETTER(solid, cp)
+MEDIUM_GETTER(solid, lambda)
+MEDIUM_GETTER(solid, rho)
+MEDIUM_GETTER(solid, delta)
+MEDIUM_GETTER(solid, temperature)
+#undef MEDIUM_GETTER
+
+static struct sdis_medium*
+create_fluid
+ (struct sdis_device* dev,
+ const struct fluid* param)
{
- CHK(data != NULL && vtx != NULL);
- return ((const struct solid*)sdis_data_cget(data))->cp;
-}
+ struct sdis_fluid_shader fluid_shader = DUMMY_FLUID_SHADER;
+ struct sdis_data* data = NULL;
+ struct sdis_medium* fluid = NULL;
-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;
-}
+ CHK(param != NULL);
-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;
-}
+ /* Copy the fluid parameters into the Stardis memory space */
+ OK(sdis_data_create
+ (dev, sizeof(struct fluid), ALIGNOF(struct fluid), NULL, &data));
+ memcpy(sdis_data_get(data), param, sizeof(struct fluid));
-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;
+ /* Setup the fluid shader */
+ fluid_shader.calorific_capacity = fluid_get_cp;
+ fluid_shader.volumic_mass = fluid_get_rho;
+ fluid_shader.temperature = fluid_get_temperature;
+
+ /* Create the fluid medium */
+ OK(sdis_fluid_create(dev, &fluid_shader, data, &fluid));
+ OK(sdis_data_ref_put(data));
+
+ return fluid;
}
-static double
-solid_get_temperature
- (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+static struct sdis_medium*
+create_solid
+ (struct sdis_device* dev,
+ const struct solid* param)
{
- CHK(data != NULL && vtx != NULL);
- return ((const struct solid*)sdis_data_cget(data))->temperature;
+ struct sdis_solid_shader solid_shader = DUMMY_SOLID_SHADER;
+ struct sdis_data* data = NULL;
+ struct sdis_medium* solid = NULL;
+
+ CHK(param != NULL);
+
+ /* Copy the solid parameters into the Stardis memory space */
+ OK(sdis_data_create
+ (dev, sizeof(struct solid), ALIGNOF(struct solid), NULL, &data));
+ memcpy(sdis_data_get(data), param, sizeof(struct solid));
+
+ /* Setup the solid shader */
+ solid_shader.calorific_capacity = solid_get_cp;
+ solid_shader.thermal_conductivity = solid_get_lambda;
+ solid_shader.volumic_mass = solid_get_rho;
+ solid_shader.delta = solid_get_delta;
+ solid_shader.temperature = solid_get_temperature;
+
+ /* Create the solid medium */
+ OK(sdis_solid_create(dev, &solid_shader, data, &solid));
+ OK(sdis_data_ref_put(data));
+
+ return solid;
}
/*******************************************************************************
- * Interface
+ * Interfaces
******************************************************************************/
struct interf {
double hc;
@@ -243,50 +355,76 @@ static const struct interf INTERF_NULL = {
0, 0, 0, SDIS_TEMPERATURE_NONE, SDIS_TEMPERATURE_NONE
};
-static double
-interface_get_convection_coef
- (const struct sdis_interface_fragment* frag, struct sdis_data* data)
-{
- CHK(data != NULL && frag != NULL);
- return ((const struct interf*)sdis_data_cget(data))->hc;
-}
+#define INTERFACE_GETTER(Prop) \
+ static double \
+ interface_get_##Prop \
+ (const struct sdis_interface_fragment* frag, \
+ struct sdis_data* data) \
+ { \
+ CHK(data && frag); \
+ return ((const struct interf*)sdis_data_cget(data))->Prop; \
+ }
+INTERFACE_GETTER(hc)
+INTERFACE_GETTER(temperature)
+INTERFACE_GETTER(reference_temperature)
+#undef INTERFACE_GETTER
+
+#define INTERFACE_GETTER(Prop) \
+ static double \
+ interface_get_##Prop \
+ (const struct sdis_interface_fragment* frag, \
+ const unsigned source_id, \
+ struct sdis_data* data) \
+ { \
+ CHK(data && frag); \
+ (void)source_id; \
+ return ((const struct interf*)sdis_data_cget(data))->Prop; \
+ }
+INTERFACE_GETTER(epsilon)
+INTERFACE_GETTER(specular_fraction)
+#undef INTERFACE_GETTER
-static double
-interface_get_emissivity
- (const struct sdis_interface_fragment* frag,
- const unsigned source_id,
- struct sdis_data* data)
+static struct sdis_interface*
+create_interface
+ (struct sdis_device* dev,
+ struct sdis_medium* mdm_front,
+ struct sdis_medium* mdm_back,
+ const struct interf* param)
{
- (void)source_id;
- CHK(data != NULL && frag != NULL);
- return ((const struct interf*)sdis_data_cget(data))->epsilon;
-}
+ struct sdis_interface_shader interface_shader = SDIS_INTERFACE_SHADER_NULL;
+ struct sdis_data* data = NULL;
+ struct sdis_interface* interf = NULL;
-static double
-interface_get_specular_fraction
- (const struct sdis_interface_fragment* frag,
- const unsigned source_id,
- struct sdis_data* data)
-{
- (void)source_id;
- CHK(data != NULL && frag != NULL);
- return ((const struct interf*)sdis_data_cget(data))->specular_fraction;
-}
+ CHK(mdm_front != NULL);
+ CHK(mdm_back != NULL);
-static double
-interface_get_temperature
- (const struct sdis_interface_fragment* frag, struct sdis_data* data)
-{
- CHK(data != NULL && frag != NULL);
- return ((const struct interf*)sdis_data_cget(data))->temperature;
-}
+ /* Copy the interface parameters into the Stardis memory space */
+ OK(sdis_data_create
+ (dev, sizeof(struct interf), ALIGNOF(struct interf), NULL, &data));
+ memcpy(sdis_data_get(data), param, sizeof(struct interf));
-static double
-interface_get_reference_temperature
- (const struct sdis_interface_fragment* frag, struct sdis_data* data)
-{
- CHK(data != NULL && frag != NULL);
- return ((const struct interf*)sdis_data_cget(data))->reference_temperature;
+ /* Setup the interface shader */
+ interface_shader.convection_coef = interface_get_hc;
+ interface_shader.front.temperature = interface_get_temperature;
+ interface_shader.back.temperature = interface_get_temperature;
+ if(sdis_medium_get_type(mdm_front) == SDIS_FLUID) {
+ interface_shader.front.emissivity = interface_get_epsilon;
+ interface_shader.front.specular_fraction = interface_get_specular_fraction;
+ interface_shader.front.reference_temperature =
+ interface_get_reference_temperature;
+ }
+ if(sdis_medium_get_type(mdm_back) == SDIS_FLUID) {
+ interface_shader.back.emissivity = interface_get_epsilon;
+ interface_shader.back.specular_fraction = interface_get_specular_fraction;
+ interface_shader.back.reference_temperature =
+ interface_get_reference_temperature;
+ }
+ /* Create the interface */
+ OK(sdis_interface_create
+ (dev, mdm_front, mdm_back, &interface_shader, data, &interf));
+ OK(sdis_data_ref_put(data));
+
+ return interf;
}
/*******************************************************************************
@@ -339,254 +477,289 @@ create_radenv
}
/*******************************************************************************
- * Helper functions
+ * Scene
******************************************************************************/
-static void
-create_solid
+static struct sdis_scene*
+create_scene
(struct sdis_device* dev,
- const struct solid* param,
- struct sdis_medium** solid)
+ struct sdis_radiative_env* radenv,
+ struct geometry* geom)
{
- struct sdis_data* data = NULL;
- struct solid* solid_args = NULL;
- struct sdis_solid_shader solid_shader = DUMMY_SOLID_SHADER;
+ struct sdis_scene_create_args scn_args = SDIS_SCENE_CREATE_ARGS_DEFAULT;
+ struct sdis_scene* scn = NULL;
+ size_t ntris = 0;
+ size_t npos = 0;
+ CHK(geom);
- CHK(param != NULL);
- CHK(solid != NULL);
+ ntris = sa_size(geom->indices) / 3; /* #primitives */
+ npos = sa_size(geom->positions) / 3; /* #positions */
- /* Copy the solid parameters into the Stardis memory space */
- OK(sdis_data_create
- (dev, sizeof(struct solid), ALIGNOF(struct solid), NULL, &data));
- solid_args = sdis_data_get(data);
- memcpy(solid_args, param, sizeof(struct solid));
+ scn_args.get_indices = geometry_get_indices;
+ scn_args.get_interface = geometry_get_interface;
+ scn_args.get_position = geometry_get_position;
+ scn_args.nprimitives = ntris;
+ scn_args.nvertices = npos;
+ scn_args.t_range[0] = 300;
+ scn_args.t_range[1] = 350;
+ scn_args.radenv = radenv;
+ scn_args.context = geom;
- /* Setup the solid shader */
- solid_shader.calorific_capacity = solid_get_calorific_capacity;
- 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_get_delta;
- solid_shader.temperature = solid_get_temperature;
+ OK(sdis_scene_create(dev, &scn_args, &scn));
+ return scn;
+}
- /* Create the solid medium */
- OK(sdis_solid_create(dev, &solid_shader, data, solid));
+/*******************************************************************************
+ * Rendering point of view
+ ******************************************************************************/
+static struct sdis_camera*
+create_camera(struct sdis_device* dev)
+{
+ const double pos[3] = {3, 3, 3};
+ const double tgt[3] = {0, 0, 0};
+ const double up[3] = {0, 0, 1};
+ struct sdis_camera* cam = NULL;
- /* Release the ownership onto the Stardis memory space storing the solid
- * parameters */
- OK(sdis_data_ref_put(data));
+ OK(sdis_camera_create(dev, &cam));
+ OK(sdis_camera_set_proj_ratio(cam, (double)IMG_WIDTH/(double)IMG_HEIGHT));
+ OK(sdis_camera_set_fov(cam, MDEG2RAD(70)));
+ OK(sdis_camera_look_at(cam, pos, tgt, up));
+ return cam;
}
+/*******************************************************************************
+ * Draw the scene
+ ******************************************************************************/
static void
-create_fluid
- (struct sdis_device* dev,
- const struct fluid* param,
- struct sdis_medium** fluid)
+check_pixel(const struct sdis_estimator* pixel)
{
- struct sdis_data* data = NULL;
- struct fluid* fluid_args = NULL;
- struct sdis_fluid_shader fluid_shader = DUMMY_FLUID_SHADER;
-
- CHK(param != NULL);
- CHK(fluid != NULL);
-
- /* Copy the fluid parameters into the Stardis memory space */
- OK(sdis_data_create
- (dev, sizeof(struct fluid), ALIGNOF(struct fluid), NULL, &data));
- fluid_args = sdis_data_get(data);
- memcpy(fluid_args, param, sizeof(struct fluid));
-
- /* Setup the fluid shader */
- fluid_shader.calorific_capacity = fluid_get_calorific_capacity;
- fluid_shader.volumic_mass = fluid_get_volumic_mass;
- fluid_shader.temperature = fluid_get_temperature;
+ struct sdis_mc T = SDIS_MC_NULL;
+ struct sdis_mc time = SDIS_MC_NULL;
+ size_t nfails = 0;
+ size_t nreals = 0;
- /* Create the fluid medium */
- OK(sdis_fluid_create(dev, &fluid_shader, data, fluid));
+ OK(sdis_estimator_get_realisation_count(pixel, &nreals));
+ OK(sdis_estimator_get_failure_count(pixel, &nfails));
+ OK(sdis_estimator_get_temperature(pixel, &T));
+ OK(sdis_estimator_get_realisation_time(pixel, &time));
- /* Release the ownership onto the Stardis memory space storing the fluid
- * parameters */
- OK(sdis_data_ref_put(data));
+ CHK(nreals + nfails == SPP);
+ CHK(T.E > 0);
+ CHK(time.E > 0);
}
static void
-create_interface
- (struct sdis_device* dev,
- struct sdis_medium* mdm_front,
- struct sdis_medium* mdm_back,
- const struct interf* param,
- struct sdis_interface** interf)
+check_image(const struct sdis_estimator_buffer* img)
{
- struct sdis_data* data = NULL;
- struct interf* interface_args = NULL;
- struct sdis_interface_shader interface_shader = SDIS_INTERFACE_SHADER_NULL;
-
- CHK(mdm_front != NULL);
- CHK(mdm_back != NULL);
- CHK(param != NULL);
- CHK(interf != NULL);
-
- /* Copy the interface parameters into the Stardis memory space */
- OK(sdis_data_create
- (dev, sizeof(struct interf), ALIGNOF(struct interf), NULL, &data));
- interface_args = sdis_data_get(data);
- memcpy(interface_args, param, sizeof(struct interf));
-
- /* Setup the interface shader */
- interface_shader.convection_coef = interface_get_convection_coef;
- interface_shader.front.temperature = interface_get_temperature;
- interface_shader.back.temperature = interface_get_temperature;
- if(sdis_medium_get_type(mdm_front) == SDIS_FLUID) {
- interface_shader.front.emissivity = interface_get_emissivity;
- interface_shader.front.specular_fraction = interface_get_specular_fraction;
- interface_shader.front.reference_temperature =
- interface_get_reference_temperature;
- }
- if(sdis_medium_get_type(mdm_back) == SDIS_FLUID) {
- interface_shader.back.emissivity = interface_get_emissivity;
- interface_shader.back.specular_fraction = interface_get_specular_fraction;
- interface_shader.back.reference_temperature =
- interface_get_reference_temperature;
+ struct sdis_mc T = SDIS_MC_NULL;
+ struct sdis_mc time = SDIS_MC_NULL;
+ const struct sdis_estimator* pixel = NULL;
+ struct ssp_rng* rng_state = NULL;
+ size_t definition[2];
+ size_t nreals, nfails;
+ size_t x, y;
+
+ BA(sdis_estimator_buffer_get_realisation_count(NULL, &nreals));
+ BA(sdis_estimator_buffer_get_realisation_count(img, NULL));
+ OK(sdis_estimator_buffer_get_realisation_count(img, &nreals));
+
+ BA(sdis_estimator_buffer_get_failure_count(NULL, &nfails));
+ BA(sdis_estimator_buffer_get_failure_count(img, NULL));
+ OK(sdis_estimator_buffer_get_failure_count(img, &nfails));
+
+ BA(sdis_estimator_buffer_get_temperature(NULL, &T));
+ BA(sdis_estimator_buffer_get_temperature(img, NULL));
+ OK(sdis_estimator_buffer_get_temperature(img, &T));
+
+ BA(sdis_estimator_buffer_get_realisation_time(NULL, &time));
+ BA(sdis_estimator_buffer_get_realisation_time(img, NULL));
+ OK(sdis_estimator_buffer_get_realisation_time(img, &time));
+
+ BA(sdis_estimator_buffer_get_rng_state(NULL, &rng_state));
+ BA(sdis_estimator_buffer_get_rng_state(img, NULL));
+ OK(sdis_estimator_buffer_get_rng_state(img, &rng_state));
+
+ CHK(nreals + nfails == IMG_WIDTH*IMG_HEIGHT*SPP);
+
+ fprintf(stderr, "Overall temperature ~ %g +/- %g\n", T.E, T.SE);
+ fprintf(stderr, "Time per realisation (in usec) ~ %g +/- %g\n", time.E, time.SE);
+ fprintf(stderr, "#failures = %lu/%lu\n",
+ (unsigned long)nfails, (unsigned long)(IMG_WIDTH*IMG_HEIGHT*SPP));
+
+ BA(sdis_estimator_buffer_get_definition(NULL, definition));
+ BA(sdis_estimator_buffer_get_definition(img, NULL));
+ OK(sdis_estimator_buffer_get_definition(img, definition));
+ CHK(definition[0] == IMG_WIDTH);
+ CHK(definition[1] == IMG_HEIGHT);
+
+ BA(sdis_estimator_buffer_at(NULL, 0, 0, &pixel));
+ BA(sdis_estimator_buffer_at(img, IMG_WIDTH+1,0 , &pixel));
+ BA(sdis_estimator_buffer_at(img, 0, IMG_HEIGHT+1, &pixel));
+ BA(sdis_estimator_buffer_at(img, 0, 0, NULL));
+
+ /* Pixels ordered by row */
+ FOR_EACH(y, 0, IMG_HEIGHT) {
+ FOR_EACH(x, 0, IMG_WIDTH) {
+ OK(sdis_estimator_buffer_at(img, x, y, &pixel));
+ check_pixel(pixel);
+ }
}
- /* Create the interface */
- OK(sdis_interface_create
- (dev, mdm_front, mdm_back, &interface_shader, data, interf));
+}
- /* Release the ownership onto the Stardis memory space storing the interface
- * parameters */
- OK(sdis_data_ref_put(data));
+/* Check that the images, although compatible from an estimation point of view,
+ * are not the same. This should be the case if different random sequences are
+ * used for each image */
+static void
+check_image_difference
+ (const struct sdis_estimator_buffer* img0,
+ const struct sdis_estimator_buffer* img1)
+{
+ struct sdis_mc T0 = SDIS_MC_NULL;
+ struct sdis_mc T1 = SDIS_MC_NULL;
+ size_t definition0[2];
+ size_t definition1[2];
+ size_t nreals0, nfails0;
+ size_t nreals1, nfails1;
+
+ OK(sdis_estimator_buffer_get_realisation_count(img0, &nreals0));
+ OK(sdis_estimator_buffer_get_realisation_count(img1, &nreals1));
+
+ OK(sdis_estimator_buffer_get_failure_count(img0, &nfails0));
+ OK(sdis_estimator_buffer_get_failure_count(img1, &nfails1));
+ CHK(nreals0 + nfails0 == nreals1 + nfails1);
+
+ OK(sdis_estimator_buffer_get_definition(img0, definition0));
+ OK(sdis_estimator_buffer_get_definition(img1, definition1));
+ CHK(definition0[0] == definition1[0]);
+ CHK(definition0[1] == definition1[1]);
+
+ OK(sdis_estimator_buffer_get_temperature(img0, &T0));
+ OK(sdis_estimator_buffer_get_temperature(img1, &T1));
+ CHK(T0.E != T1.E || (T0.SE == 0 && T1.SE == 0));
+ CHK(T0.E + 3*T0.SE >= T1.E - 3*T1.SE);
+ CHK(T0.E - 3*T0.SE <= T1.E + 3*T1.SE);
}
+/* Write an image in htrdr-image(5) format */
static void
-geometry_add_shape
- (struct geometry* geom,
- const double* positions,
- const size_t nverts,
- const size_t* indices,
- const size_t nprims,
- const double transform[12], /* May be NULL <=> no transformation */
- struct sdis_interface* interf)
+write_image(FILE* stream, struct sdis_estimator_buffer* img)
{
- struct map_interf* geom_interf = NULL;
- size_t nverts_prev = 0;
- size_t i;
+ size_t x, y;
- CHK(geom != NULL);
- CHK(positions != NULL);
- CHK(indices != NULL);
- CHK(nverts != 0);
- CHK(nprims != 0);
- CHK(interf != NULL);
+ /* Header */
+ fprintf(stream, "%d %d\n", IMG_WIDTH, IMG_HEIGHT);
- /* Save the previous number of vertices/primitives of the geometry */
- nverts_prev = sa_size(geom->positions) / 3;
+ /* Pixels ordered by row */
+ FOR_EACH(y, 0, IMG_HEIGHT) {
+ FOR_EACH(x, 0, IMG_WIDTH) {
+ struct sdis_mc T = SDIS_MC_NULL;
+ const struct sdis_estimator* pixel = NULL;
- /* Add the vertices */
- FOR_EACH(i, 0, nverts) {
- double* pos = sa_add(geom->positions, 3);
- d3_set(pos, positions + i*3);
- if(transform) {
- d33_muld3(pos, transform, pos);
- d3_add(pos, transform+9, pos);
+ OK(sdis_estimator_buffer_at(img, x, y, &pixel));
+ OK(sdis_estimator_get_temperature(pixel, &T));
+ fprintf(stream, "%g %g 0 0 0 0 0 0\n", T.E, T.SE);
}
}
+}
- /* Add the indices */
- FOR_EACH(i, 0, nprims) {
- sa_push(geom->indices, indices[i*3+0] + nverts_prev);
- sa_push(geom->indices, indices[i*3+1] + nverts_prev);
- sa_push(geom->indices, indices[i*3+2] + nverts_prev);
- }
-
- geom_interf = sa_add(geom->interfaces, 1);
- geom_interf->key = sa_size(geom->indices) / 3 - 1;
- geom_interf->interf = interf;
+static void
+invalid_draw(struct sdis_scene* scn, struct sdis_camera* cam)
+{
+ struct sdis_solve_camera_args args = SDIS_SOLVE_CAMERA_ARGS_DEFAULT;
+ struct sdis_estimator_buffer* img = NULL;
+
+ CHK(cam && scn);
+
+ args.cam = cam;
+ args.time_range[0] = INF;
+ args.time_range[1] = INF;
+ args.image_definition[0] = IMG_WIDTH;
+ args.image_definition[1] = IMG_HEIGHT;
+ args.spp = SPP;
+
+ BA(sdis_solve_camera(NULL, &args, &img));
+ BA(sdis_solve_camera(scn, NULL, &img));
+ BA(sdis_solve_camera(scn, &args, NULL));
+
+ /* Invalid camera */
+ args.cam = NULL;
+ BA(sdis_solve_camera(scn, &args, &img));
+ args.cam = cam;
+
+ /* Invald time range */
+ args.time_range[0] = args.time_range[1] = -1;
+ BA(sdis_solve_camera(scn, &args, &img));
+ args.time_range[0] = 1;
+ BA(sdis_solve_camera(scn, &args, &img));
+ args.time_range[1] = 0;
+ BA(sdis_solve_camera(scn, &args, &img));
+ args.time_range[0] = args.time_range[1] = INF;
+
+ /* Invalid image definition */
+ args.image_definition[0] = 0;
+ BA(sdis_solve_camera(scn, &args, &img));
+ args.image_definition[0] = IMG_WIDTH;
+ args.image_definition[1] = 0;
+ BA(sdis_solve_camera(scn, &args, &img));
+ args.image_definition[1] = IMG_HEIGHT;
+
+ /* Invalid number of samples per pixel */
+ args.spp = 0;
+ BA(sdis_solve_camera(scn, &args, &img));
+ args.spp = SPP;
}
static void
-dump_image(const struct sdis_estimator_buffer* buf)
+draw
+ (struct sdis_scene* scn,
+ struct sdis_camera* cam,
+ const int is_master_process)
{
- struct image img;
- double* temps = NULL;
- double Tmax = -DBL_MAX;
- double Tmin = DBL_MAX;
- double norm;
- size_t definition[2];
- size_t ix, iy;
-
- CHK(buf != NULL);
- OK(sdis_estimator_buffer_get_definition(buf, definition));
-
- temps = mem_alloc(definition[0]*definition[1]*sizeof(double));
- CHK(temps != NULL);
-
- /* Check the results validity */
- FOR_EACH(iy, 0, definition[1]) {
- FOR_EACH(ix, 0, definition[0]) {
- const struct sdis_estimator* estimator;
- struct sdis_mc T;
- struct sdis_mc time;
- size_t nreals;
- size_t nfails;
-
- BA(sdis_estimator_buffer_at(NULL, ix, iy, &estimator));
- BA(sdis_estimator_buffer_at(buf, definition[0]+1, iy, &estimator));
- BA(sdis_estimator_buffer_at(buf, ix, definition[1]+1, &estimator));
- BA(sdis_estimator_buffer_at(buf, ix, iy, NULL));
- OK(sdis_estimator_buffer_at(buf, ix, iy, &estimator));
-
- OK(sdis_estimator_get_realisation_count(estimator, &nreals));
- OK(sdis_estimator_get_failure_count(estimator, &nfails));
- OK(sdis_estimator_get_temperature(estimator, &T));
- OK(sdis_estimator_get_realisation_time(estimator, &time));
-
- CHK(nreals + nfails == SPP);
- CHK(T.E > 0);
- CHK(time.E > 0);
- }
- }
+ struct sdis_solve_camera_args args = SDIS_SOLVE_CAMERA_ARGS_DEFAULT;
+ struct sdis_estimator_buffer* img0 = NULL;
+ struct sdis_estimator_buffer* img1 = NULL;
+ struct ssp_rng* rng = NULL;
- /* Compute the per pixel temperature */
- FOR_EACH(iy, 0, definition[1]) {
- double* row = temps + iy * definition[0];
- FOR_EACH(ix, 0, definition[0]) {
- const struct sdis_estimator* estimator;
- struct sdis_mc T;
+ args.cam = cam;
+ args.time_range[0] = INF;
+ args.time_range[1] = INF;
+ args.image_definition[0] = IMG_WIDTH;
+ args.image_definition[1] = IMG_HEIGHT;
+ args.spp = SPP;
- OK(sdis_estimator_buffer_at(buf, ix, iy, &estimator));
- OK(sdis_estimator_get_temperature(estimator, &T));
+ OK(sdis_solve_camera(scn, &args, &img0));
- row[ix] = T.E;
- Tmax = MMAX(row[ix], Tmax);
- Tmin = MMIN(row[ix], Tmin);
- }
- }
- if(Tmax != Tmin) {
- norm = Tmax - Tmin;
+ if(!is_master_process) {
+ CHK(img0 == NULL);
} else {
- Tmin = 0;
- norm = 1;
+ check_image(img0);
+ write_image(stdout, img0);
}
- /* Allocate the image memory space */
- OK(image_init(NULL, &img));
- OK(image_setup(&img, IMG_WIDTH, IMG_HEIGHT, IMG_WIDTH*3, IMAGE_RGB8, NULL));
-
- FOR_EACH(iy, 0, definition[1]) {
- const double* src_row = temps + iy*definition[0];
- char* dst_row = img.pixels + iy*img.pitch;
- FOR_EACH(ix, 0, definition[0]) {
- unsigned char* pixels = (unsigned char*)
- (dst_row + ix * sizeof_image_format(img.format));
- const unsigned char T = (unsigned char)
- ((src_row[ix] - Tmin)/ norm * 255.0);
- pixels[0] = T;
- pixels[1] = T;
- pixels[2] = T;
- }
+ /* Provide a RNG state and check that the results, although compatible, are
+ * not the same */
+ OK(ssp_rng_create(NULL, SSP_RNG_THREEFRY, &rng));
+ OK(ssp_rng_discard(rng, 3141592653589)); /* Move the RNG state */
+ args.rng_state = rng;
+ args.rng_type = SSP_RNG_TYPE_NULL;
+ OK(sdis_solve_camera(scn, &args, &img1));
+ if(is_master_process) {
+ check_image_difference(img0, img1);
+ OK(sdis_estimator_buffer_ref_put(img1));
+ }
+ OK(ssp_rng_ref_put(rng));
+
+ /* Change the RNG type and check that the results, although compatible, are
+ * not the same */
+ args.rng_state = NULL;
+ args.rng_type = SDIS_SOLVE_CAMERA_ARGS_DEFAULT.rng_type == SSP_RNG_THREEFRY
+ ? SSP_RNG_MT19937_64 : SSP_RNG_THREEFRY;
+ OK(sdis_solve_camera(scn, &args, &img1));
+ if(is_master_process) {
+ check_image_difference(img0, img1);
+ OK(sdis_estimator_buffer_ref_put(img1));
}
- OK(image_write_ppm_stream(&img, 0/*binary?*/, stdout));
- image_release(&img);
- mem_rm(temps);
+
+ if(is_master_process) OK(sdis_estimator_buffer_ref_put(img0));
}
/*******************************************************************************
@@ -596,56 +769,44 @@ int
main(int argc, char** argv)
{
struct geometry geom = GEOMETRY_NULL;
- struct s3dut_mesh* msh = NULL;
- struct s3dut_mesh_data msh_data;
- struct sdis_mc T = SDIS_MC_NULL;
- struct sdis_mc T2 = SDIS_MC_NULL;
- struct sdis_mc time = SDIS_MC_NULL;
struct sdis_camera* cam = NULL;
struct sdis_device* dev = NULL;
- struct sdis_estimator_buffer* buf = NULL;
- struct sdis_estimator_buffer* buf2 = NULL;
struct sdis_medium* solid = NULL;
struct sdis_medium* fluid0 = NULL;
struct sdis_medium* fluid1 = NULL;
+ struct sdis_medium* fluid2 = NULL;
struct sdis_interface* interf0 = NULL;
struct sdis_interface* interf1 = NULL;
+ struct sdis_interface* interf2 = NULL;
struct sdis_radiative_env* radenv = NULL;
struct sdis_scene* scn = NULL;
- struct sdis_scene_create_args scn_args = SDIS_SCENE_CREATE_ARGS_DEFAULT;
- struct sdis_solve_camera_args solve_args = SDIS_SOLVE_CAMERA_ARGS_DEFAULT;
- struct ssp_rng* rng = NULL;
- struct ssp_rng* rng_state = NULL;
struct fluid fluid_args = FLUID_NULL;
struct solid solid_args = SOLID_NULL;
struct interf interface_args = INTERF_NULL;
- struct fluid* pfluid_args = NULL;
- struct radenv* pradenv_args = NULL;
- size_t ntris, npos;
- size_t nreals, nfails;
- size_t definition[2];
- double pos[3];
- double tgt[3];
- double up[3];
int is_master_process;
(void)argc, (void)argv;
create_default_device(&argc, &argv, &is_master_process, &dev);
radenv = create_radenv(dev, 300, 300);
- pradenv_args = sdis_data_get(sdis_radiative_env_get_data(radenv));
/* Create the fluid0 */
fluid_args.temperature = 350;
fluid_args.rho = 0;
fluid_args.cp = 0;
- create_fluid(dev, &fluid_args, &fluid0);
+ fluid0 = create_fluid(dev, &fluid_args);
/* Create the fluid1 */
fluid_args.temperature = 300;
fluid_args.rho = 0;
fluid_args.cp = 0;
- create_fluid(dev, &fluid_args, &fluid1);
+ fluid1 = create_fluid(dev, &fluid_args);
+
+ /* Create the fluid2 */
+ fluid_args.temperature = 280;
+ fluid_args.rho = 0;
+ fluid_args.cp = 0;
+ fluid2 = create_fluid(dev, &fluid_args);
/* Create the solid medium */
solid_args.cp = 1.0;
@@ -653,14 +814,14 @@ main(int argc, char** argv)
solid_args.rho = 1.0;
solid_args.delta = 1.0/20.0;
solid_args.temperature = SDIS_TEMPERATURE_NONE;
- create_solid(dev, &solid_args, &solid);
+ solid = create_solid(dev, &solid_args);
/* Create the fluid0/solid interface */
interface_args.hc = 1;
interface_args.epsilon = 0;
interface_args.specular_fraction = 0;
interface_args.temperature = SDIS_TEMPERATURE_NONE;
- create_interface(dev, solid, fluid0, &interface_args, &interf0);
+ interf0 = create_interface(dev, solid, fluid0, &interface_args);
/* Create the fluid1/solid interface */
interface_args.hc = 0.1;
@@ -668,178 +829,45 @@ main(int argc, char** argv)
interface_args.specular_fraction = 1;
interface_args.temperature = SDIS_TEMPERATURE_NONE;
interface_args.reference_temperature = 300;
- create_interface(dev, fluid1, solid, &interface_args, &interf1);
-
- /* Setup the cube geometry */
- OK(s3dut_create_cuboid(NULL, 2, 2, 2, &msh));
- OK(s3dut_mesh_get_data(msh, &msh_data));
- geometry_add_shape(&geom, msh_data.positions, msh_data.nvertices,
- msh_data.indices, msh_data.nprimitives, NULL, interf1);
- OK(s3dut_mesh_ref_put(msh));
+ interf1 = create_interface(dev, fluid1, solid, &interface_args);
- /* Setup the sphere geometry */
- OK(s3dut_create_sphere(NULL, 0.5, 32, 16, &msh));
- OK(s3dut_mesh_get_data(msh, &msh_data));
- geometry_add_shape(&geom, msh_data.positions, msh_data.nvertices,
- msh_data.indices, msh_data.nprimitives, NULL, interf0);
- OK(s3dut_mesh_ref_put(msh));
+ /* Create the fluid2/ground interface */
+ interface_args.hc = 0.2;
+ interface_args.epsilon = 1;
+ interface_args.specular_fraction = 1;
+ interface_args.temperature = SDIS_TEMPERATURE_NONE;
+ interface_args.reference_temperature = 300;
+ interf2 = create_interface(dev, fluid2, solid, &interface_args);
- /* Setup the scene */
- ntris = sa_size(geom.indices) / 3; /* #primitives */
- npos = sa_size(geom.positions) / 3; /* #positions */
- scn_args.get_indices = geometry_get_indices;
- scn_args.get_interface = geometry_get_interface;
- scn_args.get_position = geometry_get_position;
- scn_args.nprimitives = ntris;
- scn_args.nvertices = npos;
- scn_args.t_range[0] = 300;
- scn_args.t_range[1] = 350;
- scn_args.radenv = radenv;
- scn_args.context = &geom;
- OK(sdis_scene_create(dev, &scn_args, &scn));
+ add_cube(&geom, interf1);
+ add_sphere(&geom, interf0);
+ add_ground(&geom, interf2);
#if 0
dump_mesh(stdout, geom.positions, sa_size(geom.positions)/3, geom.indices,
sa_size(geom.indices)/3);
#endif
- /* Setup the camera */
- d3(pos, 3, 3, 3);
- d3(tgt, 0, 0, 0);
- d3(up, 0, 0, 1);
- OK(sdis_camera_create(dev, &cam));
- OK(sdis_camera_set_proj_ratio(cam, (double)IMG_WIDTH/(double)IMG_HEIGHT));
- OK(sdis_camera_set_fov(cam, MDEG2RAD(70)));
- OK(sdis_camera_look_at(cam, pos, tgt, up));
-
-#if 0
- dump_mesh(stdout, geom.positions, npos, geom.indices, ntris);
- exit(0);
-#endif
- solve_args.cam = cam;
- solve_args.time_range[0] = INF;
- solve_args.time_range[0] = INF;
- solve_args.image_definition[0] = IMG_WIDTH;
- solve_args.image_definition[1] = IMG_HEIGHT;
- solve_args.spp = SPP;
-
- BA(sdis_solve_camera(NULL, &solve_args, &buf));
- BA(sdis_solve_camera(scn, NULL, &buf));
- BA(sdis_solve_camera(scn, &solve_args, NULL));
- solve_args.cam = NULL;
- BA(sdis_solve_camera(scn, &solve_args, &buf));
- solve_args.cam = cam;
- pradenv_args->temperature = SDIS_TEMPERATURE_NONE;
- BA(sdis_solve_camera(scn, &solve_args, &buf));
- pradenv_args->temperature = 300;
- solve_args.time_range[0] = solve_args.time_range[1] = -1;
- BA(sdis_solve_camera(scn, &solve_args, &buf));
- solve_args.time_range[0] = 1;
- BA(sdis_solve_camera(scn, &solve_args, &buf));
- solve_args.time_range[1] = 0;
- BA(sdis_solve_camera(scn, &solve_args, &buf));
- solve_args.time_range[0] = solve_args.time_range[1] = INF;
-
- /* Launch the simulation */
- OK(sdis_solve_camera(scn, &solve_args, &buf));
+ scn = create_scene(dev, radenv, &geom);
+ cam = create_camera(dev);
- if(!is_master_process) {
- CHK(buf == NULL);
- } else {
- BA(sdis_estimator_buffer_get_realisation_count(NULL, &nreals));
- BA(sdis_estimator_buffer_get_realisation_count(buf, NULL));
- OK(sdis_estimator_buffer_get_realisation_count(buf, &nreals));
-
- BA(sdis_estimator_buffer_get_failure_count(NULL, &nfails));
- BA(sdis_estimator_buffer_get_failure_count(buf, NULL));
- OK(sdis_estimator_buffer_get_failure_count(buf, &nfails));
-
- BA(sdis_estimator_buffer_get_temperature(NULL, &T));
- BA(sdis_estimator_buffer_get_temperature(buf, NULL));
- OK(sdis_estimator_buffer_get_temperature(buf, &T));
-
- BA(sdis_estimator_buffer_get_realisation_time(NULL, &time));
- BA(sdis_estimator_buffer_get_realisation_time(buf, NULL));
- OK(sdis_estimator_buffer_get_realisation_time(buf, &time));
-
- BA(sdis_estimator_buffer_get_rng_state(NULL, &rng_state));
- BA(sdis_estimator_buffer_get_rng_state(buf, NULL));
- OK(sdis_estimator_buffer_get_rng_state(buf, &rng_state));
-
- CHK(nreals + nfails == IMG_WIDTH*IMG_HEIGHT*SPP);
-
- fprintf(stderr, "Overall temperature ~ %g +/- %g\n", T.E, T.SE);
- fprintf(stderr, "Time per realisation (in usec) ~ %g +/- %g\n", time.E, time.SE);
- fprintf(stderr, "#failures = %lu/%lu\n",
- (unsigned long)nfails, (unsigned long)(IMG_WIDTH*IMG_HEIGHT*SPP));
-
- BA(sdis_estimator_buffer_get_definition(NULL, definition));
- BA(sdis_estimator_buffer_get_definition(buf, NULL));
- OK(sdis_estimator_buffer_get_definition(buf, definition));
- CHK(definition[0] == IMG_WIDTH);
- CHK(definition[1] == IMG_HEIGHT);
-
- /* Write the image */
- dump_image(buf);
- OK(sdis_estimator_buffer_ref_put(buf));
- }
-
- /* Check RNG type */
- solve_args.rng_state = NULL;
- solve_args.rng_type = SSP_RNG_TYPE_NULL;
- BA(sdis_solve_camera(scn, &solve_args, &buf2));
- solve_args.rng_type =
- SDIS_SOLVE_CAMERA_ARGS_DEFAULT.rng_type == SSP_RNG_THREEFRY
- ? SSP_RNG_MT19937_64 : SSP_RNG_THREEFRY;
- OK(sdis_solve_camera(scn, &solve_args, &buf2));
- if(is_master_process) {
- OK(sdis_estimator_buffer_get_temperature(buf2, &T2));
- CHK(T.E != T2.E || (T2.SE == 0 && T.SE ==0));
- CHK(T2.E + 3*T2.SE >= T.E - 3*T.SE
- && T2.E - 3*T2.SE <= T.E + 3*T.SE);
- OK(sdis_estimator_buffer_ref_put(buf2));
- }
-
- /* Check the RNG state */
- OK(ssp_rng_create(NULL, SSP_RNG_THREEFRY, &rng));
- OK(ssp_rng_discard(rng, 3141592653589)); /* Move the RNG state */
- solve_args.rng_state = rng;
- solve_args.rng_type = SSP_RNG_TYPE_NULL;
- OK(sdis_solve_camera(scn, &solve_args, &buf2));
- OK(ssp_rng_ref_put(rng));
- if(is_master_process) {
- OK(sdis_estimator_buffer_get_temperature(buf2, &T2));
- CHK(T.E != T2.E || (T2.SE == 0 && T.SE == 0));
- CHK(T2.E + 3*T2.SE >= T.E - 3*T.SE
- && T2.E - 3*T2.SE <= T.E + 3*T.SE);
- OK(sdis_estimator_buffer_ref_put(buf2));
- }
-
- /* Restore args */
- solve_args.rng_state = SDIS_SOLVE_CAMERA_ARGS_DEFAULT.rng_state;
- solve_args.rng_type = SDIS_SOLVE_CAMERA_ARGS_DEFAULT.rng_type;
-
- pfluid_args = sdis_data_get(sdis_medium_get_data(fluid1));
- pfluid_args->temperature = SDIS_TEMPERATURE_NONE;
-
- /* Check simulation error handling */
- BA(sdis_solve_camera(scn, &solve_args, &buf));
- solve_args.register_paths = SDIS_HEAT_PATH_ALL;
- BA(sdis_solve_camera(scn, &solve_args, &buf));
+ invalid_draw(scn, cam);
+ draw(scn, cam, is_master_process);
/* Release memory */
OK(sdis_medium_ref_put(solid));
OK(sdis_medium_ref_put(fluid0));
OK(sdis_medium_ref_put(fluid1));
+ OK(sdis_medium_ref_put(fluid2));
OK(sdis_radiative_env_ref_put(radenv));
OK(sdis_scene_ref_put(scn));
OK(sdis_camera_ref_put(cam));
OK(sdis_interface_ref_put(interf0));
OK(sdis_interface_ref_put(interf1));
+ OK(sdis_interface_ref_put(interf2));
free_default_device(dev);
geometry_release(&geom);
CHK(mem_allocated_size() == 0);
return 0;
}
-
