commit cae30778d715efd903e10c8045b8a099bf1ea9d2
parent 1e912557c634f88db4c0190cec6434caa6764dac
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Thu, 14 Mar 2024 12:05:29 +0100
Merge remote-tracking branch 'origin/develop' into feature_wos
Diffstat:
27 files changed, 1475 insertions(+), 388 deletions(-)
diff --git a/Makefile b/Makefile
@@ -213,7 +213,8 @@ TEST_SRC_MPI =\
src/test_sdis_solve_boundary.c\
src/test_sdis_solve_boundary_flux.c\
src/test_sdis_solve_probe2.c\
- src/test_sdis_solve_probe_list.c
+ src/test_sdis_solve_probe_list.c\
+ src/test_sdis_solve_probe_boundary_list.c
TEST_OBJ =\
$(TEST_SRC:.c=.o)\
$(TEST_SRC_MPI:.c=.o)\
@@ -431,18 +432,21 @@ test_sdis_solve_medium_2d \
src/test_sdis_compute_power.d \
src/test_sdis_solve_camera.d \
src/test_sdis_solve_medium.d \
+src/test_sdis_solve_probe_boundary_list.d \
: config.mk sdis-local.pc
@$(CC) $(TEST_CFLAGS_MPI) $(S3DUT_CFLAGS) -MM -MT "$(@:.d=.o) $@" $(@:.d=.c) -MF $@
src/test_sdis_compute_power.o \
src/test_sdis_solve_camera.o \
src/test_sdis_solve_medium.o \
+src/test_sdis_solve_probe_boundary_list.o \
: config.mk sdis-local.pc
$(CC) $(TEST_CFLAGS_MPI) $(S3DUT_CFLAGS) -c $(@:.o=.c) -o $@
test_sdis_compute_power \
test_sdis_solve_camera \
test_sdis_solve_medium \
+test_sdis_solve_probe_boundary_list \
: config.mk sdis-local.pc $(LIBNAME) src/test_sdis_utils.o
$(CC) $(TEST_CFLAGS_MPI) $(S3DUT_CFLAGS) -o $@ src/$@.o $(TEST_LIBS_MPI) $(S3DUT_LIBS)
diff --git a/sdis.pc.in b/sdis.pc.in
@@ -2,10 +2,12 @@ prefix=@PREFIX@
includedir=${prefix}/include
libdir=${prefix}/lib
-Requires: rsys >= @RSYS_VERSION@, star-sp >= @SSP_VERSION@
-Requires.private:\
+Requires: \
+ rsys >= @RSYS_VERSION@,\
s2d >= @S2D_VERSION@,\
s3d >= @S3D_VERSION@,\
+ star-sp >= @SSP_VERSION@
+Requires.private:\
senc2d >= @SENC3D_VERSION@,\
senc3d >= @SENC3D_VERSION@ @MPI@
Name: sdis
diff --git a/src/sdis.c b/src/sdis.c
@@ -476,10 +476,10 @@ compute_process_index_range
per_process_indices = nindices / (size_t)dev->mpi_nprocs;
range[0] = per_process_indices * (size_t)dev->mpi_rank;
- range[1] = range[0] + per_process_indices; /* Upper bound is _exclusive */
+ range[1] = range[0] + per_process_indices; /* Upper bound is _exclusive_ */
ASSERT(range[0] <= range[1]);
- /* Set the remaining number of indexes that are not managed by one process */
+ /* Set the remaining number of indices that are not managed by one process */
remaining_indices =
nindices - per_process_indices * (size_t)dev->mpi_nprocs;
diff --git a/src/sdis.h b/src/sdis.h
@@ -16,6 +16,8 @@
#ifndef SDIS_H
#define SDIS_H
+#include <star/s2d.h>
+#include <star/s3d.h>
#include <star/ssp.h>
#include <rsys/hash.h>
@@ -147,25 +149,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. */
+/* 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);
+/* Type of functor used to retrieve the source's power relative to time */
+typedef double
+(*sdis_get_power_T)
+ (const double time,
+ 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] */
+ sdis_get_power_T power; /* Total power [W] */
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__;
+struct sdis_scene_find_closest_point_args {
+ double position[3]; /* Query position */
+ double radius; /* Maxium search distance around pos */
+
+ /* User defined filter function */
+ s2d_hit_filter_function_T filter_2d;
+ s3d_hit_filter_function_T filter_3d;
+ void* filter_data; /* Filter function data */
+};
+#define SDIS_SCENE_FIND_CLOSEST_POINT_ARGS_NULL__ {{0,0,0}, 0, NULL, NULL, NULL}
+static const struct sdis_scene_find_closest_point_args
+SDIS_SCENE_FIND_CLOSEST_POINT_ARGS_NULL = SDIS_SCENE_FIND_CLOSEST_POINT_ARGS_NULL__;
+
/*******************************************************************************
* Estimation data types
******************************************************************************/
@@ -566,6 +587,27 @@ static const struct sdis_solve_probe_boundary_args
SDIS_SOLVE_PROBE_BOUNDARY_ARGS_DEFAULT =
SDIS_SOLVE_PROBE_BOUNDARY_ARGS_DEFAULT__;
+/* Input arguments of the solve function that distributes the calculations of
+ * several boundary probes rather than the realizations of a probe */
+struct sdis_solve_probe_boundary_list_args {
+ struct sdis_solve_probe_boundary_args* probes; /* List of probes to compute */
+ size_t nprobes; /* Total number of probes */
+
+ /* State/type of the RNG to use for the list of probes to calculate.
+ * The state/type defines per probe is ignored */
+ struct ssp_rng* rng_state; /* Initial RNG state. May be NULL */
+ enum ssp_rng_type rng_type; /* RNG type to use if `rng_state' is NULL */
+};
+#define SDIS_SOLVE_PROBE_BOUNDARY_LIST_ARGS_DEFAULT__ { \
+ NULL, /* List of probes */ \
+ 0, /* #probes */ \
+ NULL, /* RNG state */ \
+ SSP_RNG_THREEFRY /* RNG type */ \
+}
+static const struct sdis_solve_probe_boundary_list_args
+SDIS_SOLVE_PROBE_BOUNDARY_LIST_ARGS_DEFAULT =
+ SDIS_SOLVE_PROBE_BOUNDARY_LIST_ARGS_DEFAULT__;
+
struct sdis_solve_boundary_args {
size_t nrealisations; /* #realisations */
const size_t* primitives; /* List of boundary primitives to handle */
@@ -1006,10 +1048,10 @@ SDIS_API res_T
sdis_source_ref_put
(struct sdis_source* source);
-SDIS_API res_T
+SDIS_API double
sdis_source_get_power
- (const struct sdis_source* src,
- double* power);/* [W] */
+ (struct sdis_source* source,
+ const double time); /* [s] */
/*******************************************************************************
* A scene is a collection of primitives. Each primitive is the geometric
@@ -1119,8 +1161,7 @@ sdis_scene_set_temperature_range
SDIS_API res_T
sdis_scene_find_closest_point
(const struct sdis_scene* scn,
- const double pos[], /* Query position */
- const double radius, /* Maximum search distance around pos */
+ const struct sdis_scene_find_closest_point_args* args,
size_t* iprim, /* Primitive index onto which the closest point lies */
double uv[]); /* Parametric cordinate onto the primitive */
@@ -1452,18 +1493,6 @@ sdis_solve_probe
const struct sdis_solve_probe_args* args,
struct sdis_estimator** estimator);
-/* Calculate temperature for a list of probe points. Unlike its
- * single-probe counterpart, this function parallelizes the list of
- * probes, rather than calculating a single probe. Calling this function
- * is therefore more advantageous in terms of load distribution when the
- * number of probe points to be evaluated is large compared to the cost
- * of calculating a single probe point. */
-SDIS_API res_T
-sdis_solve_probe_list
- (struct sdis_scene* scn,
- const struct sdis_solve_probe_list_args* args,
- struct sdis_estimator_buffer** buf);
-
SDIS_API res_T
sdis_solve_probe_boundary
(struct sdis_scene* scn,
@@ -1514,6 +1543,27 @@ sdis_compute_power
struct sdis_estimator** estimator);
/*******************************************************************************
+ * Solvers of a list of probes
+ *
+ * Unlike their single-probe counterpart, this function parallelizes the list of
+ * probes, rather than calculating a single probe. Calling these functions is
+ * therefore more advantageous in terms of load distribution when the number of
+ * probes to be evaluated is large compared to the cost of calculating a single
+ * probe.
+ ******************************************************************************/
+SDIS_API res_T
+sdis_solve_probe_list
+ (struct sdis_scene* scn,
+ const struct sdis_solve_probe_list_args* args,
+ struct sdis_estimator_buffer** buf);
+
+SDIS_API res_T
+sdis_solve_probe_boundary_list
+ (struct sdis_scene* scn,
+ const struct sdis_solve_probe_boundary_list_args* args,
+ struct sdis_estimator_buffer** buf);
+
+/*******************************************************************************
* Green solvers.
*
* Note that only the interfaces/media with flux/volumic power defined during
@@ -1524,7 +1574,7 @@ sdis_compute_power
*
* Also note that the green solvers assume that the interface fluxes are
* constant in time and space. The same applies to the volumic power of the
- * solid media.
+ * solid media and the power of external sources.
*
* If these assumptions are not ensured by the caller, the behavior of the
* estimated green function is undefined.
diff --git a/src/sdis_Xd_begin.h b/src/sdis_Xd_begin.h
@@ -42,6 +42,7 @@
#define FORMAT_VECX "%g, %g, %g"
#define SPLITX(V) SPLIT3(V)
+
#else
#error "Invalid dimension."
#endif
@@ -61,13 +62,16 @@
#define SXD_FLOAT2 CONCAT(CONCAT(S, DIM), D_FLOAT2)
#define SXD_FLOAT3 CONCAT(CONCAT(S, DIM), D_FLOAT3)
#define SXD_FLOATX CONCAT(CONCAT(CONCAT(S,DIM), D_FLOAT), DIM)
+#define SXD_GET_PRIMITIVE CONCAT(CONCAT(S, DIM), D_GET_PRIMITIVE)
#define SXD_SAMPLE CONCAT(CONCAT(S, DIM), D_SAMPLE)
+#define SXD_TRACE CONCAT(CONCAT(S, DIM), D_TRACE)
#define SXD_PRIMITIVE_EQ CONCAT(CONCAT(S, DIM), D_PRIMITIVE_EQ)
/* Vector macros generic to SDIS_XD_DIMENSION */
#define dX(Func) CONCAT(CONCAT(CONCAT(d, DIM), _), Func)
#define fX(Func) CONCAT(CONCAT(CONCAT(f, DIM), _), Func)
#define fX_set_dX CONCAT(CONCAT(CONCAT(f, DIM), _set_d), DIM)
+#define fXX_mulfX CONCAT(CONCAT(CONCAT(CONCAT(f, DIM), DIM), _mulf), DIM)
#define dX_set_fX CONCAT(CONCAT(CONCAT(d, DIM), _set_f), DIM)
/* Macro making generic its submitted name to SDIS_XD_DIMENSION */
@@ -82,6 +86,8 @@
#define SDIS_3D_H
#endif
+struct rwalk_context;
+
/* Current state of the random walk */
struct XD(rwalk) {
struct sdis_rwalk_vertex vtx; /* Position and time of the Random walk */
diff --git a/src/sdis_Xd_end.h b/src/sdis_Xd_end.h
@@ -31,11 +31,14 @@
#undef SXD_FLOAT2
#undef SXD_FLOAT3
#undef SXD_FLOATX
+#undef SXD_GET_PRIMITIVE
#undef SXD_SAMPLE
+#undef SXD_TRACE
#undef dX
#undef fX
#undef fX_set_dX
+#undef fXX_mulfX
#undef dX_set_fX
#undef FORMAT_VECX
diff --git a/src/sdis_estimator_buffer.c b/src/sdis_estimator_buffer.c
@@ -273,3 +273,8 @@ estimator_buffer_save_rng_state
return create_rng_from_rng_proxy(buf->dev, proxy, &buf->rng);
}
+/* Define the functions generic to the observable type */
+#define SDIS_X_OBS probe
+#include "sdis_estimator_buffer_X_obs.h"
+#define SDIS_X_OBS probe_boundary
+#include "sdis_estimator_buffer_X_obs.h"
diff --git a/src/sdis_estimator_buffer_X_obs.h b/src/sdis_estimator_buffer_X_obs.h
@@ -0,0 +1,118 @@
+/* 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/>. */
+
+/* Define functions on estimator buffer generic to the observable type. */
+
+#ifndef SDIS_ESTIMATOR_BUFFER_X_OBS_H
+#define SDIS_ESTIMATOR_BUFFER_X_OBS_H
+
+#include "sdis.h"
+#include "sdis_estimator_c.h"
+#include "sdis_estimator_buffer_c.h"
+#include "sdis_misc.h"
+
+#endif /* SDIS_ESTIMATOR_BUFFER_X_OBS_H */
+
+/* Check the generic observable parameter */
+#ifndef SDIS_X_OBS
+ #error "The SDIS_X_OBS macro must be defined."
+#endif
+
+#define X_OBS(Name) CONCAT(CONCAT(Name, _), SDIS_X_OBS)
+#define SDIS_SOLVE_X_OBS_ARGS CONCAT(CONCAT(sdis_solve_, SDIS_X_OBS),_args)
+
+res_T
+X_OBS(estimator_buffer_create_from_observable_list)
+ (struct sdis_device* dev,
+ struct ssp_rng_proxy* rng_proxy,
+ const struct SDIS_SOLVE_X_OBS_ARGS obs_list_args[],
+ const struct accum* per_obs_acc_temp,
+ const struct accum* per_obs_acc_time,
+ const size_t nobs, /* #observables */
+ struct sdis_estimator_buffer** out_estim_buffer)
+{
+ /* Accumulators throughout the buffer */
+ struct accum acc_temp = ACCUM_NULL;
+ struct accum acc_time = ACCUM_NULL;
+ size_t nrealisations = 0;
+
+ struct sdis_estimator_buffer* estim_buf = NULL;
+ size_t iobs = 0;
+ res_T res = RES_OK;
+
+ ASSERT(dev && rng_proxy);
+ ASSERT(obs_list_args || !nobs);
+ ASSERT(per_obs_acc_time && per_obs_acc_time && out_estim_buffer);
+
+ res = estimator_buffer_create(dev, nobs, 1, &estim_buf);
+ if(res != RES_OK) {
+ log_err(dev, "Unable to allocate the estimator buffer.\n");
+ goto error;
+ }
+
+ FOR_EACH(iobs, 0, nobs) {
+ const struct SDIS_SOLVE_X_OBS_ARGS* obs_args = NULL;
+ const struct accum* obs_acc_temp = NULL;
+ const struct accum* obs_acc_time = NULL;
+ struct sdis_estimator* estim = NULL;
+
+ /* Get observable data */
+ obs_args = obs_list_args + iobs;
+ obs_acc_temp = per_obs_acc_temp + iobs;
+ obs_acc_time = per_obs_acc_time + iobs;
+ ASSERT(obs_acc_temp->count == obs_acc_time->count);
+
+ /* Setup the estimator of the current observable */
+ estim = estimator_buffer_grab(estim_buf, iobs, 0);
+ estimator_setup_realisations_count
+ (estim, obs_args->nrealisations, obs_acc_temp->count);
+ estimator_setup_temperature
+ (estim, obs_acc_temp->sum, obs_acc_temp->sum2);
+ estimator_setup_realisation_time
+ (estim, obs_acc_time->sum, obs_acc_time->sum2);
+
+ /* Update global accumulators */
+ acc_temp.sum += obs_acc_temp->sum;
+ acc_temp.sum2 += obs_acc_temp->sum2;
+ acc_temp.count += obs_acc_temp->count;
+ acc_time.sum += obs_acc_time->sum;
+ acc_time.sum2 += obs_acc_time->sum2;
+ acc_time.count += obs_acc_time->count;
+ nrealisations += obs_args->nrealisations;
+ }
+
+ ASSERT(acc_temp.count == acc_time.count);
+
+ /* Setup global estimator */
+ estimator_buffer_setup_realisations_count
+ (estim_buf, nrealisations, acc_temp.count);
+ estimator_buffer_setup_temperature
+ (estim_buf, acc_temp.sum, acc_temp.sum2);
+ estimator_buffer_setup_realisation_time
+ (estim_buf, acc_time.sum, acc_time.sum2);
+
+ res = estimator_buffer_save_rng_state(estim_buf, rng_proxy);
+ if(res != RES_OK) goto error;
+
+exit:
+ *out_estim_buffer = estim_buf;
+ return res;
+error:
+ goto exit;
+}
+
+#undef X_OBS
+#undef SDIS_SOLVE_X_OBS_ARGS
+#undef SDIS_X_OBS
diff --git a/src/sdis_estimator_buffer_c.h b/src/sdis_estimator_buffer_c.h
@@ -59,5 +59,25 @@ estimator_buffer_save_rng_state
(struct sdis_estimator_buffer* buf,
const struct ssp_rng_proxy* proxy);
+extern LOCAL_SYM res_T
+estimator_buffer_create_from_observable_list_probe
+ (struct sdis_device* dev,
+ struct ssp_rng_proxy* rng_proxy,
+ const struct sdis_solve_probe_args obs_list_args[],
+ const struct accum* per_obs_acc_temp,
+ const struct accum* per_obs_acc_time,
+ const size_t nobs, /* #observables */
+ struct sdis_estimator_buffer** out_estim_buffer);
+
+extern LOCAL_SYM res_T
+estimator_buffer_create_from_observable_list_probe_boundary
+ (struct sdis_device* dev,
+ struct ssp_rng_proxy* rng_proxy,
+ const struct sdis_solve_probe_boundary_args obs_list_args[],
+ const struct accum* per_obs_acc_temp,
+ const struct accum* per_obs_acc_time,
+ const size_t nobs, /* #observables */
+ struct sdis_estimator_buffer** out_estim_buffer);
+
#endif /* SDIS_ESTIMATOR_BUFFER_C_H */
diff --git a/src/sdis_green.c b/src/sdis_green.c
@@ -453,8 +453,10 @@ green_function_solve_path
/* Compute external flux */
external_flux = 0;
if(green->scn->source) {
+ /* NOTE: The power of the source is assumed to be constant over time and is
+ * therefore recovered at steady state */
external_flux =
- path->external_flux_term * source_get_power(green->scn->source);
+ path->external_flux_term * source_get_power(green->scn->source, INF);
}
/* Compute path's end temperature */
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
@@ -505,7 +505,7 @@ XD(handle_external_net_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);
+ T->value += external_flux_term * source_get_power(scn->source, frag.time);
/* Register the external net flux term */
if(args->green_path) {
diff --git a/src/sdis_heat_path_boundary_Xd_solid_fluid_picard1.h b/src/sdis_heat_path_boundary_Xd_solid_fluid_picard1.h
@@ -49,6 +49,7 @@ XD(check_Tref)
func_name, scn->tmax, Tref, SPLITX(pos));
return RES_BAD_OP_IRRECOVERABLE;
}
+
return RES_OK;
}
diff --git a/src/sdis_scene.c b/src/sdis_scene.c
@@ -13,14 +13,6 @@
* 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_scene_Xd.h"
-
-/* Generate the Generic functions of the scene */
-#define SDIS_SCENE_DIMENSION 2
-#include "sdis_scene_Xd.h"
-#define SDIS_SCENE_DIMENSION 3
-#include "sdis_scene_Xd.h"
-
#include "sdis.h"
#include "sdis_interface_c.h"
#include "sdis_scene_c.h"
@@ -29,6 +21,12 @@
#include <float.h>
#include <limits.h>
+/* Generate the Generic functions of the scene */
+#define SDIS_XD_DIMENSION 2
+#include "sdis_scene_Xd.h"
+#define SDIS_XD_DIMENSION 3
+#include "sdis_scene_Xd.h"
+
/*******************************************************************************
* Helper function
******************************************************************************/
@@ -240,16 +238,15 @@ sdis_scene_set_temperature_range
res_T
sdis_scene_find_closest_point
(const struct sdis_scene* scn,
- const double pos[],
- const double radius,
+ const struct sdis_scene_find_closest_point_args* args,
size_t* iprim,
double uv[])
{
if(!scn) return RES_BAD_ARG;
if(scene_is_2d(scn)) {
- return scene_find_closest_point_2d(scn, pos, radius, iprim, uv);
+ return scene_find_closest_point_2d(scn, args, iprim, uv);
} else {
- return scene_find_closest_point_3d(scn, pos, radius, iprim, uv);
+ return scene_find_closest_point_3d(scn, args, iprim, uv);
}
}
diff --git a/src/sdis_scene_Xd.h b/src/sdis_scene_Xd.h
@@ -13,8 +13,6 @@
* 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_SCENE_DIMENSION
-
#ifndef SDIS_SCENE_XD_H
#define SDIS_SCENE_XD_H
@@ -123,8 +121,20 @@ check_sdis_scene_create_args(const struct sdis_scene_create_args* args)
&& args->fp_to_meter > 0;
}
+static INLINE res_T
+check_sdis_scene_find_closest_point_args
+ (const struct sdis_scene_find_closest_point_args* args)
+{
+ /* Undefined input arguments */
+ if(!args) return RES_BAD_ARG;
+
+ /* Invalid radius */
+ if(args->radius <= 0) return RES_BAD_ARG;
+
+ return RES_OK;
+}
+
#endif /* SDIS_SCENE_XD_H */
-#else /* !SDIS_SCENE_DIMENSION */
#include "sdis_device_c.h"
@@ -137,39 +147,18 @@ check_sdis_scene_create_args(const struct sdis_scene_create_args* args)
#include <limits.h>
/* Check the submitted dimension and include its specific headers */
-#define SENCXD_DIM SDIS_SCENE_DIMENSION
-#if (SDIS_SCENE_DIMENSION == 2)
+#define SENCXD_DIM SDIS_XD_DIMENSION
+#if (SDIS_XD_DIMENSION == 2)
#include <star/sencX2d.h>
#include <star/s2d.h>
-#elif (SDIS_SCENE_DIMENSION == 3)
+#elif (SDIS_XD_DIMENSION == 3)
#include <star/sencX3d.h>
#include <star/s3d.h>
#else
- #error "Invalid SDIS_SCENE_DIMENSION value."
+ #error "Invalid SDIS_XD_DIMENSION value."
#endif
-/* Syntactic sugar */
-#define DIM SDIS_SCENE_DIMENSION
-
-/* Star-XD macros generic to SDIS_SCENE_DIMENSION */
-#define sXd(Name) CONCAT(CONCAT(CONCAT(s, DIM), d_), Name)
-#define SXD CONCAT(CONCAT(S, DIM), D)
-#define SXD_VERTEX_DATA_NULL CONCAT(CONCAT(S,DIM),D_VERTEX_DATA_NULL)
-#define SXD_POSITION CONCAT(CONCAT(S, DIM), D_POSITION)
-#define SXD_TRACE CONCAT(CONCAT(S,DIM), D_TRACE)
-#define SXD_SAMPLE CONCAT(CONCAT(S,DIM), D_SAMPLE)
-#define SXD_GET_PRIMITIVE CONCAT(CONCAT(S,DIM), D_GET_PRIMITIVE)
-#define SXD_HIT_NONE CONCAT(CONCAT(S,DIM), D_HIT_NONE)
-#define SXD_PRIMITIVE_EQ CONCAT(CONCAT(S,DIM), D_PRIMITIVE_EQ)
-#define SXD_FLOATX CONCAT(CONCAT(CONCAT(S,DIM), D_FLOAT), DIM)
-
-/* Vector macros generic to SDIS_SCENE_DIMENSION */
-#define fX(Func) CONCAT(CONCAT(CONCAT(f, DIM), _), Func)
-#define fX_set_dX CONCAT(CONCAT(CONCAT(f, DIM), _set_d), DIM)
-#define fXX_mulfX CONCAT(CONCAT(CONCAT(CONCAT(f, DIM), DIM), _mulf), DIM)
-
-/* Macro making generic its subimitted name to SDIS_SCENE_DIMENSION */
-#define XD(Name) CONCAT(CONCAT(CONCAT(Name, _), DIM), d)
+#include "sdis_Xd_begin.h"
#if DIM == 2
#define HIT_ON_BOUNDARY hit_on_vertex
@@ -339,22 +328,17 @@ static int
hit_shared_edge
(const struct s3d_primitive* tri0,
const struct s3d_primitive* tri1,
+ const float uv0[2], /* Barycentric coordinates of tested position on tri0 */
+ const float uv1[2], /* Barycentric coordinates of tested position on tri1 */
const float pos0[3], /* Tested position onto the triangle 0 */
const float pos1[3]) /* Tested Position onto the triangle 1 */
{
struct s3d_attrib tri0_vertices[3]; /* Vertex positions of the triangle 0 */
struct s3d_attrib tri1_vertices[3]; /* Vertex positions of the triangle 1 */
- float E0[3], E1[3]; /* Temporary variables storing triangle edges */
- float N0[3], N1[3]; /* Temporary Normals */
- float tri0_2area, tri1_2area; /* 2*area of the submitted triangles */
- float tmp0_2area, tmp1_2area;
- float cos_normals;
int tri0_edge[2] = {-1, -1}; /* Shared edge vertex ids for the triangle 0 */
int tri1_edge[2] = {-1, -1}; /* Shared edge vertex ids for the triangle 1 */
int edge_ivertex = 0; /* Temporary variable */
int tri0_ivertex, tri1_ivertex;
- int iv0, iv1, iv2;
- int hit_edge;
ASSERT(tri0 && tri1 && pos0 && pos1);
/* Fetch the vertices of the triangle 0 */
@@ -386,62 +370,102 @@ hit_shared_edge
}}
/* The triangles do not have a common edge */
- if(edge_ivertex < 2) return 0;
-
- /* Ensure that the vertices of the shared edge are registered in the right
- * order regarding the triangle vertices, i.e. (0,1), (1,2) or (2,0) */
- if((tri0_edge[0]+1)%3 != tri0_edge[1]) SWAP(int, tri0_edge[0], tri0_edge[1]);
- if((tri1_edge[0]+1)%3 != tri1_edge[1]) SWAP(int, tri1_edge[0], tri1_edge[1]);
-
- /* Compute the shared edge normal lying in the triangle 0 plane */
- iv0 = tri0_edge[0];
- iv1 = tri0_edge[1];
- iv2 = (tri0_edge[1]+1) % 3;
- f3_sub(E0, tri0_vertices[iv1].value, tri0_vertices[iv0].value);
- f3_sub(E1, tri0_vertices[iv2].value, tri0_vertices[iv0].value);
- f3_cross(N0, E0, E1); /* Triangle 0 normal */
- tri0_2area = f3_len(N0);
- f3_cross(N0, N0, E0);
-
- /* Compute the shared edge normal lying in the triangle 1 plane */
- iv0 = tri1_edge[0];
- iv1 = tri1_edge[1];
- iv2 = (tri1_edge[1]+1) % 3;
- f3_sub(E0, tri1_vertices[iv1].value, tri1_vertices[iv0].value);
- f3_sub(E1, tri1_vertices[iv2].value, tri1_vertices[iv0].value);
- f3_cross(N1, E0, E1);
- tri1_2area = f3_len(N1);
- f3_cross(N1, N1, E0);
-
- /* Compute the cosine between the 2 edge normals */
- f3_normalize(N0, N0);
- f3_normalize(N1, N1);
- cos_normals = f3_dot(N0, N1);
-
- /* The angle formed by the 2 triangles is sharp */
- if(cos_normals > SHARP_ANGLE_COS_THRESOLD) return 0;
+ if(edge_ivertex == 0) {
+ return 0;
+
+ /* The triangles have a common vertex */
+ } else if(edge_ivertex == 1) {
+ float bcoord0, bcoord1;
+ int hit_vertex;
+
+ /* Retrieve the barycentric coordinate of the position on triangle 0
+ * corresponding to the vertex shared between the 2 triangles. */
+ switch(tri0_edge[0]) {
+ case 0: bcoord0 = uv0[0]; break;
+ case 1: bcoord0 = uv0[1]; break;
+ case 2: bcoord0 = CLAMP(1.f - uv0[0] - uv0[1], 0.f, 1.f); break;
+ default: FATAL("Unreachable code\n"); break;
+ }
- /* Compute the 2 times the area of the (pos0, shared_edge.vertex0,
- * shared_edge.vertex1) triangles */
- f3_sub(E0, tri0_vertices[tri0_edge[0]].value, pos0);
- f3_sub(E1, tri0_vertices[tri0_edge[1]].value, pos0);
- tmp0_2area = f3_len(f3_cross(N0, E0, E1));
-
- /* Compute the 2 times the area of the (pos1, shared_edge.vertex0,
- * shared_edge.vertex1) triangles */
- f3_sub(E0, tri1_vertices[tri1_edge[0]].value, pos1);
- f3_sub(E1, tri1_vertices[tri1_edge[1]].value, pos1);
- tmp1_2area = f3_len(f3_cross(N1, E0, E1));
-
- 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;
+ /* Retrieve the barycentric coordinate of the position on triangle 1
+ * corresponding to the vertex shared between the 2 triangles. */
+ switch(tri1_edge[0]) {
+ case 0: bcoord1 = uv1[0]; break;
+ case 1: bcoord1 = uv1[1]; break;
+ case 2: bcoord1 = CLAMP(1.f - uv0[0] - uv0[1], 0.f, 1.f); break;
+ default: FATAL("Unreachable code\n"); break;
+ }
+
+ /* Check that the both positions lie on the shared vertex */
+ hit_vertex = eq_epsf(1.f, bcoord0, ON_VERTEX_EPSILON)
+ && eq_epsf(1.f, bcoord1, ON_VERTEX_EPSILON);
+ return hit_vertex;
+
+ /* The triangles have a common edge */
+ } else {
+ float E0[3], E1[3]; /* Temporary variables storing triangle edges */
+ float N0[3], N1[3]; /* Temporary Normals */
+ float tri0_2area, tri1_2area; /* 2*area of the submitted triangles */
+ float tmp0_2area, tmp1_2area;
+ float cos_normals;
+ int iv0, iv1, iv2;
+ int hit_edge;
+
+ /* Ensure that the vertices of the shared edge are registered in the right
+ * order regarding the triangle vertices, i.e. (0,1), (1,2) or (2,0) */
+ if((tri0_edge[0]+1)%3 != tri0_edge[1]) SWAP(int, tri0_edge[0], tri0_edge[1]);
+ if((tri1_edge[0]+1)%3 != tri1_edge[1]) SWAP(int, tri1_edge[0], tri1_edge[1]);
+
+ /* Compute the shared edge normal lying in the triangle 0 plane */
+ iv0 = tri0_edge[0];
+ iv1 = tri0_edge[1];
+ iv2 = (tri0_edge[1]+1) % 3;
+ f3_sub(E0, tri0_vertices[iv1].value, tri0_vertices[iv0].value);
+ f3_sub(E1, tri0_vertices[iv2].value, tri0_vertices[iv0].value);
+ f3_cross(N0, E0, E1); /* Triangle 0 normal */
+ tri0_2area = f3_len(N0);
+ f3_cross(N0, N0, E0);
+
+ /* Compute the shared edge normal lying in the triangle 1 plane */
+ iv0 = tri1_edge[0];
+ iv1 = tri1_edge[1];
+ iv2 = (tri1_edge[1]+1) % 3;
+ f3_sub(E0, tri1_vertices[iv1].value, tri1_vertices[iv0].value);
+ f3_sub(E1, tri1_vertices[iv2].value, tri1_vertices[iv0].value);
+ f3_cross(N1, E0, E1);
+ tri1_2area = f3_len(N1);
+ f3_cross(N1, N1, E0);
+
+ /* Compute the cosine between the 2 edge normals */
+ f3_normalize(N0, N0);
+ f3_normalize(N1, N1);
+ cos_normals = f3_dot(N0, N1);
+
+ /* The angle formed by the 2 triangles is sharp */
+ if(cos_normals > SHARP_ANGLE_COS_THRESOLD) return 0;
+
+ /* Compute the 2 times the area of the (pos0, shared_edge.vertex0,
+ * shared_edge.vertex1) triangles */
+ f3_sub(E0, tri0_vertices[tri0_edge[0]].value, pos0);
+ f3_sub(E1, tri0_vertices[tri0_edge[1]].value, pos0);
+ tmp0_2area = f3_len(f3_cross(N0, E0, E1));
+
+ /* Compute the 2 times the area of the (pos1, shared_edge.vertex0,
+ * shared_edge.vertex1) triangles */
+ f3_sub(E0, tri1_vertices[tri1_edge[0]].value, pos1);
+ f3_sub(E1, tri1_vertices[tri1_edge[1]].value, pos1);
+ tmp1_2area = f3_len(f3_cross(N1, E0, E1));
+
+ 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 */
@@ -454,15 +478,26 @@ XD(hit_filter_function)
const float org[DIM],
const float dir[DIM],
const float range[2],
- void* ray_data,
+ void* query_data,
void* global_data)
{
- const struct hit_filter_data* filter_data = ray_data;
- const struct sXd(hit)* hit_from = &filter_data->XD(hit);
+ const struct hit_filter_data* filter_data = query_data;
+ const struct sXd(hit)* hit_from = NULL;
(void)org, (void)dir, (void)global_data, (void)range;
/* No user defined data. Do not filter */
- if(!ray_data || SXD_HIT_NONE(hit_from)) return 0;
+ if(!filter_data) return 0;
+
+ /* Call the custom filter function if it exists
+ * or perform regular filtering otherwise */
+ if(filter_data->XD(custom_filter)) {
+ return filter_data->XD(custom_filter)
+ (hit, org, dir, range, filter_data->custom_filter_data, global_data);
+ }
+
+ /* There is no intersection to discard */
+ hit_from = &filter_data->XD(hit);
+ if(SXD_HIT_NONE(hit_from)) return 0;
if(SXD_PRIMITIVE_EQ(&hit_from->prim, &hit->prim)) return 1;
@@ -471,14 +506,17 @@ XD(hit_filter_function)
if(eq_epsf(hit->distance, 0, (float)filter_data->epsilon)) {
float pos[DIM];
+ int reject_hit = 0;
fX(add)(pos, org, fX(mulf)(pos, dir, hit->distance));
/* If the targeted point is near of the origin, check that it lies on an
* edge/vertex shared by the 2 primitives */
#if DIM == 2
- return hit_shared_vertex(&hit_from->prim, &hit->prim, org, pos);
+ reject_hit = hit_shared_vertex(&hit_from->prim, &hit->prim, org, pos);
#else
- return hit_shared_edge(&hit_from->prim, &hit->prim, org, pos);
+ reject_hit = hit_shared_edge
+ (&hit_from->prim, &hit->prim, hit_from->uv, hit->uv, org, pos);
#endif
+ return reject_hit;
}
return 0;
}
@@ -988,8 +1026,7 @@ error:
static res_T
XD(scene_find_closest_point)
(const struct sdis_scene* scn,
- const double pos[3],
- const double radius,
+ const struct sdis_scene_find_closest_point_args* args,
size_t* iprim,
double uv[2])
{
@@ -998,30 +1035,37 @@ XD(scene_find_closest_point)
float query_radius;
res_T res = RES_OK;
- if(!scn || !pos || radius <= 0 || !iprim || !uv
- || scene_is_2d(scn) != (DIM == 2)) {
+ if(!scn || !iprim || !uv || scene_is_2d(scn) != (DIM == 2)) {
res = RES_BAD_ARG;
goto error;
}
+ res = check_sdis_scene_find_closest_point_args(args);
+ if(res != RES_OK) goto error;
/* Avoid a null query radius due to casting in single-precision */
- query_radius = MMAX((float)radius, FLT_MIN);
+ query_radius = MMAX((float)args->radius, FLT_MIN);
+
+ fX_set_dX(query_pos, args->position);
+
+ /* Do not filter anything */
+ if(!args->XD(filter)) {
+ res = sXd(scene_view_closest_point)
+ (scn->sXd(view), query_pos, query_radius, NULL, &hit);
+
+ /* Filter points according to user-defined filter function */
+ } else {
+ struct hit_filter_data filter_data = HIT_FILTER_DATA_NULL;
+ filter_data.XD(custom_filter) = args->XD(filter);
+ filter_data.custom_filter_data = args->filter_data;
+ res = sXd(scene_view_closest_point)
+ (scn->sXd(view), query_pos, query_radius, &filter_data, &hit);
+ }
- fX_set_dX(query_pos, pos);
- res = sXd(scene_view_closest_point)
- (scn->sXd(view), query_pos, query_radius, NULL, &hit);
if(res != RES_OK) {
-#if DIM == 2
log_err(scn->dev,
- "%s: error querying the closest position at {%g, %g} "
+ "%s: error querying the closest position at `"FORMAT_VECX"' "
"for a radius of %g -- %s.\n",
- FUNC_NAME, SPLIT2(query_pos), query_radius, res_to_cstr(res));
-#else
- log_err(scn->dev,
- "%s: error querying the closest position at {%g, %g, %g} "
- "for a radius of %g -- %s.\n",
- FUNC_NAME, SPLIT3(query_pos), query_radius, res_to_cstr(res));
-#endif
+ FUNC_NAME, SPLITX(query_pos), query_radius, res_to_cstr(res));
goto error;
}
@@ -1250,29 +1294,12 @@ error:
goto exit;
}
-#if (SDIS_SCENE_DIMENSION == 2)
-#include <star/sencX2d_undefs.h>
-#else /* SDIS_SCENE_DIMENSION == 3 */
-#include <star/sencX3d_undefs.h>
+#if (SDIS_XD_DIMENSION == 2)
+ #include <star/sencX2d_undefs.h>
+#else /* SDIS_XD_DIMENSION == 3 */
+ #include <star/sencX3d_undefs.h>
#endif
-#undef SDIS_SCENE_DIMENSION
-#undef DIM
-#undef sXd
-#undef SXD
-#undef SXD_VERTEX_DATA_NULL
-#undef SXD_POSITION
-#undef SXD_FLOAT3
-#undef SXD_TRACE
-#undef SXD_TRACE
-#undef SXD_GET_PRIMITIVE
-#undef SXD_HIT_NONE
-#undef SXD_PRIMITIVE_EQ
-#undef SXD_FLOATX
-#undef fX
-#undef fX_set_dX
-#undef fXX_mulfX
-#undef XD
#undef HIT_ON_BOUNDARY
-#endif /* !SDIS_SCENE_DIMENSION */
+#include "sdis_Xd_end.h"
diff --git a/src/sdis_scene_c.h b/src/sdis_scene_c.h
@@ -36,8 +36,15 @@ struct hit_filter_data {
struct s2d_hit hit_2d;
struct s3d_hit hit_3d;
double epsilon; /* Threshold defining roughly equal intersections */
+
+ /* Bypass the regular filter function */
+ s2d_hit_filter_function_T custom_filter_2d;
+ s3d_hit_filter_function_T custom_filter_3d;
+
+ /* Custom filter query data. It is ignored if custom_filter is NULL */
+ void* custom_filter_data;
};
-#define HIT_FILTER_DATA_NULL__ {S2D_HIT_NULL__, S3D_HIT_NULL__, 0}
+#define HIT_FILTER_DATA_NULL__ {S2D_HIT_NULL__,S3D_HIT_NULL__,0,NULL,NULL,NULL}
static const struct hit_filter_data HIT_FILTER_DATA_NULL =
HIT_FILTER_DATA_NULL__;
diff --git a/src/sdis_solve.c b/src/sdis_solve.c
@@ -59,7 +59,7 @@ sdis_solve_probe
res_T
sdis_solve_probe_list
(struct sdis_scene* scn,
- const struct sdis_solve_probe_list_args args[],
+ const struct sdis_solve_probe_list_args* args,
struct sdis_estimator_buffer** out_buf)
{
if(!scn) return RES_BAD_ARG;
@@ -99,6 +99,20 @@ sdis_solve_probe_boundary
}
res_T
+sdis_solve_probe_boundary_list
+ (struct sdis_scene* scn,
+ const struct sdis_solve_probe_boundary_list_args* args,
+ struct sdis_estimator_buffer** out_buf)
+{
+ if(!scn) return RES_BAD_ARG;
+ if(scene_is_2d(scn)) {
+ return solve_probe_boundary_list_2d(scn, args, out_buf);
+ } else {
+ return solve_probe_boundary_list_3d(scn, args, out_buf);
+ }
+}
+
+res_T
sdis_solve_probe_boundary_green_function
(struct sdis_scene* scn,
const struct sdis_solve_probe_boundary_args* args,
diff --git a/src/sdis_solve_probe_Xd.h b/src/sdis_solve_probe_Xd.h
@@ -106,84 +106,6 @@ check_solve_probe_list_args
return RES_OK;
}
-static res_T
-setup_estimator_buffer
- (struct sdis_device* dev,
- struct ssp_rng_proxy* rng_proxy,
- const struct sdis_solve_probe_list_args* solve_args,
- const struct accum* per_probe_acc_temp,
- const struct accum* per_probe_acc_time,
- struct sdis_estimator_buffer** out_estim_buffer)
-{
- /* Accumulators throughout the buffer */
- struct accum acc_temp = ACCUM_NULL;
- struct accum acc_time = ACCUM_NULL;
- size_t nrealisations = 0;
-
- struct sdis_estimator_buffer* estim_buf = NULL;
- size_t iprobe = 0;
- res_T res = RES_OK;
-
- ASSERT(dev && rng_proxy && solve_args);
- ASSERT(per_probe_acc_time && per_probe_acc_time && out_estim_buffer);
-
- res = estimator_buffer_create(dev, solve_args->nprobes, 1, &estim_buf);
- if(res != RES_OK) {
- log_err(dev, "Unable to allocate the estimator buffer.\n");
- goto error;
- }
-
- FOR_EACH(iprobe, 0, solve_args->nprobes) {
- const struct sdis_solve_probe_args* probe = NULL;
- const struct accum* probe_acc_temp = NULL;
- const struct accum* probe_acc_time = NULL;
- struct sdis_estimator* estim = NULL;
-
- /* Get probe data */
- probe = solve_args->probes + iprobe;
- probe_acc_temp = per_probe_acc_temp + iprobe;
- probe_acc_time = per_probe_acc_time + iprobe;
- ASSERT(probe_acc_temp->count == probe_acc_time->count);
-
- /* Setup probe estimator */
- estim = estimator_buffer_grab(estim_buf, iprobe, 0);
- estimator_setup_realisations_count
- (estim, probe->nrealisations, probe_acc_temp->count);
- estimator_setup_temperature
- (estim, probe_acc_temp->sum, probe_acc_temp->sum2);
- estimator_setup_realisation_time
- (estim, probe_acc_time->sum, probe_acc_time->sum2);
-
- /* Update global accumulators */
- acc_temp.sum += probe_acc_temp->sum;
- acc_temp.sum2 += probe_acc_temp->sum2;
- acc_temp.count += probe_acc_temp->count;
- acc_time.sum += probe_acc_time->sum;
- acc_time.sum2 += probe_acc_time->sum2;
- acc_time.count += probe_acc_time->count;
- nrealisations += probe->nrealisations;
- }
-
- ASSERT(acc_temp.count == acc_time.count);
-
- /* Setup global estimator */
- estimator_buffer_setup_realisations_count
- (estim_buf, nrealisations, acc_temp.count);
- estimator_buffer_setup_temperature
- (estim_buf, acc_temp.sum, acc_temp.sum2);
- estimator_buffer_setup_realisation_time
- (estim_buf, acc_time.sum, acc_time.sum2);
-
- res = estimator_buffer_save_rng_state(estim_buf, rng_proxy);
- if(res != RES_OK) goto error;
-
-exit:
- *out_estim_buffer = estim_buf;
- return res;
-error:
- goto exit;
-}
-
#endif /* SDIS_SOLVE_PROBE_XD_H */
static res_T
@@ -633,9 +555,6 @@ XD(solve_probe_list)
goto post_sync;
}
- /* Begin time registration of the computation */
- time_current(&time0);
-
/* Allocate the list of accumulators per probe. On the master process,
* allocate a complete list in which the accumulators of all processes will be
* stored. */
@@ -648,6 +567,9 @@ XD(solve_probe_list)
goto error;
}
+ /* Begin time registration of the computation */
+ time_current(&time0);
+
/* Here we go! Calculation of probe list */
omp_set_num_threads((int)scn->dev->nthreads);
#pragma omp parallel for schedule(static)
@@ -733,8 +655,9 @@ post_sync:
log_info(scn->dev, "Probes accumulator gathered in %s.\n", buf);
if(is_master_process) {
- res = setup_estimator_buffer(scn->dev, rng_proxy, args, per_probe_acc_temp,
- per_probe_acc_time, &estim_buf);
+ res = estimator_buffer_create_from_observable_list_probe
+ (scn->dev, rng_proxy, args->probes, per_probe_acc_temp,
+ per_probe_acc_time, args->nprobes, &estim_buf);
if(res != RES_OK) goto error;
}
diff --git a/src/sdis_solve_probe_boundary_Xd.h b/src/sdis_solve_probe_boundary_Xd.h
@@ -80,6 +80,40 @@ check_solve_probe_boundary_args
}
static INLINE res_T
+check_solve_probe_boundary_list_args
+ (struct sdis_device* dev,
+ const struct sdis_solve_probe_boundary_list_args* args)
+{
+ size_t iprobe = 0;
+
+ if(!args) return RES_BAD_ARG;
+
+ /* Check the list of probes */
+ if(!args->probes || !args->nprobes) {
+ return RES_BAD_ARG;
+ }
+
+ /* Check the RNG type */
+ if(!args->rng_state && args->rng_type >= SSP_RNG_TYPES_COUNT__) {
+ return RES_BAD_ARG;
+ }
+
+ FOR_EACH(iprobe, 0, args->nprobes) {
+ const res_T res = check_solve_probe_boundary_args(args->probes+iprobe);
+ if(res != RES_OK) return res;
+
+ if(args->probes[iprobe].register_paths != SDIS_HEAT_PATH_NONE) {
+ log_warn(dev,
+ "Unable to save paths for probe boundary %lu. "
+ "Saving path is not supported when solving multiple probes\n",
+ (unsigned long)iprobe);
+ }
+ }
+
+ return RES_OK;
+}
+
+static INLINE res_T
check_solve_probe_boundary_flux_args
(const struct sdis_solve_probe_boundary_flux_args* args)
{
@@ -119,6 +153,66 @@ check_solve_probe_boundary_flux_args
#endif /* SDIS_SOLVE_PROBE_BOUNDARY_XD_H */
+static res_T
+XD(solve_one_probe_boundary)
+ (struct sdis_scene* scn,
+ struct ssp_rng* rng,
+ const struct sdis_solve_probe_boundary_args* args,
+ struct accum* acc_temp,
+ struct accum* acc_time)
+{
+ size_t irealisation = 0;
+ res_T res = RES_OK;
+ ASSERT(scn && rng && check_solve_probe_boundary_args(args) == RES_OK);
+
+ *acc_temp = ACCUM_NULL;
+ *acc_time = ACCUM_NULL;
+
+ FOR_EACH(irealisation, 0, args->nrealisations) {
+ struct boundary_realisation_args realis_args = BOUNDARY_REALISATION_ARGS_NULL;
+ double w = NaN; /* MC weight */
+ double usec = 0; /* Time of a realisation */
+ double time = 0; /* Sampled observation time */
+ struct time t0, t1; /* Register the time spent solving a realisation */
+
+ /* Begin time registration of the realisation */
+ time_current(&t0);
+
+ /* Sample observation time */
+ time = sample_time(rng, args->time_range);
+
+ /* Run a realisation */
+ realis_args.rng = rng;
+ realis_args.iprim = args->iprim;
+ realis_args.time = time;
+ realis_args.picard_order = args->picard_order;
+ realis_args.side = args->side;
+ realis_args.irealisation = irealisation;
+ realis_args.uv[0] = args->uv[0];
+#if SDIS_XD_DIMENSION == 3
+ realis_args.uv[1] = args->uv[1];
+#endif
+ res = XD(boundary_realisation)(scn, &realis_args, &w);
+ if(res != RES_OK) goto error;
+
+ /* Stop time registration */
+ time_sub(&t0, time_current(&t1), &t0);
+ usec = (double)time_val(&t0, TIME_NSEC) * 0.001;
+
+ /* Update MC weights */
+ acc_temp->sum += w;
+ acc_temp->sum2 += w*w;
+ acc_temp->count += 1;
+ acc_time->sum += usec;
+ acc_time->sum2 += usec*usec;
+ acc_time->count += 1;
+ }
+exit:
+ return res;
+error:
+ goto exit;
+}
+
/*******************************************************************************
* Local functions
******************************************************************************/
@@ -422,6 +516,207 @@ error:
}
static res_T
+XD(solve_probe_boundary_list)
+ (struct sdis_scene* scn,
+ const struct sdis_solve_probe_boundary_list_args* args,
+ struct sdis_estimator_buffer** out_estim_buf)
+{
+ /* Time registration */
+ struct time time0, time1;
+ char buf[128]; /* Temporary buffer used to store formated time */
+
+ /* Device variable */
+ struct mem_allocator* allocator = NULL;
+
+ /* Stardis variables */
+ struct sdis_estimator_buffer* estim_buf = NULL;
+
+ /* Random Number generator */
+ struct ssp_rng_proxy* rng_proxy = NULL;
+ struct ssp_rng** per_thread_rng = NULL;
+
+ /* Probe variables */
+ size_t process_probes[2] = {0, 0}; /* Probes range managed by the process */
+ size_t process_nprobes = 0; /* Number of probes managed by the process */
+ size_t nprobes = 0;
+ struct accum* per_probe_acc_temp = NULL;
+ struct accum* per_probe_acc_time = NULL;
+
+ /* Miscellaneous */
+ int32_t* progress = NULL; /* Per process progress bar */
+ int is_master_process = 0;
+ int pcent_progress = 1; /* Percentage requiring progress update */
+ int64_t i = 0;
+ ATOMIC nsolved_probes = 0;
+ ATOMIC res = RES_OK;
+
+ if(!scn || !out_estim_buf) { res = RES_BAD_ARG; goto error; }
+ res = check_solve_probe_boundary_list_args(scn->dev, args);
+ if(res != RES_OK) goto error;
+ res = XD(scene_check_dimensionality)(scn);
+ if(res != RES_OK) goto error;
+
+#ifdef SDIS_ENABLE_MPI
+ is_master_process = !scn->dev->use_mpi || scn->dev->mpi_rank == 0;
+#endif
+
+ allocator = scn->dev->allocator;
+
+ /* Update the progress bar every percent if escape sequences are allowed in
+ * log messages or only every 10 percent when only plain text is allowed.
+ * This reduces the number of lines of plain text printed */
+ pcent_progress = scn->dev->no_escape_sequence ? 10 : 1;
+
+ /* Create the per threads RNGs */
+ res = create_per_thread_rng
+ (scn->dev, args->rng_state, args->rng_type, &rng_proxy, &per_thread_rng);
+ if(res != RES_OK) goto error;
+
+ /* Allocate the per process progress status */
+ res = alloc_process_progress(scn->dev, &progress);
+ if(res != RES_OK) goto error;
+
+ /* Synchronise the processes */
+ process_barrier(scn->dev);
+
+ /* Define the range of probes to manage in this process */
+ process_nprobes = compute_process_index_range
+ (scn->dev, args->nprobes, process_probes);
+
+ #define PROGRESS_MSG "Solving surface probes: "
+ print_progress(scn->dev, progress, PROGRESS_MSG);
+
+ /* If there is no work to be done on this process (i.e. no probe to
+ * calculate), simply print its completion and go straight to the
+ * synchronization barrier.*/
+ if(process_nprobes == 0) {
+ progress[0] = 100;
+ print_progress_update(scn->dev, progress, PROGRESS_MSG);
+ goto post_sync;
+ }
+
+ /* Allocate the list of accumulators per probe. On the master process,
+ * allocate a complete list in which the accumulators of all processes will be
+ * stored. */
+ nprobes = is_master_process ? args->nprobes : process_nprobes;
+ per_probe_acc_temp = MEM_CALLOC(allocator, nprobes, sizeof(*per_probe_acc_temp));
+ per_probe_acc_time = MEM_CALLOC(allocator, nprobes, sizeof(*per_probe_acc_time));
+ if(!per_probe_acc_temp || !per_probe_acc_time) {
+ log_err(scn->dev, "Unable to allocate the list of accumulators per probe.\n");
+ res = RES_MEM_ERR;
+ goto error;
+ }
+
+ /* Begin time registration of the computation */
+ time_current(&time0);
+
+ /* Calculation of probe list */
+ #pragma omp parallel for schedule(static)
+ for(i = 0; i < (int64_t)process_nprobes; ++i) {
+ /* Thread */
+ const int ithread = omp_get_thread_num(); /* Thread ID */
+ struct ssp_rng* rng = per_thread_rng[ithread]; /* RNG of the thread */
+
+ /* Probe */
+ struct accum* probe_acc_temp = NULL;
+ struct accum* probe_acc_time = NULL;
+ const struct sdis_solve_probe_boundary_args* probe_args = NULL;
+ const size_t iprobe = process_probes[0] + (size_t)i; /* Probe ID */
+
+ /* Miscellaneous */
+ size_t n = 0; /* Number of solved probes */
+ int pcent = 0; /* Current progress */
+ res_T res_local = RES_OK;
+
+ if(ATOMIC_GET(&res) != RES_OK) continue; /* An error occurred */
+
+ /* Retrieve the probe arguments */
+ probe_args = &args->probes[iprobe];
+
+ /* Retrieve the probe accumulators */
+ probe_acc_temp = &per_probe_acc_temp[i];
+ probe_acc_time = &per_probe_acc_time[i];
+
+ res_local = XD(solve_one_probe_boundary)
+ (scn, rng, probe_args, probe_acc_temp, probe_acc_time);
+ if(res_local != RES_OK) {
+ ATOMIC_SET(&res, res_local);
+ continue;
+ }
+
+ /* Update progress */
+ n = (size_t)ATOMIC_INCR(&nsolved_probes);
+ pcent = (int)((double)n * 100.0 / (double)process_nprobes + 0.5/*round*/);
+
+ #pragma omp critical
+ if(pcent/pcent_progress > progress[0]/pcent_progress) {
+ progress[0] = pcent;
+ print_progress_update(scn->dev, progress, PROGRESS_MSG);
+ }
+ }
+
+post_sync:
+ /* Synchronise processes */
+ process_barrier(scn->dev);
+
+ res = gather_res_T(scn->dev, (res_T)res);
+ if(res != RES_OK) goto error;
+
+ print_progress_completion(scn->dev, progress, PROGRESS_MSG);
+ #undef PROGRESS_MSG
+
+ /* Report computatio time */
+ time_sub(&time0, time_current(&time1), &time0);
+ time_dump(&time0, TIME_ALL, NULL, buf, sizeof(buf));
+ log_info(scn->dev, "%lu surface probes solved in %s.\n",
+ (unsigned long)args->nprobes, buf);
+
+ /* Gather the RNG proxy sequence IDs and ensure that the RNG proxy state of
+ * the master process is greater than the RNG proxy state of all other
+ * processes */
+ res = gather_rng_proxy_sequence_id(scn->dev, rng_proxy);
+ if(res != RES_OK) goto error;
+
+ time_current(&time0);
+
+ /* Gather the list of accumulators */
+ #define GATHER_ACCUMS_LIST(Msg, Acc) { \
+ res = gather_accumulators_list \
+ (scn->dev, Msg, args->nprobes, process_probes, per_probe_acc_##Acc); \
+ if(res != RES_OK) goto error; \
+ } (void)0
+ GATHER_ACCUMS_LIST(MPI_SDIS_MSG_ACCUM_TEMP, temp);
+ GATHER_ACCUMS_LIST(MPI_SDIS_MSG_ACCUM_TIME, time);
+ #undef GATHER_ACCUMS_LIST
+
+ time_sub(&time0, time_current(&time1), &time0);
+ time_dump(&time0, TIME_ALL, NULL, buf, sizeof(buf));
+ log_info(scn->dev, "Probes accumulator gathered in %s.\n", buf);
+
+ if(is_master_process) {
+ res = estimator_buffer_create_from_observable_list_probe_boundary
+ (scn->dev, rng_proxy, args->probes, per_probe_acc_temp,
+ per_probe_acc_time, args->nprobes, &estim_buf);
+ if(res != RES_OK) goto error;
+ }
+
+exit:
+ if(per_thread_rng) release_per_thread_rng(scn->dev, per_thread_rng);
+ if(rng_proxy) SSP(rng_proxy_ref_put(rng_proxy));
+ if(per_probe_acc_temp) MEM_RM(allocator, per_probe_acc_temp);
+ if(per_probe_acc_time) MEM_RM(allocator, per_probe_acc_time);
+ if(progress) free_process_progress(scn->dev, progress);
+ if(out_estim_buf) *out_estim_buf = estim_buf;
+ return (res_T)res;
+error:
+ if(estim_buf) {
+ SDIS(estimator_buffer_ref_put(estim_buf));
+ estim_buf = NULL;
+ }
+ goto exit;
+}
+
+static res_T
XD(solve_probe_boundary_flux)
(struct sdis_scene* scn,
const struct sdis_solve_probe_boundary_flux_args* args,
diff --git a/src/sdis_source.c b/src/sdis_source.c
@@ -47,15 +47,14 @@ check_spherical_source_create_args
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);
+ if(!args->power) {
+ log_err(dev, "%s: the power functor is missing.\n", func_name);
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);
+ 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;
}
@@ -126,12 +125,11 @@ sdis_source_ref_put(struct sdis_source* src)
return RES_OK;
}
-res_T
-sdis_source_get_power(const struct sdis_source* src, double* power)
+double
+sdis_source_get_power(struct sdis_source* src, const double time /* [s] */)
{
- if(!src || !power) return RES_BAD_ARG;
- *power = source_get_power(src);
- return RES_OK;
+ ASSERT(src);
+ return source_get_power(src, time);
}
/*******************************************************************************
@@ -151,14 +149,23 @@ source_sample
double cos_half_angle; /* [radians] */
double dst; /* [m] */
double radius; /* Source radius [m] */
+ double power; /* Source power [W] */
double area; /* Source area [m^2] */
res_T res = RES_OK;
ASSERT(src && rng && pos && sample);
- /* Retrieve current source position and radius */
+ /* Retrieve current source position, radius and power */
src->spherical.position(time, src_pos, src->spherical.data);
+ power = src->spherical.power(time, src->spherical.data);
radius = src->spherical.radius;
+ if(power < 0) {
+ log_err(src->dev, "%s: invalid source power '%g' W. It cannot be negative.\n",
+ FUNC_NAME, power);
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
area = 4*PI*radius*radius; /* [m^2] */
/* compute the direction of `pos' toward the center of the source */
@@ -186,7 +193,7 @@ source_sample
d3_set(sample->dir, main_dir);
sample->pdf = 1;
sample->dst = dst;
- sample->power = src->spherical.power; /* [W] */
+ sample->power = 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] */
@@ -201,7 +208,7 @@ source_sample
/* Set other sample variables */
sample->dst = dst - radius; /* From pos to source boundaries [m] */
- sample->power = src->spherical.power; /* [W] */
+ sample->power = power; /* [W] */
sample->radiance_term = 1.0 / (PI*area); /* [W/m^2/sr] */
sample->radiance = sample->power * sample->radiance_term; /* [W/m^2/sr] */
}
@@ -223,6 +230,7 @@ source_trace_to
double src_pos[3]; /* [m] */
double main_dir[3];
double radius; /* [m] */
+ double power; /* [W] */
double dst; /* Distance from pos to the source center [m] */
double half_angle; /* [radian] */
res_T res = RES_OK;
@@ -237,8 +245,16 @@ source_trace_to
goto exit;
}
- /* Retrieve current source position and radius */
+ /* Retrieve current source position and power */
src->spherical.position(time, src_pos, src->spherical.data);
+ power = src->spherical.power(time, src->spherical.data);
+
+ if(power < 0) {
+ log_err(src->dev, "%s: invalid source power '%g' W. It cannot be negative.\n",
+ FUNC_NAME, power);
+ res = RES_BAD_ARG;
+ goto error;
+ }
/* compute the direction of `pos' toward the center of the source */
d3_sub(main_dir, src_pos, pos);
@@ -274,7 +290,7 @@ source_trace_to
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->power = power; /* [W] */
sample->radiance_term = 1.0 / (PI*area); /* [W/m^2/sr] */
sample->radiance = sample->power * sample->radiance_term; /* [W/m^2/sr] */
}
@@ -286,11 +302,11 @@ error:
goto exit;
}
-double
-source_get_power(const struct sdis_source* src)
+double /* [W] */
+source_get_power(const struct sdis_source* src, const double time /* [s] */)
{
ASSERT(src);
- return src->spherical.power;
+ return src->spherical.power(time, src->spherical.data);
}
void
diff --git a/src/sdis_source_c.h b/src/sdis_source_c.h
@@ -63,7 +63,8 @@ source_trace_to
extern LOCAL_SYM double /* [W] */
source_get_power
- (const struct sdis_source* source);
+ (const struct sdis_source* source,
+ const double time); /* [s] */
extern LOCAL_SYM void
source_compute_signature
diff --git a/src/test_sdis_draw_external_flux.c b/src/test_sdis_draw_external_flux.c
@@ -14,10 +14,10 @@
* along with this program. If not, see <http://www.gnu.org/licenses/>. */
#include "sdis.h"
+#include "test_sdis_mesh.h"
#include "test_sdis_utils.h"
#include <star/s3dut.h>
-#include <rsys/stretchy_array.h>
#define IMG_WIDTH 320
#define IMG_HEIGHT 240
@@ -37,80 +37,6 @@
/*******************************************************************************
* 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)
{
@@ -125,7 +51,8 @@ mesh_add_super_shape(struct mesh* mesh)
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);
+ mesh_append(mesh, sshape_data.positions, sshape_data.nvertices,
+ sshape_data.indices, sshape_data.nprimitives, NULL);
OK(s3dut_mesh_ref_put(sshape));
}
@@ -140,7 +67,8 @@ mesh_add_ground(struct mesh* mesh)
OK(s3dut_create_cuboid(NULL, 20, 20, DEPTH, &ground));
OK(s3dut_mesh_get_data(ground, &ground_data));
- mesh_append(mesh, &ground_data, translate);
+ mesh_append(mesh, ground_data.positions, ground_data.nvertices,
+ ground_data.indices, ground_data.nprimitives, translate);
OK(s3dut_mesh_ref_put(ground));
#undef DEPTH
}
@@ -254,6 +182,13 @@ source_get_position
pos[2] = 5.0;
}
+static double
+source_get_power(const double time, struct sdis_data* data)
+{
+ (void)time, (void)data; /* Avoid the "unusued variable" warning */
+ return SOURCE_POWER; /* [W] */
+}
+
static struct sdis_source*
create_source(struct sdis_device* sdis)
{
@@ -261,9 +196,9 @@ create_source(struct sdis_device* sdis)
struct sdis_source* source = NULL;
args.position = source_get_position;
+ args.power = source_get_power;
args.data = NULL;
args.radius = 3e-1; /* [m] */
- args.power = SOURCE_POWER; /* [W] */
OK(sdis_spherical_source_create(sdis, &args, &source));
return source;
}
diff --git a/src/test_sdis_external_flux.c b/src/test_sdis_external_flux.c
@@ -261,6 +261,13 @@ source_get_position
pos[2] = 0;
}
+static double
+source_get_power(const double time, struct sdis_data* data)
+{
+ (void)time, (void)data; /* Avoid the "unusued variable" warning */
+ return 3.845e26; /* [W] */
+}
+
static struct sdis_source*
create_source(struct sdis_device* sdis)
{
@@ -268,9 +275,9 @@ create_source(struct sdis_device* sdis)
struct sdis_source* src = NULL;
args.position = source_get_position;
+ args.power = source_get_power;
args.data = NULL;
args.radius = 6.5991756e8; /* [m] */
- args.power = 3.845e26; /* [W] */
OK(sdis_spherical_source_create(sdis, &args, &src));
return src;
}
diff --git a/src/test_sdis_mesh.h b/src/test_sdis_mesh.h
@@ -0,0 +1,120 @@
+/* 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 TEST_SDIS_MESH_H
+#define TEST_SDIS_MESH_H
+
+#include <rsys/rsys.h>
+#include <rsys/stretchy_array.h>
+
+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 INLINE void
+mesh_append
+ (struct mesh* mesh,
+ const double* in_positions,
+ const size_t in_nvertices,
+ const size_t* in_indices,
+ const size_t in_ntriangles,
+ 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, in_nvertices*3);
+ indices = sa_add(mesh->indices, in_ntriangles*3);
+
+ if(in_translate) {
+ translate[0] = in_translate[0];
+ translate[1] = in_translate[1];
+ translate[2] = in_translate[2];
+ }
+
+ FOR_EACH(i, 0, in_nvertices) {
+ positions[i*3 + 0] = in_positions[i*3 + 0] + translate[0];
+ positions[i*3 + 1] = in_positions[i*3 + 1] + translate[1];
+ positions[i*3 + 2] = in_positions[i*3 + 2] + translate[2];
+ }
+
+ FOR_EACH(i, 0, in_ntriangles) {
+ indices[i*3 + 0] = in_indices[i*3 + 0] + ivert;
+ indices[i*3 + 1] = in_indices[i*3 + 1] + ivert;
+ indices[i*3 + 2] = in_indices[i*3 + 2] + ivert;
+ }
+}
+
+static INLINE void
+mesh_dump(const struct mesh* mesh, FILE* stream)
+{
+ size_t i, n;
+ CHK(mesh != NULL);
+
+ n = mesh_nvertices(mesh);
+ FOR_EACH(i, 0, n) {
+ fprintf(stream, "v %g %g %g\n", SPLIT3(mesh->positions+i*3));
+ }
+
+ n = mesh_ntriangles(mesh);
+ FOR_EACH(i, 0, n) {
+ fprintf(stream, "f %lu %lu %lu\n",
+ (unsigned long)(mesh->indices[i*3+0] + 1),
+ (unsigned long)(mesh->indices[i*3+1] + 1),
+ (unsigned long)(mesh->indices[i*3+2] + 1));
+ }
+ fflush(stream);
+}
+
+#endif /* TEST_SDIS_MESH_H */
diff --git a/src/test_sdis_scene.c b/src/test_sdis_scene.c
@@ -98,6 +98,8 @@ test_scene_3d(struct sdis_device* dev, struct sdis_interface* interf)
struct senc2d_scene* scn2d;
struct senc3d_scene* scn3d;
struct sdis_scene_create_args scn_args = SDIS_SCENE_CREATE_ARGS_DEFAULT;
+ struct sdis_scene_find_closest_point_args closest_pt_args =
+ SDIS_SCENE_FIND_CLOSEST_POINT_ARGS_NULL;
struct sdis_scene* scn = NULL;
struct sdis_device* dev2 = NULL;
size_t ntris, npos;
@@ -195,12 +197,15 @@ test_scene_3d(struct sdis_device* dev, struct sdis_interface* interf)
BA(sdis_scene_boundary_project_position(scn, 6, pos, NULL));
OK(sdis_scene_boundary_project_position(scn, 6, pos, uv1));
- BA(sdis_scene_find_closest_point(NULL, pos, INF, &iprim, uv2));
- BA(sdis_scene_find_closest_point(scn, NULL, INF, &iprim, uv2));
- BA(sdis_scene_find_closest_point(scn, pos, 0, &iprim, uv2));
- BA(sdis_scene_find_closest_point(scn, pos, INF, NULL, uv2));
- BA(sdis_scene_find_closest_point(scn, pos, INF, &iprim, NULL));
- OK(sdis_scene_find_closest_point(scn, pos, INF, &iprim, uv2));
+ closest_pt_args.position[0] = pos[0];
+ closest_pt_args.position[1] = pos[1];
+ closest_pt_args.position[2] = pos[2];
+ closest_pt_args.radius = INF;
+ BA(sdis_scene_find_closest_point(NULL, &closest_pt_args, &iprim, uv2));
+ BA(sdis_scene_find_closest_point(scn, NULL, &iprim, uv2));
+ BA(sdis_scene_find_closest_point(scn, &closest_pt_args, NULL, uv2));
+ BA(sdis_scene_find_closest_point(scn, &closest_pt_args, &iprim, NULL));
+ OK(sdis_scene_find_closest_point(scn, &closest_pt_args, &iprim, uv2));
CHK(iprim == 6);
CHK(d2_eq_eps(uv0, uv1, 1.e-6));
@@ -209,7 +214,10 @@ test_scene_3d(struct sdis_device* dev, struct sdis_interface* interf)
pos[0] = 0.5;
pos[1] = 0.1;
pos[2] = 0.25;
- OK(sdis_scene_find_closest_point(scn, pos, INF, &iprim, uv2));
+ closest_pt_args.position[0] = pos[0];
+ closest_pt_args.position[1] = pos[1];
+ closest_pt_args.position[2] = pos[2];
+ OK(sdis_scene_find_closest_point(scn, &closest_pt_args, &iprim, uv2));
CHK(iprim == 10);
OK(sdis_scene_boundary_project_position(scn, 10, pos, uv0));
@@ -219,7 +227,11 @@ test_scene_3d(struct sdis_device* dev, struct sdis_interface* interf)
dst = d3_len(d3_sub(pos1, pos, pos1));
CHK(eq_eps(dst, 0.1, 1.e-6));
- OK(sdis_scene_find_closest_point(scn, pos, 0.09, &iprim, uv2));
+ closest_pt_args.position[0] = pos[0];
+ closest_pt_args.position[1] = pos[1];
+ closest_pt_args.position[2] = pos[2];
+ closest_pt_args.radius = 0.09;
+ OK(sdis_scene_find_closest_point(scn, &closest_pt_args, &iprim, uv2));
CHK(iprim == SDIS_PRIMITIVE_NONE);
FOR_EACH(i, 0, 64) {
@@ -228,7 +240,12 @@ test_scene_3d(struct sdis_device* dev, struct sdis_interface* interf)
OK(sdis_scene_get_boundary_position(scn, 4, uv0, pos));
OK(sdis_scene_boundary_project_position(scn, 4, pos, uv1));
- OK(sdis_scene_find_closest_point(scn, pos, INF, &iprim, uv2));
+
+ closest_pt_args.position[0] = pos[0];
+ closest_pt_args.position[1] = pos[1];
+ closest_pt_args.position[2] = pos[2];
+ closest_pt_args.radius = INF;
+ OK(sdis_scene_find_closest_point(scn, &closest_pt_args, &iprim, uv2));
CHK(d2_eq_eps(uv0, uv1, 1.e-6));
CHK(d2_eq_eps(uv1, uv2, 1.e-6));
CHK(iprim == 4);
@@ -264,6 +281,8 @@ test_scene_2d(struct sdis_device* dev, struct sdis_interface* interf)
double duplicated_vertices[] = { 0, 0, 0, 0 };
struct sdis_scene* scn = NULL;
struct sdis_scene_create_args scn_args = SDIS_SCENE_CREATE_ARGS_DEFAULT;
+ struct sdis_scene_find_closest_point_args closest_pt_args =
+ SDIS_SCENE_FIND_CLOSEST_POINT_ARGS_NULL;
struct sdis_ambient_radiative_temperature trad =
SDIS_AMBIENT_RADIATIVE_TEMPERATURE_NULL;
double lower[2], upper[2];
@@ -407,12 +426,14 @@ test_scene_2d(struct sdis_device* dev, struct sdis_interface* interf)
BA(sdis_scene_boundary_project_position(scn, 1, pos, NULL));
OK(sdis_scene_boundary_project_position(scn, 1, pos, &u1));
- BA(sdis_scene_find_closest_point(NULL, pos, INF, &iprim, &u2));
- BA(sdis_scene_find_closest_point(scn, NULL, INF, &iprim, &u2));
- BA(sdis_scene_find_closest_point(scn, pos, 0, &iprim, &u2));
- BA(sdis_scene_find_closest_point(scn, pos, INF, NULL, &u2));
- BA(sdis_scene_find_closest_point(scn, pos, INF, &iprim, NULL));
- OK(sdis_scene_find_closest_point(scn, pos, INF, &iprim, &u2));
+ closest_pt_args.position[0] = pos[0];
+ closest_pt_args.position[1] = pos[1];
+ closest_pt_args.radius = INF;
+ BA(sdis_scene_find_closest_point(NULL, &closest_pt_args, &iprim, &u2));
+ BA(sdis_scene_find_closest_point(scn, NULL, &iprim, &u2));
+ BA(sdis_scene_find_closest_point(scn, &closest_pt_args, NULL, &u2));
+ BA(sdis_scene_find_closest_point(scn, &closest_pt_args, &iprim, NULL));
+ OK(sdis_scene_find_closest_point(scn, &closest_pt_args, &iprim, &u2));
CHK(eq_eps(u0, u1, 1.e-6));
CHK(eq_eps(u1, u2, 1.e-6));
@@ -420,7 +441,10 @@ test_scene_2d(struct sdis_device* dev, struct sdis_interface* interf)
pos[0] = 0.5;
pos[1] = 0.1;
- OK(sdis_scene_find_closest_point(scn, pos, INF, &iprim, &u2));
+ closest_pt_args.position[0] = pos[0];
+ closest_pt_args.position[1] = pos[1];
+ closest_pt_args.radius = INF;
+ OK(sdis_scene_find_closest_point(scn, &closest_pt_args, &iprim, &u2));
CHK(iprim == 0);
OK(sdis_scene_boundary_project_position(scn, 0, pos, &u0));
@@ -430,7 +454,10 @@ test_scene_2d(struct sdis_device* dev, struct sdis_interface* interf)
dst = d2_len(d2_sub(pos1, pos, pos1));
CHK(eq_eps(dst, 0.1, 1.e-6));
- OK(sdis_scene_find_closest_point(scn, pos, 0.09, &iprim, &u2));
+ closest_pt_args.position[0] = pos[0];
+ closest_pt_args.position[1] = pos[1];
+ closest_pt_args.radius = 0.09;
+ OK(sdis_scene_find_closest_point(scn, &closest_pt_args, &iprim, &u2));
CHK(iprim == SDIS_PRIMITIVE_NONE);
FOR_EACH(i, 0, 64) {
@@ -438,7 +465,11 @@ test_scene_2d(struct sdis_device* dev, struct sdis_interface* interf)
OK(sdis_scene_get_boundary_position(scn, 2, &u0, pos));
OK(sdis_scene_boundary_project_position(scn, 2, pos, &u1));
- OK(sdis_scene_find_closest_point(scn, pos, INF, &iprim, &u2));
+
+ closest_pt_args.position[0] = pos[0];
+ closest_pt_args.position[1] = pos[1];
+ closest_pt_args.radius = INF;
+ OK(sdis_scene_find_closest_point(scn, &closest_pt_args, &iprim, &u2));
CHK(eq_eps(u0, u1, 1.e-6));
CHK(eq_eps(u1, u2, 1.e-6));
CHK(iprim == 2);
diff --git a/src/test_sdis_solve_probe_boundary_list.c b/src/test_sdis_solve_probe_boundary_list.c
@@ -0,0 +1,493 @@
+/* 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 "test_sdis_mesh.h"
+
+#include <star/s3dut.h>
+#include <rsys/math.h>
+
+#include <math.h>
+
+/*
+ * The system is a trilinear profile of the temperature at steady state, i.e. at
+ * each point of the system we can calculate the temperature analytically. Two
+ * forms are immersed in this temperature field: a super shape and a sphere
+ * included in the super shape. On the Monte Carlo side, the temperature is
+ * unknown everywhere except on the surface of the super shape whose
+ * temperature is defined from the aformentionned trilinear profile. We will
+ * estimate the temperature at the sphere boundary at several probe points. We
+ * should find by Monte Carlo the temperature of the trilinear profile at the
+ * position of the probe. It's the test.
+ * /\ <-- T(x,y,z)
+ * ___/ \___
+ * T(z) \ __ /
+ * | T(y) \ / \ /
+ * |/ T=? *__/ \
+ * o--- T(x) /_ __ _\
+ * \/ \/
+ */
+
+/*******************************************************************************
+ * Helper functions
+ ******************************************************************************/
+static double
+trilinear_profile(const double pos[3])
+{
+ /* Range in X, Y and Z in which the trilinear profile is defined */
+ const double lower = -4;
+ const double upper = +4;
+
+ /* Upper temperature limit in X, Y and Z [K]. Lower temperature limit is
+ * implicitly 0 */
+ const double a = 333; /* Upper temperature limit in X [K] */
+ const double b = 432; /* Upper temperature limit in Y [K] */
+ const double c = 579; /* Upper temperature limit in Z [K] */
+
+ double x, y, z;
+
+ /* Check pre-conditions */
+ CHK(pos);
+ CHK(lower <= pos[0] && pos[0] <= upper);
+ CHK(lower <= pos[1] && pos[1] <= upper);
+ CHK(lower <= pos[2] && pos[2] <= upper);
+
+ x = (pos[0] - lower) / (upper - lower);
+ y = (pos[1] - lower) / (upper - lower);
+ z = (pos[2] - lower) / (upper - lower);
+ return a*x + b*y + c*z;
+}
+
+static INLINE float
+rand_canonic(void)
+{
+ return (float)rand() / (float)((unsigned)RAND_MAX+1);
+}
+
+static void
+sample_sphere(double pos[3])
+{
+ const double phi = rand_canonic() * 2 * PI;
+ const double v = rand_canonic();
+ const double cos_theta = 1 - 2 * v;
+ const double sin_theta = 2 * sqrt(v * (1 - v));
+ pos[0] = cos(phi) * sin_theta;
+ pos[1] = sin(phi) * sin_theta;
+ pos[2] = cos_theta;
+}
+
+static INLINE void
+check_intersection
+ (const double val0,
+ const double eps0,
+ const double val1,
+ const double eps1)
+{
+ double interval0[2], interval1[2];
+ double intersection[2];
+ interval0[0] = val0 - eps0;
+ interval0[1] = val0 + eps0;
+ interval1[0] = val1 - eps1;
+ interval1[1] = val1 + eps1;
+ intersection[0] = MMAX(interval0[0], interval1[0]);
+ intersection[1] = MMIN(interval0[1], interval1[1]);
+ CHK(intersection[0] <= intersection[1]);
+}
+
+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 = 2;
+ 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.positions, sshape_data.nvertices,
+ sshape_data.indices, sshape_data.nprimitives, NULL);
+ OK(s3dut_mesh_ref_put(sshape));
+}
+
+static void
+mesh_add_sphere(struct mesh* mesh)
+{
+ struct s3dut_mesh* sphere = NULL;
+ struct s3dut_mesh_data sphere_data;
+ const double radius = 1;
+ const unsigned nslices = 128;
+
+ OK(s3dut_create_sphere(NULL, radius, nslices, nslices/2, &sphere));
+ OK(s3dut_mesh_get_data(sphere, &sphere_data));
+ mesh_append(mesh, sphere_data.positions, sphere_data.nvertices,
+ sphere_data.indices, sphere_data.nprimitives, NULL);
+ OK(s3dut_mesh_ref_put(sphere));
+}
+
+/*******************************************************************************
+ * Solid, i.e. medium of super shape and sphere
+ ******************************************************************************/
+#define SOLID_PROP(Prop, Val) \
+ static double \
+ solid_get_##Prop \
+ (const struct sdis_rwalk_vertex* vtx, \
+ struct sdis_data* data) \
+ { \
+ (void)vtx, (void)data; /* Avoid the "unused variable" warning */ \
+ return Val; \
+ }
+SOLID_PROP(calorific_capacity, 500.0) /* [J/K/Kg] */
+SOLID_PROP(thermal_conductivity, 25.0) /* [W/m/K] */
+SOLID_PROP(volumic_mass, 7500.0) /* [kg/m^3] */
+SOLID_PROP(temperature, -1/*<=> unknown*/) /* [K] */
+SOLID_PROP(delta, 1.0/20.0) /* [m] */
+
+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;
+}
+
+/*******************************************************************************
+ * Dummy environment, i.e. environment surrounding the super shape. It is
+ * defined only for Stardis compliance: in Stardis, an interface must divide 2
+ * media.
+ ******************************************************************************/
+static struct sdis_medium*
+create_dummy(struct sdis_device* sdis)
+{
+ struct sdis_fluid_shader shader = SDIS_FLUID_SHADER_NULL;
+ struct sdis_medium* dummy = NULL;
+
+ shader.calorific_capacity = dummy_medium_getter;
+ shader.volumic_mass = dummy_medium_getter;
+ shader.temperature = dummy_medium_getter;
+ OK(sdis_fluid_create(sdis, &shader, NULL, &dummy));
+ return dummy;
+}
+
+/*******************************************************************************
+ * Interfaces
+ ******************************************************************************/
+static double
+interface_get_temperature
+ (const struct sdis_interface_fragment* frag,
+ struct sdis_data* data)
+{
+ (void)data; /* Avoid the "unused variable" warning */
+ return trilinear_profile(frag->P);
+}
+
+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;
+
+ /* Fluid/solid interface: fix the temperature to the temperature profile */
+ if(sdis_medium_get_type(front) != sdis_medium_get_type(back)) {
+ shader.front.temperature = interface_get_temperature;
+ shader.back.temperature = interface_get_temperature;
+ }
+
+ OK(sdis_interface_create(sdis, front, back, &shader, NULL, &interf));
+ return interf;
+}
+
+/*******************************************************************************
+ * Scene, i.e. the system to simulate
+ ******************************************************************************/
+struct scene_context {
+ const struct mesh* mesh;
+ size_t sshape_end_id; /* Last triangle index of the super shape */
+ struct sdis_interface* sshape;
+ struct sdis_interface* sphere;
+};
+#define SCENE_CONTEXT_NULL__ {NULL, 0, 0, 0}
+static const struct scene_context SCENE_CONTEXT_NULL = SCENE_CONTEXT_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));
+ /* Flip the indices to ensure that the normal points into the super shape */
+ ids[0] = (unsigned)context->mesh->indices[itri*3+0];
+ ids[1] = (unsigned)context->mesh->indices[itri*3+2];
+ ids[2] = (unsigned)context->mesh->indices[itri*3+1];
+}
+
+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));
+ if(itri < context->sshape_end_id) {
+ *interf = context->sshape;
+ } else {
+ *interf = context->sphere;
+ }
+}
+
+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, struct scene_context* ctx)
+
+{
+ struct sdis_scene* scn = NULL;
+ struct sdis_scene_create_args scn_args = SDIS_SCENE_CREATE_ARGS_DEFAULT;
+
+ 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(ctx->mesh);
+ scn_args.nvertices = mesh_nvertices(ctx->mesh);
+ scn_args.context = ctx;
+ OK(sdis_scene_create(sdis, &scn_args, &scn));
+ return scn;
+}
+
+/*******************************************************************************
+ * Validations
+ ******************************************************************************/
+static void
+check_probe_boundary_list_api
+ (struct sdis_scene* scn,
+ const struct sdis_solve_probe_boundary_list_args* in_args)
+{
+ struct sdis_solve_probe_boundary_list_args args = *in_args;
+ struct sdis_estimator_buffer* estim_buf = NULL;
+
+ /* Check API */
+ BA(sdis_solve_probe_boundary_list(NULL, &args, &estim_buf));
+ BA(sdis_solve_probe_boundary_list(scn, NULL, &estim_buf));
+ BA(sdis_solve_probe_boundary_list(scn, &args, NULL));
+ args.nprobes = 0;
+ BA(sdis_solve_probe_boundary_list(scn, &args, &estim_buf));
+ args.nprobes = in_args->nprobes;
+ args.probes = NULL;
+ BA(sdis_solve_probe_boundary_list(scn, &args, &estim_buf));
+}
+
+/* Check the estimators against the analytical solution */
+static void
+check_estimator_buffer
+ (const struct sdis_estimator_buffer* estim_buf,
+ struct sdis_scene* scn,
+ const struct sdis_solve_probe_boundary_list_args* args)
+{
+ struct sdis_mc T = SDIS_MC_NULL;
+ size_t iprobe = 0;
+
+ /* Variables used to check the global estimation results */
+ size_t total_nrealisations = 0;
+ size_t total_nrealisations_ref = 0;
+ size_t total_nfailures = 0;
+ size_t total_nfailures_ref = 0;
+ double sum = 0; /* Global sum of weights */
+ double sum2 = 0; /* Global sum of squared weights */
+ double N = 0; /* Number of (successful) realisations */
+ double E = 0; /* Expected value */
+ double V = 0; /* Variance */
+ double SE = 0; /* Standard Error */
+
+ /* Check the results */
+ FOR_EACH(iprobe, 0, args->nprobes) {
+ const struct sdis_estimator* estimator = NULL;
+ size_t probe_nrealisations = 0;
+ size_t probe_nfailures = 0;
+ double pos[3] = {0,0,0};
+ double probe_sum = 0;
+ double probe_sum2 = 0;
+ double ref = 0;
+
+ OK(sdis_scene_get_boundary_position(scn, args->probes[iprobe].iprim,
+ args->probes[iprobe].uv, pos));
+
+ /* Fetch result */
+ OK(sdis_estimator_buffer_at(estim_buf, iprobe, 0, &estimator));
+ OK(sdis_estimator_get_temperature(estimator, &T));
+ OK(sdis_estimator_get_realisation_count(estimator, &probe_nrealisations));
+ OK(sdis_estimator_get_failure_count(estimator, &probe_nfailures));
+
+ /* Check probe estimation */
+ ref = trilinear_profile(pos);
+ printf("T(%g, %g, %g) = %g ~ %g +/- %g\n", SPLIT3(pos), ref, T.E, T.SE);
+ CHK(eq_eps(ref, T.E, 3*T.SE));
+
+ /* Check miscellaneous results */
+ CHK(probe_nrealisations == args->probes[iprobe].nrealisations);
+
+ /* Compute the sum of weights and sum of squared weights of the probe */
+ probe_sum = T.E * (double)probe_nrealisations;
+ probe_sum2 = (T.V + T.E*T.E) * (double)probe_nrealisations;
+
+ /* Update the global estimation */
+ total_nrealisations_ref += probe_nrealisations;
+ total_nfailures_ref += probe_nfailures;
+ sum += probe_sum;
+ sum2 += probe_sum2;
+ }
+
+ /* Check the overall estimate of the estimate buffer. This estimate is not
+ * really a result expected by the caller since the probes are all
+ * independent. But to be consistent with the estimate buffer API, we need to
+ * provide it. And so we check its validity */
+ OK(sdis_estimator_buffer_get_temperature(estim_buf, &T));
+ OK(sdis_estimator_buffer_get_realisation_count(estim_buf, &total_nrealisations));
+ OK(sdis_estimator_buffer_get_failure_count(estim_buf, &total_nfailures));
+
+ CHK(total_nrealisations == total_nrealisations_ref);
+ CHK(total_nfailures == total_nfailures_ref);
+
+ N = (double)total_nrealisations_ref;
+ E = sum / N;
+ V = sum2 / N - E*E;
+ SE = sqrt(V/N);
+ check_intersection(E, SE, T.E, T.SE);
+}
+
+static void
+check_probe_boundary_list(struct sdis_scene* scn, const int is_master_process)
+{
+ #define NPROBES 10
+
+ /* Estimations */
+ struct sdis_estimator_buffer* estim_buf = NULL;
+
+ /* Probe variables */
+ struct sdis_solve_probe_boundary_args probes[NPROBES];
+ struct sdis_solve_probe_boundary_list_args args =
+ SDIS_SOLVE_PROBE_BOUNDARY_LIST_ARGS_DEFAULT;
+ size_t iprobe;
+
+ /* Miscellaneous */
+ struct sdis_scene_find_closest_point_args closest_pt_args =
+ SDIS_SCENE_FIND_CLOSEST_POINT_ARGS_NULL;
+
+ (void)is_master_process;
+
+ /* Setup the list of probes to calculate */
+ args.probes = probes;
+ args.nprobes = NPROBES;
+ FOR_EACH(iprobe, 0, NPROBES) {
+ sample_sphere(closest_pt_args.position);
+ closest_pt_args.radius = INF;
+
+ probes[iprobe] = SDIS_SOLVE_PROBE_BOUNDARY_ARGS_DEFAULT;
+ probes[iprobe].nrealisations = 10000;
+ probes[iprobe].side = SDIS_FRONT;
+
+ OK(sdis_scene_find_closest_point
+ (scn, &closest_pt_args, &probes[iprobe].iprim, probes[iprobe].uv));
+ }
+
+ check_probe_boundary_list_api(scn, &args);
+
+ /* Solve the probes */
+ OK(sdis_solve_probe_boundary_list(scn, &args, &estim_buf));
+ if(!is_master_process) {
+ CHK(estim_buf == NULL);
+ } else {
+ check_estimator_buffer(estim_buf, scn, &args);
+ OK(sdis_estimator_buffer_ref_put(estim_buf));
+ }
+
+ #undef NPROBES
+}
+
+/*******************************************************************************
+ * The test
+ ******************************************************************************/
+int
+main(int argc, char** argv)
+{
+ /* Stardis */
+ struct sdis_device* dev = NULL;
+ struct sdis_interface* solid_fluid = NULL;
+ struct sdis_interface* solid_solid = NULL;
+ struct sdis_medium* fluid = NULL;
+ struct sdis_medium* solid = NULL;
+ struct sdis_scene* scn = NULL;
+
+ /* Miscellaneous */
+ struct scene_context ctx = SCENE_CONTEXT_NULL;
+ struct mesh mesh = MESH_NULL;
+ size_t sshape_end_id = 0; /* Last index of the super shape */
+ int is_master_process = 1;
+ (void)argc, (void)argv;
+
+ create_default_device(&argc, &argv, &is_master_process, &dev);
+
+ /* Setup the mesh */
+ mesh_init(&mesh);
+ mesh_add_super_shape(&mesh);
+ sshape_end_id = mesh_ntriangles(&mesh);
+ mesh_add_sphere(&mesh);
+
+ /* Setup physical properties */
+ fluid = create_dummy(dev);
+ solid = create_solid(dev);
+ solid_fluid = create_interface(dev, solid, fluid);
+ solid_solid = create_interface(dev, solid, solid);
+
+ /* Create the scene */
+ ctx.mesh = &mesh;
+ ctx.sshape_end_id = sshape_end_id;
+ ctx.sshape = solid_fluid;
+ ctx.sphere = solid_solid;
+ scn = create_scene(dev, &ctx);
+
+ check_probe_boundary_list(scn, is_master_process);
+
+ mesh_release(&mesh);
+ OK(sdis_interface_ref_put(solid_fluid));
+ OK(sdis_interface_ref_put(solid_solid));
+ OK(sdis_medium_ref_put(fluid));
+ OK(sdis_medium_ref_put(solid));
+ OK(sdis_scene_ref_put(scn));
+
+ free_default_device(dev);
+
+ CHK(mem_allocated_size() == 0);
+ return 0;
+}
diff --git a/src/test_sdis_source.c b/src/test_sdis_source.c
@@ -26,7 +26,16 @@ spherical_source_get_position
struct sdis_data* data)
{
(void)time, (void)data;
- pos[0] = pos[1] = pos[2] = 1.234;
+ pos[0] = pos[1] = pos[2] = 1.234; /* [m] */
+}
+
+static double
+spherical_source_get_power
+ (const double time,
+ struct sdis_data* data)
+{
+ (void)time, (void)data;
+ return 10; /* [W] */
}
static void
@@ -36,25 +45,21 @@ check_spherical_source(struct sdis_device* dev)
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.power = spherical_source_get_power;
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);
+ CHK(sdis_source_get_power(src, INF) == 10);
BA(sdis_source_ref_get(NULL));
OK(sdis_source_ref_get(src));
@@ -67,6 +72,12 @@ check_spherical_source(struct sdis_device* dev)
args.data = NULL;
OK(sdis_spherical_source_create(dev, &args, &src));
OK(sdis_source_ref_put(src));
+
+ args.position = NULL;
+ BA(sdis_spherical_source_create(dev, &args, &src));
+ args.position = spherical_source_get_position;
+ args.power = NULL;
+ BA(sdis_spherical_source_create(dev, &args, &src));
}
/*******************************************************************************
diff --git a/src/test_sdis_utils.c b/src/test_sdis_utils.c
@@ -99,7 +99,6 @@ solve_green_path(struct sdis_green_path* path, void* ctx)
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] */
@@ -139,8 +138,9 @@ solve_green_path(struct sdis_green_path* path, void* ctx)
if(source == NULL) {
CHK(external_flux == 0);
} else {
- OK(sdis_source_get_power(source, &source_power));
- external_flux *= source_power;
+ /* NOTE: source power is assumed constant in time and is therefore retrieved
+ * at steady state*/
+ external_flux *= sdis_source_get_power(source, INF); /* [W] */
}
BA(sdis_green_path_get_end_type(NULL, NULL));
@@ -498,4 +498,3 @@ check_green_serialization
OK(sdis_estimator_ref_put(e2));
OK(sdis_green_function_ref_put(green2));
}
-
