commit 922814aeb8a1977f27ee6897a5160d2c957e71da
parent b07f9bf0e5b655f3458e2ab4fb0c81a2bec31df3
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Wed, 17 Jan 2024 11:35:55 +0100
Merge branch 'feature_external_flux' into develop
Diffstat:
21 files changed, 2357 insertions(+), 61 deletions(-)
diff --git a/Makefile b/Makefile
@@ -46,6 +46,7 @@ SRC =\
src/sdis_scene.c\
src/sdis_solve.c\
src/sdis_solve_camera.c\
+ src/sdis_source.c\
src/sdis_tile.c\
$($(DISTRIB_PARALLELISM)_SRC)
OBJ = $(SRC:.c=.o)
@@ -176,6 +177,7 @@ TEST_SRC =\
src/test_sdis_convection.c\
src/test_sdis_convection_non_uniform.c\
src/test_sdis_data.c\
+ src/test_sdis_draw_external_flux.c\
src/test_sdis_enclosure_limit_conditions.c\
src/test_sdis_flux.c\
src/test_sdis_flux2.c\
@@ -190,6 +192,7 @@ TEST_SRC =\
src/test_sdis_solve_probe_2d.c\
src/test_sdis_solve_probe2_2d.c\
src/test_sdis_solve_probe3_2d.c\
+ src/test_sdis_source.c\
src/test_sdis_transcient.c\
src/test_sdis_unstationary_atm.c\
src/test_sdis_volumic_power.c\
@@ -202,6 +205,7 @@ TEST_SRC_MPI =\
src/test_sdis.c\
src/test_sdis_compute_power.c\
src/test_sdis_device.c\
+ src/test_sdis_external_flux.c\
src/test_sdis_solve_camera.c\
src/test_sdis_solve_medium.c\
src/test_sdis_solve_medium_2d.c\
@@ -287,6 +291,7 @@ src/test_sdis_solve_probe.d \
src/test_sdis_solve_probe_2d.d \
src/test_sdis_solve_probe2_2d.d \
src/test_sdis_solve_probe3_2d \
+src/test_sdis_source.d \
src/test_sdis_transcient.d \
src/test_sdis_unstationary_atm.d \
src/test_sdis_utils.d \
@@ -317,6 +322,7 @@ src/test_sdis_solve_probe.o \
src/test_sdis_solve_probe_2d.o \
src/test_sdis_solve_probe2_2d.o \
src/test_sdis_solve_probe3_2d.o \
+src/test_sdis_source.o \
src/test_sdis_transcient.o \
src/test_sdis_unstationary_atm.o \
src/test_sdis_utils.o \
@@ -347,6 +353,7 @@ test_sdis_solve_probe \
test_sdis_solve_probe_2d \
test_sdis_solve_probe2_2d \
test_sdis_solve_probe3_2d \
+test_sdis_source \
test_sdis_transcient \
test_sdis_unstationary_atm \
test_sdis_volumic_power \
@@ -360,16 +367,19 @@ test_sdis_volumic_power4 \
################################################################################
# Tests based on Star-3DUT
################################################################################
+src/test_sdis_draw_external_flux.d \
src/test_sdis_solid_random_walk_robustness.d \
src/test_sdis_solve_probe3.d \
: config.mk sdis-local.pc
@$(CC) $(TEST_CFLAGS) $(S3DUT_CFLAGS) -MM -MT "$(@:.d=.o) $@" $(@:.d=.c) -MF $@
+src/test_sdis_draw_external_flux.o \
src/test_sdis_solid_random_walk_robustness.o \
src/test_sdis_solve_probe3.o \
: config.mk sdis-local.pc
$(CC) $(TEST_CFLAGS) $(S3DUT_CFLAGS) -c $(@:.o=.c) -o $@
+test_sdis_draw_external_flux \
test_sdis_solid_random_walk_robustness \
test_sdis_solve_probe3 \
: config.mk sdis-local.pc $(LIBNAME) src/test_sdis_utils.o
@@ -395,18 +405,21 @@ test_sdis_scene \
################################################################################
src/test_sdis.d \
src/test_sdis_device.d \
+src/test_sdis_external_flux.d \
src/test_sdis_solve_medium_2d.d \
: config.mk sdis-local.pc
@$(CC) $(TEST_CFLAGS_MPI) -MM -MT "$(@:.d=.o) $@" $(@:.d=.c) -MF $@
src/test_sdis.o \
src/test_sdis_device.o \
+src/test_sdis_external_flux.o \
src/test_sdis_solve_medium_2d.o \
: config.mk sdis-local.pc
$(CC) $(TEST_CFLAGS_MPI) -c $(@:.o=.c) -o $@
test_sdis \
test_sdis_device \
+test_sdis_external_flux \
test_sdis_solve_medium_2d \
: config.mk sdis-local.pc $(LIBNAME) src/test_sdis_utils.o
$(CC) $(TEST_CFLAGS_MPI) -o $@ src/$@.o $(TEST_LIBS_MPI)
diff --git a/src/sdis.h b/src/sdis.h
@@ -69,6 +69,7 @@ struct sdis_green_function;
struct sdis_interface;
struct sdis_medium;
struct sdis_scene;
+struct sdis_source;
/* Forward declaration of non ref counted types */
struct sdis_green_path;
@@ -113,25 +114,6 @@ struct sdis_interface_fragment {
static const struct sdis_interface_fragment SDIS_INTERFACE_FRAGMENT_NULL =
SDIS_INTERFACE_FRAGMENT_NULL__;
-/*******************************************************************************
- * Estimation data types
- ******************************************************************************/
-enum sdis_estimator_type {
- SDIS_ESTIMATOR_TEMPERATURE, /* In Kelvin */
- SDIS_ESTIMATOR_FLUX, /* In Watt/m^2 */
- SDIS_ESTIMATOR_POWER, /* In Watt */
- SDIS_ESTIMATOR_TYPES_COUNT__
-};
-
-/* Monte-Carlo estimation */
-struct sdis_mc {
- double E; /* Expected value */
- double V; /* Variance */
- double SE; /* Standard error */
-};
-#define SDIS_MC_NULL__ {0, 0, 0}
-static const struct sdis_mc SDIS_MC_NULL = SDIS_MC_NULL__;
-
/* Input arguments of the sdis_device_create function */
struct sdis_device_create_args {
struct logger* logger; /* NULL <=> default logger */
@@ -158,6 +140,44 @@ struct sdis_info {
#define SDIS_INFO_NULL__ {0}
static const struct sdis_info SDIS_INFO_NULL = SDIS_INFO_NULL__;
+/* Type of functor used to retrieve the source's position relative to time. */
+typedef void
+(*sdis_get_position_T)
+ (const double time,
+ double pos[3],
+ struct sdis_data* data);
+
+/* Input arguments of the sdis_spherical_source_create function */
+struct sdis_spherical_source_create_args {
+ sdis_get_position_T position; /* [m] */
+ struct sdis_data* data; /* Data sent to the position functor */
+ double radius; /* [m] */
+ double power; /* Total power [W] */
+};
+#define SDIS_SPHERICAL_SOURCE_CREATE_ARGS_NULL__ {NULL, NULL, 0, 0}
+static const struct sdis_spherical_source_create_args
+SDIS_SPHERICAL_SOURCE_CREATE_ARGS_NULL =
+ SDIS_SPHERICAL_SOURCE_CREATE_ARGS_NULL__;
+
+/*******************************************************************************
+ * Estimation data types
+ ******************************************************************************/
+enum sdis_estimator_type {
+ SDIS_ESTIMATOR_TEMPERATURE, /* In Kelvin */
+ SDIS_ESTIMATOR_FLUX, /* In Watt/m^2 */
+ SDIS_ESTIMATOR_POWER, /* In Watt */
+ SDIS_ESTIMATOR_TYPES_COUNT__
+};
+
+/* Monte-Carlo estimation */
+struct sdis_mc {
+ double E; /* Expected value */
+ double V; /* Variance */
+ double SE; /* Standard error */
+};
+#define SDIS_MC_NULL__ {0, 0, 0}
+static const struct sdis_mc SDIS_MC_NULL = SDIS_MC_NULL__;
+
/*******************************************************************************
* Data type used to describe physical properties
******************************************************************************/
@@ -171,15 +191,15 @@ enum sdis_medium_type {
* medium. */
typedef double
(*sdis_medium_getter_T)
- (const struct sdis_rwalk_vertex* vert,
- struct sdis_data* data);
+ (const struct sdis_rwalk_vertex* vert, /* Medium position */
+ struct sdis_data* data); /* User data */
/* Functor type used to retrieve the spatio temporal physical properties of an
* interface. */
typedef double
(*sdis_interface_getter_T)
- (const struct sdis_interface_fragment* frag,
- struct sdis_data* data);
+ (const struct sdis_interface_fragment* frag, /* Interface position */
+ struct sdis_data* data); /* User data */
/* Define the physical properties of a solid */
struct sdis_solid_shader {
@@ -198,6 +218,7 @@ struct sdis_solid_shader {
* unknown for the submitted random walk vertex.
* This getter is always called at time >= t0 (see below). */
sdis_medium_getter_T temperature;
+
/* The time until the initial condition is maintained for this solid;
* can neither be negative nor infinity, default is 0. */
double t0;
@@ -238,8 +259,12 @@ struct sdis_interface_side_shader {
/* Reference temperature used in Picard 1 */
sdis_interface_getter_T reference_temperature;
+
+ /* Define whether external sources interact with the interface, i.e. whether
+ * external fluxes should be processed or not */
+ int handle_external_flux;
};
-#define SDIS_INTERFACE_SIDE_SHADER_NULL__ { NULL, NULL, NULL, NULL, NULL }
+#define SDIS_INTERFACE_SIDE_SHADER_NULL__ { NULL, NULL, NULL, NULL, NULL, 1 }
static const struct sdis_interface_side_shader SDIS_INTERFACE_SIDE_SHADER_NULL =
SDIS_INTERFACE_SIDE_SHADER_NULL__;
@@ -414,6 +439,10 @@ struct sdis_scene_create_args {
/* Min/max temperature used to linearise the radiative temperature */
double t_range[2];
+
+ /* External source. Can be NULL <=> no external flux will be calculated on
+ * scene interfaces */
+ struct sdis_source* source;
};
#define SDIS_SCENE_CREATE_ARGS_DEFAULT__ { \
@@ -425,7 +454,8 @@ struct sdis_scene_create_args {
0, /* #vertices */ \
1.0, /* #Floating point to meter scale factor */ \
SDIS_AMBIENT_RADIATIVE_TEMPERATURE_NULL__,/* Ambient radiative temperature */\
- {0.0, -1.0} /* Temperature range */ \
+ {0.0, -1.0}, /* Temperature range */ \
+ NULL /* source */ \
}
static const struct sdis_scene_create_args SDIS_SCENE_CREATE_ARGS_DEFAULT =
SDIS_SCENE_CREATE_ARGS_DEFAULT__;
@@ -931,6 +961,29 @@ sdis_interface_get_id
(const struct sdis_interface* interf);
/*******************************************************************************
+ * External source API. When a scene has external sources, an external flux
+ * (in both its direct and diffuse parts) is imposed on the interfaces.
+ ******************************************************************************/
+SDIS_API res_T
+sdis_spherical_source_create
+ (struct sdis_device* dev,
+ struct sdis_spherical_source_create_args* args,
+ struct sdis_source** source);
+
+SDIS_API res_T
+sdis_source_ref_get
+ (struct sdis_source* source);
+
+SDIS_API res_T
+sdis_source_ref_put
+ (struct sdis_source* source);
+
+SDIS_API res_T
+sdis_source_get_power
+ (const struct sdis_source* src,
+ double* power);/* [W] */
+
+/*******************************************************************************
* A scene is a collection of primitives. Each primitive is the geometric
* support of the interface between 2 media.
******************************************************************************/
@@ -1117,6 +1170,11 @@ sdis_scene_get_device
(struct sdis_scene* scn,
struct sdis_device** device);
+SDIS_API res_T
+sdis_scene_get_source
+ (struct sdis_scene* scn,
+ struct sdis_source** src); /* The returned pointer can be NULL <=> no source */
+
/*******************************************************************************
* An estimator stores the state of a simulation
******************************************************************************/
@@ -1316,6 +1374,14 @@ sdis_green_path_for_each_flux_term
sdis_process_interface_flux_term_T func,
void* context);
+/* Return the external flux term, i.e. the relative net flux along the path from
+ * the external source. Multiply it by the power of the source to obtain its
+ * contribution to the path. */
+SDIS_API res_T
+sdis_green_path_get_external_flux_term
+ (struct sdis_green_path* path,
+ double* external_flux_term); /* [W/m^2] */
+
/*******************************************************************************
* Heat path API
******************************************************************************/
diff --git a/src/sdis_Xd_begin.h b/src/sdis_Xd_begin.h
@@ -34,11 +34,10 @@ struct rwalk_context {
double That2; /* That^2 */
double That3; /* That^3 */
- /* Maximum branchings i.e. the maximum number of times
- * XD(compute_temperature) can be called. It controls the number of
- * ramifications of the heat path and currently is correlated to the Picard
- * order used to estimate the radiative temperature. max_branchings ==
- * picard_order-1 */
+ /* Maximum branchings i.e. the maximum number of times XD(sample_coupled_path)
+ * can be called. It controls the number of ramifications of the heat path and
+ * currently is correlated to the Picard order used to estimate the radiative
+ * temperature. max_branchings == picard_order-1 */
size_t max_branchings;
/* Number of heat path branchings */
diff --git a/src/sdis_green.c b/src/sdis_green.c
@@ -21,6 +21,7 @@
#include "sdis_medium_c.h"
#include "sdis_misc.h"
#include "sdis_scene_c.h"
+#include "sdis_source_c.h"
#include <star/ssp.h>
@@ -83,6 +84,7 @@ flux_term_init(struct mem_allocator* allocator, struct flux_term* term)
struct green_path {
double elapsed_time;
+ double external_flux_term; /* [W/m^2] */
struct darray_flux_term flux_terms; /* List of flux terms */
struct darray_power_term power_terms; /* List of volumic power terms */
union {
@@ -105,6 +107,7 @@ green_path_init(struct mem_allocator* allocator, struct green_path* path)
darray_flux_term_init(allocator, &path->flux_terms);
darray_power_term_init(allocator, &path->power_terms);
path->elapsed_time = -INF;
+ path->external_flux_term = 0;
path->limit.vertex = SDIS_RWALK_VERTEX_NULL;
path->limit.fragment = SDIS_INTERFACE_FRAGMENT_NULL;
path->limit_id = UINT_MAX;
@@ -127,6 +130,7 @@ green_path_copy(struct green_path* dst, const struct green_path* src)
res_T res = RES_OK;
ASSERT(dst && src);
dst->elapsed_time = src->elapsed_time;
+ dst->external_flux_term = src->external_flux_term;
dst->limit = src->limit;
dst->limit_id = src->limit_id;
dst->end_type = src->end_type;
@@ -145,6 +149,7 @@ green_path_copy_and_clear(struct green_path* dst, struct green_path* src)
res_T res = RES_OK;
ASSERT(dst && src);
dst->elapsed_time = src->elapsed_time;
+ dst->external_flux_term = src->external_flux_term;
dst->limit = src->limit;
dst->limit_id = src->limit_id;
dst->end_type = src->end_type;
@@ -164,6 +169,7 @@ green_path_copy_and_release(struct green_path* dst, struct green_path* src)
res_T res = RES_OK;
ASSERT(dst && src);
dst->elapsed_time = src->elapsed_time;
+ dst->external_flux_term = src->external_flux_term;
dst->limit = src->limit;
dst->limit_id = src->limit_id;
dst->end_type = src->end_type;
@@ -192,6 +198,7 @@ green_path_write(const struct green_path* path, FILE* stream)
/* Write elapsed time */
WRITE(&path->elapsed_time, 1);
+ WRITE(&path->external_flux_term, 1);
/* Write the list of flux terms */
sz = darray_flux_term_size_get(&path->flux_terms);
@@ -242,6 +249,7 @@ green_path_read(struct green_path* path, FILE* stream)
/* Read elapsed time */
READ(&path->elapsed_time, 1);
+ READ(&path->external_flux_term, 1);
/* Read the list of flux terms */
READ(&sz, 1);
@@ -397,7 +405,7 @@ green_function_solve_path
const size_t ipath,
double* weight)
{
- struct sdis_ambient_radiative_temperature trad =
+ struct sdis_ambient_radiative_temperature trad =
SDIS_AMBIENT_RADIATIVE_TEMPERATURE_NULL;
const struct power_term* power_terms = NULL;
const struct flux_term* flux_terms = NULL;
@@ -409,6 +417,7 @@ green_function_solve_path
struct sdis_interface_fragment frag = SDIS_INTERFACE_FRAGMENT_NULL;
double power;
double flux;
+ double external_flux;
double end_temperature;
size_t i, n;
res_T res = RES_OK;
@@ -420,7 +429,7 @@ green_function_solve_path
goto error;
}
- /* Compute medium power terms */
+ /* Compute medium powers */
power = 0;
n = darray_power_term_size_get(&path->power_terms);
power_terms = darray_power_term_cdata_get(&path->power_terms);
@@ -441,6 +450,13 @@ green_function_solve_path
flux += flux_terms[i].term * interface_side_get_flux(interf, &frag);
}
+ /* Compute external flux */
+ external_flux = 0;
+ if(green->scn->source) {
+ external_flux =
+ path->external_flux_term * source_get_power(green->scn->source);
+ }
+
/* Compute path's end temperature */
switch(path->end_type) {
case SDIS_GREEN_PATH_END_AT_INTERFACE:
@@ -466,7 +482,7 @@ green_function_solve_path
}
/* Compute the path weight */
- *weight = power + flux + end_temperature;
+ *weight = power + flux + external_flux + end_temperature;
exit:
return res;
@@ -1032,7 +1048,7 @@ error:
res_T
sdis_green_function_get_scene
- (const struct sdis_green_function* green,
+ (const struct sdis_green_function* green,
struct sdis_scene** scn)
{
if(!green || !scn) return RES_BAD_ARG;
@@ -1352,6 +1368,34 @@ error:
goto exit;
}
+res_T
+sdis_green_path_get_external_flux_term
+ (struct sdis_green_path* path_handle,
+ double* external_flux_term)
+{
+ const struct green_path* path = NULL;
+ struct sdis_green_function* green = NULL;
+ res_T res = RES_OK;
+
+ if(!path_handle || !external_flux_term) {
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ green = path_handle->green__; (void)green;
+ ASSERT(path_handle->id__ < darray_green_path_size_get(&green->paths));
+
+ path = darray_green_path_cdata_get(&green->paths) + path_handle->id__;
+
+ *external_flux_term = path->external_flux_term;
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
+
/*******************************************************************************
* Local functions
******************************************************************************/
@@ -1710,3 +1754,13 @@ exit:
error:
goto exit;
}
+
+res_T
+green_path_add_external_flux_term
+ (struct green_path_handle* handle,
+ const double val) /* [W/m^2/sr] */
+{
+ ASSERT(handle);
+ handle->path->external_flux_term += val;
+ return RES_OK;
+}
diff --git a/src/sdis_green.h b/src/sdis_green.h
@@ -20,9 +20,9 @@
#include <rsys/rsys.h>
/* Current version the green function data structure. One should increment it
- * and perform a version management onto serialized data when the gren function
+ * and perform a version management onto serialized data when the green function
* data structure is updated. */
-static const int SDIS_GREEN_FUNCTION_VERSION = 2;
+static const int SDIS_GREEN_FUNCTION_VERSION = 3;
/* Forward declaration */
struct accum;
@@ -106,5 +106,10 @@ green_path_add_flux_term
const struct sdis_interface_fragment* fragment,
const double term);
+extern LOCAL_SYM res_T
+green_path_add_external_flux_term
+ (struct green_path_handle* handle,
+ const double term); /* [W/m^2/sr] */
+
#endif /* SDIS_GREEN_H */
diff --git a/src/sdis_heat_path_boundary.c b/src/sdis_heat_path_boundary.c
@@ -41,3 +41,7 @@
#include "sdis_heat_path_boundary_Xd.h"
#define SDIS_XD_DIMENSION 3
#include "sdis_heat_path_boundary_Xd.h"
+#define SDIS_XD_DIMENSION 2
+#include "sdis_heat_path_boundary_Xd_handle_external_net_flux.h"
+#define SDIS_XD_DIMENSION 3
+#include "sdis_heat_path_boundary_Xd_handle_external_net_flux.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
@@ -0,0 +1,522 @@
+/* Copyright (C) 2016-2023 |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_boundary_c.h"
+#include "sdis_interface_c.h"
+#include "sdis_log.h"
+#include "sdis_scene_c.h"
+#include "sdis_source_c.h"
+
+#include <rsys/cstr.h> /* res_to_cstr */
+
+#include "sdis_Xd_begin.h"
+
+/*******************************************************************************
+ * Non generic data types and function
+ ******************************************************************************/
+#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__;
+
+/* 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)
+{
+ 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) {
+ 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) {
+ 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 */
+
+/*******************************************************************************
+ * Generic helper functions
+ ******************************************************************************/
+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)
+{
+ int net_flux = 0;
+ res_T res = RES_OK;
+
+ /* Handle bugs */
+ ASSERT(scn && func_name && args);
+ ASSERT(args->interf && args->frag);
+ ASSERT(!SXD_HIT_NONE(args->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);
+
+ net_flux = interface_side_is_external_flux_handled(args->interf, args->frag);
+ net_flux = net_flux && (scn->source != NULL);
+
+ if(net_flux && args->picard_order != 1) {
+ log_err(scn->dev,
+ "%s: Impossible to process external fluxes when Picard order is not "
+ "equal to 1; Picard order is currently set to %lu.\n",
+ func_name, (unsigned long)args->picard_order);
+ res = RES_BAD_ARG;
+ return res;
+ }
+
+ if(sdis_medium_get_type(args->interf->medium_back)
+ == sdis_medium_get_type(args->interf->medium_front)) {
+ log_err(scn->dev,
+ "%s: external fluxes can only be processed on fluid/solid interfaces.\n",
+ func_name);
+ res = RES_BAD_ARG;
+ return res;
+ }
+
+ 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,
+ struct source_sample* sample,
+ const double pos[DIM],
+ const struct sXd(hit)* hit_from)
+{
+ struct sXd(hit) hit = SXD_HIT_NULL;
+ ASSERT(scn && sample && pos && hit_from);
+
+ /* Is the source hidden */
+ XD(trace_ray)(scn, pos, sample->dir, sample->dst, hit_from, &hit);
+ if(!SXD_HIT_NONE(&hit)) return 0; /* [W/m^2/sr] */
+
+ /* Note that the value returned is not the source's actual radiance, but the
+ * radiance relative to the source's power. Care must therefore be taken to
+ * multiply it by the power of the source to obtain its real contribution.
+ * This trick makes it possible to manage the external flux in the green
+ * function. */
+ 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(setup_brdf)
+ (struct sdis_device* dev,
+ struct brdf* brdf,
+ const struct sdis_interface* interf,
+ const struct sdis_interface_fragment* frag)
+{
+ double epsilon = 0;
+ double alpha = 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);
+
+ epsilon = interface_side_get_emissivity(interf, 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, 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,
+ struct ssp_rng* rng,
+ const double in_pos[DIM], /* position */
+ const double in_N[DIM], /* Surface normal. (Away from the surface) */
+ const double time,
+ const struct sXd(hit)* in_hit, /* Current intersection */
+ double* out_flux) /* [W/m^2] */
+{
+ struct sXd(hit) hit = SXD_HIT_NULL;
+ double pos[3] = {0}; /* In 3D for ray tracing ray to the source */
+ double dir[3] = {0}; /* Incident direction (toward the surface). Always 3D.*/
+ double N[3] = {0}; /* Surface normal. Always 3D */
+ double incident_diffuse_flux = 0; /* [W/m^2] */
+ res_T res = RES_OK;
+ ASSERT(in_pos && in_N && in_hit);
+
+ /* Local copy of input argument */
+ dX(set)(pos, in_pos);
+ dX(set)(N, in_N);
+ hit = *in_hit;
+
+ /* Sample a diffusive direction in 3D */
+ ssp_ran_hemisphere_cos(rng, N, dir, NULL);
+
+ for(;;) {
+ /* External sources */
+ struct source_sample src_sample = SOURCE_SAMPLE_NULL;
+
+ /* Interface */
+ struct sdis_interface_fragment frag = SDIS_INTERFACE_FRAGMENT_NULL;
+ struct sdis_interface* interf = NULL;
+
+ /* BRDF */
+ struct brdf brdf = BRDF_NULL;
+ struct brdf_sample brdf_sample = BRDF_SAMPLE_NULL;
+
+ /* Miscellaneous */
+ double L = 0; /* incident direct flux to bounce position */
+ double wi[3] = {0}; /* Incident direction (outward the surface). Always 3D */
+ double vec[DIM] = {0}; /* Temporary variable */
+
+ d3_minus(wi, dir); /* Always in 3D */
+
+ /* Find the following surface along the direction of propagation */
+ XD(trace_ray)(scn, pos, dir, INF, &hit, &hit);
+ if(SXD_HIT_NONE(&hit)) break; /* No surface */
+
+ /* Retrieve the current position and normal */
+ dX(add)(pos, pos, dX(muld)(vec, 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, dir, time, N, &hit);
+
+ /* Check that the path reaches a valid interface */
+ res = check_interface(interf, &frag);
+ if(res != RES_OK) goto error;
+
+ XD(setup_brdf)(scn->dev, &brdf, interf, &frag);
+
+ /* Check if path is absorbed */
+ if(ssp_rng_canonical(rng) < brdf.emissivity) break;
+
+ /* Sample rebound direction */
+ if(frag.side == SDIS_BACK) dX(minus)(N, N); /* Revert normal if necessary */
+ sample_brdf(&brdf, rng, wi, N, &brdf_sample);
+ d3_set(dir, brdf_sample.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(res != RES_OK) goto error;
+
+ if(!SOURCE_SAMPLE_NONE(&src_sample)) {
+ const double Ld = XD(direct_contribution)(scn, &src_sample, pos, &hit);
+ L = Ld; /* [W/m^2] */
+ }
+
+ /* Calculate the direct contribution of the rebound is diffuse */
+ } else {
+ double cos_theta = 0;
+ ASSERT(brdf_sample.cpnt == BRDF_DIFFUSE);
+
+ /* Sample an external source to handle its direct contribution at the
+ * bounce position */
+ res = source_sample(scn->source, rng, pos, time, &src_sample);
+ CHK(res == RES_OK);
+ cos_theta = d3_dot(src_sample.dir, N);
+
+ /* The source is behind the surface */
+ if(cos_theta <= 0) {
+ L = 0; /* [W/m^2] */
+
+ /* The source is above the surface */
+ } else {
+ const double Ld = XD(direct_contribution)(scn, &src_sample, pos, &hit);
+ L = Ld * cos_theta / (PI * src_sample.pdf); /* [W/m^2] */
+ }
+ }
+ incident_diffuse_flux += L;
+ }
+
+ incident_diffuse_flux *= PI;
+
+exit:
+ *out_flux = incident_diffuse_flux;
+ return res;
+error:
+ goto exit;
+}
+
+/*******************************************************************************
+ * Local functions
+ ******************************************************************************/
+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)
+{
+ /* Sampling external sources */
+ struct source_sample src_sample = SOURCE_SAMPLE_NULL;
+
+ /* External flux */
+ double incident_flux = 0; /* [W/m^2] */
+ double incident_flux_diffuse = 0; /* [W/m^2] */
+ double incident_flux_direct = 0; /* [W/m^2] */
+ double net_flux = 0; /* [W/m^2] */
+ double external_flux_term = 0; /* [W/m^2] */
+
+ /* Sampled path */
+ double N[3] = {0}; /* Normal. Always in 3D */
+
+ /* On the fluid side */
+ struct sdis_interface_fragment frag = SDIS_INTERFACE_FRAGMENT_NULL;
+
+ /* Miscellaneous */
+ double emissivity = 0; /* Emissivity */
+ double Ld = 0; /* Incident radiance [W/m^2/sr] */
+ double cos_theta = 0;
+ int handle_flux = 0;
+ res_T res = RES_OK;
+ ASSERT(scn && args && T);
+
+ res = XD(check_handle_external_net_flux_args)(scn, FUNC_NAME, args);
+ if(res != RES_OK) goto error;
+
+ /* Setup the interface fragment on flud side */
+ frag = *args->frag;
+ if(sdis_medium_get_type(args->interf->medium_front) == SDIS_FLUID) {
+ frag.side = SDIS_FRONT;
+ } else {
+ ASSERT(sdis_medium_get_type(args->interf->medium_back) == SDIS_FLUID);
+ frag.side = SDIS_BACK;
+ }
+
+ /* No external sources <=> no external fluxes. Nothing to do */
+ handle_flux = interface_side_is_external_flux_handled(args->interf, &frag);
+ handle_flux = handle_flux && (scn->source != NULL);
+ if(!handle_flux) goto exit;
+
+ /* Sample the external source */
+ res = source_sample
+ (scn->source, rng, frag.P, frag.time, &src_sample);
+ if(res != RES_OK) goto error;
+
+ /* Setup the normal to ensure that it points toward the fluid medium */
+ dX(set)(N, frag.Ng);
+ if(frag.side == SDIS_BACK) dX(minus)(N, N);
+
+ /* Calculate the incident direct flux if the external source is above the
+ * 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);
+ 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);
+ if(res != RES_OK) goto error;
+
+ /* Calculate the global incident flux */
+ incident_flux = incident_flux_direct + incident_flux_diffuse; /* [W/m^2] */
+
+ /* Calculate the net flux */
+ emissivity = interface_side_get_emissivity(args->interf, &frag);
+ res = interface_side_check_emissivity(scn->dev, emissivity, frag.P, frag.time);
+ if(res != RES_OK) goto error;
+ net_flux = incident_flux * emissivity; /* [W/m^2] */
+
+ /* Until now, the net flux was calculated in relation to the source power.
+ * What is calculated is the external flux term of the green function. This
+ * must be multiplied by the source power to obtain the actual external flux*/
+ external_flux_term = net_flux / (args->h_radi + args->h_conv + args->h_cond);
+
+ /* Update the Monte Carlo weight */
+ T->value += external_flux_term * source_get_power(scn->source);
+
+ /* Register the external net flux term */
+ if(args->green_path) {
+ res = green_path_add_external_flux_term(args->green_path, external_flux_term);
+ if(res != RES_OK) goto error;
+ }
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
+#include "sdis_Xd_end.h"
diff --git a/src/sdis_heat_path_boundary_Xd_solid_fluid_picard1.h b/src/sdis_heat_path_boundary_Xd_solid_fluid_picard1.h
@@ -122,6 +122,10 @@ XD(solid_fluid_boundary_picard1_path)
/* 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);
+
/* 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);
@@ -238,6 +242,18 @@ XD(solid_fluid_boundary_picard1_path)
res = XD(handle_net_flux)(scn, &handle_net_flux_args, T);
if(res != RES_OK) goto error;
+ /* 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.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;
+ handle_external_net_flux_args.h_conv = h_conv;
+ handle_external_net_flux_args.h_radi = h_radi_hat;
+ res = XD(handle_external_net_flux)(scn, rng, &handle_external_net_flux_args, T);
+ if(res != RES_OK) goto error;
+
/* Fetch the last registered heat path vertex */
if(ctx->heat_path) hvtx = *heat_path_get_last_vertex(ctx->heat_path);
diff --git a/src/sdis_heat_path_boundary_c.h b/src/sdis_heat_path_boundary_c.h
@@ -169,6 +169,58 @@ handle_net_flux_3d
struct temperature_3d* T);
/*******************************************************************************
+ * Handle external flux
+ ******************************************************************************/
+struct handle_external_net_flux_args_2d {
+ struct sdis_interface* interf;
+ const struct sdis_interface_fragment* frag;
+ const struct s2d_hit* hit;
+
+ struct green_path_handle* green_path; /* Store the propagator */
+ struct sdis_heat_path* heat_path; /* Save paths */
+
+ size_t picard_order;
+ double h_cond; /* Convective coefficient, i.e. lambda/delta */
+ double h_conv; /* Condutive coefficient */
+ double h_radi; /* Radiative coefficient */
+};
+
+struct handle_external_net_flux_args_3d {
+ struct sdis_interface* interf;
+ const struct sdis_interface_fragment* frag;
+ const struct s3d_hit* hit;
+
+ struct green_path_handle* green_path; /* Store the propagator */
+ struct sdis_heat_path* heat_path; /* Save paths */
+
+ size_t picard_order;
+ double h_cond; /* Convective coefficient, i.e. lambda/delta */
+ double h_conv; /* Condutive coefficient */
+ 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;
+
+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);
+
+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);
+
+/*******************************************************************************
* Boundary sub-paths
******************************************************************************/
extern LOCAL_SYM res_T
diff --git a/src/sdis_interface_c.h b/src/sdis_interface_c.h
@@ -17,6 +17,8 @@
#define SDIS_INTERFACE_C_H
#include "sdis.h"
+#include "sdis_log.h"
+
#include <rsys/free_list.h>
#include <rsys/ref_count.h>
#include <float.h>
@@ -183,5 +185,76 @@ interface_side_get_reference_temperature
? shader->reference_temperature(frag, interf->data) : -1;
}
-#endif /* SDIS_INTERFACE_C_H */
+static INLINE int
+interface_side_is_external_flux_handled
+ (const struct sdis_interface* interf,
+ const struct sdis_interface_fragment* frag)
+{
+ const struct sdis_interface_side_shader* shader;
+ ASSERT(interf && frag);
+ switch(frag->side) {
+ case SDIS_FRONT: shader = &interf->shader.front; break;
+ case SDIS_BACK: shader = &interf->shader.back; break;
+ default: FATAL("Unreachable code\n"); break;
+ }
+ return shader->handle_external_flux;
+}
+
+/*******************************************************************************
+ * Check interface properties
+ ******************************************************************************/
+#define DEFINE_INTERF_CHK_PROP_FUNC(Interf, Prop, Low, Upp, LowIsInc, UppIsInc)\
+ static INLINE res_T \
+ Interf##_check_##Prop \
+ (struct sdis_device* dev, \
+ const double val, /* Value of the property */ \
+ const double pos[3], /* Position at which the property was queried */ \
+ const double time) /* Time at which the property was queried */ \
+ { \
+ const int low_test = LowIsInc ? Low <= val : Low < val; \
+ const int upp_test = UppIsInc ? Upp >= val : Upp > val; \
+ const char low_char = LowIsInc ? '[' : ']'; \
+ const char upp_char = UppIsInc ? ']' : '['; \
+ ASSERT(dev && pos); \
+ \
+ if(!low_test || !upp_test) { \
+ log_err(dev, \
+ "invalid "STR(Interf)" "PROP_STR(Prop)" '%g': " \
+ "it must be in %c%g, %g%c -- position=%g, %g, %g; time=%g\n", \
+ val, low_char, (double)Low, (double)Upp, upp_char, SPLIT3(pos), time); \
+ return RES_BAD_ARG; \
+ } \
+ return RES_OK; \
+ }
+#define PROP_STR(Prop) CONCAT(PROP_STR_, Prop)
+#define PROP_STR_convection_coef "convection coefficient"
+#define PROP_STR_thermal_contact_resistance "thermal contact resistance"
+#define PROP_STR_convection_coef_upper_bound "convection coefficient upper bound"
+#define PROP_STR_temperature "temperature"
+#define PROP_STR_flux "net flux"
+#define PROP_STR_emissivity "emissivity"
+#define PROP_STR_specular_fraction "specular fraction"
+#define PROP_STR_reference_temperature "reference temperature"
+
+DEFINE_INTERF_CHK_PROP_FUNC(interface, convection_coef, 0, INF, 1, 1)
+DEFINE_INTERF_CHK_PROP_FUNC(interface, thermal_contact_resistance, 0, INF, 1, 1)
+DEFINE_INTERF_CHK_PROP_FUNC(interface, convection_coef_upper_bound, 0, INF, 1, 1)
+DEFINE_INTERF_CHK_PROP_FUNC(interface_side, temperature, 0, INF, 1, 1)
+DEFINE_INTERF_CHK_PROP_FUNC(interface_side, flux, -INF, INF, 1, 1)
+DEFINE_INTERF_CHK_PROP_FUNC(interface_side, emissivity, 0, 1, 1, 1)
+DEFINE_INTERF_CHK_PROP_FUNC(interface_side, specular_fraction, 0, 1, 1, 1)
+DEFINE_INTERF_CHK_PROP_FUNC(interface_side, reference_temperature, 0, INF, 1, 1)
+
+#undef DEFINE_INTERF_CHK_PROP_FUNC
+#undef PROP_STR
+#undef PROP_STR_convection_coef
+#undef PROP_STR_thermal_contact_resistance
+#undef PROP_STR_convection_coef_upper_bound
+#undef PROP_STR_temperature
+#undef PROP_STR_flux
+#undef PROP_STR_emissivity
+#undef PROP_STR_specular_fraction
+#undef PROP_STR_reference_temperature
+
+#endif /* SDIS_INTERFACE_C_H */
diff --git a/src/sdis_scene.c b/src/sdis_scene.c
@@ -24,6 +24,7 @@
#include "sdis.h"
#include "sdis_interface_c.h"
#include "sdis_scene_c.h"
+#include "sdis_source_c.h"
#include <float.h>
#include <limits.h>
@@ -108,6 +109,7 @@ scene_release(ref_T * ref)
if(scn->s3d_view) S3D(scene_view_ref_put(scn->s3d_view));
if(scn->senc2d_scn) SENC2D(scene_ref_put(scn->senc2d_scn));
if(scn->senc3d_scn) SENC3D(scene_ref_put(scn->senc3d_scn));
+ if(scn->source) SDIS(source_ref_put(scn->source));
MEM_RM(dev->allocator, scn);
SDIS(device_ref_put(dev));
}
@@ -418,6 +420,14 @@ sdis_scene_get_device(struct sdis_scene* scn, struct sdis_device** device)
return RES_OK;
}
+res_T
+sdis_scene_get_source(struct sdis_scene* scn, struct sdis_source** source)
+{
+ if(!scn || !source) return RES_BAD_ARG;
+ *source = scn->source;
+ return RES_OK;
+}
+
/*******************************************************************************
* Local miscellaneous function
******************************************************************************/
@@ -472,6 +482,13 @@ scene_compute_hash(const struct sdis_scene* scn, hash256_T hash)
SHA256_UPD(&scn->trad.reference, 1);
SHA256_UPD(&scn->tmax, 1);
SHA256_UPD(&scn->fp_to_meter, 1);
+
+ if(scn->source) {
+ hash256_T src_hash;
+ source_compute_signature(scn->source, src_hash);
+ sha256_ctx_update(&sha256_ctx, src_hash, sizeof(hash256_T));
+ }
+
FOR_EACH(iprim, 0, nprims) {
struct sdis_interface* interf = NULL;
size_t ivert;
diff --git a/src/sdis_scene_Xd.h b/src/sdis_scene_Xd.h
@@ -354,6 +354,7 @@ hit_shared_edge
int edge_ivertex = 0; /* Temporary variable */
int tri0_ivertex, tri1_ivertex;
int iv0, iv1, iv2;
+ int hit_edge;
ASSERT(tri0 && tri1 && pos0 && pos1);
/* Fetch the vertices of the triangle 0 */
@@ -432,8 +433,15 @@ hit_shared_edge
f3_sub(E1, tri1_vertices[tri1_edge[1]].value, pos1);
tmp1_2area = f3_len(f3_cross(N1, E0, E1));
- return (eq_epsf(tri0_2area, 0, 1.e-6f) || tmp0_2area/tri0_2area < ON_EDGE_EPSILON)
- && (eq_epsf(tri1_2area, 0, 1.e-6f) || tmp1_2area/tri1_2area < ON_EDGE_EPSILON);
+ hit_edge =
+ ( eq_epsf(tri0_2area, 0, 1.e-6f)
+ || eq_epsf(tmp0_2area, 0, 1.e-6f)
+ || tmp0_2area/tri0_2area < ON_EDGE_EPSILON);
+ hit_edge = hit_edge &&
+ ( eq_epsf(tri1_2area, 0, 1.e-6f)
+ || eq_epsf(tmp1_2area, 0, 1.e-6f)
+ || tmp1_2area/tri1_2area < ON_EDGE_EPSILON);
+ return hit_edge;
}
#undef ON_EDGE_EPSILON
#endif /* DIM == 2 */
@@ -929,6 +937,11 @@ XD(scene_create)
htable_enclosure_init(dev->allocator, &scn->enclosures);
htable_d_init(dev->allocator, &scn->tmp_hc_ub);
+ if(args->source) {
+ SDIS(source_ref_get(args->source));
+ scn->source = args->source;
+ }
+
res = XD(run_analyze)
(scn,
args->nprimitives,
diff --git a/src/sdis_scene_c.h b/src/sdis_scene_c.h
@@ -218,6 +218,8 @@ struct sdis_scene {
double tmin; /* Minimum temperature of the system (In Kelvin) */
double tmax; /* Maximum temperature of the system (In Kelvin) */
+ struct sdis_source* source; /* External source. May be NULL */
+
ref_T ref;
struct sdis_device* dev;
};
diff --git a/src/sdis_source.c b/src/sdis_source.c
@@ -0,0 +1,311 @@
+/* Copyright (C) 2016-2023 |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_device_c.h"
+#include "sdis_log.h"
+#include "sdis_source_c.h"
+
+#include <rsys/mem_allocator.h>
+#include <rsys/ref_count.h>
+
+#include <star/ssp.h>
+
+struct sdis_source {
+ struct sdis_spherical_source_create_args spherical;
+
+ struct sdis_device* dev;
+ ref_T ref;
+};
+
+/*******************************************************************************
+ * Helper functions
+ ******************************************************************************/
+static res_T
+check_spherical_source_create_args
+ (struct sdis_device* dev,
+ const char* func_name,
+ struct sdis_spherical_source_create_args* args)
+{
+ ASSERT(func_name);
+ if(!args) return RES_BAD_ARG;
+
+ if(!args->position) {
+ log_err(dev, "%s: the position functor is missing.\n", func_name);
+ return RES_BAD_ARG;
+ }
+
+ if(args->radius < 0) {
+ log_err(dev, "%s: invalid source radius '%g' m. It cannot be negative.\n",
+ func_name, args->radius);
+ return RES_BAD_ARG;
+ }
+
+ if(args->power < 0) {
+ log_err(dev, "%s: invalid source power '%g' W. It cannot be negative.\n",
+ func_name, args->power);
+ return RES_BAD_ARG;
+ }
+
+ return RES_OK;
+}
+
+static void
+release_source(ref_T* ref)
+{
+ struct sdis_device* dev = NULL;
+ struct sdis_source* src = CONTAINER_OF(ref, struct sdis_source, ref);
+ ASSERT(ref);
+ dev = src->dev;
+ if(src->spherical.data) SDIS(data_ref_put(src->spherical.data));
+ MEM_RM(dev->allocator, src);
+ SDIS(device_ref_put(dev));
+}
+
+/*******************************************************************************
+ * Exported symbols
+ ******************************************************************************/
+res_T
+sdis_spherical_source_create
+ (struct sdis_device* dev,
+ struct sdis_spherical_source_create_args* args,
+ struct sdis_source** out_src)
+{
+ struct sdis_source* src = NULL;
+ res_T res = RES_OK;
+
+ if(!dev || !out_src) { res = RES_BAD_ARG; goto error; }
+ res = check_spherical_source_create_args(dev, FUNC_NAME, args);
+ if(res != RES_OK) goto error;
+
+ src = MEM_CALLOC(dev->allocator, 1, sizeof(*src));
+ if(!src) {
+ log_err(dev, "%s: cannot allocate spherical source.\n", FUNC_NAME);
+ res = RES_OK;
+ goto error;
+ }
+ ref_init(&src->ref);
+ SDIS(device_ref_get(dev));
+ if(args->data) SDIS(data_ref_get(args->data));
+ src->spherical = *args;
+ src->dev = dev;
+
+exit:
+ if(out_src) *out_src = src;
+ return res;
+error:
+ if(src) { SDIS(source_ref_put(src)); src = NULL; }
+ goto exit;
+}
+
+res_T
+sdis_source_ref_get(struct sdis_source* src)
+{
+ if(!src) return RES_BAD_ARG;
+ ref_get(&src->ref);
+ return RES_OK;
+}
+
+res_T
+sdis_source_ref_put(struct sdis_source* src)
+{
+ if(!src) return RES_BAD_ARG;
+ ref_put(&src->ref, release_source);
+ return RES_OK;
+}
+
+res_T
+sdis_source_get_power(const struct sdis_source* src, double* power)
+{
+ if(!src || !power) return RES_BAD_ARG;
+ *power = source_get_power(src);
+ return RES_OK;
+}
+
+/*******************************************************************************
+ * Local functions
+ ******************************************************************************/
+res_T
+source_sample
+ (const struct sdis_source* src,
+ struct ssp_rng* rng,
+ const double pos[3],
+ const double time,
+ struct source_sample* sample)
+{
+ double src_pos[3]; /* [m] */
+ double main_dir[3];
+ double half_angle; /* [radians] */
+ double cos_half_angle; /* [radians] */
+ double dst; /* [m] */
+ double radius; /* Source radius [m] */
+ double area; /* Source area [m^2] */
+ res_T res = RES_OK;
+ ASSERT(src && rng && pos && sample);
+
+ /* Retrieve current source position and radius */
+ src->spherical.position(time, src_pos, src->spherical.data);
+ radius = src->spherical.radius;
+
+ area = 4*PI*radius*radius; /* [m^2] */
+
+ /* compute the direction of `pos' toward the center of the source */
+ d3_sub(main_dir, src_pos, pos);
+
+ /* Normalize the direction and keep the distance from `pos' to the center of
+ * the source */
+ dst = d3_normalize(main_dir, main_dir);
+ if(dst <= radius) {
+ log_err(src->dev,
+ "%s: the position from which the external source is sampled "
+ "is included in the source:\n"
+ "\tsource position = %g, %g, %g\n"
+ "\tsource radius = %g\n"
+ "\tposition = %g, %g, %g\n"
+ "\ttime = %g\n"
+ "\tdistance from position to source = %g\n",
+ FUNC_NAME, SPLIT3(src_pos), radius, SPLIT3(pos), time, dst);
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ /* Point source */
+ if(area == 0) {
+ d3_set(sample->dir, main_dir);
+ sample->pdf = 1;
+ sample->dst = dst;
+ sample->power = src->spherical.power; /* [W] */
+ sample->radiance_term = 1.0 / (4*PI*dst*dst); /* [W/m^2/sr] */
+ sample->radiance = sample->power * sample->radiance_term; /* [W/m^2/sr] */
+
+ /* Spherical source */
+ } else {
+ /* Sample the source according to its solid angle,
+ * i.e. 2*PI*(1 - cos(half_angle)) */
+ half_angle = asin(radius/dst);
+ cos_half_angle = cos(half_angle);
+ ssp_ran_sphere_cap_uniform /* pdf = 1/(2*PI*(1-cos(half_angle))) */
+ (rng, main_dir, cos_half_angle, sample->dir, &sample->pdf);
+
+ /* Set other sample variables */
+ sample->dst = dst - radius; /* From pos to source boundaries [m] */
+ sample->power = src->spherical.power; /* [W] */
+ sample->radiance_term = 1.0 / (PI*area); /* [W/m^2/sr] */
+ sample->radiance = sample->power * sample->radiance_term; /* [W/m^2/sr] */
+ }
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
+res_T
+source_trace_to
+ (const struct sdis_source* src,
+ const double pos[3], /* Ray origin */
+ const double dir[3], /* Ray direction */
+ const double time, /* Time at which ray is traced */
+ struct source_sample* sample)
+{
+ double src_pos[3]; /* [m] */
+ double main_dir[3];
+ double radius; /* [m] */
+ double dst; /* Distance from pos to the source center [m] */
+ double half_angle; /* [radian] */
+ res_T res = RES_OK;
+ ASSERT(src && pos && dir && sample);
+ ASSERT(d3_is_normalized(dir));
+
+ radius = src->spherical.radius;
+
+ /* Point sources cannot be targeted */
+ if(radius == 0) {
+ *sample = SOURCE_SAMPLE_NULL;
+ goto exit;
+ }
+
+ /* Retrieve current source position and radius */
+ src->spherical.position(time, src_pos, src->spherical.data);
+
+ /* compute the direction of `pos' toward the center of the source */
+ d3_sub(main_dir, src_pos, pos);
+
+ /* Normalize the direction and keep the distance from `pos' to the center of
+ * the source */
+ dst = d3_normalize(main_dir, main_dir);
+ if(dst <= radius) {
+ log_err(src->dev,
+ "%s: the position from which the external source is targeted "
+ "is included in the source:\n"
+ "\tsource position = %g, %g, %g\n"
+ "\tsource radius = %g\n"
+ "\tposition = %g, %g, %g\n"
+ "\ttime = %g\n"
+ "\tdistance from position to source = %g\n",
+ FUNC_NAME, SPLIT3(src_pos), radius, SPLIT3(pos), time, dst);
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ /* Compute the half angle of the source as seen from pos */
+ half_angle = asin(radius/dst);
+
+ /* The source is missed */
+ if(d3_dot(dir, main_dir) < cos(half_angle)) {
+ *sample = SOURCE_SAMPLE_NULL;
+
+ /* The source is intersected */
+ } else {
+ const double area = 4*PI*radius*radius; /* [m^2] */
+
+ d3_set(sample->dir, dir);
+ sample->pdf = 1;
+ sample->dst = dst - radius; /* From pos to source boundaries [m] */
+ sample->power = src->spherical.power; /* [W] */
+ sample->radiance_term = 1.0 / (PI*area); /* [W/m^2/sr] */
+ sample->radiance = sample->power * sample->radiance_term; /* [W/m^2/sr] */
+ }
+
+exit:
+ return res;
+error:
+ *sample = SOURCE_SAMPLE_NULL;
+ goto exit;
+}
+
+double
+source_get_power(const struct sdis_source* src)
+{
+ ASSERT(src);
+ return src->spherical.power;
+}
+
+void
+source_compute_signature(const struct sdis_source* src, hash256_T hash)
+{
+ struct sha256_ctx ctx;
+ ASSERT(src && hash);
+
+ /* Calculate the source signature. Currently, it is only the source radius.
+ * But the Source API is designed to be independent of source type. In the
+ * future, the source will not necessarily be spherical, so the data to be
+ * hashed will depend on the type of source. This function anticipate this by
+ * calculating a hash even if it is currently dispensable. */
+ sha256_ctx_init(&ctx);
+ sha256_ctx_update
+ (&ctx, (const char*)&src->spherical.radius, sizeof(src->spherical.radius));
+ sha256_ctx_finalize(&ctx, hash);
+}
diff --git a/src/sdis_source_c.h b/src/sdis_source_c.h
@@ -0,0 +1,73 @@
+/* Copyright (C) 2016-2023 |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_SOURCE_C_H
+#define SDIS_SOURCE_C_H
+
+#include <rsys/hash.h>
+#include <rsys/rsys.h>
+
+struct sdis_source;
+struct ssp_rng;
+
+struct source_sample {
+ double dir[3]; /* Direction _to_ the source */
+ double pdf; /* pdf of sampled direction */
+ double dst; /* Distance to the source [m] */
+ double power; /* [W] */
+ double radiance; /* [W/m^2/sr] */
+
+ /* Radiance relative to power, i.e. the source power is assumed to be equal to
+ * 1. It must be multiplied by the source power to obtain the actual radiance
+ * of the source. In other words, this variable defines the contribution of
+ * the source independently of its power, and can therefore be recorded in the
+ * green function */
+ double radiance_term; /* [W/m^2/sr] */
+};
+#define SOURCE_SAMPLE_NULL__ {{0,0,0}, 0, 0, 0, 0, 0}
+static const struct source_sample SOURCE_SAMPLE_NULL = SOURCE_SAMPLE_NULL__;
+
+/* Helper macro used to define whether a sample is valid or not */
+#define SOURCE_SAMPLE_NONE(Sample) ((Sample)->pdf == 0)
+
+extern LOCAL_SYM res_T
+source_sample
+ (const struct sdis_source* source,
+ struct ssp_rng* rng,
+ const double pos[3], /* Position from which the source is sampled */
+ const double time, /* Time at which the source is sampled */
+ struct source_sample* sample);
+
+/* Trace a ray toward the source. The returned sample has a pdf of 1 or 0
+ * whether the source is intersected by the ray or not, respectively. You can
+ * use the SOURCE_SAMPLE_NONE macro to check this */
+extern LOCAL_SYM res_T
+source_trace_to
+ (const struct sdis_source* source,
+ const double pos[3], /* Ray origin */
+ const double dir[3], /* Ray direction */
+ const double time, /* Time at which ray is traced */
+ struct source_sample* sample); /* pdf == 0 if no source is reached */
+
+extern LOCAL_SYM double /* [W] */
+source_get_power
+ (const struct sdis_source* source);
+
+extern LOCAL_SYM void
+source_compute_signature
+ (const struct sdis_source* source,
+ hash256_T hash);
+
+#endif /* SDIS_SOURCE_C_H */
diff --git a/src/test_sdis_draw_external_flux.c b/src/test_sdis_draw_external_flux.c
@@ -0,0 +1,445 @@
+/* Copyright (C) 2016-2023 |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 "test_sdis_utils.h"
+
+#include <star/s3dut.h>
+#include <rsys/stretchy_array.h>
+
+#define IMG_WIDTH 320
+#define IMG_HEIGHT 240
+#define IMG_SPP 64
+#define SOURCE_POWER 30e6 /* W */
+#define EMISSIVITY 0.5
+#define T_RAD 300 /* [K] */
+#define T_REF 300 /* [K] */
+
+/*
+ * The system consists of a floor (i.e. a parallelepiped), a super-form floating
+ * above it and a spherical external source. The test attempts to render the
+ * scene and thus verify the entire sampling of conductive-radiative paths
+ * coupled to an external flux from the external source.
+ */
+
+/*******************************************************************************
+ * Geometries
+ ******************************************************************************/
+struct mesh {
+ double* positions; /* List of double 3 */
+ size_t* indices; /* List of size_t 3 */
+};
+#define MESH_NULL__ {NULL, NULL}
+static const struct mesh MESH_NULL = MESH_NULL__;
+
+static INLINE void
+mesh_init(struct mesh* mesh)
+{
+ CHK(mesh);
+ *mesh = MESH_NULL;
+}
+
+static INLINE void
+mesh_release(struct mesh* mesh)
+{
+ CHK(mesh);
+ sa_release(mesh->positions);
+ sa_release(mesh->indices);
+}
+
+/* Number of vertices */
+static INLINE size_t
+mesh_nvertices(const struct mesh* mesh)
+{
+ CHK(mesh);
+ return sa_size(mesh->positions) / 3/* #coords per vertex */;
+}
+
+/* Number of triangles */
+static INLINE size_t
+mesh_ntriangles(const struct mesh* mesh)
+{
+ CHK(mesh);
+ return sa_size(mesh->indices) / 3/* #indices per triangle */;
+}
+
+static void
+mesh_append
+ (struct mesh* mesh,
+ const struct s3dut_mesh_data* mesh_data,
+ const double in_translate[3]) /* May be NULL */
+{
+ double translate[3] = {0, 0, 0};
+ double* positions = NULL;
+ size_t* indices = NULL;
+ size_t ivert = 0;
+ size_t i = 0;
+ CHK(mesh != NULL);
+
+ ivert = mesh_nvertices(mesh);
+ positions = sa_add(mesh->positions, mesh_data->nvertices*3);
+ indices = sa_add(mesh->indices, mesh_data->nprimitives*3);
+
+ if(in_translate) {
+ translate[0] = in_translate[0];
+ translate[1] = in_translate[1];
+ translate[2] = in_translate[2];
+ }
+
+ FOR_EACH(i, 0, mesh_data->nvertices) {
+ positions[i*3 + 0] = mesh_data->positions[i*3 + 0] + translate[0];
+ positions[i*3 + 1] = mesh_data->positions[i*3 + 1] + translate[1];
+ positions[i*3 + 2] = mesh_data->positions[i*3 + 2] + translate[2];
+ }
+
+ FOR_EACH(i, 0, mesh_data->nprimitives) {
+ indices[i*3 + 0] = mesh_data->indices[i*3 + 0] + ivert;
+ indices[i*3 + 1] = mesh_data->indices[i*3 + 1] + ivert;
+ indices[i*3 + 2] = mesh_data->indices[i*3 + 2] + ivert;
+ }
+}
+
+static void
+mesh_add_super_shape(struct mesh* mesh)
+{
+ struct s3dut_mesh* sshape = NULL;
+ struct s3dut_mesh_data sshape_data;
+ struct s3dut_super_formula f0 = S3DUT_SUPER_FORMULA_NULL;
+ struct s3dut_super_formula f1 = S3DUT_SUPER_FORMULA_NULL;
+ const double radius = 1;
+ const unsigned nslices = 256;
+
+ f0.A = 1.5; f0.B = 1; f0.M = 11.0; f0.N0 = 1; f0.N1 = 1; f0.N2 = 2.0;
+ f1.A = 1.0; f1.B = 2; f1.M = 3.6; f1.N0 = 1; f1.N1 = 2; f1.N2 = 0.7;
+ OK(s3dut_create_super_shape(NULL, &f0, &f1, radius, nslices, nslices/2, &sshape));
+ OK(s3dut_mesh_get_data(sshape, &sshape_data));
+ mesh_append(mesh, &sshape_data, NULL);
+ OK(s3dut_mesh_ref_put(sshape));
+}
+
+static void
+mesh_add_ground(struct mesh* mesh)
+{
+ #define DEPTH 2.0
+ const double translate[3] = {0, 0, -DEPTH-1};
+
+ struct s3dut_mesh* ground = NULL;
+ struct s3dut_mesh_data ground_data;
+
+ OK(s3dut_create_cuboid(NULL, 20, 20, DEPTH, &ground));
+ OK(s3dut_mesh_get_data(ground, &ground_data));
+ mesh_append(mesh, &ground_data, translate);
+ OK(s3dut_mesh_ref_put(ground));
+ #undef DEPTH
+}
+
+/*******************************************************************************
+ * Media & interface
+ ******************************************************************************/
+#define MEDIUM_PROP(Type, Prop, Val) \
+ static double \
+ Type##_get_##Prop \
+ (const struct sdis_rwalk_vertex* vtx, \
+ struct sdis_data* data) \
+ { \
+ (void)vtx, (void)data; /* Avoid the "unused variable" warning */ \
+ return Val; \
+ }
+MEDIUM_PROP(solid, calorific_capacity, 500.0) /* [J/K/Kg] */
+MEDIUM_PROP(solid, thermal_conductivity, 25.0) /* [W/m/K] */
+MEDIUM_PROP(solid, volumic_mass, 7500.0) /* [kg/m^3] */
+MEDIUM_PROP(solid, temperature, 310) /* [K] */
+MEDIUM_PROP(solid, delta, 1.0/20.0) /* [m] */
+MEDIUM_PROP(fluid, calorific_capacity, 2.0) /* [J/K/Kg] */
+MEDIUM_PROP(fluid, volumic_mass, 25.0) /* |kg/m^3] */
+MEDIUM_PROP(fluid, temperature, -1/*<=> unknown*/) /* [K] */
+#undef MEDIUM_PROP
+
+static struct sdis_medium*
+create_solid(struct sdis_device* sdis)
+{
+ struct sdis_solid_shader shader = SDIS_SOLID_SHADER_NULL;
+ struct sdis_medium* solid = NULL;
+
+ shader.calorific_capacity = solid_get_calorific_capacity;
+ shader.thermal_conductivity = solid_get_thermal_conductivity;
+ shader.volumic_mass = solid_get_volumic_mass;
+ shader.delta = solid_get_delta;
+ shader.temperature = solid_get_temperature;
+ OK(sdis_solid_create(sdis, &shader, NULL, &solid));
+ return solid;
+}
+
+static struct sdis_medium*
+create_fluid(struct sdis_device* sdis)
+{
+ struct sdis_fluid_shader shader = SDIS_FLUID_SHADER_NULL;
+ struct sdis_medium* fluid = NULL;
+
+ shader.calorific_capacity = fluid_get_calorific_capacity;
+ shader.volumic_mass = fluid_get_volumic_mass;
+ shader.temperature = fluid_get_temperature;
+ OK(sdis_fluid_create(sdis, &shader, NULL, &fluid));
+ return fluid;
+}
+
+static double
+interface_get_reference_temperature
+ (const struct sdis_interface_fragment* frag,
+ struct sdis_data* data)
+{
+ (void)frag, (void)data; /* Avoid the "unused variable" warning */
+ return T_REF;
+}
+
+static double
+interface_get_temperature
+ (const struct sdis_interface_fragment* frag,
+ struct sdis_data* data)
+{
+ (void)frag, (void)data;/* Avoid the "unused variable" warning */
+ return -1;
+}
+
+static double
+interface_get_emissivity
+ (const struct sdis_interface_fragment* frag,
+ struct sdis_data* data)
+{
+ (void)frag, (void)data;/* Avoid the "unused variable" warning */
+ return EMISSIVITY;
+}
+
+static struct sdis_interface*
+create_interface
+ (struct sdis_device* sdis,
+ struct sdis_medium* front,
+ struct sdis_medium* back)
+{
+ struct sdis_interface* interf = NULL;
+ struct sdis_interface_shader shader = SDIS_INTERFACE_SHADER_NULL;
+
+ shader.front.temperature = interface_get_temperature;
+ shader.front.reference_temperature = interface_get_reference_temperature;
+ shader.front.emissivity = interface_get_emissivity;
+ shader.back.temperature = interface_get_temperature;
+ OK(sdis_interface_create(sdis, front, back, &shader, NULL, &interf));
+ return interf;
+}
+
+/*******************************************************************************
+ * External source
+ ******************************************************************************/
+static void
+source_get_position
+ (const double time,
+ double pos[3],
+ struct sdis_data* data)
+{
+ (void)time, (void)data; /* Avoid the "unusued variable" warning */
+ pos[0] = 5.0;
+ pos[1] = 5.0;
+ pos[2] = 5.0;
+}
+
+static struct sdis_source*
+create_source(struct sdis_device* sdis)
+{
+ struct sdis_spherical_source_create_args args = SDIS_SPHERICAL_SOURCE_CREATE_ARGS_NULL;
+ struct sdis_source* source = NULL;
+
+ args.position = source_get_position;
+ args.data = NULL;
+ args.radius = 3e-1; /* [m] */
+ args.power = SOURCE_POWER; /* [W] */
+ OK(sdis_spherical_source_create(sdis, &args, &source));
+ return source;
+}
+
+/*******************************************************************************
+ * Scene, i.e. the system to simulate
+ ******************************************************************************/
+struct scene_context {
+ const struct mesh* mesh;
+ struct sdis_interface* interf;
+};
+static const struct scene_context SCENE_CONTEXT_NULL = {NULL, NULL};
+
+static void
+scene_get_indices(const size_t itri, size_t ids[3], void* ctx)
+{
+ struct scene_context* context = ctx;
+ CHK(ids && context && itri < mesh_ntriangles(context->mesh));
+ ids[0] = (unsigned)context->mesh->indices[itri*3+0];
+ ids[1] = (unsigned)context->mesh->indices[itri*3+1];
+ ids[2] = (unsigned)context->mesh->indices[itri*3+2];
+}
+
+static void
+scene_get_interface(const size_t itri, struct sdis_interface** interf, void* ctx)
+{
+ struct scene_context* context = ctx;
+ CHK(interf && context && itri < mesh_ntriangles(context->mesh));
+ *interf = context->interf;
+}
+
+static void
+scene_get_position(const size_t ivert, double pos[3], void* ctx)
+{
+ struct scene_context* context = ctx;
+ CHK(pos && context && ivert < mesh_nvertices(context->mesh));
+ pos[0] = context->mesh->positions[ivert*3+0];
+ pos[1] = context->mesh->positions[ivert*3+1];
+ pos[2] = context->mesh->positions[ivert*3+2];
+}
+
+static struct sdis_scene*
+create_scene
+ (struct sdis_device* sdis,
+ const struct mesh* mesh,
+ struct sdis_interface* interf,
+ struct sdis_source* source)
+{
+ struct sdis_scene* scn = NULL;
+ struct sdis_scene_create_args scn_args = SDIS_SCENE_CREATE_ARGS_DEFAULT;
+ struct scene_context context = SCENE_CONTEXT_NULL;
+
+ context.mesh = mesh;
+ context.interf = interf;
+
+ scn_args.get_indices = scene_get_indices;
+ scn_args.get_interface = scene_get_interface;
+ scn_args.get_position = scene_get_position;
+ scn_args.nprimitives = mesh_ntriangles(mesh);
+ scn_args.nvertices = mesh_nvertices(mesh);
+ scn_args.trad.temperature = T_RAD;
+ scn_args.trad.reference = T_REF;
+ scn_args.t_range[0] = MMIN(T_RAD, T_REF);
+ scn_args.t_range[1] = MMAX(T_RAD, T_REF);
+ scn_args.source = source;
+ scn_args.context = &context;
+ OK(sdis_scene_create(sdis, &scn_args, &scn));
+ return scn;
+}
+
+/*******************************************************************************
+ * Rendering point of view
+ ******************************************************************************/
+static struct sdis_camera*
+create_camera(struct sdis_device* sdis)
+{
+ const double pos[3] = {-3.81, 11.23, 5.29};
+ const double tgt[3] = {-0.46, 0, -0.32};
+ const double up[3] = {0, 0, 1};
+ struct sdis_camera* cam = NULL;
+
+ OK(sdis_camera_create(sdis, &cam));
+ OK(sdis_camera_set_proj_ratio(cam, (double)IMG_WIDTH/(double)IMG_HEIGHT));
+ OK(sdis_camera_set_fov(cam, MDEG2RAD(30)));
+ OK(sdis_camera_look_at(cam, pos, tgt, up));
+ return cam;
+}
+
+/*******************************************************************************
+ * Draw the submitted scene
+ ******************************************************************************/
+/* Write an image in htrdr-image(5) format */
+static void
+write_image(FILE* stream, struct sdis_estimator_buffer* image)
+{
+ size_t x, y;
+
+ /* Header */
+ fprintf(stream, "%d %d\n", IMG_WIDTH, IMG_HEIGHT);
+
+ /* Pixels ordered by row */
+ FOR_EACH(y, 0, IMG_HEIGHT) {
+ FOR_EACH(x, 0, IMG_WIDTH) {
+ const struct sdis_estimator* pixel = NULL;
+ struct sdis_mc temperature = SDIS_MC_NULL;
+
+ OK(sdis_estimator_buffer_at(image, x, y, &pixel));
+ OK(sdis_estimator_get_temperature(pixel, &temperature));
+ fprintf(stream, "%g %g 0 0 0 0 0 0\n", temperature.E, temperature.SE);
+ }
+ }
+}
+
+static void
+draw(struct sdis_scene* scn, struct sdis_camera* cam)
+{
+ struct sdis_solve_camera_args args = SDIS_SOLVE_CAMERA_ARGS_DEFAULT;
+ struct sdis_estimator_buffer* image = NULL;
+
+ args.cam = cam;
+ args.image_definition[0] = IMG_WIDTH;
+ args.image_definition[1] = IMG_HEIGHT;
+ args.spp = IMG_SPP;
+
+ OK(sdis_solve_camera(scn, &args, &image));
+ write_image(stdout, image);
+ OK(sdis_estimator_buffer_ref_put(image));
+}
+
+/*******************************************************************************
+ * The test
+ ******************************************************************************/
+int
+main(int argc, char** argv)
+{
+ /* Stardis */
+ struct sdis_camera* cam = NULL;
+ struct sdis_device* dev = NULL;
+ struct sdis_interface* interf = NULL;
+ struct sdis_medium* fluid = NULL;
+ struct sdis_medium* solid = NULL;
+ struct sdis_scene* scn = NULL;
+ struct sdis_source* source = NULL;
+
+ /* Miscellaneous */
+ struct mesh mesh = MESH_NULL;
+ (void)argc, (void)argv;
+
+ OK(sdis_device_create(&SDIS_DEVICE_CREATE_ARGS_DEFAULT, &dev));
+
+ mesh_init(&mesh);
+ mesh_add_super_shape(&mesh);
+ mesh_add_ground(&mesh);
+
+ solid = create_solid(dev);
+ fluid = create_fluid(dev);
+ interf = create_interface(dev, fluid, solid);
+ source = create_source(dev);
+ scn = create_scene(dev, &mesh, interf, source);
+ cam = create_camera(dev);
+
+ draw(scn, cam);
+
+ /* For debug of scene geometry */
+ /*dump_mesh(stdout,
+ mesh.positions, mesh_nvertices(&mesh),
+ mesh.indices, mesh_ntriangles(&mesh));*/
+
+ mesh_release(&mesh);
+ OK(sdis_camera_ref_put(cam));
+ OK(sdis_device_ref_put(dev));
+ OK(sdis_interface_ref_put(interf));
+ OK(sdis_medium_ref_put(fluid));
+ OK(sdis_medium_ref_put(solid));
+ OK(sdis_scene_ref_put(scn));
+ OK(sdis_source_ref_put(source));
+ CHK(mem_allocated_size() == 0);
+ return 0;
+}
diff --git a/src/test_sdis_external_flux.c b/src/test_sdis_external_flux.c
@@ -0,0 +1,527 @@
+/* Copyright (C) 2016-2023 |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 "test_sdis_utils.h"
+#include "sdis.h"
+
+#include <rsys/rsys.h>
+
+/*
+ * The system consists of 2 parallelepipeds: a vertical one called the wall, and
+ * a horizontal one representing the floor. The wall is a black body, while the
+ * floor is a perfectly reflective surface. The surrounding fluid has a fixed
+ * temperature and, finally, an external spherical source represents the sun.
+ * This test calculates the steady temperature at a position in the wall and
+ * compares it with the analytical solution given for a perfectly diffuse or
+ * specular ground.
+ *
+ * (-0.1,1500)
+ * +---+ External source
+ * | E=1 T_FLUID ##
+ * Probe x | _\ ####
+ * | | / / ##
+ * +---+ \__/
+ * (0,500)
+ *
+ * (0,0)
+ * Y *--------E=0------------- - - -
+ * | |
+ * o--X *------------------------ - - -
+ * / (0,-1)
+ * Z
+ *
+ */
+
+#define T_FLUID 300.0 /* [K] */
+#define T_REF 0 /* [K] */
+
+/*******************************************************************************
+ * Geometries
+ ******************************************************************************/
+static const double positions_2d[] = {
+ /* Ground */
+ 1.0e12, -1.0,
+ 0.0, -1.0,
+ 0.0, 0.0,
+ 1.0e12, 0.0,
+
+ /* Wall */
+ 0.0, 500.0,
+ -0.1, 500.0,
+ -0.1, 1500.0,
+ 0.0, 1500.0
+};
+static const size_t nvertices_2d = sizeof(positions_2d) / (sizeof(double)*2);
+
+static const double positions_3d[] = {
+ /* Ground */
+ 0.0, -1.0, -1.0e6,
+ 1.0e12, -1.0, -1.0e6,
+ 0.0, 1.0, -1.0e6,
+ 1.0e12, 1.0, -1.0e6,
+ 0.0, -1.0, 1.0e6,
+ 1.0e12, -1.0, 1.0e6,
+ 0.0, 1.0, 1.0e6,
+ 1.0e12, 1.0, 1.0e6,
+
+ /* Wall */
+ -0.1, 500.0, -500.0,
+ 0.0, 500.0, -500.0,
+ -0.1, 1500.0, -500.0,
+ 0.0, 1500.0, -500.0,
+ -0.1, 500.0, 500.0,
+ 0.0, 500.0, 500.0,
+ -0.1, 1500.0, 500.0,
+ 0.0, 1500.0, 500.0
+};
+static const size_t nvertices_3d = sizeof(positions_3d) / (sizeof(double)*3);
+
+static const size_t indices_2d[] = {
+ /* Ground */
+ 0, 1, /* -y */
+ 1, 2, /* -x */
+ 2, 3, /* +y */
+ 3, 0, /* +x */
+
+ /* Wall */
+ 4, 5, /* -y */
+ 5, 6, /* -x */
+ 6, 7, /* +y */
+ 7, 4 /* +x */
+};
+static const size_t nsegments = sizeof(indices_2d) / (sizeof(size_t)*2);
+
+static const size_t indices_3d[] = {
+ /* Ground */
+ 0, 2, 1, 1, 2, 3, /* -z */
+ 0, 4, 2, 2, 4, 6, /* -x */
+ 4, 5, 6, 6, 5, 7, /* +z */
+ 3, 7, 5, 5, 1, 3, /* +x */
+ 2, 6, 7, 7, 3, 2, /* +y */
+ 0, 1, 5, 5, 4, 0, /* -y */
+
+ /* Wall */
+ 8, 10, 9, 9, 10, 11, /* -z */
+ 8, 12, 10, 10, 12, 14, /* -x */
+ 12, 13, 14, 14, 13, 15, /* +z */
+ 11, 15, 13, 13, 9, 11, /* +x */
+ 10, 14, 15, 15, 11, 10, /* +y */
+ 8, 9, 13, 13, 12, 8 /* -y */
+};
+static const size_t ntriangles = sizeof(indices_3d) / (sizeof(size_t)*3);
+
+/*******************************************************************************
+ * Media
+ ******************************************************************************/
+#define MEDIUM_PROP(Type, Prop, Val) \
+ static double \
+ Type##_get_##Prop \
+ (const struct sdis_rwalk_vertex* vtx, \
+ struct sdis_data* data) \
+ { \
+ (void)vtx, (void)data; /* Avoid the "unused variable" warning */ \
+ return Val; \
+ }
+MEDIUM_PROP(medium, volumic_mass, 1700) /* [kj/m^3] */
+MEDIUM_PROP(medium, calorific_capacity, 800) /* [J/K/Kg] */
+MEDIUM_PROP(solid, thermal_conductivity, 1.15) /* [W/m/K] */
+MEDIUM_PROP(solid, delta, 0.1/20.0) /* [m] */
+MEDIUM_PROP(solid, temperature, -1/*<=> unknown*/) /* [K] */
+MEDIUM_PROP(fluid, temperature, T_FLUID) /* [K] */
+#undef MEDIUM_PROP
+
+static struct sdis_medium*
+create_solid(struct sdis_device* sdis)
+{
+ struct sdis_solid_shader shader = SDIS_SOLID_SHADER_NULL;
+ struct sdis_medium* solid = NULL;
+
+ shader.calorific_capacity = medium_get_calorific_capacity;
+ shader.thermal_conductivity = solid_get_thermal_conductivity;
+ shader.volumic_mass = medium_get_volumic_mass;
+ shader.delta = solid_get_delta;
+ shader.temperature = solid_get_temperature;
+ OK(sdis_solid_create(sdis, &shader, NULL, &solid));
+ return solid;
+}
+
+static struct sdis_medium*
+create_fluid(struct sdis_device* sdis)
+{
+ struct sdis_fluid_shader shader = SDIS_FLUID_SHADER_NULL;
+ struct sdis_medium* fluid = NULL;
+
+ shader.calorific_capacity = medium_get_calorific_capacity;
+ shader.volumic_mass = medium_get_volumic_mass;
+ shader.temperature = fluid_get_temperature;
+ OK(sdis_fluid_create(sdis, &shader, NULL, &fluid));
+ return fluid;
+}
+
+/*******************************************************************************
+ * Interfaces
+ ******************************************************************************/
+struct interface {
+ double emissivity;
+ double specular_fraction;
+ double convection_coef; /* [W/m^2/K] */
+};
+
+#define INTERF_PROP(Prop, Val) \
+ static double \
+ interface_get_##Prop \
+ (const struct sdis_interface_fragment* frag, \
+ struct sdis_data* data) \
+ { \
+ (void)frag, (void)data; /* Avoid the "unused variable" warning */ \
+ return Val; \
+ }
+INTERF_PROP(temperature, -1/*<=> unknown*/) /* [K] */
+INTERF_PROP(reference_temperature, T_REF) /* [K] */
+#undef INTERF_PROP
+
+#define INTERF_PROP(Prop) \
+ static double \
+ interface_get_##Prop \
+ (const struct sdis_interface_fragment* frag, \
+ struct sdis_data* data) \
+ { \
+ struct interface* interf_data = NULL; \
+ (void)frag; /* Avoid the "unused variable" warning */ \
+ interf_data = sdis_data_get(data); \
+ return interf_data->Prop; \
+ }
+INTERF_PROP(emissivity)
+INTERF_PROP(specular_fraction)
+INTERF_PROP(convection_coef) /* [W/m^2/K] */
+#undef INTERF_PROP
+
+static struct sdis_interface*
+create_interface
+ (struct sdis_device* sdis,
+ struct sdis_medium* front,
+ struct sdis_medium* back,
+ const double emissivity,
+ const double convection_coef,
+ struct interface** out_interf_data) /* May be NULL */
+{
+ struct sdis_data* data = NULL;
+ struct sdis_interface* interf = NULL;
+ struct sdis_interface_shader shader = SDIS_INTERFACE_SHADER_NULL;
+ struct interface* interf_data = NULL;
+
+ OK(sdis_data_create
+ (sdis, sizeof(struct interface), ALIGNOF(struct interface), NULL, &data));
+ interf_data = sdis_data_get(data);
+ interf_data->emissivity = emissivity;
+ interf_data->convection_coef = convection_coef; /* [W/m^2/K] */
+ interf_data->specular_fraction = 0; /* Diffuse */
+ if(out_interf_data) *out_interf_data = interf_data;
+
+ shader.front.temperature = interface_get_temperature;
+ shader.back.temperature = interface_get_temperature;
+ shader.back.reference_temperature = interface_get_reference_temperature;
+ shader.back.emissivity = interface_get_emissivity;
+ shader.back.specular_fraction = interface_get_specular_fraction;
+ shader.convection_coef = interface_get_convection_coef;
+ shader.convection_coef_upper_bound = convection_coef;
+ OK(sdis_interface_create(sdis, front, back, &shader, data, &interf));
+ OK(sdis_data_ref_put(data));
+
+ return interf;
+}
+
+/*******************************************************************************
+ * External source
+ ******************************************************************************/
+static void
+source_get_position
+ (const double time,
+ double pos[3],
+ struct sdis_data* data)
+{
+ const double elevation = MDEG2RAD(30); /* [radian] */
+ const double distance = 1.5e11; /* [m] */
+ (void)time, (void)data; /* Avoid the "unusued variable" warning */
+
+ pos[0] = cos(elevation) * distance;
+ pos[1] = sin(elevation) * distance;
+ pos[2] = 0;
+}
+
+static struct sdis_source*
+create_source(struct sdis_device* sdis)
+{
+ struct sdis_spherical_source_create_args args = SDIS_SPHERICAL_SOURCE_CREATE_ARGS_NULL;
+ struct sdis_source* src = NULL;
+
+ args.position = source_get_position;
+ args.data = NULL;
+ args.radius = 6.5991756e8; /* [m] */
+ args.power = 3.845e26; /* [W] */
+ OK(sdis_spherical_source_create(sdis, &args, &src));
+ return src;
+}
+
+/*******************************************************************************
+ * Scene
+ ******************************************************************************/
+struct scene_context {
+ struct sdis_interface* interf_ground;
+ struct sdis_interface* interf_wall;
+};
+static const struct scene_context SCENE_CONTEXT_NULL = {NULL, NULL};
+
+static void
+scene_get_indices_2d(const size_t iseg, size_t ids[2], void* ctx)
+{
+ struct scene_context* context = ctx;
+ CHK(ids && context && iseg < nsegments);
+ ids[0] = (unsigned)indices_2d[iseg*2+0];
+ ids[1] = (unsigned)indices_2d[iseg*2+1];
+}
+
+static void
+scene_get_indices_3d(const size_t itri, size_t ids[3], void* ctx)
+{
+ struct scene_context* context = ctx;
+ CHK(ids && context && itri < ntriangles);
+ ids[0] = (unsigned)indices_3d[itri*3+0];
+ ids[1] = (unsigned)indices_3d[itri*3+1];
+ ids[2] = (unsigned)indices_3d[itri*3+2];
+}
+
+static void
+scene_get_interface_2d(const size_t iseg, struct sdis_interface** interf, void* ctx)
+{
+ struct scene_context* context = ctx;
+ CHK(interf && context && iseg < nsegments);
+ *interf = iseg < 4 ? context->interf_ground : context->interf_wall;
+}
+
+static void
+scene_get_interface_3d(const size_t itri, struct sdis_interface** interf, void* ctx)
+{
+ struct scene_context* context = ctx;
+ CHK(interf && context && itri < ntriangles);
+ *interf = itri < 12 ? context->interf_ground : context->interf_wall;
+}
+
+static void
+scene_get_position_2d(const size_t ivert, double pos[2], void* ctx)
+{
+ struct scene_context* context = ctx;
+ CHK(pos && context && ivert < nvertices_2d);
+ pos[0] = positions_2d[ivert*2+0];
+ pos[1] = positions_2d[ivert*2+1];
+}
+
+static void
+scene_get_position_3d(const size_t ivert, double pos[3], void* ctx)
+{
+ struct scene_context* context = ctx;
+ CHK(pos && context && ivert < nvertices_3d);
+ pos[0] = positions_3d[ivert*3+0];
+ pos[1] = positions_3d[ivert*3+1];
+ pos[2] = positions_3d[ivert*3+2];
+}
+
+static struct sdis_scene*
+create_scene_2d
+ (struct sdis_device* sdis,
+ struct sdis_interface* interf_ground,
+ struct sdis_interface* interf_wall,
+ struct sdis_source* source)
+{
+ struct sdis_scene* scn = NULL;
+ struct sdis_source* src = NULL;
+ struct sdis_scene_create_args scn_args = SDIS_SCENE_CREATE_ARGS_DEFAULT;
+ struct scene_context context = SCENE_CONTEXT_NULL;
+
+ context.interf_ground = interf_ground;
+ context.interf_wall = interf_wall;
+
+ scn_args.get_indices = scene_get_indices_2d;
+ scn_args.get_interface = scene_get_interface_2d;
+ scn_args.get_position = scene_get_position_2d;
+ scn_args.nprimitives = nsegments;
+ scn_args.nvertices = nvertices_2d;
+ scn_args.trad.temperature = 0; /* [K] */
+ scn_args.trad.reference = T_REF; /* [K] */
+ scn_args.t_range[0] = 0; /* [K] */
+ scn_args.t_range[1] = 0; /* [K] */
+ scn_args.source = source;
+ scn_args.context = &context;
+ OK(sdis_scene_2d_create(sdis, &scn_args, &scn));
+
+ BA(sdis_scene_get_source(NULL, &src));
+ BA(sdis_scene_get_source(scn, NULL));
+ OK(sdis_scene_get_source(scn, &src));
+ CHK(src == source);
+
+ return scn;
+}
+
+static struct sdis_scene*
+create_scene_3d
+ (struct sdis_device* sdis,
+ struct sdis_interface* interf_ground,
+ struct sdis_interface* interf_wall,
+ struct sdis_source* source)
+{
+ struct sdis_scene* scn = NULL;
+ struct sdis_source* src = NULL;
+ struct sdis_scene_create_args scn_args = SDIS_SCENE_CREATE_ARGS_DEFAULT;
+ struct scene_context context = SCENE_CONTEXT_NULL;
+
+ context.interf_ground = interf_ground;
+ context.interf_wall = interf_wall;
+
+ scn_args.get_indices = scene_get_indices_3d;
+ scn_args.get_interface = scene_get_interface_3d;
+ scn_args.get_position = scene_get_position_3d;
+ scn_args.nprimitives = ntriangles;
+ scn_args.nvertices = nvertices_3d;
+ scn_args.trad.temperature = 0; /* [K] */
+ scn_args.trad.reference = T_REF; /* [K] */
+ scn_args.t_range[0] = 0; /* [K] */
+ scn_args.t_range[1] = 0; /* [K] */
+ scn_args.source = source;
+ scn_args.context = &context;
+ OK(sdis_scene_create(sdis, &scn_args, &scn));
+
+ BA(sdis_scene_get_source(NULL, &src));
+ BA(sdis_scene_get_source(scn, NULL));
+ OK(sdis_scene_get_source(scn, &src));
+ CHK(src == source);
+
+ return scn;
+}
+
+/*******************************************************************************
+ * Validations
+ ******************************************************************************/
+static void
+check
+ (struct sdis_scene* scn,
+ const size_t nrealisations,
+ const double analytical_ref,
+ const int is_master_process)
+{
+ struct sdis_solve_probe_args probe_args = SDIS_SOLVE_PROBE_ARGS_DEFAULT;
+ struct sdis_mc T = SDIS_MC_NULL;
+ struct sdis_estimator* estimator = NULL;
+
+ probe_args.position[0] = -0.05;
+ probe_args.position[1] = 1000;
+ probe_args.position[2] = 0;
+ probe_args.nrealisations = nrealisations;
+ OK(sdis_solve_probe(scn, &probe_args, &estimator));
+
+ if(!is_master_process) return;
+
+ OK(sdis_estimator_get_temperature(estimator, &T));
+ printf("T(%g, %g, %g) = %g ~ %g +/- %g\n",
+ SPLIT3(probe_args.position), analytical_ref, T.E, T.SE);
+ OK(sdis_estimator_ref_put(estimator));
+
+ CHK(eq_eps(analytical_ref, T.E, 3*T.SE));
+}
+
+static void
+check_green
+ (struct sdis_scene* scn,
+ const size_t nrealisations,
+ const double analytical_ref,
+ const int is_master_process)
+{
+ struct sdis_solve_probe_args probe_args = SDIS_SOLVE_PROBE_ARGS_DEFAULT;
+ struct sdis_mc T = SDIS_MC_NULL;
+ struct sdis_green_function* green = NULL;
+ struct sdis_estimator* estimator = NULL;
+
+ probe_args.position[0] = -0.05;
+ probe_args.position[1] = 1000;
+ probe_args.position[2] = 0;
+ probe_args.nrealisations = nrealisations;
+ OK(sdis_solve_probe_green_function(scn, &probe_args, &green));
+
+ if(!is_master_process) return;
+
+ OK(sdis_green_function_solve(green, &estimator));
+ check_green_function(green);
+
+ OK(sdis_estimator_get_temperature(estimator, &T));
+
+ printf("T(%g, %g, %g) = %g ~ %g +/- %g\n",
+ SPLIT3(probe_args.position), analytical_ref, T.E, T.SE);
+ OK(sdis_green_function_ref_put(green));
+ OK(sdis_estimator_ref_put(estimator));
+
+ CHK(eq_eps(analytical_ref, T.E, 3*T.SE));
+}
+
+/*******************************************************************************
+ * The test
+ ******************************************************************************/
+int
+main(int argc, char** argv)
+{
+ struct sdis_device* dev = NULL;
+ struct sdis_medium* fluid = NULL;
+ struct sdis_medium* solid = NULL;
+ struct sdis_interface* interf_ground = NULL;
+ struct sdis_interface* interf_wall = NULL;
+ struct sdis_source* src = NULL;
+ struct sdis_scene* scn_2d = NULL;
+ struct sdis_scene* scn_3d = NULL;
+
+ struct interface* ground_interf_data = NULL;
+ int is_master_process = 0;
+ (void) argc, (void)argv;
+
+ create_default_device(&argc, &argv, &is_master_process, &dev);
+
+ fluid = create_fluid(dev);
+ solid = create_solid(dev);
+ interf_ground = create_interface
+ (dev, solid, fluid, 0/*emissivity*/, 0/*h*/, &ground_interf_data);
+ interf_wall = create_interface
+ (dev, solid, fluid, 1/*emissivity*/, 10/*h*/, NULL);
+ src = create_source(dev);
+ scn_2d = create_scene_2d(dev, interf_ground, interf_wall, src);
+ scn_3d = create_scene_3d(dev, interf_ground, interf_wall, src);
+
+ ground_interf_data->specular_fraction = 0; /* Lambertian */
+ check(scn_2d, 10000, 375.88, is_master_process);
+ check_green(scn_3d, 10000, 375.88, is_master_process);
+
+ ground_interf_data->specular_fraction = 1; /* Specular */
+ check(scn_2d, 100000, 417.77, is_master_process);
+ check_green(scn_3d, 100000, 417.77, is_master_process);
+
+ OK(sdis_medium_ref_put(fluid));
+ OK(sdis_medium_ref_put(solid));
+ OK(sdis_interface_ref_put(interf_ground));
+ OK(sdis_interface_ref_put(interf_wall));
+ OK(sdis_source_ref_put(src));
+ OK(sdis_scene_ref_put(scn_2d));
+ OK(sdis_scene_ref_put(scn_3d));
+
+ free_default_device(dev);
+
+ CHK(mem_allocated_size() == 0);
+ return 0;
+}
diff --git a/src/test_sdis_solve_probe_list.c b/src/test_sdis_solve_probe_list.c
@@ -36,12 +36,12 @@
* trilinear profile. Thus, we should find by Monte Carlo the temperature
* defined by the trilinear profile.
*
- * /\ <-- T(x,y,z)
- * T(z) ___/ \___
- * | \ . T=? /
- * o-- T(x) /_ __ _\
- * / \/ \/
- * T(y)
+ * T(z) /\ <-- T(x,y,z)
+ * | T(y) ___/ \___
+ * |/ \ . T=? /
+ * o--- T(x) /_ __ _\
+ * \/ \/
+ *
*/
/*******************************************************************************
@@ -118,9 +118,9 @@ create_super_shape(void)
}
/*******************************************************************************
- * View, i.e. acceleration structure used to query geometry. In this test it is
- * used to calculate the delta parameter and to sample the probe positions in
- * the supershape.
+ * View, i.e. acceleration structure used to query geometry. In this test it is
+ * used to calculate the delta parameter and to sample the probe positions in
+ * the supershape.
******************************************************************************/
static void
view_get_indices(const unsigned itri, unsigned ids[3], void* ctx)
diff --git a/src/test_sdis_source.c b/src/test_sdis_source.c
@@ -0,0 +1,88 @@
+/* Copyright (C) 2016-2023 |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 "test_sdis_utils.h"
+
+/*******************************************************************************
+ * Helper functions
+ ******************************************************************************/
+static void
+spherical_source_get_position
+ (const double time,
+ double pos[3],
+ struct sdis_data* data)
+{
+ (void)time, (void)data;
+ pos[0] = pos[1] = pos[2] = 1.234;
+}
+
+static void
+check_spherical_source(struct sdis_device* dev)
+{
+ struct sdis_spherical_source_create_args args =
+ SDIS_SPHERICAL_SOURCE_CREATE_ARGS_NULL;
+ struct sdis_source* src = NULL;
+ struct sdis_data* data = NULL;
+ double power = 0;
+
+ /* Create a data to check its memory management */
+ OK(sdis_data_create(dev, sizeof(double[3]), ALIGNOF(double[3]), NULL, &data));
+
+ args.position = spherical_source_get_position;
+ args.data = data;
+ args.radius = 1;
+ args.power = 10;
+
+ BA(sdis_spherical_source_create(NULL, &args, &src));
+ BA(sdis_spherical_source_create(dev, NULL, &src));
+ BA(sdis_spherical_source_create(dev, &args, NULL));
+ OK(sdis_spherical_source_create(dev, &args, &src));
+
+ BA(sdis_source_get_power(NULL, &power));
+ BA(sdis_source_get_power(src, NULL));
+ OK(sdis_source_get_power(src, &power));
+ CHK(power == args.power);
+
+ BA(sdis_source_ref_get(NULL));
+ OK(sdis_source_ref_get(src));
+ BA(sdis_source_ref_put(NULL));
+ OK(sdis_source_ref_put(src));
+ OK(sdis_source_ref_put(src));
+
+ OK(sdis_data_ref_put(data));
+
+ args.data = NULL;
+ OK(sdis_spherical_source_create(dev, &args, &src));
+ OK(sdis_source_ref_put(src));
+}
+
+/*******************************************************************************
+ * The test
+ ******************************************************************************/
+int
+main(int argc, char** argv)
+{
+ struct sdis_device* dev = NULL;
+ (void)argc, (void)argv;
+
+ OK(sdis_device_create(&SDIS_DEVICE_CREATE_ARGS_DEFAULT, &dev));
+
+ check_spherical_source(dev);
+
+ OK(sdis_device_ref_put(dev));
+ CHK(mem_allocated_size() == 0);
+ return 0;
+}
diff --git a/src/test_sdis_utils.c b/src/test_sdis_utils.c
@@ -92,18 +92,33 @@ solve_green_path(struct sdis_green_path* path, void* ctx)
struct sdis_solid_shader solid = SDIS_SOLID_SHADER_NULL;
struct sdis_fluid_shader fluid = SDIS_FLUID_SHADER_NULL;
struct sdis_interface_shader interf = SDIS_INTERFACE_SHADER_NULL;
+ struct sdis_green_function* green = NULL;
+ struct sdis_scene* scn = NULL;
+ struct sdis_source* source = NULL;
struct green_accum* acc = NULL;
struct sdis_data* data = NULL;
enum sdis_medium_type type;
enum sdis_green_path_end_type end_type;
+ double source_power = 0; /* [W] */
double power = 0;
double flux = 0;
+ double external_flux = 0; /* [W/m^2] */
double time, temp = 0;
double weight = 0;
CHK(path && ctx);
acc = ctx;
+ BA(sdis_green_path_get_green_function(NULL, NULL));
+ BA(sdis_green_path_get_green_function(path, NULL));
+ BA(sdis_green_path_get_green_function(NULL, &green));
+ OK(sdis_green_path_get_green_function(path, &green));
+
+ BA(sdis_green_function_get_scene(NULL, NULL));
+ BA(sdis_green_function_get_scene(NULL, &scn));
+ BA(sdis_green_function_get_scene(green, NULL));
+ OK(sdis_green_function_get_scene(green, &scn));
+
BA(sdis_green_path_for_each_power_term(NULL, accum_power_terms, &power));
BA(sdis_green_path_for_each_power_term(path, NULL, &acc));
OK(sdis_green_path_for_each_power_term(path, accum_power_terms, &power));
@@ -117,6 +132,17 @@ solve_green_path(struct sdis_green_path* path, void* ctx)
BA(sdis_green_path_get_elapsed_time(NULL, &time));
OK(sdis_green_path_get_elapsed_time(path, &time));
+ BA(sdis_green_path_get_external_flux_term(NULL, &external_flux));
+ BA(sdis_green_path_get_external_flux_term(path, NULL));
+ OK(sdis_green_path_get_external_flux_term(path, &external_flux));
+ OK(sdis_scene_get_source(scn, &source));
+ if(source == NULL) {
+ CHK(external_flux == 0);
+ } else {
+ OK(sdis_source_get_power(source, &source_power));
+ external_flux *= source_power;
+ }
+
BA(sdis_green_path_get_end_type(NULL, NULL));
BA(sdis_green_path_get_end_type(path, NULL));
BA(sdis_green_path_get_end_type(NULL, &end_type));
@@ -127,18 +153,7 @@ solve_green_path(struct sdis_green_path* path, void* ctx)
BA(sdis_green_path_get_limit_point(path, NULL));
if(end_type == SDIS_GREEN_PATH_END_RADIATIVE) {
struct sdis_ambient_radiative_temperature trad;
- struct sdis_green_function* green;
- struct sdis_scene* scn;
BO(sdis_green_path_get_limit_point(path, &pt));
- BA(sdis_green_path_get_green_function(NULL, NULL));
- BA(sdis_green_path_get_green_function(path, NULL));
- BA(sdis_green_path_get_green_function(NULL, &green));
- OK(sdis_green_path_get_green_function(path, &green));
-
- BA(sdis_green_function_get_scene(NULL, NULL));
- BA(sdis_green_function_get_scene(NULL, &scn));
- BA(sdis_green_function_get_scene(green, NULL));
- OK(sdis_green_function_get_scene(green, &scn));
BA(sdis_scene_get_ambient_radiative_temperature(NULL, NULL));
BA(sdis_scene_get_ambient_radiative_temperature(scn, NULL));
@@ -172,7 +187,7 @@ solve_green_path(struct sdis_green_path* path, void* ctx)
}
}
- weight = temp + power + flux;
+ weight = temp + power + external_flux + flux;
acc->sum += weight;
acc->sum2 += weight*weight;
diff --git a/src/test_sdis_utils.h b/src/test_sdis_utils.h
@@ -192,7 +192,8 @@ static const struct sdis_fluid_shader DUMMY_FLUID_SHADER = {
dummy_interface_getter, /* Flux */ \
dummy_interface_getter, /* Emissivity */ \
dummy_interface_getter, /* Specular fraction */ \
- dummy_interface_getter /* Reference temperature */ \
+ dummy_interface_getter, /* Reference temperature */ \
+ 1 /* Handle external flux */ \
}
static const struct sdis_interface_shader DUMMY_INTERFACE_SHADER = {
dummy_interface_getter, /* Convection coef */
