commit e8ae1bdf25bfd2138bafb004a02f856f2292e781
parent ecc997121cc8e2cde3354f4593ffc2daaebaecf8
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Wed, 4 Jul 2018 15:18:56 +0200
Merge remote-tracking branch 'origin/develop' into feature_convection
Diffstat:
34 files changed, 3689 insertions(+), 654 deletions(-)
diff --git a/README.md b/README.md
@@ -23,13 +23,70 @@ variable the install directories of its dependencies.
## Release notes
+### Version 0.4
+
+Full rewrite of how the volumetric power is taken into account.
+
+- Change the scheme of the random walk "solid re-injection": use a 2D
+ re-injection scheme in order to handle 2D effects. On one hand, this scheme
+ drastically improves the accuracy of the temperature estimation in solid with
+ a volumetric power term. On the other hand it is more sensible to numerical
+ imprecisions. The previous 1D scheme is thus used in situations where the 2D
+ scheme exhibits too numerical issues, i.e. on sharp angles.
+- Add the missing volumetric power term on solid re-injection.
+- Add a corrective term to fix the bias on the volumetric power introduced when
+ the random walk progresses at a distance of `delta` of a boundary.
+- Add several volumetric power tests.
+- Remove the `delta_boundary` parameter of the `struct sdis_solid_shader` data
+ structure.
+
+### Version 0.3
+
+- Some interface properties become double sided: the temperature, emissivity
+ and specular fraction is defined for each side of the interface. Actually,
+ only the convection coefficient is shared by the 2 sides of the interface.
+ The per side interface properties are grouped into the new `struct
+ sdis_interface_side_shader` data structure.
+- Add the support of fixed fluxes: the flux is a per side interface property.
+ Currently, the flux is handled only for the interface sides facing a solid
+ medium.
+- Add the `sdis_scene_boundary_project_pos` function that computes the
+ parametric coordinates of a world space position projected onto a given
+ primitive with respect to its normal. If the projection lies outside the
+ primitive, its parametric coordinates are wrapped against its boundaries in
+ order to ensure that they are valid coordinates into the primitive. Actually,
+ this function was mainly added to help in the definition of the probe
+ position onto a boundary as expected by the
+ `sdis_solve_probe_boundary` function.
+- Update the default comportment of the interface shader when a function is not
+ set.
+- Rename the `SDIS_MEDIUM_<FLUID|SOLID>` constants in `SDIS_<FLUID|SOLID>`.
+- Rename the `enum sdis_side_flag` enumerate in `enum sdis_side` and update its
+ values.
+
+### Version 0.2
+
+- Add the support of volumic power to solid media: add the `volumic_power`
+ functor to the `sdis_solid_shader` data structure that, once defined, should
+ return the volumic power of the solid at a specific position and time. On
+ solve invocation, the conductive random walks take into account this
+ spatio-temporal volumic power in the computation of the solid temperature.
+- Add the `sdis_solve_probe_boundary` function: it computes the temperature at
+ a given position and time onto a geometric primitive. The probe position is
+ defined by the index of the primitive and a parametric coordinates onto it.
+- Add the `sdis_scene_get_boundary_position` function: it computes a world
+ space position from the index of a geometric primitive and a parametric
+ coordinate onto it.
+- Fix how the `sdis_solve_probe` was parallelised. The submitted `threads_hint`
+ parameter was not correctly handled.
+
### Version 0.1
- Add the support of radiative temperature.
-- Add the `sdis_camera` API : it defines a pinhole camera into the scene.
-- Add the `sdis_accum_buffer` API : it is a pool of MC accumulators, i.e. a sum
+- Add the `sdis_camera` API: it defines a pinhole camera into the scene.
+- Add the `sdis_accum_buffer` API: it is a pool of MC accumulators, i.e. a sum
of MC weights and square weights.
-- Add the `sdis_solve_camera` function : it relies on a `sdis_camera` and a
+- Add the `sdis_solve_camera` function: it relies on a `sdis_camera` and a
`sdis_accum_buffer` to compute the radiative temperature that reaches each
pixel of an image whose definition is defined by the caller. Note that
actually this function uses the same underlying MC algorithm behind the
diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt
@@ -17,9 +17,14 @@ cmake_minimum_required(VERSION 3.0)
project(stardis C)
enable_testing()
+include(CMakeDependentOption)
+
set(SDIS_SOURCE_DIR ${PROJECT_SOURCE_DIR}/../src)
option(NO_TEST "Do not build tests" OFF)
+CMAKE_DEPENDENT_OPTION(ALL_TESTS
+ "Perform basic and advanced tests" OFF "NOT NO_TEST" OFF)
+
################################################################################
# Check dependencies
################################################################################
@@ -42,7 +47,7 @@ rcmake_append_runtime_dirs(_runtime_dirs RSys Star3D StarSP StarEnc StarEnc2D)
# Configure and define targets
################################################################################
set(VERSION_MAJOR 0)
-set(VERSION_MINOR 1)
+set(VERSION_MINOR 4)
set(VERSION_PATCH 0)
set(VERSION ${VERSION_MAJOR}.${VERSION_MINOR}.${VERSION_PATCH})
@@ -78,14 +83,24 @@ rcmake_prepend_path(SDIS_FILES_INC ${SDIS_SOURCE_DIR})
rcmake_prepend_path(SDIS_FILES_INC_API ${SDIS_SOURCE_DIR})
rcmake_prepend_path(SDIS_FILES_DOC ${PROJECT_SOURCE_DIR}/../)
+if(CMAKE_COMPILER_IS_GNUCC)
+ set(MATH_LIB m)
+endif()
+
+if(MSVC)
+ ### disable verbose warnings:
+ # warning C4127: conditional expression is constant
+ set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} /wd4127")
+ # warning C4938: Floating point reduction variable may cause inconsistent
+ # results under /fp:strict or #pragma fenv_access
+ set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} /wd4938")
+endif()
+
add_library(sdis SHARED
${SDIS_FILES_SRC}
${SDIS_FILES_INC}
${SDIS_FILES_INC_API})
-target_link_libraries(sdis RSys Star2D Star3D StarEnc StarEnc2D StarSP m)
-if(CMAKE_COMPILER_IS_GNUCC)
- target_link_libraries(sdis m)
-endif()
+target_link_libraries(sdis RSys Star2D Star3D StarSP StarEnc StarEnc2D ${MATH_LIB})
set_target_properties(sdis PROPERTIES
DEFINE_SYMBOL SDIS_SHARED_BUILD
@@ -95,7 +110,7 @@ set_target_properties(sdis PROPERTIES
rcmake_copy_runtime_libraries(sdis)
if(CMAKE_COMPILER_IS_GNUCC)
- set_target_properties(sdis PROPERTIES LINK_FLAGS "${OpenMP_C_FLAGS} -lm")
+ set_target_properties(sdis PROPERTIES LINK_FLAGS "${OpenMP_C_FLAGS}")
endif()
rcmake_setup_devel(sdis Stardis ${VERSION} sdis_version.h)
@@ -128,6 +143,7 @@ if(NOT NO_TEST)
new_test(test_sdis_convection)
new_test(test_sdis_data)
new_test(test_sdis_device)
+ new_test(test_sdis_flux)
new_test(test_sdis_interface)
new_test(test_sdis_medium)
new_test(test_sdis_scene)
@@ -141,12 +157,23 @@ if(NOT NO_TEST)
new_test(test_sdis_solve_probe_boundary)
new_test(test_sdis_volumic_power)
+ # Additionnal tests
+ build_test(test_sdis_volumic_power2)
+ build_test(test_sdis_volumic_power2_2d)
+ build_test(test_sdis_volumic_power3_2d)
+ build_test(test_sdis_volumic_power4_2d)
+
+ if(ALL_TESTS)
+ add_test(test_sdis_volumic_power2 test_sdis_volumic_power2)
+ add_test(test_sdis_volumic_power2_2d test_sdis_volumic_power2_2d)
+ add_test(test_sdis_volumic_power3_2d test_sdis_volumic_power3_2d)
+ add_test(test_sdis_volumic_power4_2d test_sdis_volumic_power4_2d)
+ endif()
+
target_link_libraries(test_sdis_solve_probe3 Star3DUT)
+ target_link_libraries(test_sdis_solve_probe3_2d ${MATH_LIB})
target_link_libraries(test_sdis_solve_camera Star3DUT)
- if(CMAKE_COMPILER_IS_GNUCC)
- target_link_libraries(test_sdis_solve_probe3_2d m)
- endif()
endif()
################################################################################
diff --git a/src/sdis.h b/src/sdis.h
@@ -17,6 +17,7 @@
#define SDIS_H
#include <rsys/rsys.h>
+#include <float.h>
/* Library symbol management */
#if defined(SDIS_SHARED_BUILD)
@@ -40,6 +41,9 @@
* as CPU cores */
#define SDIS_NTHREADS_DEFAULT (~0u)
+#define SDIS_VOLUMIC_POWER_NONE DBL_MAX /* <=> No volumic power */
+#define SDIS_FLUX_NONE DBL_MAX /* <=> No flux */
+
/* Forward declaration of external opaque data types */
struct logger;
struct mem_allocator;
@@ -60,15 +64,15 @@ struct sdis_interface;
struct sdis_medium;
struct sdis_scene;
-enum sdis_side_flag {
- SDIS_FRONT = BIT(0),
- SDIS_BACK = BIT(1),
- SDIS_SIDE_NULL__ = BIT(2)
+enum sdis_side {
+ SDIS_FRONT,
+ SDIS_BACK,
+ SDIS_SIDE_NULL__
};
enum sdis_medium_type {
- SDIS_MEDIUM_FLUID,
- SDIS_MEDIUM_SOLID,
+ SDIS_FLUID,
+ SDIS_SOLID,
SDIS_MEDIUM_TYPES_COUNT__
};
@@ -82,22 +86,27 @@ struct sdis_rwalk_vertex {
static const struct sdis_rwalk_vertex SDIS_RWALK_VERTEX_NULL =
SDIS_RWALK_VERTEX_NULL__;
-/* Spatiotemporal position onto an interface */
+/* Spatiotemporal position onto an interface. As a random walk vertex, it
+ * stores the position and time of the random walk, but since it lies onto an
+ * interface, it has additionnal parameters as the normal of the interface and
+ * the parametric coordinate of the position onto the interface */
struct sdis_interface_fragment {
double P[3]; /* World space position */
double Ng[3]; /* Normalized world space geometry normal at the interface */
double uv[2]; /* Parametric coordinates of the interface */
double time; /* Current time */
+ enum sdis_side side;
};
-#define SDIS_INTERFACE_FRAGMENT_NULL__ {{0}, {0}, {0}, -1}
+#define SDIS_INTERFACE_FRAGMENT_NULL__ {{0}, {0}, {0}, -1, SDIS_SIDE_NULL__}
static const struct sdis_interface_fragment SDIS_INTERFACE_FRAGMENT_NULL =
SDIS_INTERFACE_FRAGMENT_NULL__;
/* Monte-Carlo accumulator */
struct sdis_accum {
- double sum_weights; /* Sum of Monte-Carlo weight */
+ double sum_weights; /* Sum of Monte-Carlo weights */
double sum_weights_sqr; /* Sum of Monte-Carlo square weights */
size_t nweights; /* #accumulated weights */
+ size_t nfailures; /* #failures */
};
/* Monte-Carlo estimation */
@@ -109,40 +118,46 @@ struct sdis_mc {
#define SDIS_MC_NULL__ {0, 0, 0}
static const struct sdis_mc SDIS_MC_NULL = SDIS_MC_NULL__;
-/* Functor type to retrieve the medium properties. */
+/* Functor type used to retrieve the spatio temporal physical properties of a
+ * medium. */
typedef double
(*sdis_medium_getter_T)
(const struct sdis_rwalk_vertex* vert,
struct sdis_data* data);
-/* Functor type to retrieve the interface properties. */
+/* 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);
+/* Define the physical properties of a solid */
struct sdis_solid_shader {
/* Properties */
- sdis_medium_getter_T calorific_capacity;
- sdis_medium_getter_T thermal_conductivity;
- sdis_medium_getter_T volumic_mass;
+ sdis_medium_getter_T calorific_capacity; /* In J.K^-1.kg^-1 */
+ sdis_medium_getter_T thermal_conductivity; /* In W.m^-1.K^-1 */
+ sdis_medium_getter_T volumic_mass; /* In kg.m^-3 */
sdis_medium_getter_T delta_solid;
- sdis_medium_getter_T delta_boundary;
- sdis_medium_getter_T volumic_power; /* May be NULL <=> no volumic power */
+ /* May be NULL if there is no volumic power. One can also return
+ * SDIS_VOLUMIC_POWER_NONE to define that there is no volumic power at the
+ * submitted position and time */
+ sdis_medium_getter_T volumic_power; /* In W.m^-3 */
/* Initial/limit condition. A temperature < 0 means that the temperature is
* unknown for the submitted random walk vertex. */
sdis_medium_getter_T temperature;
};
-#define SDIS_SOLID_SHADER_NULL__ {NULL, NULL, NULL, NULL, NULL, NULL, NULL}
+#define SDIS_SOLID_SHADER_NULL__ {NULL, NULL, NULL, NULL, NULL, NULL}
static const struct sdis_solid_shader SDIS_SOLID_SHADER_NULL =
SDIS_SOLID_SHADER_NULL__;
+/* Define the physical properties of a fluid */
struct sdis_fluid_shader {
/* Properties */
- sdis_medium_getter_T calorific_capacity;
- sdis_medium_getter_T volumic_mass;
+ sdis_medium_getter_T calorific_capacity; /* In J.K^-1.kg^-1 */
+ sdis_medium_getter_T volumic_mass; /* In kg.m^-3 */
/* Initial/limit condition. A temperature < 0 means that the temperature is
* unknown for the submitted position and time. */
@@ -152,15 +167,33 @@ struct sdis_fluid_shader {
static const struct sdis_fluid_shader SDIS_FLUID_SHADER_NULL =
SDIS_FLUID_SHADER_NULL__;
+/* Define the physical properties of one side of an interface. */
+struct sdis_interface_side_shader {
+ /* Fixed temperature/flux. May be NULL if the temperature/flux is unknown
+ * onto the whole interface */
+ sdis_interface_getter_T temperature; /* In Kelvin. < 0 <=> Unknown temp */
+ sdis_interface_getter_T flux; /* In W.m^-2. SDIS_FLUX_NONE <=> no flux */
+
+ /* Control the emissivity of the interface. May be NULL for solid/solid
+ * interface or if the emissivity is 0 onto the whole interface. */
+ sdis_interface_getter_T emissivity; /* Overall emissivity. */
+ sdis_interface_getter_T specular_fraction; /* Specular part in [0,1] */
+};
+#define SDIS_INTERFACE_SIDE_SHADER_NULL__ { NULL, NULL, NULL, NULL }
+static const struct sdis_interface_side_shader SDIS_INTERFACE_SIDE_SHADER_NULL =
+ SDIS_INTERFACE_SIDE_SHADER_NULL__;
+
+/* Define the physical properties of an interface between 2 media .*/
struct sdis_interface_shader {
- sdis_interface_getter_T temperature; /* Limit condition. NULL <=> Unknown */
- sdis_interface_getter_T convection_coef; /* May be NULL for solid/solid */
+ /* May be NULL for solid/solid or if the convection coefficient is 0 onto
+ * the whole interface. */
+ sdis_interface_getter_T convection_coef; /* In W.K^-1.m^-2 */
- /* Interface emssivity. May be NULL for solid/solid interface */
- sdis_interface_getter_T emissivity; /* Overall emissivity */
- sdis_interface_getter_T specular_fraction; /* Specular fraction in [0, 1] */
+ struct sdis_interface_side_shader front;
+ struct sdis_interface_side_shader back;
};
-#define SDIS_INTERFACE_SHADER_NULL__ {NULL, NULL, NULL, NULL}
+#define SDIS_INTERFACE_SHADER_NULL__ \
+ {NULL, SDIS_INTERFACE_SIDE_SHADER_NULL__, SDIS_INTERFACE_SIDE_SHADER_NULL__}
static const struct sdis_interface_shader SDIS_INTERFACE_SHADER_NULL =
SDIS_INTERFACE_SHADER_NULL__;
@@ -204,7 +237,8 @@ sdis_device_ref_put
/*******************************************************************************
* A data stores in the Stardis memory space a set of user defined data. It can
- * be seen as a ref counted memory space allocated by Stardis.
+ * be seen as a ref counted memory space allocated by Stardis. It is used to
+ * attach user data to the media and to the interfaces.
******************************************************************************/
SDIS_API res_T
sdis_data_create
@@ -265,13 +299,15 @@ sdis_camera_look_at
const double up[3]);
/*******************************************************************************
- * A buffer of accumulations
+ * A buffer of accumulations is a 2D array whose each cell stores an
+ * Monte-Carlo accumulation, i.e. a sum of MC weights, the sum of their square
+ * and the overall number of summed weights (see struct sdis_accum)
******************************************************************************/
SDIS_API res_T
sdis_accum_buffer_create
(struct sdis_device* dev,
- const size_t width,
- const size_t height,
+ const size_t width, /* #cells in X */
+ const size_t height, /* #cells in Y */
struct sdis_accum_buffer** buf);
SDIS_API res_T
@@ -287,6 +323,7 @@ sdis_accum_buffer_get_layout
(const struct sdis_accum_buffer* buf,
struct sdis_accum_buffer_layout* layout);
+/* Get a read only pointer toward the memory space of the accum buffer */
SDIS_API res_T
sdis_accum_buffer_map
(const struct sdis_accum_buffer* buf,
@@ -296,7 +333,11 @@ SDIS_API res_T
sdis_accum_buffer_unmap
(const struct sdis_accum_buffer* buf);
-/* Helper function that matches the `sdis_write_accums_T' functor type */
+/* Helper function that matches the `sdis_write_accums_T' functor type. On can
+ * send this function directly to the sdis_solve_camera function, to fill the
+ * accum buffer with the estimation of the radiative temperature that reaches
+ * each pixel of an image whose definition matches the definition of the accum
+ * buffer. */
SDIS_API res_T
sdis_accum_buffer_write
(void* buf, /* User data */
@@ -333,14 +374,18 @@ SDIS_API enum sdis_medium_type
sdis_medium_get_type
(const struct sdis_medium* medium);
+SDIS_API struct sdis_data*
+sdis_medium_get_data
+ (struct sdis_medium* medium);
+
/*******************************************************************************
- * An interface is the boundary between 2 mediums.
+ * An interface is the boundary between 2 media.
******************************************************************************/
SDIS_API res_T
sdis_interface_create
(struct sdis_device* dev,
- struct sdis_medium* front,
- struct sdis_medium* back,
+ struct sdis_medium* front, /* Medium on the front side of the geometry */
+ struct sdis_medium* back, /* Medium on the back side of the geometry */
const struct sdis_interface_shader* shader,
struct sdis_data* data, /* Data sent to the shader. May be NULL */
struct sdis_interface** interf);
@@ -357,6 +402,18 @@ sdis_interface_ref_put
* A scene is a collection of primitives. Each primitive is the geometric
* support of the interface between 2 mediums.
******************************************************************************/
+/* Create a 3D scene. The geometry of the scene is defined by an indexed
+ * triangular mesh: each triangle is composed of 3 indices where each index
+ * references an absolute 3D position. The physical properties of an interface
+ * is defined by the interface of the triangle.
+ *
+ * Note that each triangle has 2 sides: a front and a back side. By convention,
+ * the front side of a triangle is the side where its vertices are clock wise
+ * ordered. The back side of a triangle is the exact opposite: it is the side
+ * where the triangle vertices are counter-clock wise ordered. The front and
+ * back media of a triangle interface directly refer to this convention and
+ * thus one has to take care of how the triangle vertices are defined to ensure
+ * that the front and the back media are correctly defined wrt the geometry. */
SDIS_API res_T
sdis_scene_create
(struct sdis_device* dev,
@@ -371,6 +428,18 @@ sdis_scene_create
void* ctx, /* Client side data sent as input of the previous callbacks */
struct sdis_scene** scn);
+/* Create a 2D scene. The geometry of the 2D scene is defined by an indexed
+ * line segments: each segment is composed of 2 indices where each index
+ * references an absolute 2D position. The physical properties of an interface
+ * is defined by the interface of the segment.
+ *
+ * Note that each segment has 2 sides: a front and a back side. By convention,
+ * the front side of a segment is the side where its vertices are clock wise
+ * ordered. The back side of a segment is the exact opposite: it is the side
+ * where the segment vertices are counter-clock wise ordered. The front and
+ * back media of a segment interface directly refer to this convention and
+ * thus one has to take care of how the segment vertices are defined to ensure
+ * that the front and the back media are correctly defined wrt the geometry. */
SDIS_API res_T
sdis_scene_2d_create
(struct sdis_device* dev,
@@ -400,13 +469,49 @@ sdis_scene_get_aabb
double lower[3],
double upper[3]);
+/* Define the world space position of a point onto the primitive `iprim' whose
+ * parametric coordinate is uv. */
SDIS_API res_T
sdis_scene_get_boundary_position
(const struct sdis_scene* scn,
const size_t iprim, /* Primitive index */
- const double uv[2], /* Parametric coordinate onto the pimitive */
+ const double uv[2], /* Parametric coordinate onto the primitive */
double pos[3]); /* World space position */
+/* Project a world space position onto a primitive wrt its normal and compute
+ * the parametric coordinates of the projected point onto the primitive. This
+ * function may help to define the probe position onto a boundary as expected
+ * by the sdis_solve_probe_boundary function.
+ *
+ * Note that the projected point can lie outside the submitted primitive. In
+ * this case, the parametric coordinates are clamped against the primitive
+ * boundaries in order to ensure that the returned parametric coordinates are
+ * valid according to the primitive. To ensure this, in 2D, the parametric
+ * coordinate is simply clamped to [0, 1]. In 3D, the `uv' coordinates are
+ * clamped against the triangle edges. For instance, let the
+ * following triangle whose vertices are `a', `b' and `c':
+ * , ,
+ * , B ,
+ * , ,
+ * b E1
+ * E0 / \ ,P
+ * / \,*'
+ * / \
+ * ....a-------c......
+ * ' '
+ * A ' E2 ' C
+ * ' '
+ * The projected point `P' is orthogonally wrapped to the edge `ab', `bc' or
+ * `ca' if it lies in the `E0', `E1' or `E2' region, respectively. If `P' is in
+ * the `A', `B' or `C' region, then it is taken back to the `a', `b' or `c'
+ * vertex, respectively. */
+SDIS_API res_T
+sdis_scene_boundary_project_position
+ (const struct sdis_scene* scn,
+ const size_t iprim,
+ const double pos[3],
+ double uv[]);
+
/*******************************************************************************
* An estimator stores the state of a simulation
******************************************************************************/
@@ -454,6 +559,7 @@ sdis_solve_probe_boundary
const size_t iprim, /* Identifier of the primitive on which the probe lies */
const double uv[2], /* Parametric coordinates of the probe onto the primitve */
const double time, /* Observation time */
+ const enum sdis_side side, /* Side of iprim on which the probe lies */
const double fp_to_meter, /* Scale from floating point units to meters */
const double ambient_radiative_temperature, /* In Kelvin */
const double reference_temperature, /* In Kelvin */
diff --git a/src/sdis_interface.c b/src/sdis_interface.c
@@ -29,26 +29,52 @@
******************************************************************************/
static int
check_interface_shader
- (const struct sdis_interface_shader* shader,
+ (struct sdis_device* dev,
+ const char* caller_name,
+ const struct sdis_interface_shader* shader,
const struct sdis_medium* front,
const struct sdis_medium* back)
{
- enum sdis_medium_type type0;
- enum sdis_medium_type type1;
- ASSERT(shader && front && back);
+ enum sdis_medium_type type[2];
+ const struct sdis_interface_side_shader* shaders[2];
+ int i;
+ ASSERT(dev && caller_name && shader && front && back);
- type0 = sdis_medium_get_type(front);
- type1 = sdis_medium_get_type(back);
+ type[0] = sdis_medium_get_type(front);
+ type[1] = sdis_medium_get_type(back);
+ shaders[0] = &shader->front;
+ shaders[1] = &shader->back;
/* Fluid<->solid interface */
- if(type0 != type1) {
- if(shader->convection_coef == NULL
- || shader->emissivity == NULL
- || shader->specular_fraction == NULL) {
- return 0;
- }
+ if(type[0] == SDIS_SOLID
+ && type[1] == SDIS_SOLID
+ && shader->convection_coef) {
+ log_warn(dev,
+ "%s: a solid/solid interface can't have a convection coefficient. This "
+ "function of the interface shader should be NULL.\n", caller_name);
}
+ FOR_EACH(i, 0, 2) {
+ switch(type[i]) {
+ case SDIS_SOLID:
+ if(shaders[i]->emissivity || shaders[i]->specular_fraction) {
+ log_warn(dev,
+ "%s: the interface side toward a solid can't have the emissivity "
+ "and specular_fraction properties. The shader functions that return "
+ "these attributes should be NULL.\n", caller_name);
+ }
+ break;
+ case SDIS_FLUID:
+ if(shaders[i]->flux) {
+ log_warn(dev,
+ "%s: the interface side toward a fluid can't have a flux property. "
+ "The shader function that returns this attribute should be NULL.\n",
+ caller_name);
+ }
+ break;
+ default: FATAL("Unreachable code.\n"); break;
+ }
+ }
return 1;
}
@@ -88,14 +114,14 @@ sdis_interface_create
goto error;
}
- if(sdis_medium_get_type(front) == SDIS_MEDIUM_FLUID
- && sdis_medium_get_type(back) == SDIS_MEDIUM_FLUID) {
+ if(sdis_medium_get_type(front) == SDIS_FLUID
+ && sdis_medium_get_type(back) == SDIS_FLUID) {
log_err(dev, "%s: invalid fluid<->fluid interface.\n", FUNC_NAME);
res = RES_BAD_ARG;
goto error;
}
- if(!check_interface_shader(shader, front, back)) {
+ if(!check_interface_shader(dev, FUNC_NAME, shader, front, back)) {
log_err(dev, "%s: invalid interface shader.\n", FUNC_NAME);
res = RES_BAD_ARG;
goto error;
@@ -154,7 +180,7 @@ sdis_interface_ref_put(struct sdis_interface* interf)
******************************************************************************/
const struct sdis_medium*
interface_get_medium
- (const struct sdis_interface* interf, const enum sdis_side_flag side)
+ (const struct sdis_interface* interf, const enum sdis_side side)
{
struct sdis_medium* mdm = NULL;
ASSERT(interf);
@@ -177,27 +203,33 @@ void
setup_interface_fragment_2d
(struct sdis_interface_fragment* frag,
const struct sdis_rwalk_vertex* vertex,
- const struct s2d_hit* hit)
+ const struct s2d_hit* hit,
+ const enum sdis_side side)
{
ASSERT(frag && vertex && hit && !S2D_HIT_NONE(hit));
+ ASSERT(side == SDIS_FRONT || side == SDIS_BACK);
d2_set(frag->P, vertex->P);
frag->P[2] = 0;
d2_normalize(frag->Ng, d2_set_f2(frag->Ng, hit->normal));
frag->Ng[2] = 0;
frag->uv[0] = hit->u;
frag->time = vertex->time;
+ frag->side = side;
}
void
setup_interface_fragment_3d
(struct sdis_interface_fragment* frag,
const struct sdis_rwalk_vertex* vertex,
- const struct s3d_hit* hit)
+ const struct s3d_hit* hit,
+ const enum sdis_side side)
{
ASSERT(frag && vertex && hit && !S3D_HIT_NONE(hit));
+ ASSERT(side == SDIS_FRONT || side == SDIS_BACK);
d3_set(frag->P, vertex->P);
d3_normalize(frag->Ng, d3_set_f3(frag->Ng, hit->normal));
d2_set_f2(frag->uv, hit->uv);
frag->time = vertex->time;
+ frag->side = side;
}
diff --git a/src/sdis_interface_c.h b/src/sdis_interface_c.h
@@ -39,7 +39,7 @@ struct sdis_interface {
extern LOCAL_SYM const struct sdis_medium*
interface_get_medium
(const struct sdis_interface* interf,
- const enum sdis_side_flag side);
+ const enum sdis_side side);
extern LOCAL_SYM unsigned
interface_get_id
@@ -49,49 +49,85 @@ extern LOCAL_SYM void
setup_interface_fragment_2d
(struct sdis_interface_fragment* frag,
const struct sdis_rwalk_vertex* vertex,
- const struct s2d_hit* hit);
+ const struct s2d_hit* hit,
+ const enum sdis_side side);
extern LOCAL_SYM void
setup_interface_fragment_3d
(struct sdis_interface_fragment* frag,
const struct sdis_rwalk_vertex* vertex,
- const struct s3d_hit* hit);
+ const struct s3d_hit* hit,
+ const enum sdis_side side);
static INLINE double
-interface_get_temperature
+interface_get_convection_coef
(const struct sdis_interface* interf,
const struct sdis_interface_fragment* frag)
{
ASSERT(interf && frag);
- if(!interf->shader.temperature) return -DBL_MAX;
- return interf->shader.temperature(frag, interf->data);
+ return interf->shader.convection_coef
+ ? interf->shader.convection_coef(frag, interf->data) : 0;
}
static INLINE double
-interface_get_convection_coef
+interface_side_get_temperature
+ (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");
+ }
+ return shader->temperature ? shader->temperature(frag, interf->data) : -1;
+}
+
+static INLINE double
+interface_side_get_flux
(const struct sdis_interface* interf,
const struct sdis_interface_fragment* frag)
{
+ const struct sdis_interface_side_shader* shader;
ASSERT(interf && frag);
- return interf->shader.convection_coef(frag, interf->data);
+ switch(frag->side) {
+ case SDIS_FRONT: shader = &interf->shader.front; break;
+ case SDIS_BACK: shader = &interf->shader.back; break;
+ default: FATAL("Unreachable code.\n");
+ }
+ return shader->flux ? shader->flux(frag, interf->data) : SDIS_FLUX_NONE;
}
static INLINE double
-interface_get_emissivity
+interface_side_get_emissivity
(const struct sdis_interface* interf,
const struct sdis_interface_fragment* frag)
{
+ const struct sdis_interface_side_shader* shader;
ASSERT(interf && frag);
- return interf->shader.emissivity(frag, interf->data);
+ 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->emissivity ? shader->emissivity(frag, interf->data) : 0;
}
static INLINE double
-interface_get_specular_fraction
+interface_side_get_specular_fraction
(const struct sdis_interface* interf,
const struct sdis_interface_fragment* frag)
{
+ const struct sdis_interface_side_shader* shader;
ASSERT(interf && frag);
- return interf->shader.specular_fraction(frag, interf->data);
+ 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->specular_fraction
+ ? shader->specular_fraction(frag, interf->data) : 0;
}
#endif /* SDIS_INTERFACE_C_H */
diff --git a/src/sdis_medium.c b/src/sdis_medium.c
@@ -39,7 +39,6 @@ check_solid_shader(const struct sdis_solid_shader* shader)
&& shader->thermal_conductivity
&& shader->volumic_mass
&& shader->delta_solid
- && shader->delta_boundary
&& shader->temperature;
}
@@ -116,7 +115,7 @@ sdis_fluid_create
goto error;
}
- res = medium_create(dev, &medium, SDIS_MEDIUM_FLUID);
+ res = medium_create(dev, &medium, SDIS_FLUID);
if(res != RES_OK) {
log_err(dev, "%s: could not create the fluid medium.\n", FUNC_NAME);
goto error;
@@ -161,7 +160,7 @@ sdis_solid_create
goto error;
}
- res = medium_create(dev, &medium, SDIS_MEDIUM_SOLID);
+ res = medium_create(dev, &medium, SDIS_SOLID);
if(res != RES_OK) {
log_err(dev, "%s: could not create the solid medium.\n", FUNC_NAME);
goto error;
@@ -207,3 +206,9 @@ sdis_medium_get_type(const struct sdis_medium* medium)
ASSERT(medium != NULL);
return medium->type;
}
+struct sdis_data*
+sdis_medium_get_data(struct sdis_medium* medium)
+{
+ ASSERT(medium);
+ return medium->data;
+}
diff --git a/src/sdis_medium_c.h b/src/sdis_medium_c.h
@@ -46,7 +46,7 @@ static INLINE double
fluid_get_calorific_capacity
(const struct sdis_medium* mdm, const struct sdis_rwalk_vertex* vtx)
{
- ASSERT(mdm && mdm->type == SDIS_MEDIUM_FLUID);
+ ASSERT(mdm && mdm->type == SDIS_FLUID);
return mdm->shader.fluid.calorific_capacity(vtx, mdm->data);
}
@@ -54,7 +54,7 @@ static INLINE double
fluid_get_volumic_mass
(const struct sdis_medium* mdm, const struct sdis_rwalk_vertex* vtx)
{
- ASSERT(mdm && mdm->type == SDIS_MEDIUM_FLUID);
+ ASSERT(mdm && mdm->type == SDIS_FLUID);
return mdm->shader.fluid.volumic_mass(vtx, mdm->data);
}
@@ -62,7 +62,7 @@ static INLINE double
fluid_get_temperature
(const struct sdis_medium* mdm, const struct sdis_rwalk_vertex* vtx)
{
- ASSERT(mdm && mdm->type == SDIS_MEDIUM_FLUID);
+ ASSERT(mdm && mdm->type == SDIS_FLUID);
return mdm->shader.fluid.temperature(vtx, mdm->data);
}
@@ -73,7 +73,7 @@ static INLINE double
solid_get_calorific_capacity
(const struct sdis_medium* mdm, const struct sdis_rwalk_vertex* vtx)
{
- ASSERT(mdm && mdm->type == SDIS_MEDIUM_SOLID);
+ ASSERT(mdm && mdm->type == SDIS_SOLID);
return mdm->shader.solid.calorific_capacity(vtx, mdm->data);
}
@@ -81,7 +81,7 @@ static INLINE double
solid_get_thermal_conductivity
(const struct sdis_medium* mdm, const struct sdis_rwalk_vertex* vtx)
{
- ASSERT(mdm && mdm->type == SDIS_MEDIUM_SOLID);
+ ASSERT(mdm && mdm->type == SDIS_SOLID);
return mdm->shader.solid.thermal_conductivity(vtx, mdm->data);
}
@@ -89,7 +89,7 @@ static INLINE double
solid_get_volumic_mass
(const struct sdis_medium* mdm, const struct sdis_rwalk_vertex* vtx)
{
- ASSERT(mdm && mdm->type == SDIS_MEDIUM_SOLID);
+ ASSERT(mdm && mdm->type == SDIS_SOLID);
return mdm->shader.solid.volumic_mass(vtx, mdm->data);
}
@@ -97,33 +97,25 @@ static INLINE double
solid_get_delta
(const struct sdis_medium* mdm, const struct sdis_rwalk_vertex* vtx)
{
- ASSERT(mdm && mdm->type == SDIS_MEDIUM_SOLID);
+ ASSERT(mdm && mdm->type == SDIS_SOLID);
return mdm->shader.solid.delta_solid(vtx, mdm->data);
}
static INLINE double
-solid_get_delta_boundary
- (const struct sdis_medium* mdm, const struct sdis_rwalk_vertex* vtx)
-{
- ASSERT(mdm && mdm->type == SDIS_MEDIUM_SOLID);
- return mdm->shader.solid.delta_boundary(vtx, mdm->data);
-}
-
-static INLINE double
solid_get_volumic_power
(const struct sdis_medium* mdm, const struct sdis_rwalk_vertex* vtx)
{
- ASSERT(mdm && mdm->type == SDIS_MEDIUM_SOLID);
- return mdm->shader.solid.volumic_power
+ ASSERT(mdm && mdm->type == SDIS_SOLID);
+ return mdm->shader.solid.volumic_power
? mdm->shader.solid.volumic_power(vtx, mdm->data)
- : 0;
+ : SDIS_VOLUMIC_POWER_NONE;
}
static INLINE double
solid_get_temperature
(const struct sdis_medium* mdm, const struct sdis_rwalk_vertex* vtx)
{
- ASSERT(mdm && mdm->type == SDIS_MEDIUM_SOLID);
+ ASSERT(mdm && mdm->type == SDIS_SOLID);
return mdm->shader.solid.temperature(vtx, mdm->data);
}
diff --git a/src/sdis_scene.c b/src/sdis_scene.c
@@ -26,10 +26,71 @@
#include "sdis.h"
#include "sdis_scene_c.h"
+#include <limits.h>
+
/*******************************************************************************
* Helper function
******************************************************************************/
static void
+project_position
+ (const double V0[3],
+ const double E0[3],
+ const double N[3],
+ const double NxE1[3],
+ const double rcp_det,
+ const double pos[3],
+ double uvw[3])
+{
+ double T[3], Q[3], k;
+ ASSERT(V0 && E0 && N && NxE1 && pos && uvw);
+
+ /* Use Moller/Trumbore intersection test the compute the parametric
+ * coordinates of the intersection between the triangle and the ray
+ * `r = pos + N*d' */
+ d3_sub(T, pos, V0);
+ uvw[0] = d3_dot(T, NxE1) * rcp_det;
+ d3_cross(Q, T, E0);
+ uvw[1] = d3_dot(Q, N) * rcp_det;
+ uvw[2] = 1.0 - uvw[0] - uvw[1];
+
+ if(uvw[0] >= 0 && uvw[1] >= 0 && uvw[2] >= 0) {/* The ray hits the triangle */
+ ASSERT(eq_eps(uvw[0] + uvw[1] + uvw[2], 1.0, 1.e-6));
+ return;
+ }
+
+ /* Clamp barycentric coordinates to triangle edges */
+ if(uvw[0] >= 0) {
+ if(uvw[1] >= 0) {
+ k = 1.0 / (uvw[0] + uvw[1]);
+ uvw[0] *= k;
+ uvw[1] *= k;
+ uvw[2] = 0;
+ } else if( uvw[2] >= 0) {
+ k = 1.0 / (uvw[0] + uvw[2]);
+ uvw[0] *= k;
+ uvw[1] = 0;
+ uvw[2] *= k;
+ } else {
+ ASSERT(uvw[0] >= 1.f);
+ d3(uvw, 1, 0, 0);
+ }
+ } else if(uvw[1] >= 0) {
+ if(uvw[2] >= 0) {
+ k = 1.0 / (uvw[1] + uvw[2]);
+ uvw[0] = 0;
+ uvw[1] *= k;
+ uvw[2] *= k;
+ } else {
+ ASSERT(uvw[1] >= 1);
+ d3(uvw, 0, 1, 0);
+ }
+ } else {
+ ASSERT(uvw[2] >= 1);
+ d3(uvw, 0, 0, 1);
+ }
+}
+
+static void
scene_release(ref_T * ref)
{
struct sdis_device* dev = NULL;
@@ -150,6 +211,69 @@ sdis_scene_get_boundary_position
return RES_OK;
}
+res_T
+sdis_scene_boundary_project_position
+ (const struct sdis_scene* scn,
+ const size_t iprim,
+ const double pos[],
+ double uv[])
+{
+ if(!scn || !pos || !uv) return RES_BAD_ARG;
+ if(iprim >= scene_get_primitives_count(scn)) return RES_BAD_ARG;
+
+ if(scene_is_2d(scn)) {
+ struct s2d_primitive prim;
+ struct s2d_attrib a;
+ double V[2][2]; /* Vertices */
+ double E[2][3]; /* V0->V1 and V0->pos */
+ double proj;
+
+ /* Retrieve the segment vertices */
+ S2D(scene_view_get_primitive(scn->s2d_view, (unsigned int)iprim, &prim));
+ S2D(primitive_get_attrib(&prim, S2D_POSITION, 0, &a)); d2_set_f2(V[0], a.value);
+ S2D(primitive_get_attrib(&prim, S2D_POSITION, 1, &a)); d2_set_f2(V[1], a.value);
+
+ /* Compute the parametric coordinate of the project of `pos' onto the
+ * segment.*/
+ d2_sub(E[0], V[1], V[0]);
+ d2_normalize(E[0], E[0]);
+ d2_sub(E[1], pos, V[0]);
+ proj = d2_dot(E[0], E[1]);
+
+ uv[0] = CLAMP(proj, 0, 1); /* Clamp the parametric coordinate in [0, 1] */
+
+ } else {
+ struct s3d_primitive prim;
+ struct s3d_attrib a;
+ double V[3][3]; /* Vertices */
+ double E[2][3]; /* V0->V1 and V0->V2 edges */
+ double N[3]; /* Normal */
+ double NxE1[3], rcp_det; /* Muller/Trumboer triangle parameters */
+ double uvw[3];
+
+ S3D(scene_view_get_primitive(scn->s3d_view, (unsigned int)iprim, &prim));
+ S3D(triangle_get_vertex_attrib(&prim, 0, S3D_POSITION, &a)); d3_set_f3(V[0], a.value);
+ S3D(triangle_get_vertex_attrib(&prim, 1, S3D_POSITION, &a)); d3_set_f3(V[1], a.value);
+ S3D(triangle_get_vertex_attrib(&prim, 2, S3D_POSITION, &a)); d3_set_f3(V[2], a.value);
+ d3_sub(E[0], V[1], V[0]);
+ d3_sub(E[1], V[2], V[0]);
+ d3_cross(N, E[0], E[1]);
+
+ /* Muller/Trumbore triangle parameters */
+ d3_cross(NxE1, N, E[1]);
+ rcp_det = 1.0 / d3_dot(NxE1, E[0]);
+
+ /* Use the Muller/Trumbore intersection test to project `pos' onto the
+ * triangle and to retrieve the parametric coordinates of the projection
+ * point */
+ project_position(V[0], E[0], N, NxE1, rcp_det, pos, uvw);
+
+ uv[0] = uvw[2];
+ uv[1] = uvw[0];
+ }
+ return RES_OK;
+}
+
/*******************************************************************************
* Local miscellaneous function
******************************************************************************/
@@ -164,10 +288,11 @@ res_T
scene_get_medium
(const struct sdis_scene* scn,
const double pos[],
+ struct get_medium_info* info,
const struct sdis_medium** out_medium)
{
return scene_is_2d(scn)
- ? scene_get_medium_2d(scn, pos, out_medium)
- : scene_get_medium_3d(scn, pos, out_medium);
+ ? scene_get_medium_2d(scn, pos, info, out_medium)
+ : scene_get_medium_3d(scn, pos, info, out_medium);
}
diff --git a/src/sdis_scene_Xd.h b/src/sdis_scene_Xd.h
@@ -154,6 +154,12 @@ clear_properties(struct sdis_scene* scn)
/* Macro making generic its subimitted name to SDIS_SCENE_DIMENSION */
#define XD(Name) CONCAT(CONCAT(CONCAT(Name, _), DIM), d)
+#if DIM == 2
+ #define HIT_ON_BOUNDARY hit_on_vertex
+#else
+ #define HIT_ON_BOUNDARY hit_on_edge
+#endif
+
/*******************************************************************************
* Helper functions
******************************************************************************/
@@ -625,7 +631,7 @@ XD(setup_enclosures)(struct sdis_scene* scn, struct sencXd(descriptor)* desc)
ASSERT(mdm);
/* Silently discard the solid and infinite enclosures */
- if(mdm->type == SDIS_MEDIUM_FLUID && !header.is_infinite) {
+ if(mdm->type == SDIS_FLUID && !header.is_infinite) {
res = XD(setup_enclosure_geometry)(scn, enc);
if(res != RES_OK) goto error;
}
@@ -717,16 +723,34 @@ error:
static INLINE res_T
XD(scene_get_medium)
(const struct sdis_scene* scn,
- const double pos[3],
+ const double pos[2],
+ struct get_medium_info* info, /* May be NULL */
const struct sdis_medium** out_medium)
{
const struct sdis_medium* medium = NULL;
size_t iprim, nprims;
size_t nfailures = 0;
const size_t max_failures = 10;
+ /* Range of the parametric coordinate into which positions are challenged */
+#if DIM == 2
+ float st[3];
+#else
+ float st[3][2];
+#endif
+ size_t nsteps = 3;
res_T res = RES_OK;
ASSERT(scn && pos);
+#if DIM == 2
+ st[0] = 0.25f;
+ st[1] = 0.50f;
+ st[2] = 0.75f;
+#else
+ f2(st[0], 1.f/6.f, 5.f/12.f);
+ f2(st[1], 5.f/12.f, 1.f/6.f);
+ f2(st[2], 5.f/12.f, 5.f/12.f);
+#endif
+
SXD(scene_view_primitives_count(scn->sXd(view), &nprims));
FOR_EACH(iprim, 0, nprims) {
struct sXd(hit) hit;
@@ -734,42 +758,53 @@ XD(scene_get_medium)
struct sXd(primitive) prim;
const float range[2] = {0.f, FLT_MAX};
float N[DIM], P[DIM], dir[DIM], cos_N_dir;
-#if DIM == 2
- float st;
- st = 1.f / 3.f;
-#else
- float st[2];
- st[0] = st[1] = 1.f / 3.f; /* Or MSVC will issue a warning */
-#endif
-
- /* Retrieve a position onto the primitive */
- SXD(scene_view_get_primitive(scn->sXd(view), (unsigned)iprim, &prim));
- SXD(primitive_get_attrib(&prim, SXD_POSITION, st, &attr));
+ size_t istep = 0;
+
+ do {
+ /* Retrieve a position onto the primitive */
+ SXD(scene_view_get_primitive(scn->sXd(view), (unsigned)iprim, &prim));
+ SXD(primitive_get_attrib(&prim, SXD_POSITION, st[istep], &attr));
+
+ /* Trace a ray from the random walk vertex toward the retrieved primitive
+ * position */
+ fX(normalize)(dir, fX(sub)(dir, attr.value, fX_set_dX(P, pos)));
+ SXD(scene_view_trace_ray(scn->sXd(view), P, dir, range, NULL, &hit));
+
+ /* Unforeseen error. One has to intersect a primitive ! */
+ if(SXD_HIT_NONE(&hit)) {
+ ++nfailures;
+ if(nfailures < max_failures) {
+ continue;
+ } else {
+ res = RES_BAD_ARG;
+ goto error;
+ }
+ }
+ /* Discard the hit if it is on a vertex, i.e. between 2 segments, and
+ * target a new position onto the current primitive */
+ } while((SXD_HIT_NONE(&hit) || HIT_ON_BOUNDARY(&hit))
+ && ++istep < nsteps);
- /* Trace a ray from the randomw walk vertex toward the retrieved primitive
- * position */
- fX(normalize)(dir, fX(sub)(dir, attr.value, fX_set_dX(P, pos)));
- SXD(scene_view_trace_ray(scn->sXd(view), P, dir, range, NULL, &hit));
+ /* The hits of all targeted positions on the current primitive are on
+ * vertices. Challenge positions on another primitive. */
+ if(istep > nsteps) continue;
fX(normalize)(N, hit.normal);
cos_N_dir = fX(dot)(N, dir);
- /* Unforeseen error. One has to intersect a primitive ! */
- if(SXD_HIT_NONE(&hit)) {
- ++nfailures;
- if(nfailures < max_failures) {
- continue;
- } else {
- res = RES_BAD_ARG;
- goto error;
- }
- }
-
if(absf(cos_N_dir) > 1.e-1f) { /* Not roughly orthognonal */
const struct sdis_interface* interf;
interf = scene_get_interface(scn, hit.prim.prim_id);
medium = interface_get_medium
(interf, cos_N_dir < 0 ? SDIS_FRONT : SDIS_BACK);
+
+ /* Register the get_medium_info */
+ if(info) {
+ fX(set)(info->pos_tgt, attr.value);
+ fX(set)(info->ray_org, P);
+ fX(set)(info->ray_dir, dir);
+ info->XD(hit) = hit;
+ }
break;
}
}
@@ -787,7 +822,6 @@ error:
#endif
goto exit;
}
-
#undef SDIS_SCENE_DIMENSION
#undef DIM
#undef sencXd
@@ -805,5 +839,6 @@ error:
#undef fX
#undef fX_set_dX
#undef XD
+#undef HIT_ON_BOUNDARY
#endif /* !SDIS_SCENE_DIMENSION */
diff --git a/src/sdis_scene_c.h b/src/sdis_scene_c.h
@@ -31,6 +31,17 @@ struct prim_prop {
unsigned back_enclosure; /* Id of the back facing enclosure */
};
+struct get_medium_info {
+ /* Targeted position */
+ float pos_tgt[3];
+ /* Ray trace to the targeted position in order to define the current medium */
+ float ray_org[3];
+ float ray_dir[3];
+ /* Hit encouters along the ray and used to define the current medium */
+ struct s2d_hit hit_2d;
+ struct s3d_hit hit_3d;
+};
+
static INLINE void
prim_prop_init(struct mem_allocator* allocator, struct prim_prop* prim)
{
@@ -176,6 +187,7 @@ extern LOCAL_SYM res_T
scene_get_medium
(const struct sdis_scene* scene,
const double position[],
+ struct get_medium_info* info, /* May be NULL */
const struct sdis_medium** medium);
static INLINE void
diff --git a/src/sdis_solve.c b/src/sdis_solve.c
@@ -88,7 +88,9 @@ solve_pixel
sum_weights += w;
sum_weights_sqr += w*w;
++N;
- } else if(res != RES_BAD_OP) {
+ } else if(res == RES_BAD_OP) {
+ res = RES_OK;
+ } else {
goto error;
}
}
@@ -96,6 +98,7 @@ solve_pixel
accum->sum_weights = sum_weights;
accum->sum_weights_sqr = sum_weights_sqr;
accum->nweights = N;
+ accum->nfailures = nrealisations - N;
exit:
return res;
@@ -174,13 +177,14 @@ sdis_solve_probe
struct ssp_rng** rngs = NULL;
double weight = 0;
double sqr_weight = 0;
- size_t irealisation = 0;
+ const int64_t rcount = (int64_t)nrealisations;
+ int64_t irealisation = 0;
size_t N = 0; /* #realisations that do not fail */
size_t i;
ATOMIC res = RES_OK;
- if(!scn || !nrealisations || !position || time < 0 || fp_to_meter <= 0
- || Tref < 0 || !out_estimator) {
+ if(!scn || !nrealisations || nrealisations > INT64_MAX || !position
+ || time < 0 || fp_to_meter <= 0 || Tref < 0 || !out_estimator) {
res = RES_BAD_ARG;
goto error;
}
@@ -207,15 +211,15 @@ sdis_solve_probe
if(res != RES_OK) goto error;
/* Retrieve the medium in which the submitted position lies */
- res = scene_get_medium(scn, position, &medium);
+ res = scene_get_medium(scn, position, NULL, &medium);
if(res != RES_OK) goto error;
/* Here we go! Launch the Monte Carlo estimation */
omp_set_num_threads((int)scn->dev->nthreads);
#pragma omp parallel for schedule(static) reduction(+:weight,sqr_weight,N)
- for(irealisation = 0; irealisation < nrealisations; ++irealisation) {
+ for(irealisation = 0; irealisation < rcount; ++irealisation) {
res_T res_local;
- double w;
+ double w = NaN;
const int ithread = omp_get_thread_num();
struct ssp_rng* rng = rngs[ithread];
@@ -273,6 +277,7 @@ sdis_solve_probe_boundary
const size_t iprim, /* Identifier of the primitive on which the probe lies */
const double uv[2], /* Parametric coordinates of the probe onto the primitve */
const double time, /* Observation time */
+ const enum sdis_side side, /* Side of iprim on which the probe lies */
const double fp_to_meter, /* Scale from floating point units to meters */
const double Tarad, /* In Kelvin */
const double Tref, /* In Kelvin */
@@ -283,13 +288,15 @@ sdis_solve_probe_boundary
struct ssp_rng** rngs = NULL;
double weight = 0;
double sqr_weight = 0;
- size_t irealisation = 0;
+ const int64_t rcount = (int64_t)nrealisations;
+ int64_t irealisation = 0;
size_t N = 0; /* #realisations that do not fail */
size_t i;
- res_T res = RES_OK;
+ ATOMIC res = RES_OK;
- if(!scn || !nrealisations || !uv || time < 0 || fp_to_meter <= 0
- || Tref < 0 || !out_estimator) {
+ if(!scn || !nrealisations || nrealisations > INT64_MAX || !uv || time < 0
+ || fp_to_meter <= 0 || Tref < 0 || (side != SDIS_FRONT && side != SDIS_BACK)
+ || !out_estimator) {
res = RES_BAD_ARG;
goto error;
}
@@ -353,9 +360,9 @@ sdis_solve_probe_boundary
/* Here we go! Launch the Monte Carlo estimation */
omp_set_num_threads((int)scn->dev->nthreads);
#pragma omp parallel for schedule(static) reduction(+:weight,sqr_weight,N)
- for(irealisation = 0; irealisation < nrealisations; ++irealisation) {
+ for(irealisation = 0; irealisation < rcount; ++irealisation) {
res_T res_local;
- double w;
+ double w = NaN;
const int ithread = omp_get_thread_num();
struct ssp_rng* rng = rngs[ithread];
@@ -363,10 +370,10 @@ sdis_solve_probe_boundary
if(scene_is_2d(scn)) {
res_local = boundary_realisation_2d
- (scn, rng, iprim, uv, time, fp_to_meter, Tarad, Tref, &w);
+ (scn, rng, iprim, uv, time, side, fp_to_meter, Tarad, Tref, &w);
} else {
res_local = boundary_realisation_3d
- (scn, rng, iprim, uv, time, fp_to_meter, Tarad, Tref, &w);
+ (scn, rng, iprim, uv, time, side, fp_to_meter, Tarad, Tref, &w);
}
if(res_local != RES_OK) {
if(res_local != RES_BAD_OP) {
@@ -397,7 +404,7 @@ exit:
}
if(rng_proxy) SSP(rng_proxy_ref_put(rng_proxy));
if(out_estimator) *out_estimator = estimator;
- return res;
+ return (res_T)res;
error:
if(estimator) {
SDIS(estimator_ref_put(estimator));
@@ -445,10 +452,10 @@ sdis_solve_camera
}
/* Retrieve the medium in which the submitted position lies */
- res = scene_get_medium(scn, cam->position, &medium);
+ res = scene_get_medium(scn, cam->position, NULL, &medium);
if(res != RES_OK) goto error;
- if(medium->type != SDIS_MEDIUM_FLUID) {
+ if(medium->type != SDIS_FLUID) {
log_err(scn->dev, "%s: the camera position `%g %g %g' is not in a fluid.\n",
FUNC_NAME, SPLIT3(cam->position));
res = RES_BAD_ARG;
@@ -490,7 +497,7 @@ sdis_solve_camera
pix_sz[1] = 1.0 / (double)height;
omp_set_num_threads((int)scn->dev->nthreads);
- #pragma omp parallel for schedule(static, 1/*chunki size*/)
+ #pragma omp parallel for schedule(static, 1/*chunk size*/)
for(mcode = 0; mcode < (int64_t)ntiles; ++mcode) {
size_t tile_org[2] = {0, 0};
size_t tile_sz[2] = {0, 0};
@@ -510,7 +517,7 @@ sdis_solve_camera
tile_org[0] *= TILE_SIZE;
tile_org[1] *= TILE_SIZE;
tile_sz[0] = MMIN(TILE_SIZE, width - tile_org[0]);
- tile_sz[1] = MMIN(TILE_SIZE, width - tile_org[1]);
+ tile_sz[1] = MMIN(TILE_SIZE, height - tile_org[1]);
/* Fetch the accumulations buffer */
accums = darray_accum_data_get(tiles+ithread);
diff --git a/src/sdis_solve_Xd.h b/src/sdis_solve_Xd.h
@@ -22,13 +22,26 @@
#include "sdis_medium_c.h"
#include "sdis_scene_c.h"
+#include <rsys/float2.h>
+#include <rsys/float3.h>
#include <rsys/stretchy_array.h>
#include <star/ssp.h>
-/* Emperical scale factor to apply to the upper bound of the ray range in order
+/* Define a new result code from RES_BAD_OP saying that the bad operation is
+ * definitive, i.e. in the current state, the realisation will inevitably fail.
+ * It is thus unecessary to retry a specific section of the random walk */
+#define RES_BAD_OP_IRRECOVERABLE (-RES_BAD_OP)
+
+/* Empirical scale factor to apply to the upper bound of the ray range in order
* to handle numerical imprecisions */
-#define RAY_RANGE_MAX_SCALE 1.0001f
+#define RAY_RANGE_MAX_SCALE 1.001f
+
+/* Emperical scale factor applied to the challenged reinjection distance. If
+ * the distance to reinject is less than this adjusted value, the solver
+ * switches from 2D reinjection scheme to the 1D reinjection scheme in order to
+ * avoid numerical issues. */
+#define REINJECT_DST_MIN_SCALE 0.125f
#define BOLTZMANN_CONSTANT 5.6696e-8 /* W/m^2/K^4 */
@@ -40,7 +53,22 @@ struct rwalk_context {
/* Reflect the vector V wrt the normal N. By convention V points outward the
* surface. */
static INLINE float*
-reflect(float res[3], const float V[3], const float N[3])
+reflect_2d(float res[2], const float V[2], const float N[2])
+{
+ float tmp[2];
+ float cos_V_N;
+ ASSERT(res && V && N);
+ ASSERT(f2_is_normalized(V) && f2_is_normalized(N));
+ cos_V_N = f2_dot(V, N);
+ f2_mulf(tmp, N, 2*cos_V_N);
+ f2_sub(res, tmp, V);
+ return res;
+}
+
+/* Reflect the vector V wrt the normal N. By convention V points outward the
+ * surface. */
+static INLINE float*
+reflect_3d(float res[3], const float V[3], const float N[3])
{
float tmp[3];
float cos_V_N;
@@ -94,9 +122,10 @@ struct XD(rwalk) {
struct sdis_rwalk_vertex vtx; /* Position and time of the Random walk */
const struct sdis_medium* mdm; /* Medium in which the random walk lies */
struct sXd(hit) hit; /* Hit of the random walk */
+ enum sdis_side hit_side;
};
static const struct XD(rwalk) XD(RWALK_NULL) = {
- SDIS_RWALK_VERTEX_NULL__, NULL, SXD_HIT_NULL__
+ SDIS_RWALK_VERTEX_NULL__, NULL, SXD_HIT_NULL__, SDIS_SIDE_NULL__
};
struct XD(temperature) {
@@ -162,6 +191,48 @@ XD(move_pos)(double pos[DIM], const float dir[DIM], const float delta)
#endif
}
+static FINLINE void
+XD(sample_reinjection_dir)
+ (const struct XD(rwalk)* rwalk, struct ssp_rng* rng, float dir[DIM])
+{
+#if DIM == 2
+ /* The sampled directions is defined by rotating the normal around the Z axis
+ * of an angle of PI/4 or -PI/4. Let the rotation matrix defined as
+ * | cos(a) -sin(a) |
+ * | sin(a) cos(a) |
+ * with a = PI/4, dir = sqrt(2)/2 * | 1 -1 | . N
+ * | 1 1 |
+ * with a =-PI/4, dir = sqrt(2)/2 * | 1 1 | . N
+ * |-1 1 |
+ * Note that since the sampled direction is finally normalized, we can
+ * discard the sqrt(2)/2 constant. */
+ const uint64_t r = ssp_rng_uniform_uint64(rng, 0, 1);
+ ASSERT(rwalk && dir);
+ if(r) {
+ dir[0] = rwalk->hit.normal[0] - rwalk->hit.normal[1];
+ dir[1] = rwalk->hit.normal[0] + rwalk->hit.normal[1];
+ } else {
+ dir[0] = rwalk->hit.normal[0] + rwalk->hit.normal[1];
+ dir[1] =-rwalk->hit.normal[0] + rwalk->hit.normal[1];
+ }
+ f2_normalize(dir, dir);
+#else
+ /* Sample a random direction around the normal whose cosine is 1/sqrt(3). To
+ * do so we sample a position onto a cone whose height is 1/sqrt(2) and the
+ * radius of its base is 1. */
+ float frame[9];
+ ASSERT(fX(is_normalized)(rwalk->hit.normal));
+
+ ssp_ran_circle_uniform_float(rng, dir, NULL);
+ dir[2] = (float)(1.0/sqrt(2));
+
+ f33_basis(frame, rwalk->hit.normal);
+ f33_mulf3(dir, frame, dir);
+ f3_normalize(dir, dir);
+ ASSERT(eq_epsf(f3_dot(dir, rwalk->hit.normal), (float)(1.0/sqrt(3)), 1.e-4f));
+#endif
+}
+
/* Check that the interface fragment is consistent with the current state of
* the random walk */
static INLINE int
@@ -231,35 +302,41 @@ XD(trace_radiative_path)
(scn->sXd(view), pos, dir, range, &rwalk->hit, &rwalk->hit));
#endif
if(SXD_HIT_NONE(&rwalk->hit)) { /* Fetch the ambient radiative temperature */
+ rwalk->hit_side = SDIS_SIDE_NULL__;
if(ctx->Tarad >= 0) {
T->value += ctx->Tarad;
T->done = 1;
break;
} else {
log_err(scn->dev,
-"%s: the random walk reaches an invalid ambient radiative temperature of `%gK'\n"
-"at position `%g %g %g'. This may be due to numerical inaccuracies or to\n"
-"inconsistency in the simulated system (eg: unclosed geometry). For systems\n"
-"where the random walks can reach such temperature, one has to setup a valid\n"
-"ambient radiative temperature, i.e. it must be greater or equal to 0.\n",
+ "%s: the random walk reaches an invalid ambient radiative temperature "
+ "of `%gK' at position `%g %g %g'. This may be due to numerical "
+ "inaccuracies or to inconsistency in the simulated system (eg: "
+ "unclosed geometry). For systems where the random walks can reach "
+ "such temperature, one has to setup a valid ambient radiative "
+ "temperature, i.e. it must be greater or equal to 0.\n",
FUNC_NAME,
ctx->Tarad,
SPLIT3(rwalk->vtx.P));
- res = RES_BAD_ARG;
+ res = RES_BAD_OP;
goto error;
}
}
+ /* Define the hit side */
+ rwalk->hit_side = fX(dot)(dir, rwalk->hit.normal) < 0
+ ? SDIS_FRONT : SDIS_BACK;
+
/* Move the random walk to the hit position */
XD(move_pos)(rwalk->vtx.P, dir, rwalk->hit.distance);
/* Fetch the new interface and setup the hit fragment */
interf = scene_get_interface(scn, rwalk->hit.prim.prim_id);
- XD(setup_interface_fragment)(&frag, &rwalk->vtx, &rwalk->hit);
+ XD(setup_interface_fragment)(&frag, &rwalk->vtx, &rwalk->hit, rwalk->hit_side);
/* Fetch the interface emissivity */
- epsilon = interface_get_emissivity(interf, &frag);
- if(epsilon > 1 && epsilon >= 0) {
+ epsilon = interface_side_get_emissivity(interf, &frag);
+ if(epsilon > 1 || epsilon < 0) {
log_err(scn->dev,
"%s: invalid overall emissivity `%g' at position `%g %g %g'.\n",
FUNC_NAME, epsilon, SPLIT3(rwalk->vtx.P));
@@ -278,7 +355,7 @@ XD(trace_radiative_path)
/* Normalize the normal of the interface and ensure that it points toward the
* current medium */
fX(normalize)(N, rwalk->hit.normal);
- if(f3_dot(N, dir) > 0) {
+ if(rwalk->hit_side == SDIS_BACK){
chk_mdm = interf->medium_back;
fX(minus)(N, N);
} else {
@@ -288,13 +365,13 @@ XD(trace_radiative_path)
if(chk_mdm != rwalk->mdm) {
log_err(scn->dev, "%s: inconsistent medium definition at `%g %g %g'.\n",
FUNC_NAME, SPLIT3(rwalk->vtx.P));
- res = RES_BAD_ARG;
+ res = RES_BAD_OP;
goto error;
}
- alpha = interface_get_specular_fraction(interf, &frag);
+ alpha = interface_side_get_specular_fraction(interf, &frag);
r = ssp_rng_canonical(rng);
if(r < alpha) { /* Sample specular part */
- reflect(dir, f3_minus(dir, dir), N);
+ reflect_3d(dir, f3_minus(dir, dir), N);
} else { /* Sample diffuse part */
ssp_ran_hemisphere_cos_float(rng, N, dir, NULL);
}
@@ -315,8 +392,6 @@ XD(radiative_temperature)
struct ssp_rng* rng,
struct XD(temperature)* T)
{
- const struct sdis_interface* interf;
-
/* The radiative random walk is always perform in 3D. In 2D, the geometry are
* assumed to be extruded to the infinty along the Z dimension. */
float N[3] = {0, 0, 0};
@@ -327,13 +402,10 @@ XD(radiative_temperature)
ASSERT(!SXD_HIT_NONE(&rwalk->hit));
(void)fp_to_meter;
- /* Fetch the current interface */
- interf = scene_get_interface(scn, rwalk->hit.prim.prim_id);
-
/* Normalize the normal of the interface and ensure that it points toward the
* current medium */
fX(normalize(N, rwalk->hit.normal));
- if(interf->medium_back == rwalk->mdm) {
+ if(rwalk->hit_side == SDIS_BACK) {
fX(minus(N, N));
}
@@ -378,7 +450,7 @@ XD(fluid_temperature)
#endif
(void)rng, (void)fp_to_meter, (void)ctx;
ASSERT(scn && fp_to_meter > 0 && ctx && rwalk && rng && T);
- ASSERT(rwalk->mdm->type == SDIS_MEDIUM_FLUID);
+ ASSERT(rwalk->mdm->type == SDIS_FLUID);
tmp = fluid_get_temperature(rwalk->mdm, &rwalk->vtx);
if(tmp >= 0) {
@@ -398,6 +470,7 @@ XD(fluid_temperature)
/* Init the path hit field required to define the current enclosure and
* fetch the interface data */
SXD(scene_view_trace_ray(scn->sXd(view), org, dir, range, NULL, &rwalk->hit));
+ rwalk->hit_side = fX(dot)(rwalk->hit.normal, dir) < 0 ? SDIS_FRONT : SDIS_BACK;
if(SXD_HIT_NONE(&rwalk->hit)) {
log_err(scn->dev,
@@ -409,7 +482,7 @@ XD(fluid_temperature)
}
/* Setup the fragment of the interface */
- XD(setup_interface_fragment)(&frag, &rwalk->vtx, &rwalk->hit);
+ XD(setup_interface_fragment)(&frag, &rwalk->vtx, &rwalk->hit, rwalk->hit_side);
/* Fetch the current interface and its associated enclosures */
interf = scene_get_interface(scn, rwalk->hit.prim.prim_id);
@@ -501,17 +574,27 @@ XD(solid_solid_boundary_temperature)
struct ssp_rng* rng,
struct XD(temperature)* T)
{
+ struct sXd(hit) hit0, hit1, hit2, hit3;
+ struct sXd(hit)* hit;
const struct sdis_interface* interf = NULL;
const struct sdis_medium* solid_front = NULL;
const struct sdis_medium* solid_back = NULL;
+ const struct sdis_medium* mdm;
double lambda_front, lambda_back;
- double delta_front_boundary, delta_back_boundary;
- double delta_front_boundary_meter, delta_back_boundary_meter;
+ double delta_front, delta_back;
+ double delta_boundary_front, delta_boundary_back;
double delta_boundary;
+ double reinject_dst_front, reinject_dst_back;
+ double reinject_dst;
double proba;
double tmp;
double r;
- float pos[DIM], dir[DIM], range[2];
+ double power;
+ float range0[2], range1[2];
+ float dir0[DIM], dir1[DIM], dir2[DIM], dir3[DIM];
+ float* dir;
+ float pos[DIM];
+ int dim = DIM;
ASSERT(scn && fp_to_meter > 0 && ctx && frag && rwalk && rng && T);
ASSERT(XD(check_rwalk_fragment_consistency)(rwalk, frag));
(void)frag, (void)ctx;
@@ -520,44 +603,111 @@ XD(solid_solid_boundary_temperature)
interf = scene_get_interface(scn, rwalk->hit.prim.prim_id);
solid_front = interface_get_medium(interf, SDIS_FRONT);
solid_back = interface_get_medium(interf, SDIS_BACK);
- ASSERT(solid_front->type == SDIS_MEDIUM_SOLID);
- ASSERT(solid_back->type == SDIS_MEDIUM_SOLID);
+ ASSERT(solid_front->type == SDIS_SOLID);
+ ASSERT(solid_back->type == SDIS_SOLID);
/* Fetch the properties of the media */
lambda_front = solid_get_thermal_conductivity(solid_front, &rwalk->vtx);
lambda_back = solid_get_thermal_conductivity(solid_back, &rwalk->vtx);
- delta_front_boundary = solid_get_delta_boundary(solid_front, &rwalk->vtx);
- delta_back_boundary = solid_get_delta_boundary(solid_back, &rwalk->vtx);
-
- /* Convert the delta boundary from floating point units to meters */
- delta_front_boundary_meter = fp_to_meter * delta_front_boundary;
- delta_back_boundary_meter = fp_to_meter * delta_back_boundary;
-
- /* Compute the proba to switch in solid0 or in solid1 */
- tmp = lambda_front / delta_front_boundary_meter;
- proba = tmp / (tmp + lambda_back / delta_back_boundary_meter);
+ /* Note that reinjection distance is *FIXED*. It MUST ensure that the orthogonal
+ * distance from the boundary to the point to challenge is equal to delta. */
+ delta_front = solid_get_delta(solid_front, &rwalk->vtx);
+ delta_back = solid_get_delta(solid_back, &rwalk->vtx);
+ delta_boundary_front = delta_front*sqrt(DIM);
+ delta_boundary_back = delta_back *sqrt(DIM);
+
+ /* Sample a reinjection direction and reflect it around the normal. Then
+ * reflect them on the back side of the interface. */
+ XD(sample_reinjection_dir)(rwalk, rng, dir0);
+ XD(reflect)(dir2, dir0, rwalk->hit.normal);
+ fX(minus)(dir1, dir0);
+ fX(minus)(dir3, dir2);
+
+ /* Trace the sampled directions on both sides of the interface to adjust the
+ * reinjection distance of the random walk . */
+ fX_set_dX(pos, rwalk->vtx.P);
+ f2(range0, 0, (float)delta_boundary_front*RAY_RANGE_MAX_SCALE);
+ f2(range1, 0, (float)delta_boundary_back *RAY_RANGE_MAX_SCALE);
+ SXD(scene_view_trace_ray(scn->sXd(view), pos, dir0, range0, &rwalk->hit, &hit0));
+ SXD(scene_view_trace_ray(scn->sXd(view), pos, dir1, range1, &rwalk->hit, &hit1));
+ SXD(scene_view_trace_ray(scn->sXd(view), pos, dir2, range0, &rwalk->hit, &hit2));
+ SXD(scene_view_trace_ray(scn->sXd(view), pos, dir3, range1, &rwalk->hit, &hit3));
+
+ /* Adjust the reinjection distance */
+ reinject_dst_front = MMIN(MMIN(delta_boundary_front, hit0.distance), hit2.distance);
+ reinject_dst_back = MMIN(MMIN(delta_boundary_back, hit1.distance), hit3.distance);
+
+ /* Define the reinjection side. Note that the proba should be :
+ * Lf/Df' / (Lf/Df' + Lb/Db')
+ *
+ * with L<f|b> the lambda of the <front|back> side and D<f|b>' the adjusted
+ * delta of the <front|back> side, i.e. :
+ * D<f|b>' = reinject_dst_<front|back> / sqrt(DIM)
+ *
+ * Anyway, one can avoid to compute the adjusted delta by directly using the
+ * adjusted reinjection distance since the resulting proba is strictly the
+ * same; sqrt(DIM) can be simplified. */
r = ssp_rng_canonical(rng);
- fX(normalize)(dir, rwalk->hit.normal);
+ proba = (lambda_front/reinject_dst_front)
+ / (lambda_front/reinject_dst_front + lambda_back/reinject_dst_back);
if(r < proba) { /* Reinject in front */
- delta_boundary = delta_front_boundary;
- rwalk->mdm = solid_front;
+ dir = dir0;
+ hit = &hit0;
+ mdm = solid_front;
+ reinject_dst = reinject_dst_front;
+ delta_boundary = delta_boundary_front;
} else { /* Reinject in back */
- delta_boundary = delta_back_boundary;
- fX(minus)(dir, dir);
- rwalk->mdm = solid_back;
+ dir = dir1;
+ hit = &hit1;
+ mdm = solid_back;
+ reinject_dst = reinject_dst_back;
+ delta_boundary = delta_boundary_back;
}
- /* "Reinject" the path into the solid along the surface normal. */
- fX_set_dX(pos, rwalk->vtx.P);
- range[0] = 0, range[1] = (float)delta_boundary*RAY_RANGE_MAX_SCALE;
- SXD(scene_view_trace_ray
- (scn->sXd(view), pos, dir, range, &rwalk->hit, &rwalk->hit));
- if(!SXD_HIT_NONE(&rwalk->hit)) delta_boundary = rwalk->hit.distance * 0.5;
- XD(move_pos)(rwalk->vtx.P, dir, (float)delta_boundary);
-
- /* Switch in solid random walk */
- T->func = XD(solid_temperature);
+ /* Switch in 1D reinjection scheme */
+ if(reinject_dst < delta_boundary * REINJECT_DST_MIN_SCALE) {
+ if(dir == dir0) {
+ fX(set)(dir, rwalk->hit.normal);
+ } else {
+ fX(minus)(dir, rwalk->hit.normal);
+ }
+
+ f2(range0, 0, (float)delta_boundary*RAY_RANGE_MAX_SCALE);
+ SXD(scene_view_trace_ray(scn->sXd(view), pos, dir, range0, &rwalk->hit, hit));
+ reinject_dst = MMIN(delta_boundary, hit->distance),
+ dim = 1;
+
+ /* Hit something in 1D. Arbitrarily move the random walk to 0.5 of the hit
+ * distance */
+ if(!SXD_HIT_NONE(hit)) {
+ reinject_dst *= 0.5;
+ *hit = SXD_HIT_NULL;
+ }
+ }
+
+ /* Handle the volumic power */
+ power = solid_get_volumic_power(mdm, &rwalk->vtx);
+ if(power != SDIS_VOLUMIC_POWER_NONE) {
+ const double delta_in_meter = reinject_dst * fp_to_meter;
+ const double lambda = solid_get_thermal_conductivity(mdm, &rwalk->vtx);
+ tmp = power * delta_in_meter * delta_in_meter / (2.0 * dim * lambda);
+ T->value += tmp;
+ }
+
+ /* Reinject */
+ XD(move_pos)(rwalk->vtx.P, dir, (float)reinject_dst);
+ if(eq_epsf(hit->distance, (float)reinject_dst, 1.e-4f)) {
+ T->func = XD(boundary_temperature);
+ rwalk->mdm = NULL;
+ rwalk->hit = *hit;
+ rwalk->hit_side = fX(dot)(hit->normal, dir) < 0 ? SDIS_FRONT : SDIS_BACK;
+ } else {
+ T->func = XD(solid_temperature);
+ rwalk->mdm = mdm;
+ rwalk->hit = SXD_HIT_NULL;
+ rwalk->hit_side = SDIS_SIDE_NULL__;
+ }
}
static void
@@ -575,72 +725,254 @@ XD(solid_fluid_boundary_temperature)
const struct sdis_medium* mdm_back = NULL;
const struct sdis_medium* solid = NULL;
const struct sdis_medium* fluid = NULL;
+ struct sXd(hit) hit0 = SXD_HIT_NULL;
+ struct sXd(hit) hit1 = SXD_HIT_NULL;
+ struct sdis_interface_fragment frag_fluid;
double hc;
double hr;
double epsilon; /* Interface emissivity */
double lambda;
double fluid_proba;
double radia_proba;
+ double delta;
double delta_boundary;
double r;
double tmp;
- float dir[DIM], pos[DIM], range[2];
+ float pos[DIM];
+ float dir0[DIM], dir1[DIM];
+ float range[2];
+ int dim = DIM;
ASSERT(scn && fp_to_meter > 0 && rwalk && rng && T && ctx);
ASSERT(XD(check_rwalk_fragment_consistency)(rwalk, frag));
- /* Retrieve the solid and the fluid split by the boundary */
+ /* Retrieve the solid and the fluid split by the boundary */
interf = scene_get_interface(scn, rwalk->hit.prim.prim_id);
mdm_front = interface_get_medium(interf, SDIS_FRONT);
mdm_back = interface_get_medium(interf, SDIS_BACK);
ASSERT(mdm_front->type != mdm_back->type);
- if(mdm_front->type == SDIS_MEDIUM_SOLID) {
+
+ frag_fluid = *frag;
+ if(mdm_front->type == SDIS_SOLID) {
solid = mdm_front;
fluid = mdm_back;
+ frag_fluid.side = SDIS_BACK;
} else {
solid = mdm_back;
fluid = mdm_front;
+ frag_fluid.side = SDIS_FRONT;
}
/* Fetch the solid properties */
lambda = solid_get_thermal_conductivity(solid, &rwalk->vtx);
- delta_boundary = solid_get_delta_boundary(solid, &rwalk->vtx);
+ delta = solid_get_delta(solid, &rwalk->vtx);
+
+ /* Note that the reinjection distance is *FIXED*. It MUST ensure that the
+ * orthogonal distance from the boundary to the point to chalenge is equal to
+ * delta. */
+ delta_boundary = sqrt(DIM) * delta;
+
+ /* Sample a reinjection direction */
+ XD(sample_reinjection_dir)(rwalk, rng, dir0);
+
+ /* Reflect the sampled direction around the normal */
+ XD(reflect)(dir1, dir0, rwalk->hit.normal);
+
+ if(solid == mdm_back) {
+ fX(minus)(dir0, dir0);
+ fX(minus)(dir1, dir1);
+ }
+
+ /* Trace dir0/dir1 to adjust the reinjection distance */
+ fX_set_dX(pos, rwalk->vtx.P);
+ f2(range, 0, (float)delta_boundary*RAY_RANGE_MAX_SCALE);
+ SXD(scene_view_trace_ray(scn->sXd(view), pos, dir0, range, &rwalk->hit, &hit0));
+ SXD(scene_view_trace_ray(scn->sXd(view), pos, dir1, range, &rwalk->hit, &hit1));
+
+ /* Adjust the delta boundary to the hit distance */
+ tmp = MMIN(MMIN(delta_boundary, hit0.distance), hit1.distance);
+
+ if(tmp >= delta_boundary * REINJECT_DST_MIN_SCALE) {
+ delta_boundary = tmp;
+ /* Define the orthogonal dst from the reinjection pos to the interface */
+ delta = delta_boundary / sqrt(DIM);
+ } else { /* Switch in 1D reinjection scheme. */
+ fX(set)(dir0, rwalk->hit.normal);
+ if(solid == mdm_back) fX(minus)(dir0, dir0);
+ f2(range, 0, (float)delta*RAY_RANGE_MAX_SCALE);
+ SXD(scene_view_trace_ray(scn->sXd(view), pos, dir0, range, &rwalk->hit, &hit0));
+ delta_boundary = MMIN(hit0.distance, delta);
+
+ /* Hit something in 1D. Arbitrarily move the random walk to 0.5 of the hit
+ * distance in order to avoid infinite bounces for parallel plane */
+ if(!SXD_HIT_NONE(&hit0)) {
+ delta_boundary *= 0.5;
+ hit0 = SXD_HIT_NULL;
+ }
+
+ delta = delta_boundary;
+ dim = 1;
+ }
/* Fetch the boundary properties */
- epsilon = interface_get_emissivity(interf, frag);
+ epsilon = interface_side_get_emissivity(interf, &frag_fluid);
hc = interface_get_convection_coef(interf, frag);
/* Compute the radiative coefficient */
hr = 4.0 * BOLTZMANN_CONSTANT * ctx->Tref3 * epsilon;
- /* Compute the probas to switch in solid or fluid random walk */
- tmp = lambda / (delta_boundary*fp_to_meter);
+ /* Compute the probas to switch in solid, fluid or radiative random walk */
+ tmp = lambda / (delta*fp_to_meter);
fluid_proba = hc / (tmp + hr + hc);
radia_proba = hr / (tmp + hr + hc);
/*solid_proba = tmp / (tmp + hr + hc);*/
r = ssp_rng_canonical(rng);
if(r < radia_proba) { /* Switch in radiative random walk */
- rwalk->mdm = fluid;
T->func = XD(radiative_temperature);
- } else if(r < fluid_proba + radia_proba) { /* Switch to fluid random walk */
rwalk->mdm = fluid;
+ rwalk->hit_side = rwalk->mdm == mdm_front ? SDIS_FRONT : SDIS_BACK;
+ } else if(r < fluid_proba + radia_proba) { /* Switch to fluid random walk */
T->func = XD(fluid_temperature);
+ rwalk->mdm = fluid;
+ rwalk->hit_side = rwalk->mdm == mdm_front ? SDIS_FRONT : SDIS_BACK;
} else { /* Solid random walk */
- rwalk->mdm = solid;
- fX(normalize)(dir, rwalk->hit.normal);
- if(solid == mdm_back) fX(minus)(dir, dir);
+ /* Handle the volumic power */
+ const double power = solid_get_volumic_power(solid, &rwalk->vtx);
+ if(power != SDIS_VOLUMIC_POWER_NONE) {
+ const double delta_in_meter = delta_boundary * fp_to_meter;
+ tmp = power * delta_in_meter * delta_in_meter / (2.0 * dim * lambda);
+ T->value += tmp;
+ }
- /* "Reinject" the random walk into the solid */
- fX_set_dX(pos, rwalk->vtx.P);
- range[0] = 0, range[1] = (float)delta_boundary*RAY_RANGE_MAX_SCALE;
- SXD(scene_view_trace_ray
- (scn->sXd(view), pos, dir, range, &rwalk->hit, &rwalk->hit));
- if(!SXD_HIT_NONE(&rwalk->hit)) delta_boundary = rwalk->hit.distance * 0.5;
- XD(move_pos)(rwalk->vtx.P, dir, (float)delta_boundary);
+ /* Reinject */
+ XD(move_pos)(rwalk->vtx.P, dir0, (float)delta_boundary);
+ if(eq_epsf(hit0.distance, (float)delta_boundary, 1.e-4f)) {
+ T->func = XD(boundary_temperature);
+ rwalk->mdm = NULL;
+ rwalk->hit = hit0;
+ rwalk->hit_side = fX(dot)(hit0.normal, dir0) < 0 ? SDIS_FRONT : SDIS_BACK;
+ } else {
+ T->func = XD(solid_temperature);
+ rwalk->mdm = solid;
+ rwalk->hit = SXD_HIT_NULL;
+ rwalk->hit_side = SDIS_SIDE_NULL__;
+ }
+ }
+}
- /* Switch in solid random walk */
+static void
+XD(solid_boundary_with_flux_temperature)
+ (const struct sdis_scene* scn,
+ const double fp_to_meter,
+ const struct rwalk_context* ctx,
+ const struct sdis_interface_fragment* frag,
+ const double phi,
+ struct XD(rwalk)* rwalk,
+ struct ssp_rng* rng,
+ struct XD(temperature)* T)
+{
+ const struct sdis_interface* interf = NULL;
+ const struct sdis_medium* mdm = NULL;
+ double lambda;
+ double delta;
+ double delta_boundary;
+ double delta_in_meter;
+ double power;
+ double tmp;
+ struct sXd(hit) hit0;
+ struct sXd(hit) hit1;
+ float pos[DIM];
+ float dir0[DIM];
+ float dir1[DIM];
+ float range[2];
+ int dim = DIM;
+ ASSERT(frag && phi != SDIS_FLUX_NONE);
+ ASSERT(XD(check_rwalk_fragment_consistency)(rwalk, frag));
+ (void)ctx;
+
+ /* Fetch current interface */
+ interf = scene_get_interface(scn, rwalk->hit.prim.prim_id);
+ ASSERT(phi == interface_side_get_flux(interf, frag));
+
+ /* Fetch incoming solid */
+ mdm = interface_get_medium(interf, frag->side);
+ ASSERT(mdm->type == SDIS_SOLID);
+
+ /* Fetch medium properties */
+ lambda = solid_get_thermal_conductivity(mdm, &rwalk->vtx);
+ delta = solid_get_delta(mdm, &rwalk->vtx);
+
+ /* Compute the reinjection distance. It MUST ensure that the orthogonal
+ * distance from the boundary to the point to chalenge is equal to delta. */
+ delta_boundary = delta * sqrt(DIM);
+
+ /* Sample a reinjection direction */
+ XD(sample_reinjection_dir)(rwalk, rng, dir0);
+
+ /* Reflect the sampled direction around the normal */
+ XD(reflect)(dir1, dir0, rwalk->hit.normal);
+
+ if(frag->side == SDIS_BACK) {
+ fX(minus)(dir0, dir0);
+ fX(minus)(dir1, dir1);
+ }
+
+ /* Trace dir0/dir1 to adjust the reinjection distance wrt the geometry */
+ fX_set_dX(pos, rwalk->vtx.P);
+ f2(range, 0, (float)delta_boundary*RAY_RANGE_MAX_SCALE);
+ SXD(scene_view_trace_ray(scn->sXd(view), pos, dir0, range, &rwalk->hit, &hit0));
+ SXD(scene_view_trace_ray(scn->sXd(view), pos, dir1, range, &rwalk->hit, &hit1));
+
+ /* Adjust the delta boundary to the hit distance */
+ tmp = MMIN(MMIN(delta_boundary, hit0.distance), hit1.distance);
+
+ if(tmp >= delta_boundary * REINJECT_DST_MIN_SCALE) {
+ delta_boundary = tmp;
+ /* Define the orthogonal dst from the reinjection pos to the interface */
+ delta = delta_boundary / sqrt(DIM);
+ } else { /* Switch in 1D reinjection scheme. */
+ fX(set)(dir0, rwalk->hit.normal);
+ if(frag->side == SDIS_BACK) fX(minus)(dir0, dir0);
+ f2(range, 0, (float)delta*RAY_RANGE_MAX_SCALE);
+ SXD(scene_view_trace_ray(scn->sXd(view), pos, dir0, range, &rwalk->hit, &hit0));
+ delta_boundary = MMIN(hit0.distance, delta_boundary);
+
+ /* Hit something in 1D. Arbitrarily move the random walk to 0.5 of the hit
+ * distance in order to avoid infinite bounces for parallel plane */
+ if(!SXD_HIT_NONE(&hit0)) {
+ delta_boundary *= 0.5;
+ hit0 = SXD_HIT_NULL;
+ }
+
+ delta = delta_boundary;
+ dim = 1;
+ }
+
+ /* Handle the flux */
+ delta_in_meter = delta*fp_to_meter;
+ T->value += phi * delta_in_meter / lambda;
+
+ /* Handle the volumic power */
+ power = solid_get_volumic_power(mdm, &rwalk->vtx);
+ if(power != SDIS_VOLUMIC_POWER_NONE) {
+ delta_in_meter = delta_boundary * fp_to_meter;
+ tmp = power * delta_in_meter * delta_in_meter / (2.0 * dim * lambda);
+ T->value += tmp;
+ }
+
+ /* Reinject into the solid */
+ XD(move_pos)(rwalk->vtx.P, dir0, (float)delta_boundary);
+ if(eq_epsf(hit0.distance, (float)delta_boundary, 1.e-4f)) {
+ T->func = XD(boundary_temperature);
+ rwalk->mdm = NULL;
+ rwalk->hit = hit0;
+ rwalk->hit_side = fX(dot)(hit0.normal, dir0) < 0 ? SDIS_FRONT : SDIS_BACK;
+ } else {
T->func = XD(solid_temperature);
+ rwalk->mdm = mdm;
+ rwalk->hit = SXD_HIT_NULL;
+ rwalk->hit_side = SDIS_SIDE_NULL__;
}
}
@@ -657,24 +989,39 @@ XD(boundary_temperature)
const struct sdis_interface* interf = NULL;
const struct sdis_medium* mdm_front = NULL;
const struct sdis_medium* mdm_back = NULL;
+ const struct sdis_medium* mdm = NULL;
double tmp;
ASSERT(scn && fp_to_meter > 0 && ctx && rwalk && rng && T);
ASSERT(rwalk->mdm == NULL);
ASSERT(!SXD_HIT_NONE(&rwalk->hit));
- XD(setup_interface_fragment)(&frag, &rwalk->vtx, &rwalk->hit);
+ XD(setup_interface_fragment)(&frag, &rwalk->vtx, &rwalk->hit, rwalk->hit_side);
+
+ fX(normalize)(rwalk->hit.normal, rwalk->hit.normal);
/* Retrieve the current interface */
interf = scene_get_interface(scn, rwalk->hit.prim.prim_id);
- /* Check if the boundary condition is known */
- tmp = interface_get_temperature(interf, &frag);
+ /* Check if the boundary temperature is known */
+ tmp = interface_side_get_temperature(interf, &frag);
if(tmp >= 0) {
T->value += tmp;
T->done = 1;
return RES_OK;
}
+ /* Check if the boundary flux is known. Note that actually, only solid media
+ * can have a flux as limit condition */
+ mdm = interface_get_medium(interf, frag.side);
+ if(sdis_medium_get_type(mdm) == SDIS_SOLID ) {
+ const double phi = interface_side_get_flux(interf, &frag);
+ if(phi != SDIS_FLUX_NONE) {
+ XD(solid_boundary_with_flux_temperature)
+ (scn, fp_to_meter, ctx, &frag, phi, rwalk, rng, T);
+ return RES_OK;
+ }
+ }
+
mdm_front = interface_get_medium(interf, SDIS_FRONT);
mdm_back = interface_get_medium(interf, SDIS_BACK);
@@ -700,21 +1047,22 @@ XD(solid_temperature)
double position_start[DIM];
const struct sdis_medium* mdm;
ASSERT(scn && fp_to_meter > 0 && rwalk && rng && T);
- ASSERT(rwalk->mdm->type == SDIS_MEDIUM_SOLID);
+ ASSERT(rwalk->mdm->type == SDIS_SOLID);
(void)ctx;
/* Check the random walk consistency */
- CHK(scene_get_medium(scn, rwalk->vtx.P, &mdm) == RES_OK);
+ CHK(scene_get_medium(scn, rwalk->vtx.P, NULL, &mdm) == RES_OK);
if(mdm != rwalk->mdm) {
log_err(scn->dev, "%s: invalid solid random walk. "
"Unexpected medium at {%g, %g, %g}.\n",
FUNC_NAME, SPLIT3(rwalk->vtx.P));
- return RES_BAD_OP;
+ return RES_BAD_OP_IRRECOVERABLE;
}
/* Save the submitted position */
dX(set)(position_start, rwalk->vtx.P);
do { /* Solid random walk */
+ struct get_medium_info info;
struct sXd(hit) hit0, hit1;
double lambda; /* Thermal conductivity */
double rho; /* Volumic mass */
@@ -740,6 +1088,7 @@ XD(solid_temperature)
lambda = solid_get_thermal_conductivity(mdm, &rwalk->vtx);
rho = solid_get_volumic_mass(mdm, &rwalk->vtx);
cp = solid_get_calorific_capacity(mdm, &rwalk->vtx);
+ power = solid_get_volumic_power(mdm, &rwalk->vtx);
#if (SDIS_SOLVE_DIMENSION == 2)
/* Sample a direction around 2PI */
@@ -761,9 +1110,57 @@ XD(solid_temperature)
if(SXD_HIT_NONE(&hit0) && SXD_HIT_NONE(&hit1)) {
/* Hit nothing: move along dir0 of the original delta */
delta = delta_solid;
+
+ /* Add the volumic power density to the measured temperature */
+ if(power != SDIS_VOLUMIC_POWER_NONE) {
+ const double delta_in_meter = delta * fp_to_meter;
+ tmp = power * delta_in_meter * delta_in_meter / (2.0 * DIM * lambda);
+ T->value += tmp;
+ }
} else {
/* Hit something: move along dir0 of the minimum hit distance */
delta = MMIN(hit0.distance, hit1.distance);
+
+ /* Add the volumic power density to the measured temperature */
+ if(power != SDIS_VOLUMIC_POWER_NONE) {
+ const double delta_s_in_meter = delta_solid * fp_to_meter;
+ double h;
+ double h_in_meter;
+ double cos_U_N;
+ float N[DIM];
+
+ if(delta == hit0.distance) {
+ fX(normalize)(N, hit0.normal);
+ cos_U_N = fX(dot)(dir0, N);
+ } else {
+ ASSERT(delta == hit1.distance);
+ fX(normalize)(N, hit1.normal);
+ cos_U_N = fX(dot)(dir1, N);
+ }
+
+ h = delta * fabs(cos_U_N);
+ h_in_meter = h * fp_to_meter;
+
+ /* The regular power term at wall */
+ tmp = power * h_in_meter * h_in_meter / (2.0 * lambda);
+
+ /* Add the power corrective term */
+ if(h < delta_solid) {
+ const double sin_a = h / delta_solid;
+#if DIM==2
+ /* tmp1 = sin(2a) / (PI - 2*a) */
+ const double tmp1 = sin_a * sqrt(1 - sin_a*sin_a)/acos(sin_a);
+ tmp += -(power*delta_s_in_meter*delta_s_in_meter)/(4.0*lambda) * tmp1;
+#else
+ const double tmp1 = (sin_a*sin_a*sin_a - sin_a)/ (1-sin_a);
+ tmp += (power*delta_s_in_meter*delta_s_in_meter)/(6*lambda) * tmp1;
+#endif
+
+ } else if (h == delta_solid) {
+ tmp += -(delta_s_in_meter*delta_s_in_meter*power)/(2.0*DIM*lambda);
+ }
+ T->value += tmp;
+ }
}
/* Sample the time */
@@ -788,44 +1185,53 @@ XD(solid_temperature)
return RES_BAD_OP;
}
- /* Add the volumic power density to the measured temperature */
- power = solid_get_volumic_power(mdm, &rwalk->vtx);
- if(power > 0) {
- const double delta_in_meter = delta * fp_to_meter;
- tmp = power * delta_in_meter * delta_in_meter / (2.0 * DIM * lambda);
- T->value += tmp;
- }
-
/* Define if the random walk hits something along dir0 */
- rwalk->hit = hit0.distance > delta ? SXD_HIT_NULL : hit0;
+ if(hit0.distance > delta) {
+ rwalk->hit = SXD_HIT_NULL;
+ rwalk->hit_side = SDIS_SIDE_NULL__;
+ } else {
+ rwalk->hit = hit0;
+ rwalk->hit_side = fX(dot)(hit0.normal, dir0) < 0 ? SDIS_FRONT : SDIS_BACK;
+ }
/* Update the random walk position */
XD(move_pos)(rwalk->vtx.P, dir0, delta);
/* Fetch the current medium */
if(SXD_HIT_NONE(&rwalk->hit)) {
- CHK(scene_get_medium(scn, rwalk->vtx.P, &mdm) == RES_OK);
+ CHK(scene_get_medium(scn, rwalk->vtx.P, &info, &mdm) == RES_OK);
} else {
const struct sdis_interface* interf;
interf = scene_get_interface(scn, rwalk->hit.prim.prim_id);
- mdm = interface_get_medium
- (interf,
- fX(dot(rwalk->hit.normal, dir0)) < 0 ? SDIS_FRONT : SDIS_BACK);
+ mdm = interface_get_medium(interf, rwalk->hit_side);
}
/* Check random walk consistency */
if(mdm != rwalk->mdm) {
log_err(scn->dev,
"%s: inconsistent medium during the solid random walk.\n", FUNC_NAME);
- if(DIM == 2) {
- log_err(scn->dev,
- " start position: %g %g; current position: %g %g\n",
- SPLIT2(position_start), SPLIT2(rwalk->vtx.P));
- } else {
- log_err(scn->dev,
- " start position: %g %g %g; current position: %g %g %g\n",
- SPLIT3(position_start), SPLIT3(rwalk->vtx.P));
+#if DIM == 2
+ #define VEC_STR "%g %g"
+ #define VEC_SPLIT SPLIT2
+#else
+ #define VEC_STR "%g %g %g"
+ #define VEC_SPLIT SPLIT3
+#endif
+ log_err(scn->dev,
+ " start position: " VEC_STR "; current position: " VEC_STR "\n",
+ VEC_SPLIT(position_start), VEC_SPLIT(rwalk->vtx.P));
+ if(SXD_HIT_NONE(&rwalk->hit)) {
+ float hit_pos[DIM];
+ fX(mulf)(hit_pos, info.ray_dir, info.XD(hit).distance);
+ fX(add)(hit_pos, info.ray_org, hit_pos);
+ log_err(scn->dev, " ray org: " VEC_STR "; ray dir: " VEC_STR "\n",
+ VEC_SPLIT(info.ray_org), VEC_SPLIT(info.ray_dir));
+ log_err(scn->dev, " targeted point: " VEC_STR "\n",
+ VEC_SPLIT(info.pos_tgt));
+ log_err(scn->dev, " hit pos: " VEC_STR "\n", VEC_SPLIT(hit_pos));
}
+#undef VEC_STR
+#undef VEC_SPLIT
return RES_BAD_OP;
}
@@ -847,27 +1253,46 @@ XD(compute_temperature)
struct XD(temperature)* T)
{
#ifndef NDEBUG
- struct XD(temperature)* stack = NULL;
+ struct entry {
+ struct XD(temperature) temperature;
+ struct XD(rwalk) rwalk;
+ }* stack = NULL;
size_t istack = 0;
#endif
+ const size_t max_fails = 10;
res_T res = RES_OK;
ASSERT(scn && fp_to_meter > 0 && ctx && rwalk && rng && T);
do {
- res = T->func(scn, fp_to_meter, ctx, rwalk, rng, T);
- if(res != RES_OK) goto error;
+ /* Save the current random walk state */
+ const struct XD(rwalk) rwalk_bkp = *rwalk;
+ const struct XD(temperature) T_bkp = *T;
+
+ size_t nfails = 0;
#ifndef NDEBUG
- sa_push(stack, *T);
+ struct entry e;
+ e.temperature = *T;
+ e.rwalk = *rwalk;
+ sa_push(stack, e);
++istack;
#endif
+
+ /* Reject the current step if a BAD_OP occurs and retry up to "max_fails"
+ * times */
+ do {
+ res = T->func(scn, fp_to_meter, ctx, rwalk, rng, T);
+ if(res == RES_BAD_OP) { *rwalk = rwalk_bkp; *T = T_bkp; }
+ } while(res == RES_BAD_OP && ++nfails < max_fails);
+ if(res != RES_OK) goto error;
+
} while(!T->done);
exit:
#ifndef NDEBUG
sa_release(stack);
#endif
- return res;
+ return res == RES_BAD_OP_IRRECOVERABLE ? RES_BAD_OP : res;
error:
goto exit;
}
@@ -891,8 +1316,8 @@ XD(probe_realisation)
ASSERT(medium && position && fp_to_meter > 0 && weight && time >= 0);
switch(medium->type) {
- case SDIS_MEDIUM_FLUID: T.func = XD(fluid_temperature); break;
- case SDIS_MEDIUM_SOLID: T.func = XD(solid_temperature); break;
+ case SDIS_FLUID: T.func = XD(fluid_temperature); break;
+ case SDIS_SOLID: T.func = XD(solid_temperature); break;
default: FATAL("Unreachable code\n"); break;
}
@@ -921,6 +1346,7 @@ XD(boundary_realisation)
const size_t iprim,
const double uv[DIM],
const double time,
+ const enum sdis_side side,
const double fp_to_meter,
const double Tarad,
const double Tref,
@@ -940,6 +1366,7 @@ XD(boundary_realisation)
T.func = XD(boundary_temperature);
+ rwalk.hit_side = side;
rwalk.hit.distance = 0;
rwalk.vtx.time = time;
rwalk.mdm = NULL; /* The random walk is at an interface between 2 media */
@@ -998,11 +1425,12 @@ XD(ray_realisation)
float dir[3];
res_T res = RES_OK;
ASSERT(scn && position && direction && time>=0 && fp_to_meter>0 && weight);
- ASSERT(medium && medium->type == SDIS_MEDIUM_FLUID);
+ ASSERT(medium && medium->type == SDIS_FLUID);
dX(set)(rwalk.vtx.P, position);
rwalk.vtx.time = time;
rwalk.hit = SXD_HIT_NULL;
+ rwalk.hit_side = SDIS_SIDE_NULL__;
rwalk.mdm = medium;
ctx.Tarad = Tarad;
diff --git a/src/test_sdis_accum_buffer.c b/src/test_sdis_accum_buffer.c
@@ -60,7 +60,7 @@ main(int argc, char** argv)
accums_tmp = MEM_CALLOC
(&allocator, layout.width*layout.height, sizeof(struct sdis_accum));
CHK(accums_tmp != NULL);
- memcpy(accums_tmp, accums_tmp,
+ memmove(accums_tmp, accums_tmp,
layout.width*layout.height*sizeof(struct sdis_accum));
MEM_RM(&allocator, accums_tmp);
diff --git a/src/test_sdis_conducto_radiative.c b/src/test_sdis_conducto_radiative.c
@@ -168,18 +168,10 @@ solid_get_delta
return 1.0/10.0;
}
-static double
-solid_get_delta_boundary
- (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
-{
- CHK(vtx != NULL); (void)data;
- return 2.1/10.0;
-}
-
/*******************************************************************************
* Interface
******************************************************************************/
-struct interface {
+struct interfac {
double temperature;
double convection_coef;
double emissivity;
@@ -191,7 +183,7 @@ interface_get_temperature
(const struct sdis_interface_fragment* frag, struct sdis_data* data)
{
CHK(data != NULL && frag != NULL);
- return ((const struct interface*)sdis_data_cget(data))->temperature;
+ return ((const struct interfac*)sdis_data_cget(data))->temperature;
}
static double
@@ -199,7 +191,7 @@ interface_get_convection_coef
(const struct sdis_interface_fragment* frag, struct sdis_data* data)
{
CHK(data != NULL && frag != NULL);
- return ((const struct interface*)sdis_data_cget(data))->convection_coef;
+ return ((const struct interfac*)sdis_data_cget(data))->convection_coef;
}
static double
@@ -207,7 +199,7 @@ interface_get_emissivity
(const struct sdis_interface_fragment* frag, struct sdis_data* data)
{
CHK(data != NULL && frag != NULL);
- return ((const struct interface*)sdis_data_cget(data))->emissivity;
+ return ((const struct interfac*)sdis_data_cget(data))->emissivity;
}
static double
@@ -215,7 +207,7 @@ interface_get_specular_fraction
(const struct sdis_interface_fragment* frag, struct sdis_data* data)
{
CHK(data != NULL && frag != NULL);
- return ((const struct interface*)sdis_data_cget(data))->specular_fraction;
+ return ((const struct interfac*)sdis_data_cget(data))->specular_fraction;
}
/*******************************************************************************
@@ -226,22 +218,31 @@ create_interface
(struct sdis_device* dev,
struct sdis_medium* front,
struct sdis_medium* back,
- const struct interface* interf,
+ const struct interfac* interf,
struct sdis_interface** out_interf)
{
- struct sdis_interface_shader shader = DUMMY_INTERFACE_SHADER;
+ struct sdis_interface_shader shader = SDIS_INTERFACE_SHADER_NULL;
struct sdis_data* data = NULL;
CHK(interf != NULL);
- shader.temperature = interface_get_temperature;
- shader.convection_coef = interface_get_convection_coef;
- shader.emissivity = interface_get_emissivity;
- shader.specular_fraction = interface_get_specular_fraction;
+ shader.front.temperature = interface_get_temperature;
+ shader.back.temperature = interface_get_temperature;
+ if(sdis_medium_get_type(front) != sdis_medium_get_type(back)) {
+ shader.convection_coef = interface_get_convection_coef;
+ }
+ if(sdis_medium_get_type(front) == SDIS_FLUID) {
+ shader.front.emissivity = interface_get_emissivity;
+ shader.front.specular_fraction = interface_get_specular_fraction;
+ }
+ if(sdis_medium_get_type(back) == SDIS_FLUID) {
+ shader.back.emissivity = interface_get_emissivity;
+ shader.back.specular_fraction = interface_get_specular_fraction;
+ }
- CHK(sdis_data_create(dev, sizeof(struct interface), ALIGNOF(struct interface),
+ CHK(sdis_data_create(dev, sizeof(struct interfac), ALIGNOF(struct interfac),
NULL, &data) == RES_OK);
- *((struct interface*)sdis_data_get(data)) = *interf;
+ *((struct interfac*)sdis_data_get(data)) = *interf;
CHK(sdis_interface_create(dev, front, back, &shader, data, out_interf) == RES_OK);
CHK(sdis_data_ref_put(data) == RES_OK);
@@ -254,7 +255,7 @@ int
main(int argc, char** argv)
{
struct mem_allocator allocator;
- struct interface interf;
+ struct interfac interf;
struct geometry geom;
struct sdis_data* data = NULL;
struct sdis_device* dev = NULL;
@@ -294,7 +295,6 @@ main(int argc, char** argv)
solid_shader.thermal_conductivity = solid_get_thermal_conductivity;
solid_shader.volumic_mass = solid_get_volumic_mass;
solid_shader.delta_solid = solid_get_delta;
- solid_shader.delta_boundary = solid_get_delta_boundary;
solid_shader.temperature = temperature_unknown;
CHK(sdis_solid_create(dev, &solid_shader, data, &solid) == RES_OK);
CHK(sdis_data_ref_put(data) == RES_OK);
@@ -307,7 +307,6 @@ main(int argc, char** argv)
solid_shader.thermal_conductivity = solid_get_thermal_conductivity;
solid_shader.volumic_mass = solid_get_volumic_mass;
solid_shader.delta_solid = solid_get_delta;
- solid_shader.delta_boundary = solid_get_delta_boundary;
solid_shader.temperature = temperature_unknown;
CHK(sdis_solid_create(dev, &solid_shader, data, &solid2) == RES_OK);
CHK(sdis_data_ref_put(data) == RES_OK);
@@ -390,12 +389,13 @@ main(int argc, char** argv)
double ref, u;
size_t nreals = 0;
size_t nfails = 0;
+ const size_t N = 10000;
pos[0] = ssp_rng_uniform_double(rng, -0.9, 0.9);
pos[1] = ssp_rng_uniform_double(rng, -0.9, 0.9);
pos[2] = ssp_rng_uniform_double(rng, -0.9, 0.9);
- CHK(sdis_solve_probe(scn, 10000, pos, INF, 1, -1, Tref, &estimator) == RES_OK);
+ CHK(sdis_solve_probe(scn, N, pos, INF, 1, -1, Tref, &estimator) == RES_OK);
CHK(sdis_estimator_get_realisation_count(estimator, &nreals) == RES_OK);
CHK(sdis_estimator_get_failure_count(estimator, &nfails) == RES_OK);
CHK(sdis_estimator_get_temperature(estimator, &T) == RES_OK);
@@ -404,7 +404,10 @@ main(int argc, char** argv)
ref = u * Ts1 + (1-u) * Ts0;
printf("Temperature at (%g, %g, %g) = %g ~ %g +/- %g\n",
SPLIT3(pos), ref, T.E, T.SE);
+ printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N);
+ CHK(nfails + nreals == N);
+ CHK(nfails < N/1000);
CHK(eq_eps(T.E, ref, 2*T.SE) == 1);
CHK(sdis_estimator_ref_put(estimator) == RES_OK);
diff --git a/src/test_sdis_conducto_radiative_2d.c b/src/test_sdis_conducto_radiative_2d.c
@@ -149,54 +149,78 @@ solid_get_delta
return 1.0/10.0;
}
-static double
-solid_get_delta_boundary
- (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
-{
- CHK(vtx != NULL); (void)data;
- return 2.1/10.0;
-}
-
/*******************************************************************************
* Interface
******************************************************************************/
-struct interface {
- double temperature;
+struct interfac {
double convection_coef;
- double emissivity;
- double specular_fraction;
+ struct {
+ double temperature;
+ double emissivity;
+ double specular_fraction;
+ } front, back;
+};
+
+static const struct interfac INTERFACE_NULL = {
+ 0, {-1, -1, -1}, {-1, -1, -1}
};
static double
interface_get_temperature
(const struct sdis_interface_fragment* frag, struct sdis_data* data)
{
+ const struct interfac* interf;
+ double T = -1;
CHK(data != NULL && frag != NULL);
- return ((const struct interface*)sdis_data_cget(data))->temperature;
+ interf = sdis_data_cget(data);
+ switch(frag->side) {
+ case SDIS_FRONT: T = interf->front.temperature; break;
+ case SDIS_BACK: T = interf->back.temperature; break;
+ default: FATAL("Unreachable code.\n"); break;
+ }
+ return T;
}
static double
interface_get_convection_coef
(const struct sdis_interface_fragment* frag, struct sdis_data* data)
{
+ const struct interfac* interf;
CHK(data != NULL && frag != NULL);
- return ((const struct interface*)sdis_data_cget(data))->convection_coef;
+ interf = sdis_data_cget(data);
+ return interf->convection_coef;
}
static double
interface_get_emissivity
(const struct sdis_interface_fragment* frag, struct sdis_data* data)
{
+ const struct interfac* interf;
+ double e = -1;
CHK(data != NULL && frag != NULL);
- return ((const struct interface*)sdis_data_cget(data))->emissivity;
+ interf = sdis_data_cget(data);
+ switch(frag->side) {
+ case SDIS_FRONT: e = interf->front.emissivity; break;
+ case SDIS_BACK: e = interf->back.emissivity; break;
+ default: FATAL("Unreachable code.\n"); break;
+ }
+ return e;
}
static double
interface_get_specular_fraction
(const struct sdis_interface_fragment* frag, struct sdis_data* data)
{
+ const struct interfac* interf;
+ double f = -1;
CHK(data != NULL && frag != NULL);
- return ((const struct interface*)sdis_data_cget(data))->specular_fraction;
+ interf = sdis_data_cget(data);
+ switch(frag->side) {
+ case SDIS_FRONT: f = interf->front.specular_fraction; break;
+ case SDIS_BACK: f = interf->back.specular_fraction; break;
+ default: FATAL("Unreachable code.\n"); break;
+ }
+ return f;
}
/*******************************************************************************
@@ -207,28 +231,38 @@ create_interface
(struct sdis_device* dev,
struct sdis_medium* front,
struct sdis_medium* back,
- const struct interface* interf,
+ const struct interfac* interf,
struct sdis_interface** out_interf)
{
- struct sdis_interface_shader shader = DUMMY_INTERFACE_SHADER;
+ struct sdis_interface_shader shader = SDIS_INTERFACE_SHADER_NULL;
struct sdis_data* data = NULL;
+ const enum sdis_medium_type type_f = sdis_medium_get_type(front);
+ const enum sdis_medium_type type_b = sdis_medium_get_type(back);
CHK(interf != NULL);
- shader.temperature = interface_get_temperature;
- shader.convection_coef = interface_get_convection_coef;
- shader.emissivity = interface_get_emissivity;
- shader.specular_fraction = interface_get_specular_fraction;
+ shader.back.temperature = interface_get_temperature;
+ shader.front.temperature = interface_get_temperature;
- CHK(sdis_data_create(dev, sizeof(struct interface), ALIGNOF(struct interface),
+ if(type_f != type_b) {
+ shader.convection_coef = interface_get_convection_coef;
+ }
+ if(type_f == SDIS_FLUID) {
+ shader.front.emissivity = interface_get_emissivity;
+ shader.front.specular_fraction = interface_get_specular_fraction;
+ }
+ if(type_b == SDIS_FLUID) {
+ shader.back.emissivity = interface_get_emissivity;
+ shader.back.specular_fraction = interface_get_specular_fraction;
+ }
+ CHK(sdis_data_create(dev, sizeof(struct interfac), ALIGNOF(struct interfac),
NULL, &data) == RES_OK);
- *((struct interface*)sdis_data_get(data)) = *interf;
+ *((struct interfac*)sdis_data_get(data)) = *interf;
CHK(sdis_interface_create(dev, front, back, &shader, data, out_interf) == RES_OK);
CHK(sdis_data_ref_put(data) == RES_OK);
}
-
/*******************************************************************************
* Test
******************************************************************************/
@@ -236,7 +270,7 @@ int
main(int argc, char** argv)
{
struct mem_allocator allocator;
- struct interface interf;
+ struct interfac interf;
struct geometry geom;
struct ssp_rng* rng = NULL;
struct sdis_scene* scn = NULL;
@@ -275,7 +309,6 @@ main(int argc, char** argv)
solid_shader.thermal_conductivity = solid_get_thermal_conductivity;
solid_shader.volumic_mass = solid_get_volumic_mass;
solid_shader.delta_solid = solid_get_delta;
- solid_shader.delta_boundary = solid_get_delta_boundary;
solid_shader.temperature = temperature_unknown;
CHK(sdis_solid_create(dev, &solid_shader, data, &solid) == RES_OK);
CHK(sdis_data_ref_put(data) == RES_OK);
@@ -288,43 +321,39 @@ main(int argc, char** argv)
solid_shader.thermal_conductivity = solid_get_thermal_conductivity;
solid_shader.volumic_mass = solid_get_volumic_mass;
solid_shader.delta_solid = solid_get_delta;
- solid_shader.delta_boundary = solid_get_delta_boundary;
CHK(sdis_solid_create(dev, &solid_shader, data, &solid2) == RES_OK);
CHK(sdis_data_ref_put(data) == RES_OK);
/* Create the interface that forces to keep in conduction */
- interf.temperature = UNKNOWN_TEMPERATURE;
- interf.convection_coef = -1;
- interf.emissivity = -1;
- interf.specular_fraction = -1;
+ interf = INTERFACE_NULL;
create_interface(dev, solid, solid2, &interf, interfaces+0);
/* Create the interface that emits radiative heat from the solid */
- interf.temperature = UNKNOWN_TEMPERATURE;
- interf.convection_coef = 0;
- interf.emissivity = emissivity;
- interf.specular_fraction = -1;
+ interf = INTERFACE_NULL;
+ interf.back.temperature = UNKNOWN_TEMPERATURE;
+ interf.back.emissivity = emissivity;
+ interf.back.specular_fraction = -1; /* Should not be fetched */
create_interface(dev, solid, fluid, &interf, interfaces+1);
/* Create the interface that forces the radiative heat to bounce */
- interf.temperature = UNKNOWN_TEMPERATURE;
- interf.convection_coef = 0;
- interf.emissivity = 0;
- interf.specular_fraction = 1;
+ interf = INTERFACE_NULL;
+ interf.front.temperature = UNKNOWN_TEMPERATURE;
+ interf.front.emissivity = 0;
+ interf.front.specular_fraction = 1;
create_interface(dev, fluid, solid2, &interf, interfaces+2);
/* Create the interface with a limit condition of T0 Kelvin */
- interf.temperature = T0;
- interf.convection_coef = 0;
- interf.emissivity = 1;
- interf.specular_fraction = 1;
+ interf = INTERFACE_NULL;
+ interf.front.temperature = T0;
+ interf.front.emissivity = 1;
+ interf.front.specular_fraction = 1;
create_interface(dev, fluid, solid2, &interf, interfaces+3);
/* Create the interface with a limit condition of T1 Kelvin */
- interf.temperature = T1;
- interf.convection_coef = 0;
- interf.emissivity = 1;
- interf.specular_fraction = 1;
+ interf = INTERFACE_NULL;
+ interf.front.temperature = T1;
+ interf.front.emissivity = 1;
+ interf.front.specular_fraction = 1;
create_interface(dev, fluid, solid2, &interf, interfaces+4);
/* Setup the per primitive interface of the solid medium */
@@ -364,6 +393,7 @@ main(int argc, char** argv)
double ref, u;
size_t nreals = 0;
size_t nfails = 0;
+ const size_t N = 10000;
pos[0] = ssp_rng_uniform_double(rng, -0.9, 0.9);
pos[1] = ssp_rng_uniform_double(rng, -0.9, 0.9);
@@ -377,8 +407,11 @@ main(int argc, char** argv)
ref = u * Ts1 + (1-u) * Ts0;
printf("Temperature at (%g, %g) = %g ~ %g +/- %g\n",
SPLIT2(pos), ref, T.E, T.SE);
+ printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N);
- CHK(eq_eps(T.E, ref, 2*T.SE) == 1);
+ CHK(nfails + nreals == N);
+ CHK(nfails < N/1000);
+ CHK(eq_eps(T.E, ref, 3*T.SE) == 1);
CHK(sdis_estimator_ref_put(estimator) == RES_OK);
}
diff --git a/src/test_sdis_convection.c b/src/test_sdis_convection.c
@@ -204,10 +204,10 @@ main(int argc, char** argv)
CHK(sdis_solid_create(dev, &solid_shader, NULL, &solid) == RES_OK);
/* Setup the interface shader */
- interf_shader.temperature = interface_get_temperature;
interf_shader.convection_coef = interface_get_convection_coef;
- interf_shader.emissivity = interface_get_emissivity;
- interf_shader.specular_fraction = interface_get_specular_fraction;
+ interf_shader.front.temperature = interface_get_temperature;
+ interf_shader.front.emissivity = interface_get_emissivity;
+ interf_shader.front.specular_fraction = interface_get_specular_fraction;
/* Create the interfaces */
interf_T0 = create_interface(dev, fluid, solid, &interf_shader, T0);
diff --git a/src/test_sdis_flux.c b/src/test_sdis_flux.c
@@ -0,0 +1,272 @@
+/* Copyright (C) 2016-2018 |Meso|Star> (contact@meso-star.com)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. */
+
+#include "sdis.h"
+#include "test_sdis_utils.h"
+
+#include <rsys/double3.h>
+
+/*
+ * The scene is composed of a solid cube/square whose temperature is unknown.
+ * The temperature is fixed at T0 on the +X face. The Flux of the -X face is
+ * fixed to PHI. The flux on the other faces is null (i.e. adiabatic). This
+ * test computes the temperature of a probe position pos into the solid and
+ * check that it is equal to:
+ *
+ * T(pos) = T0 + (A-pos) * PHI/LAMBDA
+ *
+ * with LAMBDA the conductivity of the solid and A the size of cube/square.
+ *
+ * 3D 2D
+ *
+ * ///// (1,1,1) ///// (1,1)
+ * +-------+ +-------+
+ * /' /| | |
+ * +-------+ T0 PHI T0
+ * PHI +.....|.+ | |
+ * |, |/ +-------+
+ * +-------+ (0,0) /////
+ * (0,0,0) /////
+ */
+
+#define UNKNOWN_TEMPERATURE -1
+#define N 10000
+
+#define PHI 10.0
+#define T0 320.0
+#define LAMBDA 0.1
+
+/*******************************************************************************
+ * Media
+ ******************************************************************************/
+static double
+solid_get_calorific_capacity
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ (void)data;
+ CHK(vtx != NULL);
+ return 2.0;
+}
+
+static double
+solid_get_thermal_conductivity
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ (void)data;
+ CHK(vtx != NULL);
+ return LAMBDA;
+}
+
+static double
+solid_get_volumic_mass
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ (void)data;
+ CHK(vtx != NULL);
+ return 25.0;
+}
+
+static double
+solid_get_delta
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ (void)data;
+ CHK(vtx != NULL);
+ return 1.0/20.0;
+}
+
+static double
+solid_get_temperature
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ (void)data;
+ CHK(vtx != NULL);
+ return UNKNOWN_TEMPERATURE;
+}
+
+/*******************************************************************************
+ * Interfaces
+ ******************************************************************************/
+struct interf {
+ double temperature;
+ double phi;
+};
+
+static double
+interface_get_temperature
+ (const struct sdis_interface_fragment* frag, struct sdis_data* data)
+{
+ const struct interf* interf = sdis_data_cget(data);
+ CHK(frag && data);
+ return interf->temperature;
+}
+
+static double
+interface_get_flux
+ (const struct sdis_interface_fragment* frag, struct sdis_data* data)
+{
+ const struct interf* interf = sdis_data_cget(data);
+ CHK(frag && data);
+ return interf->phi;
+}
+
+/*******************************************************************************
+ * Test
+ ******************************************************************************/
+int
+main(int argc, char** argv)
+{
+ struct mem_allocator allocator;
+ struct sdis_mc T = SDIS_MC_NULL;
+ struct sdis_data* data = NULL;
+ struct sdis_device* dev = NULL;
+ struct sdis_medium* fluid = NULL;
+ struct sdis_medium* solid = NULL;
+ struct sdis_interface* interf_adiabatic = NULL;
+ struct sdis_interface* interf_T0 = NULL;
+ struct sdis_interface* interf_phi = NULL;
+ struct sdis_scene* box_scn = NULL;
+ struct sdis_scene* square_scn = NULL;
+ struct sdis_estimator* estimator = NULL;
+ struct sdis_fluid_shader fluid_shader = DUMMY_FLUID_SHADER;
+ struct sdis_solid_shader solid_shader = DUMMY_SOLID_SHADER;
+ struct sdis_interface_shader interf_shader = SDIS_INTERFACE_SHADER_NULL;
+ struct sdis_interface* box_interfaces[12 /*#triangles*/];
+ struct sdis_interface* square_interfaces[4/*#segments*/];
+ struct interf* interf_props = NULL;
+ double pos[3];
+ double ref;
+ size_t nreals;
+ size_t nfails;
+ (void)argc, (void)argv;
+
+ CHK(mem_init_proxy_allocator(&allocator, &mem_default_allocator) == RES_OK);
+ CHK(sdis_device_create
+ (NULL, &allocator, SDIS_NTHREADS_DEFAULT, 0, &dev) == RES_OK);
+
+ /* Create the dummy fluid medium */
+ CHK(sdis_fluid_create(dev, &fluid_shader, NULL, &fluid) == RES_OK);
+
+ /* Create the solid_medium */
+ solid_shader.calorific_capacity = solid_get_calorific_capacity;
+ solid_shader.thermal_conductivity = solid_get_thermal_conductivity;
+ solid_shader.volumic_mass = solid_get_volumic_mass;
+ solid_shader.delta_solid = solid_get_delta;
+ solid_shader.temperature = solid_get_temperature;
+ CHK(sdis_solid_create(dev, &solid_shader, NULL, &solid) == RES_OK);
+
+ /* Setup the interface shader */
+ interf_shader.front.temperature = interface_get_temperature;
+ interf_shader.front.flux = interface_get_flux;
+
+ /* Create the adiabatic interface */
+ CHK(sdis_data_create(dev, sizeof(struct interf), 16, NULL, &data) == RES_OK);
+ interf_props = sdis_data_get(data);
+ interf_props->temperature = UNKNOWN_TEMPERATURE;
+ interf_props->phi = 0;
+ CHK(sdis_interface_create
+ (dev, solid, fluid, &interf_shader, data, &interf_adiabatic) == RES_OK);
+ CHK(sdis_data_ref_put(data) == RES_OK);
+
+ /* Create the T0 interface */
+ CHK(sdis_data_create(dev, sizeof(struct interf), 16, NULL, &data) == RES_OK);
+ interf_props = sdis_data_get(data);
+ interf_props->temperature = T0;
+ interf_props->phi = 0; /* Unused */
+ CHK(sdis_interface_create
+ (dev, solid, fluid, &interf_shader, data, &interf_T0) == RES_OK);
+ CHK(sdis_data_ref_put(data) == RES_OK);
+
+ /* Create the PHI interface */
+ CHK(sdis_data_create(dev, sizeof(struct interf), 16, NULL, &data) == RES_OK);
+ interf_props = sdis_data_get(data);
+ interf_props->temperature = UNKNOWN_TEMPERATURE;
+ interf_props->phi = PHI;
+ CHK(sdis_interface_create
+ (dev, solid, fluid, &interf_shader, data, &interf_phi) == RES_OK);
+ CHK(sdis_data_ref_put(data) == RES_OK);
+
+ /* Release the media */
+ CHK(sdis_medium_ref_put(solid) == RES_OK);
+ CHK(sdis_medium_ref_put(fluid) == RES_OK);
+
+ /* Map the interfaces to their box triangles */
+ box_interfaces[0] = box_interfaces[1] = interf_adiabatic; /* Front */
+ box_interfaces[2] = box_interfaces[3] = interf_phi; /* Left */
+ box_interfaces[4] = box_interfaces[5] = interf_adiabatic; /* Back */
+ box_interfaces[6] = box_interfaces[7] = interf_T0; /* Right */
+ box_interfaces[8] = box_interfaces[9] = interf_adiabatic; /* Top */
+ box_interfaces[10]= box_interfaces[11]= interf_adiabatic; /* Bottom */
+
+ /* Map the interfaces to their square segments */
+ square_interfaces[0] = interf_adiabatic; /* Bottom */
+ square_interfaces[1] = interf_phi; /* Left */
+ square_interfaces[2] = interf_adiabatic; /* Top */
+ square_interfaces[3] = interf_T0; /* Right */
+
+ /* Create the box scene */
+ CHK(sdis_scene_create(dev, box_ntriangles, box_get_indices,
+ box_get_interface, box_nvertices, box_get_position, box_interfaces,
+ &box_scn) == RES_OK);
+
+ /* Create the square scene */
+ CHK(sdis_scene_2d_create(dev, square_nsegments, square_get_indices,
+ square_get_interface, square_nvertices, square_get_position,
+ square_interfaces, &square_scn) == RES_OK);
+
+ /* Release the interfaces */
+ CHK(sdis_interface_ref_put(interf_adiabatic) == RES_OK);
+ CHK(sdis_interface_ref_put(interf_T0) == RES_OK);
+ CHK(sdis_interface_ref_put(interf_phi) == RES_OK);
+
+ d3_splat(pos, 0.25);
+ ref = T0 + (1 - pos[0]) * PHI/LAMBDA;
+
+ /* Solve in 3D */
+ CHK(sdis_solve_probe(box_scn, N, pos, INF, 1.0, 0, 0, &estimator) == RES_OK);
+ CHK(sdis_estimator_get_realisation_count(estimator, &nreals) == RES_OK);
+ CHK(sdis_estimator_get_failure_count(estimator, &nfails) == RES_OK);
+ CHK(sdis_estimator_get_temperature(estimator, &T) == RES_OK);
+ CHK(sdis_estimator_ref_put(estimator) == RES_OK);
+ printf("Temperature of the box at (%g %g %g) = %g ~ %g +/- %g\n",
+ SPLIT3(pos), ref, T.E, T.SE);
+ printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N);
+ CHK(nfails + nreals == N);
+ CHK(nfails < N/1000);
+ CHK(eq_eps(T.E, ref, T.SE*3));
+
+ /* Solve in 2D */
+ CHK(sdis_solve_probe(square_scn, N, pos, INF, 1.0, 0, 0, &estimator) == RES_OK);
+ CHK(sdis_estimator_get_realisation_count(estimator, &nreals) == RES_OK);
+ CHK(sdis_estimator_get_failure_count(estimator, &nfails) == RES_OK);
+ CHK(sdis_estimator_get_temperature(estimator, &T) == RES_OK);
+ CHK(sdis_estimator_ref_put(estimator) == RES_OK);
+ printf("Temperature of the square at (%g %g) = %g ~ %g +/- %g\n",
+ SPLIT2(pos), ref, T.E, T.SE);
+ printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N);
+ CHK(nfails + nreals == N);
+ CHK(nfails < N/1000);
+ CHK(eq_eps(T.E, ref, T.SE*3));
+
+ CHK(sdis_scene_ref_put(box_scn) == RES_OK);
+ CHK(sdis_scene_ref_put(square_scn) == RES_OK);
+ CHK(sdis_device_ref_put(dev) == RES_OK);
+
+ check_memory_allocator(&allocator);
+ mem_shutdown_proxy_allocator(&allocator);
+ CHK(mem_allocated_size() == 0);
+ return 0;
+
+}
diff --git a/src/test_sdis_interface.c b/src/test_sdis_interface.c
@@ -31,11 +31,14 @@ main(int argc, char** argv)
CHK(mem_init_proxy_allocator(&allocator, &mem_default_allocator) == RES_OK);
CHK(sdis_device_create
- (NULL, &allocator, SDIS_NTHREADS_DEFAULT, 0, &dev) == RES_OK);
+ (NULL, &allocator, SDIS_NTHREADS_DEFAULT, 1, &dev) == RES_OK);
CHK(sdis_fluid_create(dev, &fluid_shader, NULL, &fluid) == RES_OK);
CHK(sdis_solid_create(dev, &solid_shader, NULL, &solid) == RES_OK);
+
+ shader = SDIS_INTERFACE_SHADER_NULL;
+
#define CREATE sdis_interface_create
CHK(CREATE(NULL, NULL, NULL, NULL, NULL, NULL) == RES_BAD_ARG);
CHK(CREATE(dev, NULL, NULL, NULL, NULL, NULL) == RES_BAD_ARG);
@@ -78,24 +81,28 @@ main(int argc, char** argv)
CHK(CREATE(dev, solid, solid, &shader, NULL, &interf) == RES_OK);
CHK(sdis_interface_ref_put(interf) == RES_OK);
- shader.convection_coef = NULL;
- shader.specular_fraction = NULL;
- shader.emissivity = NULL;
+ shader = SDIS_INTERFACE_SHADER_NULL;
CHK(CREATE(dev, solid, solid, &shader, NULL, &interf) == RES_OK);
CHK(sdis_interface_ref_put(interf) == RES_OK);
- shader.temperature = NULL;
+ shader.front.temperature = dummy_interface_getter;
CHK(CREATE(dev, solid, solid, &shader, NULL, &interf) == RES_OK);
CHK(sdis_interface_ref_put(interf) == RES_OK);
- CHK(CREATE(dev, solid, fluid, &shader, NULL, &interf) == RES_BAD_ARG);
- shader.convection_coef = DUMMY_INTERFACE_SHADER.convection_coef;
- CHK(CREATE(dev, solid, fluid, &shader, NULL, &interf) == RES_BAD_ARG);
- shader.emissivity = DUMMY_INTERFACE_SHADER.emissivity;
- CHK(CREATE(dev, solid, fluid, &shader, NULL, &interf) == RES_BAD_ARG);
- shader.specular_fraction = DUMMY_INTERFACE_SHADER.specular_fraction;
+ shader.back.emissivity = dummy_interface_getter;
+ CHK(CREATE(dev, solid, fluid, &shader, NULL, &interf) == RES_OK);
+ CHK(sdis_interface_ref_put(interf) == RES_OK);
+ shader.back.specular_fraction = dummy_interface_getter;
CHK(CREATE(dev, solid, fluid, &shader, NULL, &interf) == RES_OK);
CHK(sdis_interface_ref_put(interf) == RES_OK);
+ shader.back = SDIS_INTERFACE_SIDE_SHADER_NULL;
+ shader.front.emissivity = dummy_interface_getter;
+ CHK(CREATE(dev, solid, fluid, &shader, NULL, &interf) == RES_OK); /* Warning */
+ CHK(sdis_interface_ref_put(interf) == RES_OK);
+ shader.front.emissivity = NULL;
+ shader.front.specular_fraction = dummy_interface_getter;
+ CHK(CREATE(dev, solid, fluid, &shader, NULL, &interf) == RES_OK); /* Warning */
+ CHK(sdis_interface_ref_put(interf) == RES_OK);
#undef CREATE
CHK(sdis_device_ref_put(dev) == RES_OK);
diff --git a/src/test_sdis_medium.c b/src/test_sdis_medium.c
@@ -20,6 +20,7 @@ int
main(int argc, char** argv)
{
struct mem_allocator allocator;
+ struct sdis_data* data = NULL;
struct sdis_device* dev = NULL;
struct sdis_medium* fluid = NULL;
struct sdis_medium* solid = NULL;
@@ -40,6 +41,9 @@ main(int argc, char** argv)
CHK(sdis_fluid_create(NULL, &fluid_shader, NULL, &fluid) == RES_BAD_ARG);
CHK(sdis_fluid_create(dev, &fluid_shader, NULL, &fluid) == RES_OK);
+ CHK(sdis_medium_get_type(fluid) == SDIS_FLUID);
+ CHK(sdis_medium_get_data(fluid) == NULL);
+
CHK(sdis_medium_ref_get(NULL) == RES_BAD_ARG);
CHK(sdis_medium_ref_get(fluid) == RES_OK);
CHK(sdis_medium_ref_put(NULL) == RES_BAD_ARG);
@@ -61,15 +65,19 @@ main(int argc, char** argv)
CHK(sdis_fluid_create
(dev, &SDIS_FLUID_SHADER_NULL, NULL, &fluid) == RES_BAD_ARG);
- CHK(sdis_solid_create(NULL, NULL, NULL, NULL) == RES_BAD_ARG);
- CHK(sdis_solid_create(dev, NULL, NULL, NULL) == RES_BAD_ARG);
- CHK(sdis_solid_create(NULL, &solid_shader, NULL, NULL) == RES_BAD_ARG);
- CHK(sdis_solid_create(dev, &solid_shader, NULL, NULL) == RES_BAD_ARG);
- CHK(sdis_solid_create(NULL, NULL, NULL, &solid) == RES_BAD_ARG);
- CHK(sdis_solid_create(dev, NULL, NULL, &solid) == RES_BAD_ARG);
- CHK(sdis_solid_create(NULL, &solid_shader, NULL, &solid) == RES_BAD_ARG);
- CHK(sdis_solid_create(dev, &solid_shader, NULL, &solid) == RES_OK);
+ CHK(sdis_data_create(dev, 4, 16, NULL, &data) == RES_OK);
+ CHK(sdis_solid_create(NULL, NULL, data, NULL) == RES_BAD_ARG);
+ CHK(sdis_solid_create(dev, NULL, data, NULL) == RES_BAD_ARG);
+ CHK(sdis_solid_create(NULL, &solid_shader, data, NULL) == RES_BAD_ARG);
+ CHK(sdis_solid_create(dev, &solid_shader, data, NULL) == RES_BAD_ARG);
+ CHK(sdis_solid_create(NULL, NULL, data, &solid) == RES_BAD_ARG);
+ CHK(sdis_solid_create(dev, NULL, data, &solid) == RES_BAD_ARG);
+ CHK(sdis_solid_create(NULL, &solid_shader, data, &solid) == RES_BAD_ARG);
+ CHK(sdis_solid_create(dev, &solid_shader, data, &solid) == RES_OK);
+ CHK(sdis_medium_get_type(solid) == SDIS_SOLID);
+ CHK(sdis_medium_get_data(solid) == data);
CHK(sdis_medium_ref_put(solid) == RES_OK);
+ CHK(sdis_data_ref_put(data) == RES_OK);
solid_shader.calorific_capacity = NULL;
CHK(sdis_solid_create(dev, &solid_shader, NULL, &solid) == RES_BAD_ARG);
@@ -87,10 +95,6 @@ main(int argc, char** argv)
CHK(sdis_solid_create(dev, &solid_shader, NULL, &solid) == RES_BAD_ARG);
solid_shader.delta_solid = DUMMY_SOLID_SHADER.delta_solid;
- solid_shader.delta_boundary = NULL;
- CHK(sdis_solid_create(dev, &solid_shader, NULL, &solid) == RES_BAD_ARG);
- solid_shader.delta_boundary = DUMMY_SOLID_SHADER.delta_boundary;
-
solid_shader.temperature = NULL;
CHK(sdis_solid_create(dev, &solid_shader, NULL, &solid) == RES_BAD_ARG);
solid_shader.temperature = DUMMY_SOLID_SHADER.temperature;
diff --git a/src/test_sdis_scene.c b/src/test_sdis_scene.c
@@ -15,6 +15,9 @@
#include "sdis.h"
#include "test_sdis_utils.h"
+
+#include <rsys/double2.h>
+#include <rsys/double3.h>
#include <rsys/math.h>
struct context {
@@ -23,8 +26,14 @@ struct context {
struct sdis_interface* interf;
};
+static INLINE double
+rand_canonic()
+{
+ return (double)rand()/(double)((double)RAND_MAX+1);
+}
+
static INLINE void
-get_indices(const size_t itri, size_t ids[3], void* context)
+get_indices_3d(const size_t itri, size_t ids[3], void* context)
{
struct context* ctx = context;
CHK(ctx != NULL);
@@ -35,7 +44,17 @@ get_indices(const size_t itri, size_t ids[3], void* context)
}
static INLINE void
-get_position(const size_t ivert, double pos[3], void* context)
+get_indices_2d(const size_t iseg, size_t ids[2], void* context)
+{
+ struct context* ctx = context;
+ CHK(ctx != NULL);
+ CHK(iseg < square_nsegments);
+ ids[0] = ctx->indices[iseg*2+0];
+ ids[1] = ctx->indices[iseg*2+1];
+}
+
+static INLINE void
+get_position_3d(const size_t ivert, double pos[3], void* context)
{
struct context* ctx = context;
CHK(ctx != NULL);
@@ -46,6 +65,16 @@ get_position(const size_t ivert, double pos[3], void* context)
}
static INLINE void
+get_position_2d(const size_t ivert, double pos[2], void* context)
+{
+ struct context* ctx = context;
+ CHK(ctx != NULL);
+ CHK(ivert < square_nvertices);
+ pos[0] = ctx->positions[ivert*2+0];
+ pos[1] = ctx->positions[ivert*2+1];
+}
+
+static INLINE void
get_interface(const size_t itri, struct sdis_interface** bound, void* context)
{
struct context* ctx = context;
@@ -55,30 +84,15 @@ get_interface(const size_t itri, struct sdis_interface** bound, void* context)
*bound = ctx->interf;
}
-int
-main(int argc, char** argv)
+static void
+test_scene_3d(struct sdis_device* dev, struct sdis_interface* interf)
{
- struct mem_allocator allocator;
- struct sdis_device* dev = NULL;
- struct sdis_medium* solid = NULL;
- struct sdis_medium* fluid = NULL;
- struct sdis_interface* interf = NULL;
struct sdis_scene* scn = NULL;
- struct sdis_fluid_shader fluid_shader = DUMMY_FLUID_SHADER;
- struct sdis_solid_shader solid_shader = DUMMY_SOLID_SHADER;
- struct sdis_interface_shader interface_shader = DUMMY_INTERFACE_SHADER;
double lower[3], upper[3];
+ double uv0[2], uv1[2], pos[3], pos1[3];
struct context ctx;
size_t ntris, npos;
- (void)argc, (void)argv;
-
- CHK(mem_init_proxy_allocator(&allocator, &mem_default_allocator) == RES_OK);
- CHK(sdis_device_create(NULL, &allocator, 1, 0, &dev) == RES_OK);
-
- CHK(sdis_fluid_create(dev, &fluid_shader, NULL, &fluid) == RES_OK);
- CHK(sdis_solid_create(dev, &solid_shader, NULL, &solid) == RES_OK);
- CHK(sdis_interface_create
- (dev, solid, fluid, &interface_shader, NULL, &interf) == RES_OK);
+ size_t i;
ctx.positions = box_vertices;
ctx.indices = box_indices;
@@ -87,8 +101,8 @@ main(int argc, char** argv)
npos = box_nvertices;
#define CREATE sdis_scene_create
- #define IDS get_indices
- #define POS get_position
+ #define IDS get_indices_3d
+ #define POS get_position_3d
#define IFA get_interface
CHK(CREATE(NULL, 0, NULL, NULL, 0, NULL, &ctx, NULL) == RES_BAD_ARG);
@@ -115,17 +129,160 @@ main(int argc, char** argv)
CHK(eq_eps(upper[1], 1, 1.e-6));
CHK(eq_eps(upper[2], 1, 1.e-6));
+ uv0[0] = 0.3;
+ uv0[1] = 0.3;
+
+ CHK(sdis_scene_get_boundary_position(NULL, 6, uv0, pos) == RES_BAD_ARG);
+ CHK(sdis_scene_get_boundary_position(scn, 12, uv0, pos) == RES_BAD_ARG);
+ CHK(sdis_scene_get_boundary_position(scn, 6, NULL, pos) == RES_BAD_ARG);
+ CHK(sdis_scene_get_boundary_position(scn, 6, uv0, NULL) == RES_BAD_ARG);
+ CHK(sdis_scene_get_boundary_position(scn, 6, uv0, pos) == RES_OK);
+
+ CHK(sdis_scene_boundary_project_position(NULL, 6, pos, uv1) == RES_BAD_ARG);
+ CHK(sdis_scene_boundary_project_position(scn, 12, pos, uv1) == RES_BAD_ARG);
+ CHK(sdis_scene_boundary_project_position(scn, 6, NULL, uv1) == RES_BAD_ARG);
+ CHK(sdis_scene_boundary_project_position(scn, 6, pos, NULL) == RES_BAD_ARG);
+ CHK(sdis_scene_boundary_project_position(scn, 6, pos, uv1) == RES_OK);
+
+ CHK(d2_eq_eps(uv0, uv1, 1.e-6));
+
+ FOR_EACH(i, 0, 64) {
+ uv0[0] = rand_canonic();
+ uv0[1] = rand_canonic() * (1 - uv0[0]);
+
+ CHK(sdis_scene_get_boundary_position(scn, 4, uv0, pos) == RES_OK);
+ CHK(sdis_scene_boundary_project_position(scn, 4, pos, uv1) == RES_OK);
+ CHK(d2_eq_eps(uv0, uv1, 1.e-6));
+ }
+
+ pos[0] = 10;
+ pos[1] = 0.1;
+ pos[2] = 0.5;
+ CHK(sdis_scene_boundary_project_position(scn, 6, pos, uv1) == RES_OK);
+ CHK(sdis_scene_get_boundary_position(scn, 6, uv1, pos1) == RES_OK);
+ CHK(!d3_eq_eps(pos1, pos, 1.e-6));
+
CHK(sdis_scene_ref_get(NULL) == RES_BAD_ARG);
CHK(sdis_scene_ref_get(scn) == RES_OK);
CHK(sdis_scene_ref_put(NULL) == RES_BAD_ARG);
CHK(sdis_scene_ref_put(scn) == RES_OK);
CHK(sdis_scene_ref_put(scn) == RES_OK);
+}
+
+static void
+test_scene_2d(struct sdis_device* dev, struct sdis_interface* interf)
+{
+ struct sdis_scene* scn = NULL;
+ double lower[2], upper[2];
+ double u0, u1, pos[2];
+ struct context ctx;
+ size_t nsegs, npos;
+ size_t i;
+
+ ctx.positions = square_vertices;
+ ctx.indices = square_indices;
+ ctx.interf = interf;
+ nsegs = square_nsegments;
+ npos = square_nvertices;
+
+ #define CREATE sdis_scene_2d_create
+ #define IDS get_indices_2d
+ #define POS get_position_2d
+ #define IFA get_interface
+
+ CHK(CREATE(NULL, 0, NULL, NULL, 0, NULL, &ctx, NULL) == RES_BAD_ARG);
+ CHK(CREATE(dev, 0, IDS, IFA, npos, POS, &ctx, &scn) == RES_BAD_ARG);
+ CHK(CREATE(dev, nsegs, NULL, IFA, npos, POS, &ctx, &scn) == RES_BAD_ARG);
+ CHK(CREATE(dev, nsegs, IDS, NULL, npos, POS, &ctx, &scn) == RES_BAD_ARG);
+ CHK(CREATE(dev, nsegs, IDS, IFA, 0, POS, &ctx, &scn) == RES_BAD_ARG);
+ CHK(CREATE(dev, nsegs, IDS, IFA, npos, NULL, &ctx, &scn) == RES_BAD_ARG);
+ CHK(CREATE(dev, nsegs, IDS, IFA, npos, POS, &ctx, &scn) == RES_OK);
+
+ #undef CREATE
+ #undef IDS
+ #undef POS
+ #undef IFA
+
+ CHK(sdis_scene_get_aabb(NULL, lower, upper) == RES_BAD_ARG);
+ CHK(sdis_scene_get_aabb(scn, NULL, upper) == RES_BAD_ARG);
+ CHK(sdis_scene_get_aabb(scn, lower, NULL) == RES_BAD_ARG);
+ CHK(sdis_scene_get_aabb(scn, lower, upper) == RES_OK);
+ CHK(eq_eps(lower[0], 0, 1.e-6));
+ CHK(eq_eps(lower[1], 0, 1.e-6));
+ CHK(eq_eps(upper[0], 1, 1.e-6));
+ CHK(eq_eps(upper[1], 1, 1.e-6));
+
+ u0 = 0.5;
+
+ CHK(sdis_scene_get_boundary_position(NULL, 1, &u0, pos) == RES_BAD_ARG);
+ CHK(sdis_scene_get_boundary_position(scn, 4, &u0, pos) == RES_BAD_ARG);
+ CHK(sdis_scene_get_boundary_position(scn, 1, NULL, pos) == RES_BAD_ARG);
+ CHK(sdis_scene_get_boundary_position(scn, 1, &u0, NULL) == RES_BAD_ARG);
+ CHK(sdis_scene_get_boundary_position(scn, 1, &u0, pos) == RES_OK);
+
+ CHK(sdis_scene_boundary_project_position(NULL, 1, pos, &u1) == RES_BAD_ARG);
+ CHK(sdis_scene_boundary_project_position(scn, 4, pos, &u1) == RES_BAD_ARG);
+ CHK(sdis_scene_boundary_project_position(scn, 1, NULL, &u1) == RES_BAD_ARG);
+ CHK(sdis_scene_boundary_project_position(scn, 1, pos, NULL) == RES_BAD_ARG);
+ CHK(sdis_scene_boundary_project_position(scn, 1, pos, &u1) == RES_OK);
+
+ CHK(eq_eps(u0, u1, 1.e-6));
+
+ FOR_EACH(i, 0, 64) {
+ u0 = rand_canonic();
+
+ CHK(sdis_scene_get_boundary_position(scn, 2, &u0, pos) == RES_OK);
+ CHK(sdis_scene_boundary_project_position(scn, 2, pos, &u1) == RES_OK);
+ CHK(eq_eps(u0, u1, 1.e-6));
+ }
+
+ d2(pos, 5, 0.5);
+ CHK(sdis_scene_boundary_project_position(scn, 3, pos, &u0) == RES_OK);
+ CHK(eq_eps(u0, 0.5, 1.e-6));
+
+ d2(pos, 1, 2);
+ CHK(sdis_scene_boundary_project_position(scn, 3, pos, &u0) == RES_OK);
+ CHK(eq_eps(u0, 0, 1.e-6));
+
+ d2(pos, 1, -1);
+ CHK(sdis_scene_boundary_project_position(scn, 3, pos, &u0) == RES_OK);
+ CHK(eq_eps(u0, 1, 1.e-6));
+
+ CHK(sdis_scene_ref_put(scn) == RES_OK);
+}
+
+int
+main(int argc, char** argv)
+{
+ struct mem_allocator allocator;
+ struct sdis_device* dev = NULL;
+ struct sdis_medium* solid = NULL;
+ struct sdis_medium* fluid = NULL;
+ struct sdis_interface* interf = NULL;
+ struct sdis_fluid_shader fluid_shader = DUMMY_FLUID_SHADER;
+ struct sdis_solid_shader solid_shader = DUMMY_SOLID_SHADER;
+ struct sdis_interface_shader interface_shader = SDIS_INTERFACE_SHADER_NULL;
+ (void)argc, (void)argv;
+
+ interface_shader.convection_coef = DUMMY_INTERFACE_SHADER.convection_coef;
+
+ CHK(mem_init_proxy_allocator(&allocator, &mem_default_allocator) == RES_OK);
+ CHK(sdis_device_create(NULL, &allocator, 1, 0, &dev) == RES_OK);
+
+ CHK(sdis_fluid_create(dev, &fluid_shader, NULL, &fluid) == RES_OK);
+ CHK(sdis_solid_create(dev, &solid_shader, NULL, &solid) == RES_OK);
+ CHK(sdis_interface_create
+ (dev, solid, fluid, &interface_shader, NULL, &interf) == RES_OK);
- CHK(sdis_device_ref_put(dev) == RES_OK);
- CHK(sdis_interface_ref_put(interf) == RES_OK);
CHK(sdis_medium_ref_put(solid) == RES_OK);
CHK(sdis_medium_ref_put(fluid) == RES_OK);
+ test_scene_3d(dev, interf);
+ test_scene_2d(dev, interf);
+
+ CHK(sdis_device_ref_put(dev) == RES_OK);
+ CHK(sdis_interface_ref_put(interf) == RES_OK);
+
check_memory_allocator(&allocator);
mem_shutdown_proxy_allocator(&allocator);
CHK(mem_allocated_size() == 0);
diff --git a/src/test_sdis_solve_camera.c b/src/test_sdis_solve_camera.c
@@ -30,7 +30,7 @@
#define UNKOWN_TEMPERATURE -1
#define IMG_WIDTH 640
#define IMG_HEIGHT 480
-#define SPP 4 /* #Samples per pixel, i.e. #realisations per pixel */
+#define SPP 4 /* #Samples per pixel, i.e. #realisations per pixel */
/*
* The scene is composed of a solid cube whose temperature is unknown. The
@@ -222,14 +222,6 @@ solid_get_delta
}
static double
-solid_get_delta_boundary
- (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
-{
- CHK(data != NULL && vtx != NULL);
- return ((const struct solid*)sdis_data_cget(data))->delta * 2.1;
-}
-
-static double
solid_get_temperature
(const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
{
@@ -308,7 +300,6 @@ create_solid
solid_shader.thermal_conductivity = solid_get_thermal_conductivity;
solid_shader.volumic_mass = solid_get_volumic_mass;
solid_shader.delta_solid = solid_get_delta;
- solid_shader.delta_boundary = solid_get_delta_boundary;
solid_shader.temperature = solid_get_temperature;
/* Create the solid medium */
@@ -361,7 +352,7 @@ create_interface
{
struct sdis_data* data = NULL;
struct interf* interface_param = NULL;
- struct sdis_interface_shader interface_shader = DUMMY_INTERFACE_SHADER;
+ struct sdis_interface_shader interface_shader = SDIS_INTERFACE_SHADER_NULL;
CHK(mdm_front != NULL);
CHK(mdm_back != NULL);
@@ -376,10 +367,16 @@ create_interface
/* Setup the interface shader */
interface_shader.convection_coef = interface_get_convection_coef;
- interface_shader.temperature = interface_get_temperature;
- interface_shader.emissivity = interface_get_emissivity;
- interface_shader.specular_fraction = interface_get_specular_fraction;
-
+ interface_shader.front.temperature = interface_get_temperature;
+ interface_shader.back.temperature = interface_get_temperature;
+ if(sdis_medium_get_type(mdm_front) == SDIS_FLUID) {
+ interface_shader.front.emissivity = interface_get_emissivity;
+ interface_shader.front.specular_fraction = interface_get_specular_fraction;
+ }
+ if(sdis_medium_get_type(mdm_back) == SDIS_FLUID) {
+ interface_shader.back.emissivity = interface_get_emissivity;
+ interface_shader.back.specular_fraction = interface_get_specular_fraction;
+ }
/* Create the interface */
CHK(sdis_interface_create
(dev, mdm_front, mdm_back, &interface_shader, data, interf) == RES_OK);
@@ -445,7 +442,7 @@ dump_image(const struct sdis_accum_buffer* buf)
double Tmax = -DBL_MAX;
double Tmin = DBL_MAX;
double norm;
- size_t ix, iy;
+ size_t i, ix, iy;
CHK(buf != NULL);
CHK(sdis_accum_buffer_get_layout(buf, &layout) == RES_OK);
@@ -453,8 +450,17 @@ dump_image(const struct sdis_accum_buffer* buf)
temps = mem_alloc(layout.width*layout.height*sizeof(double));
CHK(temps != NULL);
- /* Compute the per pixel temperature */
CHK(sdis_accum_buffer_map(buf, &accums) == RES_OK);
+
+ /* Check the results validity */
+ FOR_EACH(i, 0, layout.height * layout.width) {
+ CHK(accums[i].nweights + accums[i].nfailures == SPP);
+ CHK(accums[i].nfailures <= SPP/100);
+ CHK(accums[i].sum_weights >= 0);
+ CHK(accums[i].sum_weights_sqr >= 0);
+ }
+
+ /* Compute the per pixel temperature */
FOR_EACH(iy, 0, layout.height) {
const struct sdis_accum* row_accums = accums + iy * layout.width;
double* row = temps + iy * layout.width;
@@ -599,6 +605,11 @@ main(int argc, char** argv)
/* Create the accum buffer */
CHK(sdis_accum_buffer_create(dev, IMG_WIDTH, IMG_HEIGHT, &buf) == RES_OK);
+#if 0
+ dump_mesh(stdout, geom.positions, npos, geom.indices, ntris);
+ exit(0);
+#endif
+
/* Launch the simulation */
CHK(sdis_solve_camera(scn, cam, INF, 1, 300, 300, IMG_WIDTH, IMG_HEIGHT, SPP,
sdis_accum_buffer_write, buf) == RES_OK);
diff --git a/src/test_sdis_solve_probe.c b/src/test_sdis_solve_probe.c
@@ -126,14 +126,6 @@ solid_get_delta
}
static double
-solid_get_delta_boundary
- (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
-{
- CHK(data != NULL && vtx != NULL);
- return ((const struct solid*)sdis_data_cget(data))->delta * 2.1;
-}
-
-static double
solid_get_temperature
(const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
{
@@ -191,7 +183,7 @@ main(int argc, char** argv)
struct sdis_estimator* estimator = NULL;
struct sdis_fluid_shader fluid_shader = DUMMY_FLUID_SHADER;
struct sdis_solid_shader solid_shader = DUMMY_SOLID_SHADER;
- struct sdis_interface_shader interface_shader = DUMMY_INTERFACE_SHADER;
+ struct sdis_interface_shader interface_shader = SDIS_INTERFACE_SHADER_NULL;
struct context ctx;
struct fluid* fluid_param;
struct solid* solid_param;
@@ -206,7 +198,7 @@ main(int argc, char** argv)
CHK(mem_init_proxy_allocator(&allocator, &mem_default_allocator) == RES_OK);
CHK(sdis_device_create
- (NULL, &allocator, SDIS_NTHREADS_DEFAULT, 0, &dev) == RES_OK);
+ (NULL, &allocator, SDIS_NTHREADS_DEFAULT, 1, &dev) == RES_OK);
/* Create the fluid medium */
CHK(sdis_data_create
@@ -230,7 +222,6 @@ main(int argc, char** argv)
solid_shader.thermal_conductivity = solid_get_thermal_conductivity;
solid_shader.volumic_mass = solid_get_volumic_mass;
solid_shader.delta_solid = solid_get_delta;
- solid_shader.delta_boundary = solid_get_delta_boundary;
solid_shader.temperature = solid_get_temperature;
CHK(sdis_solid_create(dev, &solid_shader, data, &solid) == RES_OK);
CHK(sdis_data_ref_put(data) == RES_OK);
@@ -243,9 +234,10 @@ main(int argc, char** argv)
interface_param->epsilon = 0;
interface_param->specular_fraction = 0;
interface_shader.convection_coef = interface_get_convection_coef;
- interface_shader.temperature = NULL;
- interface_shader.emissivity = interface_get_emissivity;
- interface_shader.specular_fraction = interface_get_specular_fraction;
+ interface_shader.front = SDIS_INTERFACE_SIDE_SHADER_NULL;
+ interface_shader.back.temperature = NULL;
+ interface_shader.back.emissivity = interface_get_emissivity;
+ interface_shader.back.specular_fraction = interface_get_specular_fraction;
CHK(sdis_interface_create
(dev, solid, fluid, &interface_shader, data, &interf) == RES_OK);
CHK(sdis_data_ref_put(data) == RES_OK);
@@ -284,8 +276,6 @@ main(int argc, char** argv)
CHK(sdis_estimator_get_failure_count(NULL, &nfails) == RES_BAD_ARG);
CHK(sdis_estimator_get_failure_count(estimator, &nfails) == RES_OK);
- CHK(nfails + nreals == N);
-
CHK(sdis_estimator_get_temperature(estimator, NULL) == RES_BAD_ARG);
CHK(sdis_estimator_get_temperature(NULL, &T) == RES_BAD_ARG);
CHK(sdis_estimator_get_temperature(estimator, &T) == RES_OK);
@@ -293,6 +283,10 @@ main(int argc, char** argv)
ref = 300;
printf("Temperature at (%g, %g, %g) = %g ~ %g +/- %g\n",
SPLIT3(pos), ref, T.E, T.SE);
+ printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N);
+
+ CHK(nfails + nreals == N);
+ CHK(nfails < N/1000);
CHK(eq_eps(T.E, ref, T.SE));
CHK(sdis_estimator_ref_get(NULL) == RES_BAD_ARG);
diff --git a/src/test_sdis_solve_probe2.c b/src/test_sdis_solve_probe2.c
@@ -114,15 +114,6 @@ solid_get_delta
return 1.0/20.0;
}
-static double
-solid_get_delta_boundary
- (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
-{
- (void)data;
- CHK(vtx != NULL);
- return 2.1/20.0;
-}
-
/*******************************************************************************
* Interface
******************************************************************************/
@@ -180,7 +171,7 @@ main(int argc, char** argv)
CHK(mem_init_proxy_allocator(&allocator, &mem_default_allocator) == RES_OK);
CHK(sdis_device_create
- (NULL, &allocator, SDIS_NTHREADS_DEFAULT, 0, &dev) == RES_OK);
+ (NULL, &allocator, SDIS_NTHREADS_DEFAULT, 1, &dev) == RES_OK);
/* Create the fluid medium */
fluid_shader.temperature = temperature_unknown;
@@ -191,15 +182,13 @@ main(int argc, char** argv)
solid_shader.thermal_conductivity = solid_get_thermal_conductivity;
solid_shader.volumic_mass = solid_get_volumic_mass;
solid_shader.delta_solid = solid_get_delta;
- solid_shader.delta_boundary = solid_get_delta_boundary;
solid_shader.temperature = temperature_unknown;
CHK(sdis_solid_create(dev, &solid_shader, NULL, &solid) == RES_OK);
/* Create the fluid/solid interface with no limit conidition */
interface_shader.convection_coef = null_interface_value;
- interface_shader.temperature = NULL;
- interface_shader.emissivity = null_interface_value;
- interface_shader.specular_fraction = null_interface_value;
+ interface_shader.front = SDIS_INTERFACE_SIDE_SHADER_NULL;
+ interface_shader.back = SDIS_INTERFACE_SIDE_SHADER_NULL;
CHK(sdis_interface_create
(dev, solid, fluid, &interface_shader, NULL, &Tnone) == RES_OK);
@@ -208,10 +197,7 @@ main(int argc, char** argv)
ALIGNOF(struct interf), NULL, &data) == RES_OK);
interface_param = sdis_data_get(data);
interface_param->temperature = 300;
- interface_shader.convection_coef = null_interface_value;
- interface_shader.temperature = interface_get_temperature;
- interface_shader.emissivity = null_interface_value;
- interface_shader.specular_fraction = null_interface_value;
+ interface_shader.front.temperature = interface_get_temperature;
CHK(sdis_interface_create
(dev, solid, fluid, &interface_shader, data, &T300) == RES_OK);
CHK(sdis_data_ref_put(data) == RES_OK);
@@ -221,10 +207,6 @@ main(int argc, char** argv)
ALIGNOF(struct interf), NULL, &data) == RES_OK);
interface_param = sdis_data_get(data);
interface_param->temperature = 350;
- interface_shader.convection_coef = null_interface_value;
- interface_shader.temperature = interface_get_temperature;
- interface_shader.emissivity = null_interface_value;
- interface_shader.specular_fraction = null_interface_value;
CHK(sdis_interface_create
(dev, solid, fluid, &interface_shader, data, &T350) == RES_OK);
CHK(sdis_data_ref_put(data) == RES_OK);
@@ -268,12 +250,12 @@ main(int argc, char** argv)
ref = 350 * pos[2] + (1-pos[2]) * 300;
printf("Temperature at (%g, %g, %g) = %g ~ %g +/- %g\n",
SPLIT3(pos), ref, T.E, T.SE);
- printf("#realisations: %lu; #failures: %lu\n",
- (unsigned long)nreals, (unsigned long)nfails);
+ printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N);
/* Check the results */
CHK(nfails + nreals == N);
- CHK(eq_eps(T.E, ref, T.SE));
+ CHK(nfails < N/1000);
+ CHK(eq_eps(T.E, ref, 3*T.SE));
/* Release data */
CHK(sdis_estimator_ref_put(estimator) == RES_OK);
diff --git a/src/test_sdis_solve_probe2_2d.c b/src/test_sdis_solve_probe2_2d.c
@@ -111,15 +111,6 @@ solid_get_delta
return 1.0/20.0;
}
-static double
-solid_get_delta_boundary
- (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
-{
- (void)data;
- CHK(vtx != NULL);
- return 2.1/20.0;
-}
-
/*******************************************************************************
* Interface
******************************************************************************/
@@ -177,7 +168,7 @@ main(int argc, char** argv)
CHK(mem_init_proxy_allocator(&allocator, &mem_default_allocator) == RES_OK);
CHK(sdis_device_create
- (NULL, &allocator, SDIS_NTHREADS_DEFAULT, 0, &dev) == RES_OK);
+ (NULL, &allocator, SDIS_NTHREADS_DEFAULT, 1, &dev) == RES_OK);
/* Create the fluid medium */
fluid_shader.temperature = temperature_unknown;
@@ -188,27 +179,26 @@ main(int argc, char** argv)
solid_shader.thermal_conductivity = solid_get_thermal_conductivity;
solid_shader.volumic_mass = solid_get_volumic_mass;
solid_shader.delta_solid = solid_get_delta;
- solid_shader.delta_boundary = solid_get_delta_boundary;
solid_shader.temperature = temperature_unknown;
CHK(sdis_solid_create(dev, &solid_shader, NULL, &solid) == RES_OK);
/* Create the fluid/solid interface with no limit conidition */
interface_shader.convection_coef = null_interface_value;
- interface_shader.temperature = NULL;
- interface_shader.emissivity = null_interface_value;
- interface_shader.specular_fraction = null_interface_value;
+ interface_shader.front = SDIS_INTERFACE_SIDE_SHADER_NULL;
+ interface_shader.back = SDIS_INTERFACE_SIDE_SHADER_NULL;
CHK(sdis_interface_create
(dev, solid, fluid, &interface_shader, NULL, &Tnone) == RES_OK);
+ interface_shader.convection_coef = null_interface_value;
+ interface_shader.front = SDIS_INTERFACE_SIDE_SHADER_NULL;
+ interface_shader.back = SDIS_INTERFACE_SIDE_SHADER_NULL;
+ interface_shader.front.temperature = interface_get_temperature;
+
/* Create the fluid/solid interface with a fixed temperature of 300K */
CHK(sdis_data_create(dev, sizeof(struct interf),
ALIGNOF(struct interf), NULL, &data) == RES_OK);
interface_param = sdis_data_get(data);
interface_param->temperature = 300;
- interface_shader.convection_coef = null_interface_value;
- interface_shader.temperature = interface_get_temperature;
- interface_shader.emissivity = null_interface_value;
- interface_shader.specular_fraction = null_interface_value;
CHK(sdis_interface_create
(dev, solid, fluid, &interface_shader, data, &T300) == RES_OK);
CHK(sdis_data_ref_put(data) == RES_OK);
@@ -218,10 +208,6 @@ main(int argc, char** argv)
ALIGNOF(struct interf), NULL, &data) == RES_OK);
interface_param = sdis_data_get(data);
interface_param->temperature = 350;
- interface_shader.convection_coef = null_interface_value;
- interface_shader.temperature = interface_get_temperature;
- interface_shader.emissivity = null_interface_value;
- interface_shader.specular_fraction = null_interface_value;
CHK(sdis_interface_create
(dev, solid, fluid, &interface_shader, data, &T350) == RES_OK);
CHK(sdis_data_ref_put(data) == RES_OK);
@@ -262,12 +248,12 @@ main(int argc, char** argv)
ref = 350 * pos[0] + (1-pos[0]) * 300;
printf("Temperature at (%g, %g) = %g ~ %g +/- %g\n",
SPLIT2(pos), ref, T.E, T.SE);
- printf("#realisations: %lu; #failures: %lu\n",
- (unsigned long)nreals, (unsigned long)nfails);
+ printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N);
/* Check the results */
CHK(nfails + nreals == N);
- CHK(eq_eps(T.E, ref, T.SE));
+ CHK(nfails < N/1000);
+ CHK(eq_eps(T.E, ref, T.SE*2));
/* Release data */
CHK(sdis_estimator_ref_put(estimator) == RES_OK);
diff --git a/src/test_sdis_solve_probe3.c b/src/test_sdis_solve_probe3.c
@@ -136,15 +136,6 @@ solid_get_delta
return 1.0/20.0;
}
-static double
-solid_get_delta_boundary
- (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
-{
- (void)data;
- CHK(vtx != NULL);
- return 2.1/20.0;
-}
-
/*******************************************************************************
* Interface
******************************************************************************/
@@ -207,7 +198,7 @@ main(int argc, char** argv)
CHK(mem_init_proxy_allocator(&allocator, &mem_default_allocator) == RES_OK);
CHK(sdis_device_create
- (NULL, &allocator, SDIS_NTHREADS_DEFAULT, 0, &dev) == RES_OK);
+ (NULL, &allocator, SDIS_NTHREADS_DEFAULT, 1, &dev) == RES_OK);
/* Create the fluid medium */
fluid_shader.temperature = temperature_unknown;
@@ -218,15 +209,13 @@ main(int argc, char** argv)
solid_shader.thermal_conductivity = solid_get_thermal_conductivity;
solid_shader.volumic_mass = solid_get_volumic_mass;
solid_shader.delta_solid = solid_get_delta;
- solid_shader.delta_boundary = solid_get_delta_boundary;
solid_shader.temperature = temperature_unknown;
CHK(sdis_solid_create(dev, &solid_shader, NULL, &solid) == RES_OK);
/* Create the fluid/solid interface with no limit conidition */
interface_shader.convection_coef = null_interface_value;
- interface_shader.temperature = NULL;
- interface_shader.emissivity = null_interface_value;
- interface_shader.specular_fraction = null_interface_value;
+ interface_shader.front = SDIS_INTERFACE_SIDE_SHADER_NULL;
+ interface_shader.back = SDIS_INTERFACE_SIDE_SHADER_NULL;
CHK(sdis_interface_create
(dev, solid, fluid, &interface_shader, NULL, &Tnone) == RES_OK);
@@ -235,10 +224,7 @@ main(int argc, char** argv)
ALIGNOF(struct interf), NULL, &data) == RES_OK);
interface_param = sdis_data_get(data);
interface_param->temperature = 300;
- interface_shader.convection_coef = null_interface_value;
- interface_shader.temperature = interface_get_temperature;
- interface_shader.emissivity = null_interface_value;
- interface_shader.specular_fraction = null_interface_value;
+ interface_shader.front.temperature = interface_get_temperature;
CHK(sdis_interface_create
(dev, solid, fluid, &interface_shader, data, &T300) == RES_OK);
CHK(sdis_data_ref_put(data) == RES_OK);
@@ -248,19 +234,12 @@ main(int argc, char** argv)
ALIGNOF(struct interf), NULL, &data) == RES_OK);
interface_param = sdis_data_get(data);
interface_param->temperature = 350;
- interface_shader.convection_coef = null_interface_value;
- interface_shader.temperature = interface_get_temperature;
- interface_shader.emissivity = null_interface_value;
- interface_shader.specular_fraction = null_interface_value;
CHK(sdis_interface_create
(dev, solid, fluid, &interface_shader, data, &T350) == RES_OK);
CHK(sdis_data_ref_put(data) == RES_OK);
/* Create the solid/solid interface */
- interface_shader.convection_coef = NULL;
- interface_shader.temperature = NULL;
- interface_shader.specular_fraction = NULL;
- interface_shader.emissivity = NULL;
+ interface_shader = SDIS_INTERFACE_SHADER_NULL;
CHK(sdis_interface_create
(dev, solid, solid, &interface_shader, NULL, &solid_solid) == RES_OK);
@@ -327,11 +306,11 @@ main(int argc, char** argv)
ref = 350 * pos[2] + (1-pos[2]) * 300;
printf("Temperature at (%g, %g, %g) = %g ~ %g +/- %g\n",
SPLIT3(pos), ref, T.E, T.SE);
- printf("#realisations: %lu; #failures: %lu\n",
- (unsigned long)nreals, (unsigned long)nfails);
+ printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N);
/* Check the results */
CHK(nfails + nreals == N);
+ CHK(nfails < N/1000);
CHK(eq_eps(T.E, ref, 2*T.SE));
/* Release data */
diff --git a/src/test_sdis_solve_probe3_2d.c b/src/test_sdis_solve_probe3_2d.c
@@ -133,15 +133,6 @@ solid_get_delta
return 1.0/20.0;
}
-static double
-solid_get_delta_boundary
- (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
-{
- (void)data;
- CHK(vtx != NULL);
- return 2.1/20.0;
-}
-
/*******************************************************************************
* Interface
******************************************************************************/
@@ -202,7 +193,7 @@ main(int argc, char** argv)
CHK(mem_init_proxy_allocator(&allocator, &mem_default_allocator) == RES_OK);
CHK(sdis_device_create
- (NULL, &allocator, SDIS_NTHREADS_DEFAULT, 0, &dev) == RES_OK);
+ (NULL, &allocator, SDIS_NTHREADS_DEFAULT, 1, &dev) == RES_OK);
/* Create the fluid medium */
fluid_shader.temperature = temperature_unknown;
@@ -213,15 +204,11 @@ main(int argc, char** argv)
solid_shader.thermal_conductivity = solid_get_thermal_conductivity;
solid_shader.volumic_mass = solid_get_volumic_mass;
solid_shader.delta_solid = solid_get_delta;
- solid_shader.delta_boundary = solid_get_delta_boundary;
solid_shader.temperature = temperature_unknown;
CHK(sdis_solid_create(dev, &solid_shader, NULL, &solid) == RES_OK);
/* Create the fluid/solid interface with no limit conidition */
- interface_shader.convection_coef = null_interface_value;
- interface_shader.temperature = NULL;
- interface_shader.emissivity = null_interface_value;
- interface_shader.specular_fraction = null_interface_value;
+ interface_shader = SDIS_INTERFACE_SHADER_NULL;
CHK(sdis_interface_create
(dev, solid, fluid, &interface_shader, NULL, &Tnone) == RES_OK);
@@ -231,9 +218,8 @@ main(int argc, char** argv)
interface_param = sdis_data_get(data);
interface_param->temperature = 300;
interface_shader.convection_coef = null_interface_value;
- interface_shader.temperature = interface_get_temperature;
- interface_shader.emissivity = null_interface_value;
- interface_shader.specular_fraction = null_interface_value;
+ interface_shader.front.temperature = interface_get_temperature;
+ interface_shader.back = SDIS_INTERFACE_SIDE_SHADER_NULL;
CHK(sdis_interface_create
(dev, solid, fluid, &interface_shader, data, &T300) == RES_OK);
CHK(sdis_data_ref_put(data) == RES_OK);
@@ -244,18 +230,14 @@ main(int argc, char** argv)
interface_param = sdis_data_get(data);
interface_param->temperature = 350;
interface_shader.convection_coef = null_interface_value;
- interface_shader.temperature = interface_get_temperature;
- interface_shader.emissivity = null_interface_value;
- interface_shader.specular_fraction = null_interface_value;
+ interface_shader.front.temperature = interface_get_temperature;
+ interface_shader.back = SDIS_INTERFACE_SIDE_SHADER_NULL;
CHK(sdis_interface_create
(dev, solid, fluid, &interface_shader, data, &T350) == RES_OK);
CHK(sdis_data_ref_put(data) == RES_OK);
/* Create the solid/solid interface */
- interface_shader.convection_coef = NULL;
- interface_shader.temperature = NULL;
- interface_shader.specular_fraction = NULL;
- interface_shader.emissivity = NULL;
+ interface_shader = SDIS_INTERFACE_SHADER_NULL;
CHK(sdis_interface_create
(dev, solid, solid, &interface_shader, NULL, &solid_solid) == RES_OK);
@@ -317,11 +299,11 @@ main(int argc, char** argv)
ref = 350 * pos[0] + (1-pos[0]) * 300;
printf("Temperature at (%g, %g) = %g ~ %g +/- %g\n",
SPLIT2(pos), ref, T.E, T.SE);
- printf("#realisations: %lu; #failures: %lu\n",
- (unsigned long)nreals, (unsigned long)nfails);
+ printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N);
/* Check the results */
CHK(nfails + nreals == N);
+ CHK(nfails < N/1000);
CHK(eq_eps(T.E, ref, 3*T.SE));
/* Release data */
diff --git a/src/test_sdis_solve_probe_2d.c b/src/test_sdis_solve_probe_2d.c
@@ -108,14 +108,6 @@ solid_get_delta
}
static double
-solid_get_delta_boundary
- (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
-{
- (void)vtx, (void)data;
- return 2.1/20.0;
-}
-
-static double
solid_get_temperature
(const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
{
@@ -131,14 +123,6 @@ interface_get_convection_coef
return 0.5;
}
-static double
-interface_null_reflectivity
- (const struct sdis_interface_fragment* frag, struct sdis_data* data)
-{
- (void)frag, (void)data;
- return 0;
-}
-
/*******************************************************************************
* Main test
******************************************************************************/
@@ -167,7 +151,7 @@ main(int argc, char** argv)
CHK(mem_init_proxy_allocator(&allocator, &mem_default_allocator) == RES_OK);
CHK(sdis_device_create
- (NULL, &allocator, SDIS_NTHREADS_DEFAULT, 0, &dev) == RES_OK);
+ (NULL, &allocator, SDIS_NTHREADS_DEFAULT, 1, &dev) == RES_OK);
/* Create the fluid medium */
fluid_shader.temperature = fluid_get_temperature;
@@ -178,15 +162,13 @@ main(int argc, char** argv)
solid_shader.thermal_conductivity = solid_get_thermal_conductivity;
solid_shader.volumic_mass = solid_get_volumic_mass;
solid_shader.delta_solid = solid_get_delta;
- solid_shader.delta_boundary = solid_get_delta_boundary;
solid_shader.temperature = solid_get_temperature;
CHK(sdis_solid_create(dev, &solid_shader, NULL, &solid) == RES_OK);
/* Create the solid/fluid interface */
interface_shader.convection_coef = interface_get_convection_coef;
- interface_shader.temperature = NULL;
- interface_shader.emissivity = interface_null_reflectivity;
- interface_shader.specular_fraction = interface_null_reflectivity;
+ interface_shader.front = SDIS_INTERFACE_SIDE_SHADER_NULL;
+ interface_shader.back = SDIS_INTERFACE_SIDE_SHADER_NULL;
CHK(sdis_interface_create
(dev, solid, fluid, &interface_shader, NULL, &interf) == RES_OK);
@@ -211,13 +193,15 @@ main(int argc, char** argv)
CHK(sdis_estimator_get_realisation_count(estimator, &nreals) == RES_OK);
CHK(sdis_estimator_get_failure_count(estimator, &nfails) == RES_OK);
- CHK(nfails + nreals == N);
-
CHK(sdis_estimator_get_temperature(estimator, &T) == RES_OK);
ref = 300;
printf("Temperature at (%g, %g) = %g ~ %g +/- %g\n",
SPLIT2(pos), ref, T.E, T.SE);
+ printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N);
+
+ CHK(nfails + nreals == N);
+ CHK(nfails < N/1000);
CHK(eq_eps(T.E, ref, T.SE));
CHK(sdis_estimator_ref_put(estimator) == RES_OK);
diff --git a/src/test_sdis_solve_probe_boundary.c b/src/test_sdis_solve_probe_boundary.c
@@ -99,15 +99,6 @@ solid_get_delta
}
static double
-solid_get_delta_boundary
- (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
-{
- (void)data;
- CHK(vtx != NULL);
- return 2.1/20.0;
-}
-
-static double
solid_get_temperature
(const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
{
@@ -142,22 +133,6 @@ interface_get_convection_coef
return interf->hc;
}
-static double
-interface_get_emissivity
- (const struct sdis_interface_fragment* frag, struct sdis_data* data)
-{
- CHK(frag && data);
- return 0;
-}
-
-static double
-interface_get_specular_fraction
- (const struct sdis_interface_fragment* frag, struct sdis_data* data)
-{
- CHK(frag && data);
- return 0;
-}
-
/*******************************************************************************
* Test
******************************************************************************/
@@ -178,7 +153,7 @@ main(int argc, char** argv)
struct sdis_estimator* estimator = NULL;
struct sdis_fluid_shader fluid_shader = DUMMY_FLUID_SHADER;
struct sdis_solid_shader solid_shader = DUMMY_SOLID_SHADER;
- struct sdis_interface_shader interf_shader = DUMMY_INTERFACE_SHADER;
+ struct sdis_interface_shader interf_shader = SDIS_INTERFACE_SHADER_NULL;
struct sdis_interface* box_interfaces[12 /*#triangles*/];
struct sdis_interface* square_interfaces[4/*#segments*/];
struct interf* interf_props = NULL;
@@ -192,7 +167,7 @@ main(int argc, char** argv)
CHK(mem_init_proxy_allocator(&allocator, &mem_default_allocator) == RES_OK);
CHK(sdis_device_create
- (NULL, &allocator, SDIS_NTHREADS_DEFAULT, 0, &dev) == RES_OK);
+ (NULL, &allocator, SDIS_NTHREADS_DEFAULT, 1, &dev) == RES_OK);
/* Create the fluid medium */
fluid_shader.temperature = fluid_get_temperature;
@@ -203,15 +178,15 @@ main(int argc, char** argv)
solid_shader.thermal_conductivity = solid_get_thermal_conductivity;
solid_shader.volumic_mass = solid_get_volumic_mass;
solid_shader.delta_solid = solid_get_delta;
- solid_shader.delta_boundary = solid_get_delta_boundary;
solid_shader.temperature = solid_get_temperature;
CHK(sdis_solid_create(dev, &solid_shader, NULL, &solid) == RES_OK);
/* Setup the interface shader */
- interf_shader.temperature = interface_get_temperature;
interf_shader.convection_coef = interface_get_convection_coef;
- interf_shader.emissivity = interface_get_emissivity;
- interf_shader.specular_fraction = interface_get_specular_fraction;
+ interf_shader.front.temperature = interface_get_temperature;
+ interf_shader.front.emissivity = NULL;
+ interf_shader.front.specular_fraction = NULL;
+ interf_shader.back = SDIS_INTERFACE_SIDE_SHADER_NULL;
/* Create the adiabatic interface */
CHK(sdis_data_create(dev, sizeof(struct interf), 16, NULL, &data) == RES_OK);
@@ -278,50 +253,44 @@ main(int argc, char** argv)
iprim = 6;
#define SOLVE sdis_solve_probe_boundary
- CHK(SOLVE(NULL, N, iprim, uv, INF, 1.0, 0, 0, &estimator) == RES_BAD_ARG);
- CHK(SOLVE(box_scn, 0, iprim, uv, INF, 1.0, 0, 0, &estimator) == RES_BAD_ARG);
- CHK(SOLVE(box_scn, N, 12, uv, INF, 1.0, 0, 0, &estimator) == RES_BAD_ARG);
- CHK(SOLVE(box_scn, N, iprim, NULL, INF, 1.0, 0, 0, &estimator) == RES_BAD_ARG);
- CHK(SOLVE(box_scn, N, iprim, uv, -1, 1.0, 0, 0, &estimator) == RES_BAD_ARG);
- CHK(SOLVE(box_scn, N, iprim, uv, INF, 1.0, 0, 0, NULL) == RES_BAD_ARG);
- CHK(SOLVE(box_scn, N, iprim, uv, INF, 1.0, 0, 0, &estimator) == RES_OK);
+ #define F SDIS_FRONT
+ CHK(SOLVE(NULL, N, iprim, uv, INF, F, 1.0, 0, 0, &estimator) == RES_BAD_ARG);
+ CHK(SOLVE(box_scn, 0, iprim, uv, INF, F, 1.0, 0, 0, &estimator) == RES_BAD_ARG);
+ CHK(SOLVE(box_scn, N, 12, uv, INF, F, 1.0, 0, 0, &estimator) == RES_BAD_ARG);
+ CHK(SOLVE(box_scn, N, iprim, NULL, INF, F, 1.0, 0, 0, &estimator) == RES_BAD_ARG);
+ CHK(SOLVE(box_scn, N, iprim, uv, -1, F, 1.0, 0, 0, &estimator) == RES_BAD_ARG);
+ CHK(SOLVE(box_scn, N, iprim, uv, INF, -1, 1.0, 0, 0, &estimator) == RES_BAD_ARG);
+ CHK(SOLVE(box_scn, N, iprim, uv, INF, F, 1.0, 0, 0, NULL) == RES_BAD_ARG);
+ CHK(SOLVE(box_scn, N, iprim, uv, INF, F, 1.0, 0, 0, &estimator) == RES_OK);
+
+ ref = (H*Tf + LAMBDA * Tb) / (H + LAMBDA);
CHK(sdis_estimator_get_realisation_count(estimator, &nreals) == RES_OK);
CHK(sdis_estimator_get_failure_count(estimator, &nfails) == RES_OK);
- CHK(nfails + nreals == N);
-
CHK(sdis_estimator_get_temperature(estimator, &T) == RES_OK);
CHK(sdis_estimator_ref_put(estimator) == RES_OK);
-
- CHK(sdis_scene_get_boundary_position(NULL, iprim, uv, pos) == RES_BAD_ARG);
- CHK(sdis_scene_get_boundary_position(box_scn, 12, uv, pos) == RES_BAD_ARG);
- CHK(sdis_scene_get_boundary_position(box_scn, iprim, NULL, pos) == RES_BAD_ARG);
- CHK(sdis_scene_get_boundary_position(box_scn, iprim, uv, NULL) == RES_BAD_ARG);
CHK(sdis_scene_get_boundary_position(box_scn, iprim, uv, pos) == RES_OK);
-
- ref = (H*Tf + LAMBDA * Tb) / (H + LAMBDA);
-
printf("Boundary temperature of the box at (%g %g %g) = %g ~ %g +/- %g\n",
SPLIT3(pos), ref, T.E, T.SE);
printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N);
- CHK(eq_eps(T.E, ref, T.SE));
+ CHK(nfails + nreals == N);
+ CHK(nfails < N/1000);
+ CHK(eq_eps(T.E, ref, 3*T.SE));
uv[0] = 0.5;
iprim = 3;
- CHK(SOLVE(square_scn, N, iprim, uv, INF, 1.0, 0, 0, &estimator) == RES_OK);
+ CHK(SOLVE(square_scn, N, iprim, uv, INF, F, 1.0, 0, 0, &estimator) == RES_OK);
CHK(sdis_estimator_get_realisation_count(estimator, &nreals) == RES_OK);
CHK(sdis_estimator_get_failure_count(estimator, &nfails) == RES_OK);
- CHK(nfails + nreals == N);
-
CHK(sdis_estimator_get_temperature(estimator, &T) == RES_OK);
CHK(sdis_estimator_ref_put(estimator) == RES_OK);
-
CHK(sdis_scene_get_boundary_position(square_scn, iprim, uv, pos) == RES_OK);
-
printf("Boundary temperature of the square at (%g %g) = %g ~ %g +/- %g\n",
SPLIT2(pos), ref, T.E, T.SE);
printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N);
- CHK(eq_eps(T.E, ref, T.SE));
+ CHK(nfails + nreals == N);
+ CHK(nfails < N/1000);
+ CHK(eq_eps(T.E, ref, 3*T.SE));
#undef SOLVE
CHK(sdis_scene_ref_put(box_scn) == RES_OK);
diff --git a/src/test_sdis_utils.h b/src/test_sdis_utils.h
@@ -100,7 +100,7 @@ static const double square_vertices[4/*#vertices*/*2/*#coords per vertex*/] = {
};
static const size_t square_nvertices = sizeof(square_vertices)/sizeof(double[2]);
-static const size_t square_indices[4/*#triangles*/*2/*#indices per segment*/]= {
+static const size_t square_indices[4/*#segments*/*2/*#indices per segment*/]= {
0, 1, /* Bottom */
1, 2, /* Left */
2, 3, /* Top */
@@ -165,7 +165,6 @@ static const struct sdis_solid_shader DUMMY_SOLID_SHADER = {
dummy_medium_getter,
dummy_medium_getter,
dummy_medium_getter,
- dummy_medium_getter,
dummy_medium_getter
};
@@ -175,11 +174,17 @@ static const struct sdis_fluid_shader DUMMY_FLUID_SHADER = {
dummy_medium_getter
};
+
+#define DUMMY_INTERFACE_SIDE_SHADER__ { \
+ dummy_interface_getter, \
+ dummy_interface_getter, \
+ dummy_interface_getter, \
+ dummy_interface_getter \
+}
static const struct sdis_interface_shader DUMMY_INTERFACE_SHADER = {
dummy_interface_getter,
- dummy_interface_getter,
- dummy_interface_getter,
- dummy_interface_getter
+ DUMMY_INTERFACE_SIDE_SHADER__,
+ DUMMY_INTERFACE_SIDE_SHADER__
};
/*******************************************************************************
diff --git a/src/test_sdis_volumic_power.c b/src/test_sdis_volumic_power.c
@@ -28,7 +28,7 @@
*
* T(pos) = P0 / (2*LAMBDA) * (A^2/4 - (pos-0.5)^2) + T0
*
- * with LAMBDA the conductivity of the cube and A the size of the cube/square,
+ * with LAMBDA the conductivity of the solid and A the size of the cube/square,
* i.e. 1.
*
* 3D 2D
@@ -49,10 +49,18 @@
#define T0 320
#define LAMBDA 0.1
#define P0 10
+#define DELTA 1.0/20.0
/*******************************************************************************
* Media
******************************************************************************/
+struct solid {
+ double lambda;
+ double rho;
+ double cp;
+ double delta;
+};
+
static double
fluid_get_temperature
(const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
@@ -66,45 +74,32 @@ static double
solid_get_calorific_capacity
(const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
{
- (void)data;
CHK(vtx != NULL);
- return 2.0;
+ return ((struct solid*)sdis_data_cget(data))->cp;
}
static double
solid_get_thermal_conductivity
(const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
{
- (void)data;
CHK(vtx != NULL);
- return LAMBDA;
+ return ((struct solid*)sdis_data_cget(data))->lambda;
}
static double
solid_get_volumic_mass
(const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
{
- (void)data;
CHK(vtx != NULL);
- return 25.0;
+ return ((struct solid*)sdis_data_cget(data))->rho;
}
static double
solid_get_delta
(const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
{
- (void)data;
CHK(vtx != NULL);
- return 1.0/20.0;
-}
-
-static double
-solid_get_delta_boundary
- (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
-{
- (void)data;
- CHK(vtx != NULL);
- return 2.1/20.0;
+ return ((struct solid*)sdis_data_cget(data))->delta;
}
static double
@@ -149,22 +144,6 @@ interface_get_convection_coef
return 0;
}
-static double
-interface_get_emissivity
- (const struct sdis_interface_fragment* frag, struct sdis_data* data)
-{
- CHK(frag && data);
- return 0;
-}
-
-static double
-interface_get_specular_fraction
- (const struct sdis_interface_fragment* frag, struct sdis_data* data)
-{
- CHK(frag && data);
- return 0;
-}
-
/*******************************************************************************
* Test
******************************************************************************/
@@ -177,6 +156,7 @@ main(int argc, char** argv)
struct sdis_device* dev = NULL;
struct sdis_medium* fluid = NULL;
struct sdis_medium* solid = NULL;
+ struct sdis_medium* solid2 = NULL; /* For debug */
struct sdis_interface* interf_adiabatic = NULL;
struct sdis_interface* interf_T0 = NULL;
struct sdis_scene* box_scn = NULL;
@@ -184,10 +164,11 @@ main(int argc, char** argv)
struct sdis_estimator* estimator = NULL;
struct sdis_fluid_shader fluid_shader = DUMMY_FLUID_SHADER;
struct sdis_solid_shader solid_shader = DUMMY_SOLID_SHADER;
- struct sdis_interface_shader interf_shader = DUMMY_INTERFACE_SHADER;
+ struct sdis_interface_shader interf_shader = SDIS_INTERFACE_SHADER_NULL;
struct sdis_interface* box_interfaces[12 /*#triangles*/];
struct sdis_interface* square_interfaces[4/*#segments*/];
struct interf* interf_props = NULL;
+ struct solid* solid_props = NULL;
double pos[3];
double x;
double ref;
@@ -197,27 +178,41 @@ main(int argc, char** argv)
CHK(mem_init_proxy_allocator(&allocator, &mem_default_allocator) == RES_OK);
CHK(sdis_device_create
- (NULL, &allocator, SDIS_NTHREADS_DEFAULT, 0, &dev) == RES_OK);
+ (NULL, &allocator, SDIS_NTHREADS_DEFAULT, 1, &dev) == RES_OK);
- /* Create the fluid medium */
fluid_shader.temperature = fluid_get_temperature;
CHK(sdis_fluid_create(dev, &fluid_shader, NULL, &fluid) == RES_OK);
- /* Create the solid_medium */
+ /* Setup the solid shader */
solid_shader.calorific_capacity = solid_get_calorific_capacity;
solid_shader.thermal_conductivity = solid_get_thermal_conductivity;
solid_shader.volumic_mass = solid_get_volumic_mass;
solid_shader.delta_solid = solid_get_delta;
- solid_shader.delta_boundary = solid_get_delta_boundary;
solid_shader.temperature = solid_get_temperature;
solid_shader.volumic_power = solid_get_volumic_power;
- CHK(sdis_solid_create(dev, &solid_shader, NULL, &solid) == RES_OK);
+
+ /* Create the solid medium */
+ CHK(sdis_data_create(dev, sizeof(struct solid), 16, NULL, &data) == RES_OK);
+ solid_props = sdis_data_get(data);
+ solid_props->lambda = LAMBDA;
+ solid_props->cp = 2;
+ solid_props->rho = 25;
+ solid_props->delta = DELTA;
+ CHK(sdis_solid_create(dev, &solid_shader, data, &solid) == RES_OK);
+ CHK(sdis_data_ref_put(data) == RES_OK);
+
+ CHK(sdis_data_create(dev, sizeof(struct solid), 16, NULL, &data) == RES_OK);
+ solid_props = sdis_data_get(data);
+ solid_props->lambda = 0;
+ solid_props->cp = 0;
+ solid_props->rho = 0;
+ solid_props->delta = DELTA/4;
+ CHK(sdis_solid_create(dev, &solid_shader, data, &solid2) == RES_OK);
+ CHK(sdis_data_ref_put(data) == RES_OK);
/* Setup the interface shader */
- interf_shader.temperature = interface_get_temperature;
interf_shader.convection_coef = interface_get_convection_coef;
- interf_shader.emissivity = interface_get_emissivity;
- interf_shader.specular_fraction = interface_get_specular_fraction;
+ interf_shader.front.temperature = interface_get_temperature;
/* Create the adiabatic interface */
CHK(sdis_data_create(dev, sizeof(struct interf), 16, NULL, &data) == RES_OK);
@@ -237,6 +232,7 @@ main(int argc, char** argv)
/* Release the media */
CHK(sdis_medium_ref_put(solid) == RES_OK);
+ CHK(sdis_medium_ref_put(solid2) == RES_OK);
CHK(sdis_medium_ref_put(fluid) == RES_OK);
/* Map the interfaces to their box triangles */
@@ -249,7 +245,7 @@ main(int argc, char** argv)
/* Map the interfaces to their square segments */
square_interfaces[0] = interf_adiabatic; /* Bottom */
- square_interfaces[1] = interf_T0; /* Lef */
+ square_interfaces[1] = interf_T0; /* Left */
square_interfaces[2] = interf_adiabatic; /* Top */
square_interfaces[3] = interf_T0; /* Right */
@@ -281,19 +277,22 @@ main(int argc, char** argv)
printf("Temperature of the box at (%g %g %g) = %g ~ %g +/- %g\n",
SPLIT3(pos), ref, T.E, T.SE);
printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N);
- CHK(eq_eps(T.E, ref, T.SE*2));
+ CHK(nfails + nreals == N);
+ CHK(nfails < N/1000);
+ CHK(eq_eps(T.E, ref, 3*T.SE));
/* Solve in 2D */
CHK(sdis_solve_probe(square_scn, N, pos, INF, 1.0, 0, 0, &estimator) == RES_OK);
CHK(sdis_estimator_get_realisation_count(estimator, &nreals) == RES_OK);
CHK(sdis_estimator_get_failure_count(estimator, &nfails) == RES_OK);
- CHK(nfails + nreals == N);
CHK(sdis_estimator_get_temperature(estimator, &T) == RES_OK);
CHK(sdis_estimator_ref_put(estimator) == RES_OK);
printf("Temperature of the square at (%g %g) = %g ~ %g +/- %g\n",
SPLIT2(pos), ref, T.E, T.SE);
printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N);
- CHK(eq_eps(T.E, ref, T.SE*2.0));
+ CHK(nfails + nreals == N);
+ CHK(nfails < N/1000);
+ CHK(eq_eps(T.E, ref, 3*T.SE));
CHK(sdis_scene_ref_put(box_scn) == RES_OK);
CHK(sdis_scene_ref_put(square_scn) == RES_OK);
diff --git a/src/test_sdis_volumic_power2.c b/src/test_sdis_volumic_power2.c
@@ -0,0 +1,468 @@
+/* Copyright (C) 2016-2018 |Meso|Star> (contact@meso-star.com)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. */
+
+#include "sdis.h"
+#include "test_sdis_utils.h"
+#include <rsys/math.h>
+
+#define N 10000 /* #realisations */
+#define Pw 10000 /* Volumic power */
+#define NONE -1
+#define DELTA 0.01
+#define DELTA_PSQUARE 0.01
+
+struct reference {
+ double pos[3];
+ double temperature_2d; /* In celcius */
+ double temperature_3d; /* In celcius */
+};
+
+/* Temperature in Celcius. The reference is computed by EDF with Syrthes
+ * #realisations: 100000
+ *
+ * >>> Check 1
+ * 0.85 0 = 190.29 ~ 190.198 +/- 0.572596; #failures: 46
+ * 0.65 0 = 259.95 ~ 259.730 +/- 0.678251; #failures: 73
+ * 0.45 0 = 286.33 ~ 285.287 +/- 0.693572; #failures: 74
+ * 0.25 0 = 235.44 ~ 235.672 +/- 0.710927; #failures: 61
+ * 0.05 0 = 192.33 ~ 192.464 +/- 0.693148; #failures: 70
+ *-0.15 0 = 156.82 ~ 157.526 +/- 0.668902; #failures: 43
+ *-0.35 0 = 123.26 ~ 124.234 +/- 0.634061; #failures: 31
+ *-0.55 0 = 90.250 ~ 91.0285 +/- 0.566423; #failures: 32
+ *
+ * >>> Check 2
+ * 0.85 0 = 678.170 ~ 671.302 +/- 4.03424; #failures: 186
+ * 0.65 0 = 1520.84 ~ 1523.42 +/- 5.38182; #failures: 442
+ * 0.45 0 = 1794.57 ~ 1790.60 +/- 5.44808; #failures: 528
+ * 0.25 0 = 1429.74 ~ 1419.80 +/- 5.33467; #failures: 406 */
+
+static const double vertices[16/*#vertices*/*3/*#coords per vertex*/] = {
+ -0.5,-1.0,-0.5,
+ -0.5, 1.0,-0.5,
+ 0.5, 1.0,-0.5,
+ 0.5,-1.0,-0.5,
+ -0.5,-1.0, 0.5,
+ -0.5, 1.0, 0.5,
+ 0.5, 1.0, 0.5,
+ 0.5,-1.0, 0.5,
+ -0.1, 0.4,-0.5,
+ -0.1, 0.6,-0.5,
+ 0.1, 0.6,-0.5,
+ 0.1, 0.4,-0.5,
+ -0.1, 0.4, 0.5,
+ -0.1, 0.6, 0.5,
+ 0.1, 0.6, 0.5,
+ 0.1, 0.4, 0.5
+};
+static const size_t nvertices = sizeof(vertices)/sizeof(double[3]);
+
+static const size_t indices[36/*#triangles*/*3/*#indices per triangle*/]= {
+ 0, 4, 5, 5, 1, 0, /* Cuboid left */
+ 1, 5, 6, 6, 2, 1, /* Cuboid top */
+ 6, 7, 3, 3, 2, 6, /* Cuboid right */
+ 0, 3, 7, 7, 4, 0, /* Cuboid bottom */
+ /* Cuboid back */
+ 0, 1, 9, 9, 8, 0,
+ 1, 2, 10, 10, 9, 1,
+ 2, 3, 11, 11, 10, 2,
+ 3, 0, 8, 8, 11, 3,
+ /* Cuboid front */
+ 5, 4, 12, 12, 13, 5,
+ 5, 13, 14, 14, 6, 5,
+ 6, 14, 15, 15, 7, 6,
+ 7, 15, 12, 12, 4, 7,
+ 8, 12, 13, 13, 9, 8, /* Cube left */
+ 9, 13, 14, 14, 10, 9, /* Cube top */
+ 14, 15, 11, 11, 10, 14, /* Cube right */
+ 8, 11, 15, 15, 12, 8, /* Cube bottom */
+ 8, 9, 10, 10, 11, 8, /* Cube back */
+ 12, 15, 14, 14, 13, 12 /* Cube front */
+};
+static const size_t ntriangles = sizeof(indices)/sizeof(size_t[3]);
+
+/*******************************************************************************
+ * Geometry
+ ******************************************************************************/
+static void
+get_indices(const size_t itri, size_t ids[3], void* context)
+{
+ (void)context;
+ CHK(ids);
+ ids[0] = indices[itri*3+0];
+ ids[1] = indices[itri*3+1];
+ ids[2] = indices[itri*3+2];
+}
+
+static void
+get_position(const size_t ivert, double pos[3], void* context)
+{
+ (void)context;
+ CHK(pos);
+ pos[0] = vertices[ivert*3+0];
+ pos[1] = vertices[ivert*3+1];
+ pos[2] = vertices[ivert*3+2];
+}
+
+static void
+get_interface(const size_t itri, struct sdis_interface** bound, void* context)
+{
+ struct sdis_interface** interfaces = context;
+ CHK(context && bound);
+ *bound = interfaces[itri/2];
+}
+
+/*******************************************************************************
+ * Solid medium
+ ******************************************************************************/
+struct solid {
+ double cp;
+ double lambda;
+ double rho;
+ double delta;
+ double P;
+ double T;
+};
+
+static double
+solid_get_calorific_capacity
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ CHK(data != NULL && vtx != NULL);
+ return ((const struct solid*)sdis_data_cget(data))->cp;
+}
+
+static double
+solid_get_thermal_conductivity
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ CHK(data != NULL && vtx != NULL);
+ return ((const struct solid*)sdis_data_cget(data))->lambda;
+}
+
+static double
+solid_get_volumic_mass
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ CHK(data != NULL && vtx != NULL);
+ return ((const struct solid*)sdis_data_cget(data))->rho;
+}
+
+static double
+solid_get_delta
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ CHK(data != NULL && vtx != NULL);
+ return ((const struct solid*)sdis_data_cget(data))->delta;
+}
+
+static double
+solid_get_temperature
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ CHK(data != NULL && vtx != NULL);
+ return ((const struct solid*)sdis_data_cget(data))->T;
+}
+
+static double
+solid_get_volumic_power
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ CHK(data != NULL && vtx != NULL);
+ return ((const struct solid*)sdis_data_cget(data))->P;
+}
+
+/*******************************************************************************
+ * Fluid medium
+ ******************************************************************************/
+struct fluid {
+ double temperature;
+};
+
+static double
+fluid_get_temperature
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ CHK(data != NULL && vtx != NULL);
+ return ((const struct fluid*)sdis_data_cget(data))->temperature;
+}
+
+/*******************************************************************************
+ * Interfaces
+ ******************************************************************************/
+struct interf {
+ double h;
+ double temperature;
+};
+
+static double
+interface_get_convection_coef
+ (const struct sdis_interface_fragment* frag, struct sdis_data* data)
+{
+ CHK(frag && data);
+ return ((const struct interf*)sdis_data_cget(data))->h;
+}
+
+static double
+interface_get_temperature
+ (const struct sdis_interface_fragment* frag, struct sdis_data* data)
+{
+ CHK(frag && data);
+ return ((const struct interf*)sdis_data_cget(data))->temperature;
+}
+
+/*******************************************************************************
+ * Helper functions
+ ******************************************************************************/
+static void
+check(struct sdis_scene* scn, const struct reference refs[], const size_t nrefs)
+{
+ struct sdis_estimator* estimator = NULL;
+ struct sdis_mc T = SDIS_MC_NULL;
+ size_t nreals;
+ size_t nfails;
+ double pos[3] = {0,0};
+ size_t i;
+
+ FOR_EACH(i, 0, nrefs) {
+ double Tc;
+ pos[0] = refs[i].pos[0];
+ pos[1] = refs[i].pos[1];
+ pos[2] = refs[i].pos[2];
+
+ CHK(sdis_solve_probe(scn, N, pos, INF, 1.f, -1, 0, &estimator) == RES_OK);
+ CHK(sdis_estimator_get_temperature(estimator, &T) == RES_OK);
+ CHK(sdis_estimator_get_realisation_count(estimator, &nreals) == RES_OK);
+ CHK(sdis_estimator_get_failure_count(estimator, &nfails) == RES_OK);
+ Tc = T.E - 273.15; /* Convert in Celcius */
+ printf("Temperature at (%g %g %g) = %g ~ %g +/- %g [%g, %g]\n",
+ SPLIT3(pos), refs[i].temperature_2d, Tc, T.SE, Tc-3*T.SE, Tc+3*T.SE);
+ printf("#realisations: %lu; #failures: %lu\n",
+ (unsigned long)nreals, (unsigned long)nfails);
+ /*CHK(eq_eps(Tc, refs[i].temperature, T.SE*3));*/
+ CHK(sdis_estimator_ref_put(estimator) == RES_OK);
+ }
+}
+
+
+int
+main(int argc, char** argv)
+{
+ struct mem_allocator allocator;
+ struct solid* solid_param = NULL;
+ struct fluid* fluid_param = NULL;
+ struct interf* interf_param = NULL;
+ struct sdis_device* dev = NULL;
+ struct sdis_data* data = NULL;
+ struct sdis_medium* fluid1 = NULL;
+ struct sdis_medium* fluid2 = NULL;
+ struct sdis_medium* solid1 = NULL;
+ struct sdis_medium* solid2 = NULL;
+ struct sdis_scene* scn = NULL;
+ struct sdis_fluid_shader fluid_shader = SDIS_FLUID_SHADER_NULL;
+ struct sdis_solid_shader solid_shader = SDIS_SOLID_SHADER_NULL;
+ struct sdis_interface_shader interf_shader = SDIS_INTERFACE_SHADER_NULL;
+ struct sdis_interface* interf_solid1_adiabatic = NULL;
+ struct sdis_interface* interf_solid2_adiabatic = NULL;
+ struct sdis_interface* interf_solid1_solid2 = NULL;
+ struct sdis_interface* interf_solid1_fluid1 = NULL;
+ struct sdis_interface* interf_solid1_fluid2 = NULL;
+ struct sdis_interface* interfaces[18 /*#rectangles*/];
+ /* In celcius. Computed by EDF with Syrthes */
+ const struct reference refs1[] = { /* Lambda1=1, Lambda2=10, Pw = 10000 */
+ {{0, 0.85, 0}, 190.29, 189.13},
+ {{0, 0.65, 0}, 259.95, 247.09},
+ {{0, 0.45, 0}, 286.33, 308.42},
+ {{0, 0.25, 0}, 235.44, 233.55},
+ {{0, 0.05, 0}, 192.33, 192.30},
+ {{0,-0.15, 0}, 156.82, 156.98},
+ {{0,-0.35, 0}, 123.26, 123.43},
+ {{0,-0.55, 0}, 90.250, 90.040}
+ };
+ const struct reference refs2[] = { /* Lambda1=0.1, Lambda2=10, Pw=10000 */
+ {{0, 0.85}, 678.170, -1},
+ {{0, 0.65}, 1520.84, -1},
+ {{0, 0.45}, 1794.57, -1},
+ {{0, 0.25}, 1429.74, -1}
+ };
+ (void)argc, (void)argv;
+
+ CHK(mem_init_proxy_allocator(&allocator, &mem_default_allocator) == RES_OK);
+ CHK(sdis_device_create
+ (NULL, &allocator, SDIS_NTHREADS_DEFAULT, 1, &dev) == RES_OK);
+
+ /* Setup the fluid shader */
+ fluid_shader.temperature = fluid_get_temperature;
+ fluid_shader.calorific_capacity = dummy_medium_getter;
+ fluid_shader.volumic_mass = dummy_medium_getter;
+
+ /* Create the fluid1 medium */
+ CHK(sdis_data_create
+ (dev, sizeof(struct fluid), ALIGNOF(struct fluid), NULL, &data) == RES_OK);
+ fluid_param = sdis_data_get(data);
+ fluid_param->temperature = 373.15;
+ CHK(sdis_fluid_create(dev, &fluid_shader, data, &fluid1) == RES_OK);
+ CHK(sdis_data_ref_put(data) == RES_OK);
+
+ /* Create the fluid2 medium */
+ CHK(sdis_data_create
+ (dev, sizeof(struct fluid), ALIGNOF(struct fluid), NULL, &data) == RES_OK);
+ fluid_param = sdis_data_get(data);
+ fluid_param->temperature = 273.15;
+ CHK(sdis_fluid_create(dev, &fluid_shader, data, &fluid2) == RES_OK);
+ CHK(sdis_data_ref_put(data) == RES_OK);
+
+ /* Setup the solid shader */
+ solid_shader.calorific_capacity = solid_get_calorific_capacity;
+ solid_shader.thermal_conductivity = solid_get_thermal_conductivity;
+ solid_shader.volumic_mass = solid_get_volumic_mass;
+ solid_shader.delta_solid = solid_get_delta;
+ solid_shader.temperature = solid_get_temperature;
+ solid_shader.volumic_power = solid_get_volumic_power;
+
+ /* Create the solid1 medium */
+ CHK(sdis_data_create
+ (dev, sizeof(struct solid), ALIGNOF(struct solid), NULL, &data) == RES_OK);
+ solid_param = sdis_data_get(data);
+ solid_param->cp = 500000;
+ solid_param->rho = 1000;
+ solid_param->lambda = 1;
+ solid_param->delta = DELTA;
+ solid_param->P = SDIS_VOLUMIC_POWER_NONE;
+ solid_param->T = -1;
+ CHK(sdis_solid_create(dev, &solid_shader, data, &solid1) == RES_OK);
+ CHK(sdis_data_ref_put(data) == RES_OK);
+
+ /* Create the solid2 medium */
+ CHK(sdis_data_create
+ (dev, sizeof(struct solid), ALIGNOF(struct solid), NULL, &data) == RES_OK);
+ solid_param = sdis_data_get(data);
+ solid_param->cp = 500000;
+ solid_param->rho = 1000;
+ solid_param->lambda = 10;
+ solid_param->delta = DELTA_PSQUARE;
+ solid_param->P = Pw;
+ solid_param->T = -1;
+ CHK(sdis_solid_create(dev, &solid_shader, data, &solid2) == RES_OK);
+ CHK(sdis_data_ref_put(data) == RES_OK);
+
+ /* Create the solid1/solid2 interface */
+ CHK(sdis_data_create (dev, sizeof(struct interf), ALIGNOF(struct interf),
+ NULL, &data) == RES_OK);
+ CHK(sdis_interface_create(dev, solid2, solid1, &SDIS_INTERFACE_SHADER_NULL,
+ NULL, &interf_solid1_solid2) == RES_OK);
+ CHK(sdis_data_ref_put(data) == RES_OK);
+
+ /* Setup the interface shader */
+ interf_shader.convection_coef = interface_get_convection_coef;
+
+ /* Create the adiabatic interfaces */
+ CHK(sdis_data_create (dev, sizeof(struct interf), ALIGNOF(struct interf),
+ NULL, &data) == RES_OK);
+ interf_param = sdis_data_get(data);
+ interf_param->h = 0;
+ CHK(sdis_interface_create(dev, solid1, fluid1, &interf_shader, data,
+ &interf_solid1_adiabatic) == RES_OK);
+ CHK(sdis_interface_create(dev, solid2, fluid1, &interf_shader, data,
+ &interf_solid2_adiabatic) == RES_OK);
+ CHK(sdis_data_ref_put(data) == RES_OK);
+
+ /* Setup the interface shader */
+ interf_shader.front.temperature = interface_get_temperature;
+
+ /* Create the solid1/fluid1 interface */
+ CHK(sdis_data_create (dev, sizeof(struct interf), ALIGNOF(struct interf),
+ NULL, &data) == RES_OK);
+ interf_param = sdis_data_get(data);
+ interf_param->h = 5;
+ interf_param->temperature = NONE;
+ CHK(sdis_interface_create(dev, solid1, fluid1, &interf_shader, data,
+ &interf_solid1_fluid1) == RES_OK);
+ CHK(sdis_data_ref_put(data) == RES_OK);
+
+ /* Create the solid1/fluid2 interace */
+ CHK(sdis_data_create (dev, sizeof(struct interf), ALIGNOF(struct interf),
+ NULL, &data) == RES_OK);
+ interf_param = sdis_data_get(data);
+ interf_param->h = 10;
+ interf_param->temperature = NONE;
+ CHK(sdis_interface_create(dev, solid1, fluid2, &interf_shader, data,
+ &interf_solid1_fluid2) == RES_OK);
+ CHK(sdis_data_ref_put(data) == RES_OK);
+
+ /* Map the interfaces to their faces */
+ interfaces[0] = interf_solid1_adiabatic;
+ interfaces[1] = interf_solid1_fluid1;
+ interfaces[2] = interf_solid1_adiabatic;
+ interfaces[3] = interf_solid1_fluid2;
+ interfaces[4] = interf_solid1_adiabatic;
+ interfaces[5] = interf_solid1_adiabatic;
+ interfaces[6] = interf_solid1_adiabatic;
+ interfaces[7] = interf_solid1_adiabatic;
+ interfaces[8] = interf_solid1_adiabatic;
+ interfaces[9] = interf_solid1_adiabatic;
+ interfaces[10] = interf_solid1_adiabatic;
+ interfaces[11] = interf_solid1_adiabatic;
+ interfaces[12] = interf_solid1_solid2;
+ interfaces[13] = interf_solid1_solid2;
+ interfaces[14] = interf_solid1_solid2;
+ interfaces[15] = interf_solid1_solid2;
+ interfaces[16] = interf_solid2_adiabatic;
+ interfaces[17] = interf_solid2_adiabatic;
+
+ /* Create the scene */
+ CHK(sdis_scene_create(dev, ntriangles, get_indices, get_interface,
+ nvertices, get_position, interfaces, &scn) == RES_OK);
+
+#if 0
+ dump_mesh(stdout, vertices, nvertices, indices, ntriangles);
+ exit(0);
+#endif
+
+ printf(">>> Check 1\n");
+ check(scn, refs1, sizeof(refs1)/sizeof(struct reference));
+
+ /* Update the scene */
+ CHK(sdis_scene_ref_put(scn) == RES_OK);
+ data = sdis_medium_get_data(solid1);
+ solid_param = sdis_data_get(data);
+ solid_param->lambda = 0.1;
+ CHK(sdis_scene_create(dev, ntriangles, get_indices, get_interface,
+ nvertices, get_position, interfaces, &scn) == RES_OK);
+
+ printf("\n>>> Check 2\n");
+ check(scn, refs2, sizeof(refs2)/sizeof(struct reference));
+
+ /* Release the interfaces */
+ CHK(sdis_interface_ref_put(interf_solid1_adiabatic) == RES_OK);
+ CHK(sdis_interface_ref_put(interf_solid2_adiabatic) == RES_OK);
+ CHK(sdis_interface_ref_put(interf_solid1_fluid1) == RES_OK);
+ CHK(sdis_interface_ref_put(interf_solid1_fluid2) == RES_OK);
+ CHK(sdis_interface_ref_put(interf_solid1_solid2) == RES_OK);
+
+ /* Release the media */
+ CHK(sdis_medium_ref_put(fluid1) == RES_OK);
+ CHK(sdis_medium_ref_put(fluid2) == RES_OK);
+ CHK(sdis_medium_ref_put(solid1) == RES_OK);
+ CHK(sdis_medium_ref_put(solid2) == RES_OK);
+
+ CHK(sdis_scene_ref_put(scn) == RES_OK);
+ CHK(sdis_device_ref_put(dev) == RES_OK);
+
+ check_memory_allocator(&allocator);
+ mem_shutdown_proxy_allocator(&allocator);
+ CHK(mem_allocated_size() == 0);
+ return 0;
+}
+
diff --git a/src/test_sdis_volumic_power2_2d.c b/src/test_sdis_volumic_power2_2d.c
@@ -0,0 +1,502 @@
+/* Copyright (C) 2016-2018 |Meso|Star> (contact@meso-star.com)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. */
+
+#include "sdis.h"
+#include "test_sdis_utils.h"
+#include <rsys/math.h>
+
+#define N 10000 /* #realisations */
+#define Pw 10000 /* Volumic power */
+#define NONE -1
+
+/* H delta T */
+#define Tboundary1 NONE
+#define Tboundary2 NONE
+#define DELTA 0.01
+#define Tref 286.83 /* In Celsius. Computed with Syrthes at the position 0.5 */
+
+/* Dirichlets */
+/*#define Tboundary1 373.15*/
+/*#define Tboundary2 273.15*/
+/*#define DELTA 0.01*/
+/*#define Tref 246.93*/ /* In Celsius. Computed with Syrthes at the position 0.5 */
+
+/* Temperature in Celcius. The reference is computed by EDF with Syrthes
+ * #realisations: 100000
+ *
+ * >>> Check1
+ * 0.85 = 190.29 ~ 189.322 +/- 0.566717; #failures: 51
+ * 0.65 = 259.95 ~ 259.995 +/- 0.674453; #failures: 82
+ * 0.45 = 286.33 ~ 285.928 +/- 0.691044; #failures: 76
+ * 0.25 = 235.44 ~ 234.672 +/- 0.700354; #failures: 80
+ * 0.05 = 192.33 ~ 191.977 +/- 0.690793; #failures: 64
+ *-0.15 = 156.82 ~ 155.765 +/- 0.660722; #failures: 40
+ *-0.35 = 123.26 ~ 122.973 +/- 0.621093; #failures: 29
+ *-0.55 = 90.250 ~ 90.3501 +/- 0.561255; #failures: 27
+ *
+ * >>> Check 2
+ * 0.85 = 678.170 ~ 662.616 +/- 3.97997; #failures: 221
+ * 0.65 = 1520.84 ~ 1486.35 +/- 5.25785; #failures: 474
+ * 0.45 = 1794.57 ~ 1767.21 +/- 5.36318; #failures: 584
+ * 0.25 = 1429.74 ~ 1401.39 +/- 5.25579; #failures: 465
+ *
+ * >>> Check 3
+ * 0.85 = 83.99 ~ 84.0098 +/- 0.115932; #failures: 51
+ * 0.65 = 73.90 ~ 73.9596 +/- 0.138835; #failures: 82
+ * 0.45 = 68.43 ~ 70.0292 +/- 0.144928; #failures: 76
+ * 0.25 = 60.61 ~ 61.4412 +/- 0.153980; #failures: 80
+ * 0.05 = 52.09 ~ 51.9452 +/- 0.158045; #failures: 64
+ *-0.15 = 42.75 ~ 42.9072 +/- 0.156546; #failures: 40
+ *-0.35 = 33.04 ~ 33.9338 +/- 0.149751; #failures: 29
+ *-0.55 = 24.58 ~ 24.7237 +/- 0.136441; #failures: 27 */
+
+/*
+ * _\ T1
+ * / /
+ * \__/
+ * ///+-----H1-------+///
+ * ///| |///
+ * ///| +------+ |///
+ * ///| |LAMBDA| |///
+ * ///| | Pw | |///
+ * ///| +------+ |///
+ * ///| |///
+ * ///| |///
+ * ///| LAMBDA1 |///
+ * ///| |///
+ * ///| |///
+ * ///| |///
+ * ///+-----H2-------+///
+ * _\ T2
+ * / /
+ * \__/
+ */
+
+struct reference {
+ double pos[2];
+ double temperature; /* In celcius */
+};
+
+static const double vertices[8/*#vertices*/*2/*#coords per vertex*/] = {
+ -0.5,-1.0,
+ -0.5, 1.0,
+ 0.5, 1.0,
+ 0.5,-1.0,
+ -0.1, 0.4,
+ -0.1, 0.6,
+ 0.1, 0.6,
+ 0.1, 0.4
+};
+static const size_t nvertices = sizeof(vertices)/sizeof(double[2]);
+
+static const size_t indices[8/*#segments*/*2/*#indices per segment*/]= {
+ 0, 1, /* Rectangle left */
+ 1, 2, /* Rectangle top */
+ 2, 3, /* Rectangle right */
+ 3, 0, /* Rectangle bottom */
+ 4, 5, /* Square left */
+ 5, 6, /* Square top */
+ 6, 7, /* Square right */
+ 7, 4 /* Square bottom */
+};
+static const size_t nsegments = sizeof(indices)/sizeof(size_t[2]);
+
+/*******************************************************************************
+ * Geometry
+ ******************************************************************************/
+static void
+get_indices(const size_t iseg, size_t ids[2], void* context)
+{
+ (void)context;
+ CHK(ids);
+ ids[0] = indices[iseg*2+0];
+ ids[1] = indices[iseg*2+1];
+}
+
+static void
+get_position(const size_t ivert, double pos[2], void* context)
+{
+ (void)context;
+ CHK(pos);
+ pos[0] = vertices[ivert*2+0];
+ pos[1] = vertices[ivert*2+1];
+}
+
+static void
+get_interface(const size_t iseg, struct sdis_interface** bound, void* context)
+{
+ struct sdis_interface** interfaces = context;
+ CHK(context && bound);
+ *bound = interfaces[iseg];
+}
+
+/*******************************************************************************
+ * Solid medium
+ ******************************************************************************/
+struct solid {
+ double cp;
+ double lambda;
+ double rho;
+ double delta;
+ double P;
+ double T;
+};
+
+static double
+solid_get_calorific_capacity
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ CHK(data != NULL && vtx != NULL);
+ return ((const struct solid*)sdis_data_cget(data))->cp;
+}
+
+static double
+solid_get_thermal_conductivity
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ CHK(data != NULL && vtx != NULL);
+ return ((const struct solid*)sdis_data_cget(data))->lambda;
+}
+
+static double
+solid_get_volumic_mass
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ CHK(data != NULL && vtx != NULL);
+ return ((const struct solid*)sdis_data_cget(data))->rho;
+}
+
+static double
+solid_get_delta
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ CHK(data != NULL && vtx != NULL);
+ return ((const struct solid*)sdis_data_cget(data))->delta;
+}
+
+static double
+solid_get_temperature
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ CHK(data != NULL && vtx != NULL);
+ return ((const struct solid*)sdis_data_cget(data))->T;
+}
+
+static double
+solid_get_volumic_power
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ CHK(data != NULL && vtx != NULL);
+ return ((const struct solid*)sdis_data_cget(data))->P;
+}
+
+/*******************************************************************************
+ * Fluid medium
+ ******************************************************************************/
+struct fluid {
+ double temperature;
+};
+
+static double
+fluid_get_temperature
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ CHK(data != NULL && vtx != NULL);
+ return ((const struct fluid*)sdis_data_cget(data))->temperature;
+}
+
+
+/*******************************************************************************
+ * Interfaces
+ ******************************************************************************/
+struct interf {
+ double h;
+ double temperature;
+};
+
+static double
+interface_get_convection_coef
+ (const struct sdis_interface_fragment* frag, struct sdis_data* data)
+{
+ CHK(frag && data);
+ return ((const struct interf*)sdis_data_cget(data))->h;
+}
+
+static double
+interface_get_temperature
+ (const struct sdis_interface_fragment* frag, struct sdis_data* data)
+{
+ CHK(frag && data);
+ return ((const struct interf*)sdis_data_cget(data))->temperature;
+}
+
+/*******************************************************************************
+ * Helper functions
+ ******************************************************************************/
+static void
+check(struct sdis_scene* scn, const struct reference refs[], const size_t nrefs)
+{
+ struct sdis_estimator* estimator = NULL;
+ struct sdis_mc T = SDIS_MC_NULL;
+ size_t nreals;
+ size_t nfails;
+ double pos[2] = {0,0};
+ size_t i;
+
+ FOR_EACH(i, 0, nrefs) {
+ double Tc;
+ pos[0] = refs[i].pos[0];
+ pos[1] = refs[i].pos[1];
+
+ CHK(sdis_solve_probe(scn, N, pos, INF, 1.f, -1, 0, &estimator) == RES_OK);
+ CHK(sdis_estimator_get_temperature(estimator, &T) == RES_OK);
+ CHK(sdis_estimator_get_realisation_count(estimator, &nreals) == RES_OK);
+ CHK(sdis_estimator_get_failure_count(estimator, &nfails) == RES_OK);
+ Tc = T.E - 273.15; /* Convert in Celcius */
+ printf("Temperature at (%g %g) = %g ~ %g +/- %g [%g, %g]\n",
+ SPLIT2(pos), refs[i].temperature, Tc, T.SE, Tc-3*T.SE, Tc+3*T.SE);
+ printf("#realisations: %lu; #failures: %lu\n",
+ (unsigned long)nreals, (unsigned long)nfails);
+ /*CHK(eq_eps(Tc, refs[i].temperature, T.SE*3));*/
+ CHK(sdis_estimator_ref_put(estimator) == RES_OK);
+ }
+}
+
+/*******************************************************************************
+ * Test
+ ******************************************************************************/
+int
+main(int argc, char** argv)
+{
+ struct mem_allocator allocator;
+ struct solid* solid_param = NULL;
+ struct fluid* fluid_param = NULL;
+ struct interf* interf_param = NULL;
+ struct sdis_device* dev = NULL;
+ struct sdis_data* data = NULL;
+ struct sdis_medium* fluid1 = NULL;
+ struct sdis_medium* fluid2 = NULL;
+ struct sdis_medium* solid1 = NULL;
+ struct sdis_medium* solid2 = NULL;
+ struct sdis_scene* scn = NULL;
+ struct sdis_fluid_shader fluid_shader = SDIS_FLUID_SHADER_NULL;
+ struct sdis_solid_shader solid_shader = SDIS_SOLID_SHADER_NULL;
+ struct sdis_interface_shader interf_shader = SDIS_INTERFACE_SHADER_NULL;
+ struct sdis_interface* interf_adiabatic = NULL;
+ struct sdis_interface* interf_solid1_solid2 = NULL;
+ struct sdis_interface* interf_solid1_fluid1 = NULL;
+ struct sdis_interface* interf_solid1_fluid2 = NULL;
+ struct sdis_interface* interfaces[8 /*#segment*/];
+
+ /* In celcius. Computed by EDF with Syrthes */
+ const struct reference refs1[] = { /* Lambda1=1, Lambda2=10, Pw = 10000 */
+ {{0, 0.85}, 190.29},
+ {{0, 0.65}, 259.95},
+ {{0, 0.45}, 286.33},
+ {{0, 0.25}, 235.44},
+ {{0, 0.05}, 192.33},
+ {{0,-0.15}, 156.82},
+ {{0,-0.35}, 123.26},
+ {{0,-0.55}, 90.250}
+ };
+ const struct reference refs2[] = { /* Lambda1=0.1, Lambda2=10, Pw=10000 */
+ {{0, 0.85}, 678.17},
+ {{0, 0.65}, 1520.84},
+ {{0, 0.45}, 1794.57},
+ {{0, 0.25}, 1429.74}
+ };
+ const struct reference refs3[] = { /* Lambda1=1, Lambda2=10, Pw=NONE */
+ {{0, 0.85}, 83.99},
+ {{0, 0.65}, 73.90},
+ {{0, 0.45}, 68.43},
+ {{0, 0.25}, 60.61},
+ {{0, 0.05}, 52.09},
+ {{0,-0.15}, 42.75},
+ {{0,-0.35}, 33.04},
+ {{0,-0.55}, 24.58}
+ };
+ (void)argc, (void)argv;
+
+ CHK(mem_init_proxy_allocator(&allocator, &mem_default_allocator) == RES_OK);
+ CHK(sdis_device_create
+ (NULL, &allocator, SDIS_NTHREADS_DEFAULT, 1, &dev) == RES_OK);
+
+ /* Setup the fluid shader */
+ fluid_shader.temperature = fluid_get_temperature;
+ fluid_shader.calorific_capacity = dummy_medium_getter;
+ fluid_shader.volumic_mass = dummy_medium_getter;
+
+ /* Create the fluid1 medium */
+ CHK(sdis_data_create
+ (dev, sizeof(struct fluid), ALIGNOF(struct fluid), NULL, &data) == RES_OK);
+ fluid_param = sdis_data_get(data);
+ fluid_param->temperature = 373.15;
+ CHK(sdis_fluid_create(dev, &fluid_shader, data, &fluid1) == RES_OK);
+ CHK(sdis_data_ref_put(data) == RES_OK);
+
+ /* Create the fluid2 medium */
+ CHK(sdis_data_create
+ (dev, sizeof(struct fluid), ALIGNOF(struct fluid), NULL, &data) == RES_OK);
+ fluid_param = sdis_data_get(data);
+ fluid_param->temperature = 273.15;
+ CHK(sdis_fluid_create(dev, &fluid_shader, data, &fluid2) == RES_OK);
+ CHK(sdis_data_ref_put(data) == RES_OK);
+
+ /* Setup the solid shader */
+ solid_shader.calorific_capacity = solid_get_calorific_capacity;
+ solid_shader.thermal_conductivity = solid_get_thermal_conductivity;
+ solid_shader.volumic_mass = solid_get_volumic_mass;
+ solid_shader.delta_solid = solid_get_delta;
+ solid_shader.temperature = solid_get_temperature;
+ solid_shader.volumic_power = solid_get_volumic_power;
+
+ /* Create the solid1 medium */
+ CHK(sdis_data_create
+ (dev, sizeof(struct solid), ALIGNOF(struct solid), NULL, &data) == RES_OK);
+ solid_param = sdis_data_get(data);
+ solid_param->cp = 500000;
+ solid_param->rho = 1000;
+ solid_param->lambda = 1;
+ solid_param->delta = DELTA;
+ solid_param->P = SDIS_VOLUMIC_POWER_NONE;
+ solid_param->T = -1;
+ CHK(sdis_solid_create(dev, &solid_shader, data, &solid1) == RES_OK);
+ CHK(sdis_data_ref_put(data) == RES_OK);
+
+ /* Create the solid2 medium */
+ CHK(sdis_data_create
+ (dev, sizeof(struct solid), ALIGNOF(struct solid), NULL, &data) == RES_OK);
+ solid_param = sdis_data_get(data);
+ solid_param->cp = 500000;
+ solid_param->rho = 1000;
+ solid_param->lambda = 10;
+ solid_param->delta = DELTA;
+ solid_param->P = Pw;
+ solid_param->T = -1;
+ CHK(sdis_solid_create(dev, &solid_shader, data, &solid2) == RES_OK);
+ CHK(sdis_data_ref_put(data) == RES_OK);
+
+ /* Create the solid1/solid2 interface */
+ CHK(sdis_data_create (dev, sizeof(struct interf), ALIGNOF(struct interf),
+ NULL, &data) == RES_OK);
+ CHK(sdis_interface_create(dev, solid2, solid1, &SDIS_INTERFACE_SHADER_NULL,
+ NULL, &interf_solid1_solid2) == RES_OK);
+ CHK(sdis_data_ref_put(data) == RES_OK);
+
+ /* Setup the interface shader */
+ interf_shader.convection_coef = interface_get_convection_coef;
+
+ /* Create the adiabatic interface */
+ CHK(sdis_data_create (dev, sizeof(struct interf), ALIGNOF(struct interf),
+ NULL, &data) == RES_OK);
+ interf_param = sdis_data_get(data);
+ interf_param->h = 0;
+ CHK(sdis_interface_create(dev, solid1, fluid1, &interf_shader, data,
+ &interf_adiabatic) == RES_OK);
+ CHK(sdis_data_ref_put(data) == RES_OK);
+
+ /* Setup the interface shader */
+ interf_shader.front.temperature = interface_get_temperature;
+
+ /* Create the solid1/fluid1 interface */
+ CHK(sdis_data_create (dev, sizeof(struct interf), ALIGNOF(struct interf),
+ NULL, &data) == RES_OK);
+ interf_param = sdis_data_get(data);
+ interf_param->h = 5;
+ interf_param->temperature = Tboundary1;
+ CHK(sdis_interface_create(dev, solid1, fluid1, &interf_shader, data,
+ &interf_solid1_fluid1) == RES_OK);
+ CHK(sdis_data_ref_put(data) == RES_OK);
+
+ /* Create the solid1/fluid2 interace */
+ CHK(sdis_data_create (dev, sizeof(struct interf), ALIGNOF(struct interf),
+ NULL, &data) == RES_OK);
+ interf_param = sdis_data_get(data);
+ interf_param->h = 10;
+ interf_param->temperature = Tboundary2;
+ CHK(sdis_interface_create(dev, solid1, fluid2, &interf_shader, data,
+ &interf_solid1_fluid2) == RES_OK);
+ CHK(sdis_data_ref_put(data) == RES_OK);
+
+
+ /* Map the interfaces to their square segments */
+ interfaces[0] = interf_adiabatic;
+ interfaces[1] = interf_solid1_fluid1;
+ interfaces[2] = interf_adiabatic;
+ interfaces[3] = interf_solid1_fluid2;
+ interfaces[4] = interf_solid1_solid2;
+ interfaces[5] = interf_solid1_solid2;
+ interfaces[6] = interf_solid1_solid2;
+ interfaces[7] = interf_solid1_solid2;
+
+ /* Create the scene */
+ CHK(sdis_scene_2d_create(dev, nsegments, get_indices, get_interface,
+ nvertices, get_position, interfaces, &scn) == RES_OK);
+
+ printf(">>> Check 1\n");
+ check(scn, refs1, sizeof(refs1)/sizeof(struct reference));
+
+ /* Update the scene */
+ CHK(sdis_scene_ref_put(scn) == RES_OK);
+ data = sdis_medium_get_data(solid1);
+ solid_param = sdis_data_get(data);
+ solid_param->lambda = 0.1;
+ CHK(sdis_scene_2d_create(dev, nsegments, get_indices, get_interface,
+ nvertices, get_position, interfaces, &scn) == RES_OK);
+
+ printf("\n>>> Check 2\n");
+ check(scn, refs2, sizeof(refs2)/sizeof(struct reference));
+
+ /* Update the scene */
+ CHK(sdis_scene_ref_put(scn) == RES_OK);
+ data = sdis_medium_get_data(solid1);
+ solid_param = sdis_data_get(data);
+ solid_param->lambda = 1;
+ data = sdis_medium_get_data(solid2);
+ solid_param = sdis_data_get(data);
+ solid_param->lambda = 10;
+ solid_param->P = SDIS_VOLUMIC_POWER_NONE;
+ CHK(sdis_scene_2d_create(dev, nsegments, get_indices, get_interface,
+ nvertices, get_position, interfaces, &scn) == RES_OK);
+
+ printf("\n>>> Check 3\n");
+ check(scn, refs3, sizeof(refs3)/sizeof(struct reference));
+
+#if 0
+ dump_segments(stdout, vertices, nvertices, indices, nsegments);
+ exit(0);
+#endif
+
+ /* Release the interfaces */
+ CHK(sdis_interface_ref_put(interf_adiabatic) == RES_OK);
+ CHK(sdis_interface_ref_put(interf_solid1_fluid1) == RES_OK);
+ CHK(sdis_interface_ref_put(interf_solid1_fluid2) == RES_OK);
+ CHK(sdis_interface_ref_put(interf_solid1_solid2) == RES_OK);
+
+ /* Release the media */
+ CHK(sdis_medium_ref_put(fluid1) == RES_OK);
+ CHK(sdis_medium_ref_put(fluid2) == RES_OK);
+ CHK(sdis_medium_ref_put(solid1) == RES_OK);
+ CHK(sdis_medium_ref_put(solid2) == RES_OK);
+
+ CHK(sdis_scene_ref_put(scn) == RES_OK);
+ CHK(sdis_device_ref_put(dev) == RES_OK);
+
+ check_memory_allocator(&allocator);
+ mem_shutdown_proxy_allocator(&allocator);
+ CHK(mem_allocated_size() == 0);
+ return 0;
+}
+
diff --git a/src/test_sdis_volumic_power3_2d.c b/src/test_sdis_volumic_power3_2d.c
@@ -0,0 +1,465 @@
+/* Copyright (C) 2016-2018 |Meso|Star> (contact@meso-star.com)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. */
+
+#include "sdis.h"
+#include "test_sdis_utils.h"
+#include <rsys/math.h>
+
+#define Pw 10000.0 /* Volumic power */
+#define LAMBDA 10.0 /* Lambda of the middle slab */
+#define LAMBDA1 1.0 /* Lambda of the upper slab */
+#define LAMBDA2 LAMBDA1 /* Lambda of the lower slab */
+#define T1 373.15 /* Temperature of the upper fluid */
+#define T2 273.15 /* Temperature of the lower fluid */
+#define H1 5.0 /* Convection coef between the upper slab and the fluid */
+#define H2 10.0 /* Convection coef between the lower slab and the fluid */
+#define DELTA 0.01 /* Delta of the middle slab */
+#define DELTA1 0.02 /* Delta of the upper slab */
+#define DELTA2 0.07 /* Delta of the lower slab */
+#define L 0.2 /* Size of the middle slab */
+#define L1 0.4 /* Size of the upper slab */
+#define L2 1.4 /* Size of the lower slab */
+
+#define N 10000 /* #realisations */
+
+/* Analitically computed temperatures wrt the previous parameters.*/
+#define Tp1 648.6217
+#define Tp2 335.4141
+#define Ta 1199.5651
+#define Tb 1207.1122
+
+/* Fixed temperatures */
+#define UNKNOWN_TEMPERATURE -1
+#define Tsolid1_fluid UNKNOWN_TEMPERATURE /*Tp1*/
+#define Tsolid2_fluid UNKNOWN_TEMPERATURE /*Tp2*/
+#define Tsolid_solid1 UNKNOWN_TEMPERATURE /*Ta*/
+#define Tsolid_solid2 UNKNOWN_TEMPERATURE /*Tb*/
+
+#define PROBE_POS 1.8
+
+/*
+ * The 2D scene is composed of 3 stacked solid slabs whose middle slab has a
+ * volumic power. The +/-X sides of the slabs are stretched far away to
+ * simulate a 1D case. The upper and lower bounds of the "sandwich" has a
+ * convective exchange with the surrounding fluid whose temperature is known.
+ *
+ * _\ T1
+ * / /
+ * \__/
+ * ... -----H1------ ... Tp1
+ * LAMBDA1
+ *
+ * ... ------------- ... Ta
+ * LAMBDA, Pw
+ * ... ------------- ... Tb
+ *
+ * LAMBDA2
+ *
+ *
+ * ... -----H2------ ... Tp2
+ * _\ T2
+ * / /
+ * \__/
+ */
+
+static const double vertices[8/*#vertices*/*2/*#coords per vertex*/] = {
+ -100000.5, 0.0, /* 0 */
+ -100000.5, 1.4, /* 1 */
+ -100000.5, 1.6, /* 2 */
+ -100000.5, 2.0, /* 3 */
+ 100000.5, 2.0, /* 4 */
+ 100000.5, 1.6, /* 5 */
+ 100000.5, 1.4, /* 6 */
+ 100000.5, 0.0 /* 7 */
+};
+static const size_t nvertices = sizeof(vertices)/sizeof(double[2]);
+
+static const size_t indices[10/*#segments*/*2/*#indices per segment*/]= {
+ 0, 1,
+ 1, 2,
+ 2, 3,
+ 3, 4,
+ 4, 5,
+ 5, 6,
+ 6, 7,
+ 7, 0,
+ 6, 1,
+ 2, 5
+};
+static const size_t nsegments = sizeof(indices)/sizeof(size_t[2]);
+
+/*******************************************************************************
+ * Geometry
+ ******************************************************************************/
+static void
+get_indices(const size_t iseg, size_t ids[2], void* context)
+{
+ (void)context;
+ CHK(ids);
+ ids[0] = indices[iseg*2+0];
+ ids[1] = indices[iseg*2+1];
+}
+
+static void
+get_position(const size_t ivert, double pos[2], void* context)
+{
+ (void)context;
+ CHK(pos);
+ pos[0] = vertices[ivert*2+0];
+ pos[1] = vertices[ivert*2+1];
+}
+
+static void
+get_interface(const size_t iseg, struct sdis_interface** bound, void* context)
+{
+ struct sdis_interface** interfaces = context;
+ CHK(context && bound);
+ *bound = interfaces[iseg];
+}
+
+/*******************************************************************************
+ * Solid medium
+ ******************************************************************************/
+struct solid {
+ double cp;
+ double lambda;
+ double rho;
+ double delta;
+ double volumic_power;
+ double temperature;
+};
+
+static double
+solid_get_calorific_capacity
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ CHK(data != NULL && vtx != NULL);
+ return ((const struct solid*)sdis_data_cget(data))->cp;
+}
+
+static double
+solid_get_thermal_conductivity
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ CHK(data != NULL && vtx != NULL);
+ return ((const struct solid*)sdis_data_cget(data))->lambda;
+}
+
+static double
+solid_get_volumic_mass
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ CHK(data != NULL && vtx != NULL);
+ return ((const struct solid*)sdis_data_cget(data))->rho;
+}
+
+static double
+solid_get_delta
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ CHK(data != NULL && vtx != NULL);
+ return ((const struct solid*)sdis_data_cget(data))->delta;
+}
+
+static double
+solid_get_temperature
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ CHK(data != NULL && vtx != NULL);
+ return ((const struct solid*)sdis_data_cget(data))->temperature;
+}
+
+static double
+solid_get_volumic_power
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ CHK(data != NULL && vtx != NULL);
+ return ((const struct solid*)sdis_data_cget(data))->volumic_power;
+}
+
+/*******************************************************************************
+ * Fluid medium
+ ******************************************************************************/
+struct fluid {
+ double temperature_lower;
+ double temperature_upper;
+};
+
+static double
+fluid_get_temperature
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ const struct fluid* fluid;
+ CHK(data != NULL && vtx != NULL);
+ fluid = sdis_data_cget(data);
+ return vtx->P[1] < 1 ? fluid->temperature_lower : fluid->temperature_upper;
+}
+
+
+/*******************************************************************************
+ * Interfaces
+ ******************************************************************************/
+struct interf {
+ double h;
+ double temperature;
+};
+
+static double
+interface_get_convection_coef
+ (const struct sdis_interface_fragment* frag, struct sdis_data* data)
+{
+ CHK(frag && data);
+ return ((const struct interf*)sdis_data_cget(data))->h;
+}
+
+static double
+interface_get_temperature
+ (const struct sdis_interface_fragment* frag, struct sdis_data* data)
+{
+ CHK(frag && data);
+ return ((const struct interf*)sdis_data_cget(data))->temperature;
+}
+
+/*******************************************************************************
+ * Test
+ ******************************************************************************/
+int
+main(int argc, char** argv)
+{
+ struct mem_allocator allocator;
+ struct solid* solid_param = NULL;
+ struct fluid* fluid_param = NULL;
+ struct interf* interf_param = NULL;
+ struct sdis_device* dev = NULL;
+ struct sdis_data* data = NULL;
+ struct sdis_medium* fluid = NULL;
+ struct sdis_medium* solid = NULL;
+ struct sdis_medium* solid1 = NULL;
+ struct sdis_medium* solid2 = NULL;
+ struct sdis_scene* scn = NULL;
+ struct sdis_estimator* estimator = NULL;
+ struct sdis_fluid_shader fluid_shader = SDIS_FLUID_SHADER_NULL;
+ struct sdis_solid_shader solid_shader = SDIS_SOLID_SHADER_NULL;
+ struct sdis_interface_shader interf_shader = SDIS_INTERFACE_SHADER_NULL;
+ struct sdis_interface* interf_solid_adiabatic = NULL;
+ struct sdis_interface* interf_solid1_adiabatic = NULL;
+ struct sdis_interface* interf_solid2_adiabatic = NULL;
+ struct sdis_interface* interf_solid_solid1 = NULL;
+ struct sdis_interface* interf_solid_solid2 = NULL;
+ struct sdis_interface* interf_solid1_fluid = NULL;
+ struct sdis_interface* interf_solid2_fluid = NULL;
+ struct sdis_interface* interfaces[10/*#segment*/];
+ struct sdis_mc T = SDIS_MC_NULL;
+ double Tref;
+ double pos[2];
+ size_t nfails;
+ size_t nreals;
+ (void)argc, (void)argv;
+
+ CHK(mem_init_proxy_allocator(&allocator, &mem_default_allocator) == RES_OK);
+ CHK(sdis_device_create
+ (NULL, &allocator, SDIS_NTHREADS_DEFAULT, 1, &dev) == RES_OK);
+
+ /* Create the fluid medium */
+ fluid_shader.temperature = fluid_get_temperature;
+ fluid_shader.calorific_capacity = dummy_medium_getter;
+ fluid_shader.volumic_mass = dummy_medium_getter;
+ CHK(sdis_data_create
+ (dev, sizeof(struct fluid), ALIGNOF(struct fluid), NULL, &data) == RES_OK);
+ fluid_param = sdis_data_get(data);
+ fluid_param->temperature_upper = T1;
+ fluid_param->temperature_lower = T2;
+ CHK(sdis_fluid_create(dev, &fluid_shader, data, &fluid) == RES_OK);
+ CHK(sdis_data_ref_put(data) == RES_OK);
+
+ /* Setup the solid shader */
+ solid_shader.calorific_capacity = solid_get_calorific_capacity;
+ solid_shader.thermal_conductivity = solid_get_thermal_conductivity;
+ solid_shader.volumic_mass = solid_get_volumic_mass;
+ solid_shader.delta_solid = solid_get_delta;
+ solid_shader.temperature = solid_get_temperature;
+ solid_shader.volumic_power = solid_get_volumic_power;
+
+ /* Create the medium of the upper slab */
+ CHK(sdis_data_create
+ (dev, sizeof(struct solid), ALIGNOF(struct solid), NULL, &data) == RES_OK);
+ solid_param = sdis_data_get(data);
+ solid_param->cp = 500000;
+ solid_param->rho = 1000;
+ solid_param->lambda = LAMBDA1;
+ solid_param->delta = DELTA1;
+ solid_param->volumic_power = SDIS_VOLUMIC_POWER_NONE;
+ solid_param->temperature = -1;
+ CHK(sdis_solid_create(dev, &solid_shader, data, &solid1) == RES_OK);
+ CHK(sdis_data_ref_put(data) == RES_OK);
+
+ /* Create the medium of the lower slab */
+ CHK(sdis_data_create
+ (dev, sizeof(struct solid), ALIGNOF(struct solid), NULL, &data) == RES_OK);
+ solid_param = sdis_data_get(data);
+ solid_param->cp = 500000;
+ solid_param->rho = 1000;
+ solid_param->lambda = LAMBDA2;
+ solid_param->delta = DELTA2;
+ solid_param->volumic_power = SDIS_VOLUMIC_POWER_NONE;
+ solid_param->temperature = -1;
+ CHK(sdis_solid_create(dev, &solid_shader, data, &solid2) == RES_OK);
+ CHK(sdis_data_ref_put(data) == RES_OK);
+
+ /* Create the medium of the middle slab */
+ CHK(sdis_data_create
+ (dev, sizeof(struct solid), ALIGNOF(struct solid), NULL, &data) == RES_OK);
+ solid_param = sdis_data_get(data);
+ solid_param->cp = 500000;
+ solid_param->rho = 1000;
+ solid_param->lambda = LAMBDA;
+ solid_param->delta = DELTA;
+ solid_param->volumic_power = Pw;
+ solid_param->temperature = -1;
+ CHK(sdis_solid_create(dev, &solid_shader, data, &solid) == RES_OK);
+ CHK(sdis_data_ref_put(data) == RES_OK);
+
+ interf_shader.front.temperature = interface_get_temperature;
+
+ /* Create the solid/solid1 interface */
+ CHK(sdis_data_create (dev, sizeof(struct interf), ALIGNOF(struct interf),
+ NULL, &data) == RES_OK);
+ interf_param = sdis_data_get(data);
+ interf_param->temperature = Tsolid_solid1;
+ CHK(sdis_interface_create(dev, solid, solid1, &interf_shader,
+ data, &interf_solid_solid1) == RES_OK);
+ CHK(sdis_data_ref_put(data) == RES_OK);
+
+ /* Create the solid/solid2 interface */
+ CHK(sdis_data_create (dev, sizeof(struct interf), ALIGNOF(struct interf),
+ NULL, &data) == RES_OK);
+ interf_param = sdis_data_get(data);
+ interf_param->temperature = Tsolid_solid2;
+ CHK(sdis_interface_create(dev, solid, solid2, &interf_shader,
+ data, &interf_solid_solid2) == RES_OK);
+ CHK(sdis_data_ref_put(data) == RES_OK);
+
+ /* Setup the interface shader */
+ interf_shader.convection_coef = interface_get_convection_coef;
+ interf_shader.front.temperature = interface_get_temperature;
+
+ /* Create the adiabatic interfaces */
+ CHK(sdis_data_create (dev, sizeof(struct interf), ALIGNOF(struct interf),
+ NULL, &data) == RES_OK);
+ interf_param = sdis_data_get(data);
+ interf_param->h = 0;
+ interf_param->temperature = -1;
+ CHK(sdis_interface_create(dev, solid, fluid, &interf_shader, data,
+ &interf_solid_adiabatic) == RES_OK);
+ CHK(sdis_interface_create(dev, solid1, fluid, &interf_shader, data,
+ &interf_solid1_adiabatic) == RES_OK);
+ CHK(sdis_interface_create(dev, solid2, fluid, &interf_shader, data,
+ &interf_solid2_adiabatic) == RES_OK);
+ CHK(sdis_data_ref_put(data) == RES_OK);
+
+ /* Create the solid1 fluid interace */
+ CHK(sdis_data_create (dev, sizeof(struct interf), ALIGNOF(struct interf),
+ NULL, &data) == RES_OK);
+ interf_param = sdis_data_get(data);
+ interf_param->h = H1;
+ interf_param->temperature = Tsolid1_fluid;
+ CHK(sdis_interface_create(dev, solid1, fluid, &interf_shader, data,
+ &interf_solid1_fluid) == RES_OK);
+ CHK(sdis_data_ref_put(data) == RES_OK);
+
+ /* Create the solid2 fluid interface */
+ CHK(sdis_data_create (dev, sizeof(struct interf), ALIGNOF(struct interf),
+ NULL, &data) == RES_OK);
+ interf_param = sdis_data_get(data);
+ interf_param->h = H2;
+ interf_param->temperature = Tsolid2_fluid;
+ CHK(sdis_interface_create(dev, solid2, fluid, &interf_shader, data,
+ &interf_solid2_fluid) == RES_OK);
+ CHK(sdis_data_ref_put(data) == RES_OK);
+
+ /* Release the media */
+ CHK(sdis_medium_ref_put(fluid) == RES_OK);
+ CHK(sdis_medium_ref_put(solid) == RES_OK);
+ CHK(sdis_medium_ref_put(solid1) == RES_OK);
+ CHK(sdis_medium_ref_put(solid2) == RES_OK);
+
+ /* Map the interfaces to their square segments */
+ interfaces[0] = interf_solid2_adiabatic;
+ interfaces[1] = interf_solid_adiabatic;
+ interfaces[2] = interf_solid1_adiabatic;
+ interfaces[3] = interf_solid1_fluid;
+ interfaces[4] = interf_solid1_adiabatic;
+ interfaces[5] = interf_solid_adiabatic;
+ interfaces[6] = interf_solid2_adiabatic;
+ interfaces[7] = interf_solid2_fluid;
+ interfaces[8] = interf_solid_solid2;
+ interfaces[9] = interf_solid_solid1;
+
+#if 0
+ dump_segments(stdout, vertices, nvertices, indices, nsegments);
+ exit(0);
+#endif
+
+ /* Create the scene */
+ CHK(sdis_scene_2d_create(dev, nsegments, get_indices, get_interface,
+ nvertices, get_position, interfaces, &scn) == RES_OK);
+
+ /* Release the interfaces */
+ CHK(sdis_interface_ref_put(interf_solid_adiabatic) == RES_OK);
+ CHK(sdis_interface_ref_put(interf_solid1_adiabatic) == RES_OK);
+ CHK(sdis_interface_ref_put(interf_solid2_adiabatic) == RES_OK);
+ CHK(sdis_interface_ref_put(interf_solid_solid1) == RES_OK);
+ CHK(sdis_interface_ref_put(interf_solid_solid2) == RES_OK);
+ CHK(sdis_interface_ref_put(interf_solid1_fluid) == RES_OK);
+ CHK(sdis_interface_ref_put(interf_solid2_fluid) == RES_OK);
+
+ pos[0] = 0;
+ pos[1] = PROBE_POS;
+
+ if(pos[1] > 0 && pos[1] < L2) { /* Lower slab */
+ Tref = Tp2 + (Tb - Tp2) * pos[1] / L2;
+ } else if(pos[1] > L2 && pos[1] < L2 + L) { /* Middle slab */
+ Tref =
+ (Ta + Tb) / 2
+ + (Ta - Tb)/L * (pos[1] - (L2+L/2))
+ + Pw * (L*L/4.0 - pow((pos[1] - (L2+L/2)), 2)) / (2*LAMBDA);
+ } else if(pos[1] > L2 + L && pos[1] < L2 + L1 + L) {
+ Tref = Ta + (Tp1 - Ta) / L1 * (pos[1] - (L+L2));
+ } else {
+ FATAL("Unreachable code.\n");
+ }
+
+ CHK(sdis_solve_probe(scn, N, pos, INF, 1.f, -1, 0, &estimator) == RES_OK);
+ CHK(sdis_estimator_get_temperature(estimator, &T) == RES_OK);
+ CHK(sdis_estimator_get_failure_count(estimator, &nfails) == RES_OK);
+ CHK(sdis_estimator_get_realisation_count(estimator, &nreals) == RES_OK);
+ printf("Temperature at (%g %g) = %g ~ %g +/- %g [%g, %g]\n",
+ SPLIT2(pos), Tref, T.E, T.SE, T.E-3*T.SE, T.E+3*T.SE);
+ printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N);
+ CHK(sdis_estimator_ref_put(estimator) == RES_OK);
+
+ CHK(nfails + nreals == N);
+ CHK(nfails < N/1000);
+ CHK(eq_eps(T.E, Tref, T.SE*3));
+
+ CHK(sdis_scene_ref_put(scn) == RES_OK);
+ CHK(sdis_device_ref_put(dev) == RES_OK);
+
+ check_memory_allocator(&allocator);
+ mem_shutdown_proxy_allocator(&allocator);
+ CHK(mem_allocated_size() == 0);
+ return 0;
+}
+
diff --git a/src/test_sdis_volumic_power4_2d.c b/src/test_sdis_volumic_power4_2d.c
@@ -0,0 +1,371 @@
+/* Copyright (C) 2016-2018 |Meso|Star> (contact@meso-star.com)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. */
+
+#include "sdis.h"
+#include "test_sdis_utils.h"
+#include <rsys/math.h>
+
+#define Tf1 0
+#define Tf2 100
+#define Power 0 /*10000*/
+#define H1 50
+#define H2 50
+#define LAMBDA 100.0
+#define DELTA (1.0/20.0)
+#define N 10000
+
+/*
+ * The 2D scene is a solid slabs stretched along the X dimension to simulate a
+ * 1D case. The slab has a volumic power and has a convective exchange with the
+ * surrounding fluid whose temperature is fixed to Tfluid.
+ *
+ *
+ * _\ TFluid
+ * / /
+ * \__/
+ *
+ * ... -----Hboundary----- ...
+ *
+ * Lambda, Power
+ *
+ * ... -----Hboundary----- ...
+ *
+ * _\ TFluid
+ * / /
+ * \__/
+ *
+ */
+
+static const double vertices[4/*#vertices*/*2/*#coords per vertex*/] = {
+ -10000.5,-0.5,
+ -10000.5, 0.5,
+ 10000.5, 0.5,
+ 10000.5,-0.5
+};
+static const size_t nvertices = sizeof(vertices)/sizeof(double[2]);
+
+static const size_t indices[4/*#segments*/*2/*#indices per segment*/]= {
+ 0, 1,
+ 1, 2,
+ 2, 3,
+ 3, 0
+};
+static const size_t nsegments = sizeof(indices)/sizeof(size_t[2]);
+
+/*******************************************************************************
+ * Geometry
+ ******************************************************************************/
+static void
+get_indices(const size_t iseg, size_t ids[2], void* context)
+{
+ (void)context;
+ CHK(ids);
+ ids[0] = indices[iseg*2+0];
+ ids[1] = indices[iseg*2+1];
+}
+
+static void
+get_position(const size_t ivert, double pos[2], void* context)
+{
+ (void)context;
+ CHK(pos);
+ pos[0] = vertices[ivert*2+0];
+ pos[1] = vertices[ivert*2+1];
+}
+
+static void
+get_interface(const size_t iseg, struct sdis_interface** bound, void* context)
+{
+ struct sdis_interface** interfaces = context;
+ CHK(context && bound);
+ *bound = interfaces[iseg];
+}
+
+/*******************************************************************************
+ * Solid medium
+ ******************************************************************************/
+struct solid {
+ double cp;
+ double lambda;
+ double rho;
+ double delta;
+ double volumic_power;
+ double temperature;
+};
+
+static double
+solid_get_calorific_capacity
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ CHK(data != NULL && vtx != NULL);
+ return ((const struct solid*)sdis_data_cget(data))->cp;
+}
+
+static double
+solid_get_thermal_conductivity
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ CHK(data != NULL && vtx != NULL);
+ return ((const struct solid*)sdis_data_cget(data))->lambda;
+}
+
+static double
+solid_get_volumic_mass
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ CHK(data != NULL && vtx != NULL);
+ return ((const struct solid*)sdis_data_cget(data))->rho;
+}
+
+static double
+solid_get_delta
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ CHK(data != NULL && vtx != NULL);
+ return ((const struct solid*)sdis_data_cget(data))->delta;
+}
+
+static double
+solid_get_temperature
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ CHK(data != NULL && vtx != NULL);
+ return ((const struct solid*)sdis_data_cget(data))->temperature;
+}
+
+static double
+solid_get_volumic_power
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ CHK(data != NULL && vtx != NULL);
+ return ((const struct solid*)sdis_data_cget(data))->volumic_power;
+}
+
+/*******************************************************************************
+ * Fluid medium
+ ******************************************************************************/
+struct fluid {
+ double temperature;
+};
+
+static double
+fluid_get_temperature
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ const struct fluid* fluid;
+ CHK(data != NULL && vtx != NULL);
+ fluid = sdis_data_cget(data);
+ return fluid->temperature;
+}
+
+/*******************************************************************************
+ * Interfaces
+ ******************************************************************************/
+struct interf {
+ double h;
+ double temperature;
+};
+
+static double
+interface_get_convection_coef
+ (const struct sdis_interface_fragment* frag, struct sdis_data* data)
+{
+ CHK(frag && data);
+ return ((const struct interf*)sdis_data_cget(data))->h;
+}
+
+static double
+interface_get_temperature
+ (const struct sdis_interface_fragment* frag, struct sdis_data* data)
+{
+ CHK(frag && data);
+ return ((const struct interf*)sdis_data_cget(data))->temperature;
+}
+
+/*******************************************************************************
+ * Test
+ ******************************************************************************/
+int
+main(int argc, char** argv)
+{
+ struct mem_allocator allocator;
+ struct solid* solid_param = NULL;
+ struct fluid* fluid_param = NULL;
+ struct interf* interf_param = NULL;
+ struct sdis_device* dev = NULL;
+ struct sdis_data* data = NULL;
+ struct sdis_medium* fluid1 = NULL;
+ struct sdis_medium* fluid2 = NULL;
+ struct sdis_medium* solid = NULL;
+ struct sdis_scene* scn = NULL;
+ struct sdis_estimator* estimator = NULL;
+ struct sdis_fluid_shader fluid_shader = SDIS_FLUID_SHADER_NULL;
+ struct sdis_solid_shader solid_shader = SDIS_SOLID_SHADER_NULL;
+ struct sdis_interface_shader interf_shader = SDIS_INTERFACE_SHADER_NULL;
+ struct sdis_interface* interf_adiabatic = NULL;
+ struct sdis_interface* interf_solid_fluid1 = NULL;
+ struct sdis_interface* interf_solid_fluid2 = NULL;
+ struct sdis_interface* interfaces[4/*#segment*/];
+ struct sdis_mc T = SDIS_MC_NULL;
+ size_t nreals, nfails;
+ double pos[2];
+ double Tref;
+ double a, b, x;
+ double L;
+ (void)argc, (void)argv;
+
+ CHK(mem_init_proxy_allocator(&allocator, &mem_default_allocator) == RES_OK);
+ CHK(sdis_device_create
+ (NULL, &allocator, SDIS_NTHREADS_DEFAULT, 1, &dev) == RES_OK);
+
+ /* Create the fluid medium */
+ fluid_shader.temperature = fluid_get_temperature;
+ fluid_shader.calorific_capacity = dummy_medium_getter;
+ fluid_shader.volumic_mass = dummy_medium_getter;
+
+ CHK(sdis_data_create
+ (dev, sizeof(struct fluid), ALIGNOF(struct fluid), NULL, &data) == RES_OK);
+ fluid_param = sdis_data_get(data);
+ fluid_param->temperature = Tf1;
+ CHK(sdis_fluid_create(dev, &fluid_shader, data, &fluid1) == RES_OK);
+ CHK(sdis_data_ref_put(data) == RES_OK);
+
+ CHK(sdis_data_create
+ (dev, sizeof(struct fluid), ALIGNOF(struct fluid), NULL, &data) == RES_OK);
+ fluid_param = sdis_data_get(data);
+ fluid_param->temperature = Tf2;
+ CHK(sdis_fluid_create(dev, &fluid_shader, data, &fluid2) == RES_OK);
+ CHK(sdis_data_ref_put(data) == RES_OK);
+
+ /* Setup the solid shader */
+ solid_shader.calorific_capacity = solid_get_calorific_capacity;
+ solid_shader.thermal_conductivity = solid_get_thermal_conductivity;
+ solid_shader.volumic_mass = solid_get_volumic_mass;
+ solid_shader.delta_solid = solid_get_delta;
+ solid_shader.temperature = solid_get_temperature;
+ solid_shader.volumic_power = solid_get_volumic_power;
+
+ /* Create the solid medium */
+ CHK(sdis_data_create
+ (dev, sizeof(struct solid), ALIGNOF(struct solid), NULL, &data) == RES_OK);
+ solid_param = sdis_data_get(data);
+ solid_param->cp = 500000;
+ solid_param->rho = 1000;
+ solid_param->lambda = LAMBDA;
+ solid_param->delta = DELTA;
+ solid_param->volumic_power = Power;
+ solid_param->temperature = -1;
+ CHK(sdis_solid_create(dev, &solid_shader, data, &solid) == RES_OK);
+ CHK(sdis_data_ref_put(data) == RES_OK);
+
+ /* Setup the interface shader */
+ interf_shader.convection_coef = interface_get_convection_coef;
+ interf_shader.front.temperature = interface_get_temperature;
+
+ /* Create the adiabatic interface */
+ CHK(sdis_data_create (dev, sizeof(struct interf), ALIGNOF(struct interf),
+ NULL, &data) == RES_OK);
+ interf_param = sdis_data_get(data);
+ interf_param->h = 0;
+ interf_param->temperature = -1;
+ CHK(sdis_interface_create(dev, solid, fluid1, &interf_shader, data,
+ &interf_adiabatic) == RES_OK);
+ CHK(sdis_data_ref_put(data) == RES_OK);
+
+ /* Create the solid fluid1 interface */
+ CHK(sdis_data_create (dev, sizeof(struct interf), ALIGNOF(struct interf),
+ NULL, &data) == RES_OK);
+ interf_param = sdis_data_get(data);
+ interf_param->h = H1;
+ interf_param->temperature = -1;
+ CHK(sdis_interface_create(dev, solid, fluid1, &interf_shader, data,
+ &interf_solid_fluid1) == RES_OK);
+ CHK(sdis_data_ref_put(data) == RES_OK);
+
+ /* Create the solid fluid2 interface */
+ CHK(sdis_data_create (dev, sizeof(struct interf), ALIGNOF(struct interf),
+ NULL, &data) == RES_OK);
+ interf_param = sdis_data_get(data);
+ interf_param->h = H2;
+ interf_param->temperature = -1;
+ CHK(sdis_interface_create(dev, solid, fluid2, &interf_shader, data,
+ &interf_solid_fluid2) == RES_OK);
+ CHK(sdis_data_ref_put(data) == RES_OK);
+
+ /* Release the media */
+ CHK(sdis_medium_ref_put(fluid1) == RES_OK);
+ CHK(sdis_medium_ref_put(fluid2) == RES_OK);
+ CHK(sdis_medium_ref_put(solid) == RES_OK);
+
+ /* Map the interfaces to their square segments */
+ interfaces[0] = interf_adiabatic;
+ interfaces[1] = interf_solid_fluid1;
+ interfaces[2] = interf_adiabatic;
+ interfaces[3] = interf_solid_fluid2;
+
+#if 0
+ dump_segments(stdout, vertices, nvertices, indices, nsegments);
+ exit(0);
+#endif
+
+ /* Create the scene */
+ CHK(sdis_scene_2d_create(dev, nsegments, get_indices, get_interface,
+ nvertices, get_position, interfaces, &scn) == RES_OK);
+
+ /* Release the interfaces */
+ CHK(sdis_interface_ref_put(interf_adiabatic) == RES_OK);
+ CHK(sdis_interface_ref_put(interf_solid_fluid1) == RES_OK);
+ CHK(sdis_interface_ref_put(interf_solid_fluid2) == RES_OK);
+
+ pos[0] = 0;
+ pos[1] = 0.25;
+
+ L = vertices[3] - vertices[1];
+#if 1
+ x = pos[1] + vertices[3];
+ a = (H2*Power*L + H1*H2*(Tf1 - Tf2) + H1*H2*Power*L*L/(2*LAMBDA))
+ / (LAMBDA * (H1 + H2) + H1*H2*L);
+ b = Tf2 + a * LAMBDA / H2;
+ Tref = -Power / (2*LAMBDA) * x*x + a * x + b;
+#else
+ tmp = LAMBDA / L;
+ T1 = H1 * (H2+tmp) / (tmp*(H1+H2) + H1*H2) * Tf1
+ + H2 * tmp / (tmp*(H1+H2) + H1*H2) * Tf2;
+ T2 = H1 * tmp / (tmp*(H1+H2) + H1*H2) * Tf1
+ + H2 * (H1+tmp) / (tmp*(H1+H2) + H1*H2) * Tf2;
+ Tref = T2 + (T1-T2)/L * (pos[1] + vertices[3]);
+#endif
+
+ CHK(sdis_solve_probe(scn, N, pos, INF, 1.f, -1, 0, &estimator) == RES_OK);
+ CHK(sdis_estimator_get_temperature(estimator, &T) == RES_OK);
+ CHK(sdis_estimator_get_realisation_count(estimator, &nreals) == RES_OK);
+ CHK(sdis_estimator_get_failure_count(estimator, &nfails) == RES_OK);
+ printf("Temperature at (%g %g) = %g ~ %g +/- %g [%g %g]\n",
+ SPLIT2(pos), Tref, T.E, T.SE, T.E-3*T.SE, T.E+3*T.SE);
+ printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N);
+ CHK(sdis_estimator_ref_put(estimator) == RES_OK);
+ CHK(nfails + nreals == N);
+ CHK(nfails < N/1000);
+ CHK(eq_eps(T.E, Tref, T.SE*3));
+
+ CHK(sdis_scene_ref_put(scn) == RES_OK);
+ CHK(sdis_device_ref_put(dev) == RES_OK);
+
+ check_memory_allocator(&allocator);
+ mem_shutdown_proxy_allocator(&allocator);
+ CHK(mem_allocated_size() == 0);
+ return 0;
+}
+
