commit 094c34a2d8b5cca1d48e7a3da5156b7149f754ae
parent 55901e8717c25ddb035749c86e8119b1084c92ba
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Thu, 17 Apr 2025 16:25:36 +0200
Improve the robustness of the sampling of radiative paths
An corner position can lead to the unexpected entry of a radiative
trajectory into a solid. This is in fact the same problem corrected in
commit 8d0142a for the calculation of the incident diffuse flux and
therefore the same solution: the problem is mitigated by moving the
position of the angle slightly away when this unexpected behaviour
occurs.
This commit refactors the calculation of the external net flux and the
sampling of the radiative path in order to share functions, in
particular those that manage the problems mentioned above.
Diffstat:
8 files changed, 661 insertions(+), 530 deletions(-)
diff --git a/Makefile b/Makefile
@@ -30,6 +30,7 @@ MPI_DEF = -DSDIS_ENABLE_MPI
MPI_SRC = src/sdis_mpi.c
SRC =\
src/sdis.c\
+ src/sdis_brdf.c\
src/sdis_camera.c\
src/sdis_data.c\
src/sdis_device.c\
diff --git a/src/sdis_brdf.c b/src/sdis_brdf.c
@@ -0,0 +1,127 @@
+/* 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.h"
+#include "sdis_brdf.h"
+#include "sdis_interface_c.h"
+
+#include <star/ssp.h>
+#include <rsys/double3.h>
+
+/*******************************************************************************
+ * Helper functions
+ ******************************************************************************/
+static INLINE res_T
+check_brdf_setup_args(const struct brdf_setup_args* args)
+{
+ const struct sdis_medium* mdm = NULL;
+
+ if(!args) return RES_BAD_ARG;
+
+ if(args->interf == NULL || args->frag == NULL)
+ return RES_BAD_ARG;
+
+ switch(args->frag->side) {
+ case SDIS_FRONT: mdm = args->interf->medium_front; break;
+ case SDIS_BACK: mdm = args->interf->medium_back; break;
+ default: FATAL("Unreachable code\n"); break;
+ }
+
+ if(sdis_medium_get_type(mdm) != SDIS_FLUID)
+ return RES_BAD_ARG;
+
+ return RES_OK;
+}
+
+/* Reflect the V wrt the normal N. By convention V points outward the surface.
+ * In fact, this function is a double-precision version of the reflect_3d
+ * function. TODO Clean this "repeat" */
+static FINLINE double*
+reflect(double res[3], const double V[3], const double N[3])
+{
+ double tmp[3];
+ double cos_V_N;
+ ASSERT(res && V && N);
+ ASSERT(d3_is_normalized(V) && d3_is_normalized(N));
+ cos_V_N = d3_dot(V, N);
+ d3_muld(tmp, N, 2*cos_V_N);
+ d3_sub(res, tmp, V);
+ return res;
+}
+
+/*******************************************************************************
+ * Local functions
+ ******************************************************************************/
+res_T
+brdf_setup
+ (struct sdis_device* dev,
+ const struct brdf_setup_args* args,
+ struct brdf* brdf)
+{
+ res_T res = RES_OK;
+
+ ASSERT(check_brdf_setup_args(args) == RES_OK);
+
+ #define GET(Attr) { \
+ brdf->Attr = interface_side_get_##Attr \
+ (args->interf, args->source_id, args->frag); \
+ \
+ res = interface_side_check_##Attr \
+ (dev, brdf->Attr, args->frag->P, args->frag->time); \
+ if(res != RES_OK) goto error; \
+ } (void)0
+
+ GET(emissivity);
+ GET(specular_fraction);
+
+ #undef GET
+
+exit:
+ return res;
+error:
+ *brdf = BRDF_NULL;
+ goto exit;
+}
+
+void
+brdf_sample
+ (const struct brdf* brdf,
+ struct ssp_rng* rng,
+ const double wi[3], /* Incident direction. Point away from the surface */
+ const double N[3], /* Surface normal */
+ struct brdf_sample* sample)
+{
+ double r = 0; /* Random number */
+
+ /* Preconditions */
+ ASSERT(brdf && rng && wi && N && sample);
+ ASSERT(d3_is_normalized(wi) && d3_is_normalized(N));
+ ASSERT(d3_dot(wi, N) > 0);
+
+ r = ssp_rng_canonical(rng);
+
+ /* Sample the specular part */
+ if(r < brdf->specular_fraction) {
+ reflect(sample->dir, wi, N);
+ sample->pdf = 1;
+ sample->cpnt = BRDF_SPECULAR;
+
+ /* Sample the diffuse part */
+ } else {
+ ssp_ran_hemisphere_cos(rng, N, sample->dir, NULL);
+ sample->pdf = 1.0/PI;
+ sample->cpnt = BRDF_DIFFUSE;
+ }
+}
diff --git a/src/sdis_brdf.h b/src/sdis_brdf.h
@@ -0,0 +1,71 @@
+/* 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/>. */
+
+#ifndef SDIS_BRDF_H
+#define SDIS_BRDF_H
+
+#include <rsys/rsys.h>
+
+enum brdf_component {
+ BRDF_SPECULAR,
+ BRDF_DIFFUSE,
+ BRDF_NONE
+};
+
+struct brdf_sample {
+ double dir[3];
+ double pdf;
+ enum brdf_component cpnt;
+};
+#define BRDF_SAMPLE_NULL__ {{0}, 0, BRDF_NONE}
+static const struct brdf_sample BRDF_SAMPLE_NULL = BRDF_SAMPLE_NULL__;
+
+struct brdf {
+ double emissivity;
+ double specular_fraction;
+};
+#define BRDF_NULL__ {0, 0}
+static const struct brdf BRDF_NULL = BRDF_NULL__;
+
+/* forward declarations */
+struct sdis_device;
+struct sdis_interface;
+struct sdis_interface_fragment;
+struct ssp_rng;
+
+struct brdf_setup_args {
+ struct sdis_interface* interf;
+ struct sdis_interface_fragment* frag;
+ unsigned source_id;
+};
+#define BRDF_SETUP_ARGS_NULL__ {NULL, NULL, SDIS_INTERN_SOURCE_ID}
+static const struct brdf_setup_args BRDF_SETUP_ARGS_NULL =
+ BRDF_SETUP_ARGS_NULL__;
+
+extern LOCAL_SYM res_T
+brdf_setup
+ (struct sdis_device* dev,
+ const struct brdf_setup_args* args,
+ struct brdf* brdf);
+
+extern LOCAL_SYM void
+brdf_sample
+ (const struct brdf* brdf,
+ struct ssp_rng* rng,
+ const double wi[3], /* Incident direction. Point away from the surface */
+ const double N[3], /* Surface normal */
+ struct brdf_sample* sample);
+
+#endif /* SDIS_BRDF_H */
diff --git a/src/sdis_heat_path.h b/src/sdis_heat_path.h
@@ -327,6 +327,47 @@ radiative_path_3d
struct ssp_rng* rng,
struct temperature* temperature);
+extern LOCAL_SYM void
+trace_ray_2d
+ (const struct sdis_scene* scn,
+ const double pos[2],
+ const double dir[3], /* Always in 3D */
+ const double distance,
+ const struct s2d_hit* hit_from,
+ struct s2d_hit* hit);
+
+extern LOCAL_SYM void
+trace_ray_3d
+ (const struct sdis_scene* scn,
+ const double pos[3],
+ const double dir[3], /* Always in 3D */
+ const double distance,
+ const struct s3d_hit* hit_from,
+ struct s3d_hit* hit);
+
+/* Trace a ray and setup the fragment at the intersection found, if any. */
+extern LOCAL_SYM res_T
+find_next_fragment_2d
+ (const struct sdis_scene* scn,
+ const double in_pos[2],
+ const double in_dir[3], /* Always in 3D */
+ const struct s2d_hit* in_hit,
+ const double time,
+ struct s2d_hit* out_hit,
+ struct sdis_interface** out_interf,
+ struct sdis_interface_fragment* out_frag);
+
+extern LOCAL_SYM res_T
+find_next_fragment_3d
+ (const struct sdis_scene* scn,
+ const double in_pos[3],
+ const double in_dir[3], /* Always in 3D */
+ const struct s3d_hit* in_hit,
+ const double time,
+ struct s3d_hit* out_hit,
+ struct sdis_interface** out_interf,
+ struct sdis_interface_fragment* out_frag);
+
/*******************************************************************************
* Convective path
******************************************************************************/
diff --git a/src/sdis_heat_path_boundary_Xd_c.h b/src/sdis_heat_path_boundary_Xd_c.h
@@ -246,130 +246,6 @@ XD(sample_reinjection_dir)
#endif
}
-
-#if DIM == 2
-static void
-XD(move_away_primitive_boundaries)
- (const struct sXd(hit)* hit,
- const double delta,
- double position[DIM]) /* Position to move */
-{
- struct sXd(attrib) attr;
- float pos[DIM];
- float dir[DIM];
- float len;
- const float st = 0.5f;
- ASSERT(!SXD_HIT_NONE(hit) && delta > 0);
-
- SXD(primitive_get_attrib(&hit->prim, SXD_POSITION, st, &attr));
-
- fX_set_dX(pos, position);
- fX(sub)(dir, attr.value, pos);
- len = fX(normalize)(dir, dir);
- len = MMIN(len, (float)(delta*0.1));
-
- XD(move_pos)(position, dir, len);
-}
-#else
-/* Move the submitted position away from the primitive boundaries to avoid
- * numerical issues leading to inconsistent random walks. */
-static void
-XD(move_away_primitive_boundaries)
- (const struct sXd(hit)* hit,
- const double delta,
- double position[DIM])
-{
- struct s3d_attrib v0, v1, v2; /* Triangle vertices */
- float E[3][4]; /* 3D edge equations */
- float dst[3]; /* Distance from current position to edge equation */
- float N[3]; /* Triangle normal */
- float P[3]; /* Random walk position */
- float tmp[3];
- float min_dst, max_dst;
- float cos_a1, cos_a2;
- float len;
- int imax = 0;
- int imin = 0;
- int imid = 0;
- int i;
- ASSERT(delta > 0 && !S3D_HIT_NONE(hit));
-
- fX_set_dX(P, position);
-
- /* Fetch triangle vertices */
- S3D(triangle_get_vertex_attrib(&hit->prim, 0, S3D_POSITION, &v0));
- S3D(triangle_get_vertex_attrib(&hit->prim, 1, S3D_POSITION, &v1));
- S3D(triangle_get_vertex_attrib(&hit->prim, 2, S3D_POSITION, &v2));
-
- /* Compute the edge vector */
- f3_sub(E[0], v1.value, v0.value);
- f3_sub(E[1], v2.value, v1.value);
- f3_sub(E[2], v0.value, v2.value);
-
- /* Compute the triangle normal */
- f3_cross(N, E[1], E[0]);
-
- /* Compute the 3D edge equation */
- f3_normalize(E[0], f3_cross(E[0], E[0], N));
- f3_normalize(E[1], f3_cross(E[1], E[1], N));
- f3_normalize(E[2], f3_cross(E[2], E[2], N));
- E[0][3] = -f3_dot(E[0], v0.value);
- E[1][3] = -f3_dot(E[1], v1.value);
- E[2][3] = -f3_dot(E[2], v2.value);
-
- /* Compute the distance from current position to the edges */
- dst[0] = f3_dot(E[0], P) + E[0][3];
- dst[1] = f3_dot(E[1], P) + E[1][3];
- dst[2] = f3_dot(E[2], P) + E[2][3];
-
- /* Retrieve the min and max distance from random walk position to triangle
- * edges */
- min_dst = MMIN(MMIN(dst[0], dst[1]), dst[2]);
- max_dst = MMAX(MMAX(dst[0], dst[1]), dst[2]);
-
- /* Sort the edges with respect to their distance to the random walk position */
- FOR_EACH(i, 0, 3) {
- if(dst[i] == min_dst) {
- imin = i;
- } else if(dst[i] == max_dst) {
- imax = i;
- } else {
- imid = i;
- }
- }
- (void)imax;
-
- /* TODO if the current position is near a vertex, one should move toward the
- * farthest edge along its normal to avoid too small displacement */
-
- /* Compute the distance `dst' from the current position to the edges to move
- * to, along the normal of the edge from which the random walk is the nearest
- *
- * +. cos(a) = d / dst => dst = d / cos_a
- * / `*.
- * / | `*.
- * / dst| a /`*.
- * / | / `*.
- * / | / d `*.
- * / |/ `*.
- * +---------o----------------+ */
- cos_a1 = f3_dot(E[imin], f3_minus(tmp, E[imid]));
- cos_a2 = f3_dot(E[imin], f3_minus(tmp, E[imax]));
- dst[imid] = cos_a1 > 0 ? dst[imid] / cos_a1 : FLT_MAX;
- dst[imax] = cos_a2 > 0 ? dst[imax] / cos_a2 : FLT_MAX;
-
- /* Compute the maximum displacement distance into the triangle along the
- * normal of the edge from which the random walk is the nearest */
- len = MMIN(dst[imid], dst[imax]);
- ASSERT(len != FLT_MAX);
-
- /* Define the displacement distance as the minimum between 10 percent of
- * delta and len / 2. */
- len = MMIN(len*0.5f, (float)(delta*0.1));
- XD(move_pos)(position, E[imin], len);
-}
-#endif
-
static res_T
XD(find_reinjection_ray)
(struct sdis_scene* scn,
@@ -421,7 +297,13 @@ XD(find_reinjection_ray)
do {
fX_set_dX(org, ray->org);
filter_data.XD(hit) = args->rwalk->XD(hit);
+
+ /* Limit the epsilon to 1.e-6, as Star-3D's single-precision floating-point
+ * representation will inevitably present numerical accuracy problems below
+ * this threshold. There's no point in going any lower */
+ /*filter_data.epsilon = MMAX(args->distance * 0.01, 1e-6);*/
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
@@ -1183,4 +1065,127 @@ error:
goto exit;
}
+#if DIM == 2
+void
+XD(move_away_primitive_boundaries)
+ (const struct sXd(hit)* hit,
+ const double delta,
+ double position[DIM]) /* Position to move */
+{
+ struct sXd(attrib) attr;
+ float pos[DIM];
+ float dir[DIM];
+ float len;
+ const float st = 0.5f;
+ ASSERT(!SXD_HIT_NONE(hit) && delta > 0);
+
+ SXD(primitive_get_attrib(&hit->prim, SXD_POSITION, st, &attr));
+
+ fX_set_dX(pos, position);
+ fX(sub)(dir, attr.value, pos);
+ len = fX(normalize)(dir, dir);
+ len = MMIN(len, (float)(delta*0.1));
+
+ XD(move_pos)(position, dir, len);
+}
+#else
+/* Move the submitted position away from the primitive boundaries to avoid
+ * numerical issues leading to inconsistent random walks. */
+void
+XD(move_away_primitive_boundaries)
+ (const struct sXd(hit)* hit,
+ const double delta,
+ double position[DIM])
+{
+ struct s3d_attrib v0, v1, v2; /* Triangle vertices */
+ float E[3][4]; /* 3D edge equations */
+ float dst[3]; /* Distance from current position to edge equation */
+ float N[3]; /* Triangle normal */
+ float P[3]; /* Random walk position */
+ float tmp[3];
+ float min_dst, max_dst;
+ float cos_a1, cos_a2;
+ float len;
+ int imax = 0;
+ int imin = 0;
+ int imid = 0;
+ int i;
+ ASSERT(delta > 0 && !S3D_HIT_NONE(hit));
+
+ fX_set_dX(P, position);
+
+ /* Fetch triangle vertices */
+ S3D(triangle_get_vertex_attrib(&hit->prim, 0, S3D_POSITION, &v0));
+ S3D(triangle_get_vertex_attrib(&hit->prim, 1, S3D_POSITION, &v1));
+ S3D(triangle_get_vertex_attrib(&hit->prim, 2, S3D_POSITION, &v2));
+
+ /* Compute the edge vector */
+ f3_sub(E[0], v1.value, v0.value);
+ f3_sub(E[1], v2.value, v1.value);
+ f3_sub(E[2], v0.value, v2.value);
+
+ /* Compute the triangle normal */
+ f3_cross(N, E[1], E[0]);
+
+ /* Compute the 3D edge equation */
+ f3_normalize(E[0], f3_cross(E[0], E[0], N));
+ f3_normalize(E[1], f3_cross(E[1], E[1], N));
+ f3_normalize(E[2], f3_cross(E[2], E[2], N));
+ E[0][3] = -f3_dot(E[0], v0.value);
+ E[1][3] = -f3_dot(E[1], v1.value);
+ E[2][3] = -f3_dot(E[2], v2.value);
+
+ /* Compute the distance from current position to the edges */
+ dst[0] = f3_dot(E[0], P) + E[0][3];
+ dst[1] = f3_dot(E[1], P) + E[1][3];
+ dst[2] = f3_dot(E[2], P) + E[2][3];
+
+ /* Retrieve the min and max distance from random walk position to triangle
+ * edges */
+ min_dst = MMIN(MMIN(dst[0], dst[1]), dst[2]);
+ max_dst = MMAX(MMAX(dst[0], dst[1]), dst[2]);
+
+ /* Sort the edges with respect to their distance to the random walk position */
+ FOR_EACH(i, 0, 3) {
+ if(dst[i] == min_dst) {
+ imin = i;
+ } else if(dst[i] == max_dst) {
+ imax = i;
+ } else {
+ imid = i;
+ }
+ }
+ (void)imax;
+
+ /* TODO if the current position is near a vertex, one should move toward the
+ * farthest edge along its normal to avoid too small displacement */
+
+ /* Compute the distance `dst' from the current position to the edges to move
+ * to, along the normal of the edge from which the random walk is the nearest
+ *
+ * +. cos(a) = d / dst => dst = d / cos_a
+ * / `*.
+ * / | `*.
+ * / dst| a /`*.
+ * / | / `*.
+ * / | / d `*.
+ * / |/ `*.
+ * +---------o----------------+ */
+ cos_a1 = f3_dot(E[imin], f3_minus(tmp, E[imid]));
+ cos_a2 = f3_dot(E[imin], f3_minus(tmp, E[imax]));
+ dst[imid] = cos_a1 > 0 ? dst[imid] / cos_a1 : FLT_MAX;
+ dst[imax] = cos_a2 > 0 ? dst[imax] / cos_a2 : FLT_MAX;
+
+ /* Compute the maximum displacement distance into the triangle along the
+ * normal of the edge from which the random walk is the nearest */
+ len = MMIN(dst[imid], dst[imax]);
+ ASSERT(len != FLT_MAX);
+
+ /* Define the displacement distance as the minimum between 10 percent of
+ * delta and len / 2. */
+ len = MMIN(len*0.5f, (float)(delta*0.1));
+ XD(move_pos)(position, E[imin], len);
+}
+#endif
+
#include "sdis_Xd_end.h"
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
@@ -13,6 +13,7 @@
* 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_brdf.h"
#include "sdis_heat_path_boundary_c.h"
#include "sdis_interface_c.h"
#include "sdis_log.h"
@@ -29,27 +30,6 @@
#ifndef SDIS_HEAT_PATH_BOUNDARY_XD_HANDLE_EXTERNAL_NET_FLUX_H
#define SDIS_HEAT_PATH_BOUNDARY_XD_HANDLE_EXTERNAL_NET_FLUX_H
-enum brdf_component {
- BRDF_SPECULAR,
- BRDF_DIFFUSE,
- BRDF_NONE
-};
-
-struct brdf_sample {
- double dir[3];
- double pdf;
- enum brdf_component cpnt;
-};
-#define BRDF_SAMPLE_NULL__ {{0}, 0, BRDF_NONE}
-static const struct brdf_sample BRDF_SAMPLE_NULL = BRDF_SAMPLE_NULL__;
-
-struct brdf {
- double emissivity;
- double specular_fraction;
-};
-#define BRDF_NULL__ {0, 0}
-static const struct brdf BRDF_NULL = BRDF_NULL__;
-
/* Incident diffuse flux is made up of two components. One corresponds to the
* diffuse flux due to the reflection of the source on surfaces. The other is
* the diffuse flux due to the source's radiation scattering at least once in
@@ -63,112 +43,6 @@ struct incident_diffuse_flux {
static const struct incident_diffuse_flux INCIDENT_DIFFUSE_FLUX_NULL =
INCIDENT_DIFFUSE_FLUX_NULL__;
-/* Reflect the V wrt the normal N. By convention V points outward the surface.
- * In fact, this function is a double-precision version of the reflect_3d
- * function. TODO Clean this "repeat" */
-static FINLINE double*
-reflect(double res[3], const double V[3], const double N[3])
-{
- double tmp[3];
- double cos_V_N;
- ASSERT(res && V && N);
- ASSERT(d3_is_normalized(V) && d3_is_normalized(N));
- cos_V_N = d3_dot(V, N);
- d3_muld(tmp, N, 2*cos_V_N);
- d3_sub(res, tmp, V);
- return res;
-}
-
-static void
-sample_brdf
- (const struct brdf* brdf,
- struct ssp_rng* rng,
- const double wi[3], /* Incident direction. Point away from the surface */
- const double N[3], /* Surface normal */
- struct brdf_sample* sample)
-{
- double r = 0; /* Random number */
-
- /* Preconditions */
- ASSERT(brdf && rng && wi && N && sample);
- ASSERT(d3_is_normalized(wi) && d3_is_normalized(N));
- ASSERT(d3_dot(wi, N) > 0);
-
- r = ssp_rng_canonical(rng);
-
- /* Sample the specular part */
- if(r < brdf->specular_fraction) {
- reflect(sample->dir, wi, N);
- sample->pdf = 1;
- sample->cpnt = BRDF_SPECULAR;
-
- /* Sample the diffuse part */
- } else {
- ssp_ran_hemisphere_cos(rng, N, sample->dir, NULL);
- sample->pdf = 1.0/PI;
- sample->cpnt = BRDF_DIFFUSE;
- }
-}
-
-/* Check that the trajectory reaches a valid interface, i.e. that it is on a
- * fluid/solid interface and has reached it from the fluid */
-static res_T
-check_interface
- (const struct sdis_interface* interf,
- const struct sdis_interface_fragment* frag,
- const int verbose) /* Control the verbosity of the function */
-{
- enum sdis_medium_type mdm_frt_type = SDIS_MEDIUM_TYPES_COUNT__;
- enum sdis_medium_type mdm_bck_type = SDIS_MEDIUM_TYPES_COUNT__;
- enum sdis_side fluid_side = SDIS_SIDE_NULL__;
- res_T res = RES_OK;
-
- mdm_frt_type = sdis_medium_get_type(interf->medium_front);
- mdm_bck_type = sdis_medium_get_type(interf->medium_back);
-
- /* Semi-transparent materials are not supported. This means that a solid/solid
- * interface must not be intersected when tracing radiative paths */
- if(mdm_frt_type == SDIS_SOLID && mdm_bck_type == SDIS_SOLID) {
- if(verbose) {
- log_err(interf->dev,
- "Error when sampling the trajectory to calculate the incident diffuse "
- "flux. The trajectory reaches a solid/solid interface, whereas this is "
- "supposed to be impossible (path position: %g, %g, %g).\n",
- SPLIT3(frag->P));
- }
- res = RES_BAD_OP;
- goto error;
- }
-
- /* Find out which side of the interface the fluid is on */
- if(mdm_frt_type == SDIS_FLUID) {
- fluid_side = SDIS_FRONT;
- } else if(mdm_bck_type == SDIS_FLUID) {
- fluid_side = SDIS_BACK;
- } else {
- FATAL("Unreachable code\n");
- }
-
- /* Check that the current position is on the correct side of the interface */
- if(frag->side != fluid_side) {
- if(verbose) {
- log_err(interf->dev,
- "Inconsistent intersection when sampling the trajectory to calculate "
- "the incident diffuse flux. The radiative path reaches an interface on "
- "its solid side, whereas this is supposed to be impossible "
- "(path position: %g, %g, %g).\n",
- SPLIT3(frag->P));
- }
- res = RES_BAD_OP;
- goto error;
- }
-
-exit:
- return res;
-error:
- goto exit;
-}
-
#endif /* SDIS_HEAT_PATH_BOUNDARY_XD_HANDLE_EXTERNAL_NET_FLUX_H */
/*******************************************************************************
@@ -214,36 +88,6 @@ XD(check_handle_external_net_flux_args)
return RES_OK;
}
-static INLINE void
-XD(trace_ray)
- (const struct sdis_scene* scn,
- const double pos[DIM],
- const double dir[3],
- const double distance,
- const struct sXd(hit)* hit_from,
- struct sXd(hit)* hit)
-{
- struct hit_filter_data filter_data = HIT_FILTER_DATA_NULL;
- float ray_org[DIM] = {0};
- float ray_dir[3] = {0};
- float ray_range[2] = {0};
- ASSERT(scn && pos && dir && distance >= 0 && hit_from && hit);
-
- fX_set_dX(ray_org, pos);
- f3_set_d3(ray_dir, dir);
- ray_range[0] = 0;
- ray_range[1] = (float)distance;
- filter_data.XD(hit) = *hit_from;
- filter_data.epsilon = 1.e-4;
-#if DIM == 2
- SXD(scene_view_trace_ray_3d
- (scn->sXd(view), ray_org, ray_dir, ray_range, &filter_data, hit));
-#else
- SXD(scene_view_trace_ray
- (scn->sXd(view), ray_org, ray_dir, ray_range, &filter_data, hit));
-#endif
-}
-
static INLINE double /* [W/m^2/sr] */
XD(direct_contribution)
(const struct sdis_scene* scn,
@@ -266,156 +110,6 @@ XD(direct_contribution)
return sample->radiance_term; /* [W/m^2/sr] */
}
-static INLINE void
-XD(setup_fragment)
- (struct sdis_interface_fragment* frag,
- const double pos[DIM],
- const double dir[DIM], /* Direction _toward_ the hit position */
- const double time, /* Current time */
- const double N[DIM],/* Surface normal */
- const struct sXd(hit)* hit)
-{
- struct sdis_rwalk_vertex vtx = SDIS_RWALK_VERTEX_NULL;
- enum sdis_side side = SDIS_SIDE_NULL__;
- ASSERT(frag && pos && dir && N);
- ASSERT(dX(is_normalized)(N));
-
- /* Setup the interface fragment at the intersection position */
- dX(set)(vtx.P, pos);
- vtx.time = time;
- side = dX(dot)(dir, N) < 0 ? SDIS_FRONT : SDIS_BACK;
- XD(setup_interface_fragment)(frag, &vtx, hit, side);
-}
-
-static INLINE res_T
-XD(find_next_fragment)
- (const struct sdis_scene* scn,
- const double in_pos[3],
- const double in_dir[3], /* Always in 3D */
- const struct sXd(hit)* in_hit,
- const double time,
- struct sXd(hit)* out_hit,
- struct sdis_interface** out_interf,
- struct sdis_interface_fragment* out_frag)
-{
- const int NATTEMPTS_MAX = 10;
- int nattempts = 0;
-
- /* Stardis */
- struct sdis_interface_fragment frag = SDIS_INTERFACE_FRAGMENT_NULL;
- struct sdis_interface* interf = NULL;
-
- struct sXd(hit) hit = SXD_HIT_NULL;
- double rt_pos[DIM] = {0};
- res_T res = RES_OK;
-
- ASSERT(scn && in_pos && in_dir && in_hit && !SXD_HIT_NONE(in_hit));
- ASSERT(out_hit && out_interf && out_frag);
-
- dX(set)(rt_pos, in_pos);
-
- do {
- struct sdis_rwalk_vertex vtx = SDIS_RWALK_VERTEX_NULL;
- struct sdis_medium* solid = NULL;
- double pos[3] = {0};
- double vec[3] = {0};
- double N[3] = {0};
- double delta = 0;
-
- /* Reset result code. It may have been modified during a previous attempt */
- res = RES_OK;
-
- /* Find the following surface along the direction of propagation */
- XD(trace_ray)(scn, rt_pos, in_dir, INF, in_hit, &hit);
- if(SXD_HIT_NONE(&hit)) break;
-
- /* Retrieve the current position and normal */
- dX(add)(pos, rt_pos, dX(muld)(vec, in_dir, hit.distance));
- dX_set_fX(N, hit.normal);
- dX(normalize(N, N));
-
- /* Retrieve the current interface properties */
- interf = scene_get_interface(scn, hit.prim.prim_id);
- XD(setup_fragment)(&frag, pos, in_dir, time, N, &hit);
-
- /* Check that the path reaches a valid interface.
- * An invalid fragment may mean that the ray position is in a corner and the
- * traced ray has missed the surface of that corner. To correct this, the
- * ray position is moved slightly away from the corner before a ray is drawn
- * in the same direction. This fallback solution is executed a number of
- * times, after which, if the fragment is still invalid, it is considered
- * that the numerical error cannot be mitigated. */
- res = check_interface(interf, &frag, nattempts == NATTEMPTS_MAX);
- if(res != RES_OK && nattempts == NATTEMPTS_MAX) goto error;
- ++nattempts;
-
- if(res != RES_OK) { /* Mitigate numerical error (see above) */
- if(sdis_medium_get_type(interf->medium_front) == SDIS_SOLID) {
- solid = interf->medium_front;
- } else {
- ASSERT(sdis_medium_get_type(interf->medium_back) == SDIS_SOLID);
- solid = interf->medium_back;
- }
-
- /* Retrieves the delta of the solid that surrounds the boundary, as it is
- * actually the only numerical parameter that says something about the
- * system. */
- vtx.P[0] = pos[0];
- vtx.P[1] = pos[1];
- vtx.P[2] = pos[2];
- vtx.time = time;
- delta = solid_get_delta(solid, &vtx);
-
- XD(move_away_primitive_boundaries)(in_hit, delta, rt_pos);
- }
- } while(res != RES_OK);
-
-exit:
- *out_hit = hit;
- *out_interf = interf;
- *out_frag = frag;
- return res;
-error:
- goto exit;
-}
-
-static INLINE res_T
-XD(setup_brdf)
- (struct sdis_device* dev,
- const struct sdis_source* src,
- struct brdf* brdf,
- const struct sdis_interface* interf,
- const struct sdis_interface_fragment* frag)
-{
- double epsilon = 0;
- double alpha = 0;
- unsigned src_id = 0;
- res_T res = RES_OK;
- ASSERT(brdf && frag);
- ASSERT((frag->side == SDIS_FRONT
- && sdis_medium_get_type(interf->medium_front) == SDIS_FLUID)
- || sdis_medium_get_type(interf->medium_back) == SDIS_FLUID);
-
- src_id = sdis_source_get_id(src);
-
- epsilon = interface_side_get_emissivity(interf, src_id, frag);
- res = interface_side_check_emissivity(dev, epsilon, frag->P, frag->time);
- if(res != RES_OK) goto error;
-
- alpha = interface_side_get_specular_fraction(interf, src_id, frag);
- res = interface_side_check_specular_fraction(dev, alpha, frag->P, frag->time);
- if(res != RES_OK) goto error;
-
- brdf->emissivity = epsilon;
- brdf->specular_fraction = alpha;
-
-exit:
- return res;
-error:
- *brdf = BRDF_NULL;
- goto exit;
-}
-
static res_T
XD(compute_incident_diffuse_flux)
(const struct sdis_scene* scn,
@@ -454,7 +148,8 @@ XD(compute_incident_diffuse_flux)
/* BRDF */
struct brdf brdf = BRDF_NULL;
- struct brdf_sample brdf_sample = BRDF_SAMPLE_NULL;
+ struct brdf_sample bounce = BRDF_SAMPLE_NULL;
+ struct brdf_setup_args brdf_setup_args = BRDF_SETUP_ARGS_NULL;
/* Miscellaneous */
double L = 0; /* incident flux to bounce position */
@@ -482,7 +177,13 @@ XD(compute_incident_diffuse_flux)
}
d3_set(pos, frag.P);
- XD(setup_brdf)(scn->dev, scn->source, &brdf, interf, &frag);
+
+ /* Retrieve BRDF at current interface position */
+ brdf_setup_args.interf = interf;
+ brdf_setup_args.frag = &frag;
+ brdf_setup_args.source_id = sdis_source_get_id(scn->source);
+ res = brdf_setup(scn->dev, &brdf_setup_args, &brdf);
+ if(res != RES_OK) goto error;
/* Check if path is absorbed */
if(ssp_rng_canonical(rng) < brdf.emissivity) break;
@@ -493,12 +194,12 @@ XD(compute_incident_diffuse_flux)
case SDIS_BACK: dX(minus)(N, frag.Ng); break;
default: FATAL("Unreachable code\n");
}
- sample_brdf(&brdf, rng, wi, N, &brdf_sample);
- d3_set(dir, brdf_sample.dir); /* Always in 3D */
+ brdf_sample(&brdf, rng, wi, N, &bounce);
+ d3_set(dir, bounce.dir); /* Always in 3D */
/* Calculate the direct contribution if the rebound is specular */
- if(brdf_sample.cpnt == BRDF_SPECULAR) {
- res = source_trace_to(scn->source, pos, brdf_sample.dir, time, &src_sample);
+ if(bounce.cpnt == BRDF_SPECULAR) {
+ res = source_trace_to(scn->source, pos, bounce.dir, time, &src_sample);
if(res != RES_OK) goto error;
if(!SOURCE_SAMPLE_NONE(&src_sample)) {
@@ -509,7 +210,7 @@ XD(compute_incident_diffuse_flux)
/* Calculate the direct contribution of the rebound is diffuse */
} else {
double cos_theta = 0;
- ASSERT(brdf_sample.cpnt == BRDF_DIFFUSE);
+ ASSERT(bounce.cpnt == BRDF_DIFFUSE);
/* Sample an external source to handle its direct contribution at the
* bounce position */
diff --git a/src/sdis_heat_path_boundary_c.h b/src/sdis_heat_path_boundary_c.h
@@ -222,6 +222,20 @@ query_medium_temperature_from_boundary_3d
struct rwalk* rwalk,
struct temperature* T);
+/* Move the submitted position away from the primitive boundaries to avoid
+ * numerical issues leading to inconsistent random walks. */
+extern LOCAL_SYM void
+move_away_primitive_boundaries_2d
+ (const struct s2d_hit* hit,
+ const double delta,
+ double position[2]); /* Position to move */
+
+extern LOCAL_SYM void
+move_away_primitive_boundaries_3d
+ (const struct s3d_hit* hit,
+ const double delta,
+ double position[3]); /* Position to move */
+
/*******************************************************************************
* Boundary sub-paths
******************************************************************************/
diff --git a/src/sdis_heat_path_radiative_Xd.h b/src/sdis_heat_path_radiative_Xd.h
@@ -13,9 +13,11 @@
* 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_brdf.h"
#include "sdis_device_c.h"
#include "sdis_green.h"
#include "sdis_heat_path.h"
+#include "sdis_heat_path_boundary_c.h"
#include "sdis_interface_c.h"
#include "sdis_log.h"
#include "sdis_medium_c.h"
@@ -28,7 +30,7 @@
#include "sdis_Xd_begin.h"
/*******************************************************************************
- * Non generic local functions
+ * Non generic helper functions
******************************************************************************/
#ifndef SDIS_HEAT_PATH_RADIATIVE_XD_H
#define SDIS_HEAT_PATH_RADIATIVE_XD_H
@@ -39,7 +41,7 @@ set_limit_radiative_temperature
struct rwalk_context* ctx,
struct rwalk* rwalk,
/* Direction along which the random walk reached the radiative environment */
- const float dir[3],
+ const double dir[3],
const int branch_id,
struct temperature* T)
{
@@ -52,7 +54,7 @@ set_limit_radiative_temperature
ASSERT(SXD_HIT_NONE(&rwalk->XD(hit)));
rwalk->hit_side = SDIS_SIDE_NULL__;
- d3_set_f3(rwalk->dir, dir);
+ d3_set(rwalk->dir, dir);
d3_normalize(rwalk->dir, rwalk->dir);
d3_set(ray.dir, rwalk->dir);
ray.time = rwalk->vtx.time;
@@ -86,7 +88,7 @@ set_limit_radiative_temperature
* distance of 0.1 meters from the surface along the direction reaching the
* radiative environment */
if(ctx->heat_path) {
- const float empirical_dst = 0.1f * (float)scn->fp_to_meter;
+ const double empirical_dst = 0.1 * (float)scn->fp_to_meter;
struct sdis_rwalk_vertex vtx = SDIS_RWALK_VERTEX_NULL;
vtx = rwalk->vtx;
@@ -104,9 +106,91 @@ error:
goto exit;
}
+/* Check that the trajectory reaches a valid interface, i.e. that it is on a
+ * fluid/solid interface and has reached it from the fluid */
+static res_T
+check_interface
+ (const struct sdis_interface* interf,
+ const struct sdis_interface_fragment* frag,
+ const int verbose) /* Control the verbosity of the function */
+{
+ enum sdis_medium_type mdm_frt_type = SDIS_MEDIUM_TYPES_COUNT__;
+ enum sdis_medium_type mdm_bck_type = SDIS_MEDIUM_TYPES_COUNT__;
+ enum sdis_side fluid_side = SDIS_SIDE_NULL__;
+ res_T res = RES_OK;
+
+ mdm_frt_type = sdis_medium_get_type(interf->medium_front);
+ mdm_bck_type = sdis_medium_get_type(interf->medium_back);
+
+ /* Semi-transparent materials are not supported. This means that a solid/solid
+ * interface must not be intersected when tracing radiative paths */
+ if(mdm_frt_type == SDIS_SOLID && mdm_bck_type == SDIS_SOLID) {
+ if(verbose) {
+ log_err(interf->dev,
+ "Error when sampling the radiatve path. The trajectory reaches a "
+ "solid/solid interface, whereas this is supposed to be impossible "
+ "(path position: %g, %g, %g).\n",
+ SPLIT3(frag->P));
+ }
+ res = RES_BAD_OP;
+ goto error;
+ }
+
+ /* Find out which side of the interface the fluid is on */
+ if(mdm_frt_type == SDIS_FLUID) {
+ fluid_side = SDIS_FRONT;
+ } else if(mdm_bck_type == SDIS_FLUID) {
+ fluid_side = SDIS_BACK;
+ } else {
+ FATAL("Unreachable code\n");
+ }
+
+ /* Check that the current position is on the correct side of the interface */
+ if(frag->side != fluid_side) {
+ if(verbose) {
+ log_err(interf->dev,
+ "Inconsistent intersection when sampling the radiative path. "
+ "The path reaches an interface on its solid side, whereas this is "
+ "supposed to be impossible (path position: %g, %g, %g).\n",
+ SPLIT3(frag->P));
+ }
+ res = RES_BAD_OP;
+ goto error;
+ }
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
#endif /* SDIS_HEAT_PATH_RADIATIVE_XD_H */
/*******************************************************************************
+ * Generic helper functions
+ ******************************************************************************/
+static INLINE void
+XD(setup_fragment)
+ (struct sdis_interface_fragment* frag,
+ const double pos[DIM],
+ const double dir[DIM], /* Direction _toward_ the hit position */
+ const double time, /* Current time */
+ const double N[DIM],/* Surface normal */
+ const struct sXd(hit)* hit)
+{
+ struct sdis_rwalk_vertex vtx = SDIS_RWALK_VERTEX_NULL;
+ enum sdis_side side = SDIS_SIDE_NULL__;
+ ASSERT(frag && pos && dir && N);
+ ASSERT(dX(is_normalized)(N));
+
+ /* Setup the interface fragment at the intersection position */
+ dX(set)(vtx.P, pos);
+ vtx.time = time;
+ side = dX(dot)(dir, N) < 0 ? SDIS_FRONT : SDIS_BACK;
+ XD(setup_interface_fragment)(frag, &vtx, hit, side);
+}
+
+/*******************************************************************************
* Local functions
******************************************************************************/
res_T
@@ -120,46 +204,40 @@ XD(trace_radiative_path)
{
/* 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. */
- float N[3] = {0, 0, 0};
- float dir[3] = {0, 0, 0};
+ double N[3] = {0,0,0};
+ double dir[3] = {0,0,0};
+ double pos[3] = {0,0,0};
int branch_id;
res_T res = RES_OK;
ASSERT(scn && ray_dir && ctx && rwalk && rng && T);
- f3_set(dir, ray_dir);
+ d3_set_f3(dir, ray_dir);
+ d3_normalize(dir, dir);
/* (int)ctx->nbranchings < 0 <=> Beginning of the realisation */
branch_id = MMAX((int)ctx->nbranchings, 0);
/* 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;
+ /* BRDF */
+ struct brdf brdf = BRDF_NULL;
+ struct brdf_sample bounce = BRDF_SAMPLE_NULL;
+ struct brdf_setup_args brdf_setup_args = BRDF_SETUP_ARGS_NULL;
+
+ /* Miscellaneous */
struct sdis_interface_fragment frag = SDIS_INTERFACE_FRAGMENT_NULL;
- unsigned enc_ids[2] = {ENCLOSURE_ID_NULL, ENCLOSURE_ID_NULL};
- unsigned chk_enc_id = ENCLOSURE_ID_NULL;
- double alpha;
- double epsilon;
- double r;
- float pos[DIM];
- const float range[2] = { 0, FLT_MAX };
-
- fX_set_dX(pos, rwalk->vtx.P);
-
- /* Trace the radiative ray */
- filter_data.XD(hit) = rwalk->XD(hit);
- filter_data.epsilon = 1.e-6;
- filter_data.scn = scn; /* Enable the filtering wrt the enclosure id */
- filter_data.enc_id = rwalk->enc_id;
-#if (SDIS_XD_DIMENSION == 2)
- SXD(scene_view_trace_ray_3d
- (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->XD(hit)));
-#endif
+ struct sdis_interface* interf = NULL;
+ struct sdis_medium* chk_mdm = NULL;
+ double wi[3] = {0,0,0};
+ float dirf[3] = {0,0,0};
+
+ d3_set(pos, rwalk->vtx.P);
+ d3_minus(wi, dir);
+
+ res = XD(find_next_fragment)(scn, pos, dir, &rwalk->XD(hit),
+ rwalk->vtx.time, &rwalk->XD(hit), &interf, &frag);
+ if(res != RES_OK) goto error;
/* The path reaches the radiative environment */
if(SXD_HIT_NONE(&rwalk->XD(hit))) {
@@ -170,60 +248,9 @@ XD(trace_radiative_path)
break; /* Stop the radiative path */
}
- /* Define the hit side */
- 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->XD(hit).distance);
-
- /* Register the random walk vertex against the heat path */
- res = register_heat_vertex(ctx->heat_path, &rwalk->vtx, T->value,
- SDIS_HEAT_VERTEX_RADIATIVE, branch_id);
- if(res != RES_OK) goto error;
-
- /* Fetch the new interface and setup the hit fragment */
- 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);
- if(epsilon > 1 || epsilon < 0) {
- log_err(scn->dev,
- "%s: invalid overall emissivity `%g' at position `%g %g %g'.\n",
- FUNC_NAME, epsilon, SPLIT3(rwalk->vtx.P));
- res = RES_BAD_OP;
- goto error;
- }
-
- /* Switch in boundary temperature ? */
- r = ssp_rng_canonical(rng);
- if(r < epsilon) {
- T->func = XD(boundary_path);
- 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->XD(hit).normal);
- if(rwalk->hit_side == SDIS_BACK) fX(minus)(N, N);
-
- /* Check that the radiative path is still within the same enclosure. Note
- * that this may not be the case, even if the filtering of intersections
- * relative to the current enclosure is enabled. This filtering is only
- * performed for intersections on a boundary between primitives. As a
- * consequence, a threshold effect on how "intersections on a boundary" are
- * detected could lead to this situation */
- 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;
- }
+ XD(move_pos)(rwalk->vtx.P, fX_set_dX(dirf, dir), rwalk->XD(hit).distance);
+ rwalk->hit_side = frag.side;
/* Verify that the intersection, although in the same enclosure, touches the
* interface of a fluid. We verify this by interface, since a radiative path
@@ -236,7 +263,8 @@ XD(trace_radiative_path)
* 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;
+ ? interf->medium_front
+ : interf->medium_back;
if(sdis_medium_get_type(chk_mdm) == SDIS_SOLID) {
log_err(scn->dev,
"%s: a radiative path cannot evolve in a solid -- pos=(%g, %g, %g)\n",
@@ -245,13 +273,34 @@ XD(trace_radiative_path)
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 */
- reflect_3d(dir, f3_minus(dir, dir), N);
- } else { /* Sample diffuse part */
- ssp_ran_hemisphere_cos_float(rng, N, dir, NULL);
+ /* Register the random walk vertex against the heat path */
+ res = register_heat_vertex(ctx->heat_path, &rwalk->vtx, T->value,
+ SDIS_HEAT_VERTEX_RADIATIVE, branch_id);
+ if(res != RES_OK) goto error;
+
+ /* Retrieve BRDF at current interface position */
+ brdf_setup_args.interf = interf;
+ brdf_setup_args.frag = &frag;
+ brdf_setup_args.source_id = SDIS_INTERN_SOURCE_ID;
+ res = brdf_setup(scn->dev, &brdf_setup_args, &brdf);
+ if(res != RES_OK) goto error;
+
+ /* Switch in boundary temperature? */
+ if(ssp_rng_canonical(rng) < brdf.emissivity) {
+ T->func = XD(boundary_path);
+ 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 */
+ switch(frag.side) {
+ case SDIS_FRONT: d3_set(N, frag.Ng); break;
+ case SDIS_BACK: d3_minus(N, frag.Ng); break;
+ default: FATAL("Unreachable code\n"); break;
+ }
+ brdf_sample(&brdf, rng, wi, N, &bounce);
+ d3_set(dir, bounce.dir); /* Always in 3D */
}
exit:
@@ -297,4 +346,126 @@ error:
goto exit;
}
+void
+XD(trace_ray)
+ (const struct sdis_scene* scn,
+ const double pos[DIM],
+ const double dir[3],
+ const double distance,
+ const struct sXd(hit)* hit_from,
+ struct sXd(hit)* hit)
+{
+ struct hit_filter_data filter_data = HIT_FILTER_DATA_NULL;
+ float ray_org[DIM] = {0};
+ float ray_dir[3] = {0};
+ float ray_range[2] = {0};
+ ASSERT(scn && pos && dir && distance >= 0 && hit_from && hit);
+
+ fX_set_dX(ray_org, pos);
+ f3_set_d3(ray_dir, dir);
+ ray_range[0] = 0;
+ ray_range[1] = (float)distance;
+ filter_data.XD(hit) = *hit_from;
+ filter_data.epsilon = 1.e-4;
+#if DIM == 2
+ SXD(scene_view_trace_ray_3d
+ (scn->sXd(view), ray_org, ray_dir, ray_range, &filter_data, hit));
+#else
+ SXD(scene_view_trace_ray
+ (scn->sXd(view), ray_org, ray_dir, ray_range, &filter_data, hit));
+#endif
+}
+
+res_T
+XD(find_next_fragment)
+ (const struct sdis_scene* scn,
+ const double in_pos[DIM],
+ const double in_dir[3], /* Always in 3D */
+ const struct sXd(hit)* in_hit,
+ const double time,
+ struct sXd(hit)* out_hit,
+ struct sdis_interface** out_interf,
+ struct sdis_interface_fragment* out_frag)
+{
+ const int NATTEMPTS_MAX = 10;
+ int nattempts = 0;
+
+ /* Stardis */
+ struct sdis_interface_fragment frag = SDIS_INTERFACE_FRAGMENT_NULL;
+ struct sdis_interface* interf = NULL;
+
+ struct sXd(hit) hit = SXD_HIT_NULL;
+ double rt_pos[DIM] = {0};
+ res_T res = RES_OK;
+
+ ASSERT(scn && in_pos && in_dir && in_hit);
+ ASSERT(out_hit && out_interf && out_frag);
+
+ dX(set)(rt_pos, in_pos);
+
+ do {
+ struct sdis_rwalk_vertex vtx = SDIS_RWALK_VERTEX_NULL;
+ struct sdis_medium* solid = NULL;
+ double pos[3] = {0};
+ double vec[3] = {0};
+ double N[3] = {0};
+ double delta = 0;
+
+ /* Reset result code. It may have been modified during a previous attempt */
+ res = RES_OK;
+
+ /* Find the following surface along the direction of propagation */
+ XD(trace_ray)(scn, rt_pos, in_dir, INF, in_hit, &hit);
+ if(SXD_HIT_NONE(&hit)) break;
+
+ /* Retrieve the current position and normal */
+ dX(add)(pos, rt_pos, dX(muld)(vec, in_dir, hit.distance));
+ dX_set_fX(N, hit.normal);
+ dX(normalize(N, N));
+
+ /* Retrieve the current interface properties */
+ interf = scene_get_interface(scn, hit.prim.prim_id);
+ XD(setup_fragment)(&frag, pos, in_dir, time, N, &hit);
+
+ /* Check that the path reaches a valid interface.
+ * An invalid fragment may mean that the ray position is in a corner and the
+ * traced ray has missed the surface of that corner. To correct this, the
+ * ray position is moved slightly away from the corner before a ray is drawn
+ * in the same direction. This fallback solution is executed a number of
+ * times, after which, if the fragment is still invalid, it is considered
+ * that the numerical error cannot be mitigated. */
+ res = check_interface(interf, &frag, nattempts == NATTEMPTS_MAX);
+ if(res != RES_OK && nattempts == NATTEMPTS_MAX) goto error;
+ ++nattempts;
+
+ if(res != RES_OK) { /* Mitigate numerical error (see above) */
+ if(sdis_medium_get_type(interf->medium_front) == SDIS_SOLID) {
+ solid = interf->medium_front;
+ } else {
+ ASSERT(sdis_medium_get_type(interf->medium_back) == SDIS_SOLID);
+ solid = interf->medium_back;
+ }
+
+ /* Retrieves the delta of the solid that surrounds the boundary, as it is
+ * actually the only numerical parameter that says something about the
+ * system. */
+ vtx.P[0] = pos[0];
+ vtx.P[1] = pos[1];
+ vtx.P[2] = pos[2];
+ vtx.time = time;
+ delta = solid_get_delta(solid, &vtx);
+
+ XD(move_away_primitive_boundaries)(in_hit, delta, rt_pos);
+ }
+ } while(res != RES_OK);
+
+exit:
+ *out_hit = hit;
+ *out_interf = interf;
+ *out_frag = frag;
+ return res;
+error:
+ goto exit;
+}
+
#include "sdis_Xd_end.h"
