commit ef013176835b8e33141e87b7cb4c437454629824
parent db52f0b369426b18d4335b8a5c80c34ab7fb0704
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Sat, 31 Mar 2018 12:05:53 +0200
Merge branch 'release_0.1'
Diffstat:
30 files changed, 3548 insertions(+), 758 deletions(-)
diff --git a/.gitignore b/.gitignore
@@ -1,3 +1,4 @@
+compile_commands.json
.gitignore
CMakeCache.txt
CMakeFiles
diff --git a/README.md b/README.md
@@ -1,7 +1,7 @@
# Stardis
The purpose of this library is to solve coupled convecto - conducto - radiative
-thermal problems.
+thermal problems in 2D and 3D environments.
## How to build
@@ -23,6 +23,18 @@ variable the install directories of its dependencies.
## Release notes
+### 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
+ of MC weights and square weights.
+- 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
+ `sdis_solve_probe` function.
+
### Version 0.0
First version and implementation of the Stardis solver API.
diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt
@@ -40,28 +40,32 @@ rcmake_append_runtime_dirs(_runtime_dirs RSys Star3D StarSP)
# Configure and define targets
################################################################################
set(VERSION_MAJOR 0)
-set(VERSION_MINOR 0)
+set(VERSION_MINOR 1)
set(VERSION_PATCH 0)
set(VERSION ${VERSION_MAJOR}.${VERSION_MINOR}.${VERSION_PATCH})
set(SDIS_FILES_SRC
+ sdis_accum_buffer.c
+ sdis_camera.c
sdis_data.c
sdis_device.c
sdis_estimator.c
sdis_interface.c
sdis_medium.c
sdis_scene.c
- sdis_solve_probe.c)
+ sdis_solve.c)
set(SDIS_FILES_INC_API
sdis.h)
set(SDIS_FILES_INC
+ sdis_camera.h
sdis_device_c.h
sdis_estimator_c.h
sdis_interface_c.h
sdis_medium_c.h
- sdis_scene_c.h)
+ sdis_scene_c.h
+ sdis_solve_Xd.h)
set(SDIS_FILES_DOC COPYING README.md)
@@ -75,7 +79,10 @@ add_library(sdis SHARED
${SDIS_FILES_SRC}
${SDIS_FILES_INC}
${SDIS_FILES_INC_API})
-target_link_libraries(sdis RSys Star2D Star3D StarSP)
+target_link_libraries(sdis RSys Star2D Star3D StarSP m)
+if(CMAKE_COMPILER_IS_GNUCC)
+ target_link_libraries(sdis m)
+endif()
set_target_properties(sdis PROPERTIES
DEFINE_SYMBOL SDIS_SHARED_BUILD
@@ -111,11 +118,16 @@ if(NOT NO_TEST)
register_test(${_name} ${_name})
endfunction()
+ new_test(test_sdis_accum_buffer)
+ new_test(test_sdis_camera)
+ new_test(test_sdis_conducto_radiative)
+ new_test(test_sdis_conducto_radiative_2d)
new_test(test_sdis_data)
new_test(test_sdis_device)
new_test(test_sdis_interface)
new_test(test_sdis_medium)
new_test(test_sdis_scene)
+ new_test(test_sdis_solve_camera)
new_test(test_sdis_solve_probe)
new_test(test_sdis_solve_probe2)
new_test(test_sdis_solve_probe3)
@@ -124,6 +136,8 @@ if(NOT NO_TEST)
new_test(test_sdis_solve_probe3_2d)
target_link_libraries(test_sdis_solve_probe3 Star3DUT)
+ target_link_libraries(test_sdis_solve_camera Star3DUT)
+
if(CMAKE_COMPILER_IS_GNUCC)
target_link_libraries(test_sdis_solve_probe3_2d m)
endif()
diff --git a/src/sdis.h b/src/sdis.h
@@ -51,6 +51,8 @@ struct mem_allocator;
* a reference on the data, i.e. they increment or decrement the reference
* counter, respectively. When this counter reaches 0, the object is silently
* destroyed and cannot be used anymore. */
+struct sdis_accum_buffer;
+struct sdis_camera;
struct sdis_data;
struct sdis_device;
struct sdis_estimator;
@@ -91,6 +93,13 @@ struct sdis_interface_fragment {
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_sqr; /* Sum of Monte-Carlo square weights */
+ size_t nweights; /* #accumulated weights */
+};
+
/* Monte-Carlo estimation */
struct sdis_mc {
double E; /* Expected value */
@@ -124,7 +133,7 @@ struct sdis_solid_shader {
* unknown for the submitted random walk vertex. */
sdis_medium_getter_T temperature;
};
-#define SDIS_SOLID_SHADER_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__;
@@ -137,18 +146,38 @@ struct sdis_fluid_shader {
* unknown for the submitted position and time. */
sdis_medium_getter_T temperature;
};
-#define SDIS_FLUID_SHADER_NULL__ {NULL}
+#define SDIS_FLUID_SHADER_NULL__ {NULL, NULL, NULL}
static const struct sdis_fluid_shader SDIS_FLUID_SHADER_NULL =
SDIS_FLUID_SHADER_NULL__;
struct sdis_interface_shader {
sdis_interface_getter_T temperature; /* Limit condition. NULL <=> Unknown */
- sdis_interface_getter_T convection_coef; /* NULL <=> Solid/Solid interface */
+ sdis_interface_getter_T convection_coef; /* May be NULL for solid/solid */
+
+ /* 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] */
};
-#define SDIS_INTERFACE_SHADER_NULL__ {NULL}
+#define SDIS_INTERFACE_SHADER_NULL__ {NULL, NULL, NULL, NULL}
static const struct sdis_interface_shader SDIS_INTERFACE_SHADER_NULL =
SDIS_INTERFACE_SHADER_NULL__;
+struct sdis_accum_buffer_layout {
+ size_t width;
+ size_t height;
+};
+#define SDIS_ACCUM_BUFFER_LAYOUT_NULL__ {0, 0}
+static const struct sdis_accum_buffer_layout SDIS_ACCUM_BUFFER_LAYOUT_NULL =
+ SDIS_ACCUM_BUFFER_LAYOUT_NULL__;
+
+/* Functor use to write accumulations performed by sdis_solve_camera */
+typedef res_T
+(*sdis_write_accums_T)
+ (void* context, /* User data */
+ const size_t origin[2], /* Coordinates of the 1st accumulation in image plane */
+ const size_t naccums[2], /* #accumulations in X and Y */
+ const struct sdis_accum* accums); /* List of row ordered accumulations */
+
BEGIN_DECLS
/*******************************************************************************
@@ -200,6 +229,80 @@ sdis_data_cget
(const struct sdis_data* data);
/*******************************************************************************
+ * A camera describes a point of view
+ ******************************************************************************/
+SDIS_API res_T
+sdis_camera_create
+ (struct sdis_device* dev,
+ struct sdis_camera** cam);
+
+SDIS_API res_T
+sdis_camera_ref_get
+ (struct sdis_camera* cam);
+
+SDIS_API res_T
+sdis_camera_ref_put
+ (struct sdis_camera* cam);
+
+/* Width/height projection ratio */
+SDIS_API res_T
+sdis_camera_set_proj_ratio
+ (struct sdis_camera* cam,
+ const double proj_ratio);
+
+SDIS_API res_T
+sdis_camera_set_fov /* Horizontal field of view */
+ (struct sdis_camera* cam,
+ const double fov); /* In radian */
+
+SDIS_API res_T
+sdis_camera_look_at
+ (struct sdis_camera* cam,
+ const double position[3],
+ const double target[3],
+ const double up[3]);
+
+/*******************************************************************************
+ * A buffer of accumulations
+ ******************************************************************************/
+SDIS_API res_T
+sdis_accum_buffer_create
+ (struct sdis_device* dev,
+ const size_t width,
+ const size_t height,
+ struct sdis_accum_buffer** buf);
+
+SDIS_API res_T
+sdis_accum_buffer_ref_get
+ (struct sdis_accum_buffer* buf);
+
+SDIS_API res_T
+sdis_accum_buffer_ref_put
+ (struct sdis_accum_buffer* buf);
+
+SDIS_API res_T
+sdis_accum_buffer_get_layout
+ (const struct sdis_accum_buffer* buf,
+ struct sdis_accum_buffer_layout* layout);
+
+SDIS_API res_T
+sdis_accum_buffer_map
+ (const struct sdis_accum_buffer* buf,
+ const struct sdis_accum** accums);
+
+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 */
+SDIS_API res_T
+sdis_accum_buffer_write
+ (void* buf, /* User data */
+ const size_t origin[2], /* Coordinates of the 1st accum in image plane */
+ const size_t naccum[2], /* #accum in X and Y */
+ const struct sdis_accum* accums); /* List of row ordered accum */
+
+/*******************************************************************************
* A medium encapsulates the properties of either a fluid or a solid.
******************************************************************************/
SDIS_API res_T
@@ -327,12 +430,28 @@ sdis_estimator_get_temperature
SDIS_API res_T
sdis_solve_probe
(struct sdis_scene* scn,
- const size_t nrealisations,
- const double position[3],
- const double time,
- const double fp_to_meter,/* Scale from floating point units to meters */
+ const size_t nrealisations, /* #realisations */
+ const double position[3], /* Probe position */
+ const double time, /* Observation time */
+ 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 */
struct sdis_estimator** estimator);
+SDIS_API res_T
+sdis_solve_camera
+ (struct sdis_scene* scn,
+ const struct sdis_camera* cam, /* Point of view */
+ const double time, /* Observation time */
+ 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 */
+ const size_t width, /* Image definition in in X */
+ const size_t height, /* Image definition in Y */
+ const size_t spp, /* #samples per pixel */
+ sdis_write_accums_T writer,
+ void* writer_data);
+
END_DECLS
#endif /* SDIS_H */
diff --git a/src/sdis_accum_buffer.c b/src/sdis_accum_buffer.c
@@ -0,0 +1,160 @@
+/* Copyright (C) |Meso|Star> 2016-2018 (contact@meso-star.com)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. */
+
+#include "sdis.h"
+#include "sdis_device_c.h"
+
+struct sdis_accum_buffer {
+ struct sdis_accum* accums;
+ size_t width;
+ size_t height;
+
+ ref_T ref;
+ struct sdis_device* dev;
+};
+
+/*******************************************************************************
+ * Helper functions
+ ******************************************************************************/
+static void
+accum_buffer_release(ref_T* ref)
+{
+ struct sdis_accum_buffer* buf;
+ struct sdis_device* dev;
+ ASSERT(ref);
+ buf = CONTAINER_OF(ref, struct sdis_accum_buffer, ref);
+ dev = buf->dev;
+ if(buf->accums) MEM_RM(dev->allocator, buf->accums);
+ MEM_RM(dev->allocator, buf);
+ SDIS(device_ref_put(dev));
+}
+
+/*******************************************************************************
+ * Export functions
+ ******************************************************************************/
+res_T
+sdis_accum_buffer_create
+ (struct sdis_device* dev,
+ const size_t width,
+ const size_t height,
+ struct sdis_accum_buffer** out_buf)
+{
+ struct sdis_accum_buffer* buf = NULL;
+ size_t sz;
+ res_T res = RES_OK;
+
+ if(!dev || !width || !height || !out_buf) {
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ buf = MEM_CALLOC(dev->allocator, 1, sizeof(struct sdis_accum_buffer));
+ if(!buf) {
+ res = RES_MEM_ERR;
+ goto error;
+ }
+ ref_init(&buf->ref);
+ SDIS(device_ref_get(dev));
+ buf->dev = dev;
+ buf->width = width;
+ buf->height = height;
+
+ sz = buf->width * buf->height * sizeof(struct sdis_accum);
+ buf->accums = MEM_ALLOC_ALIGNED(dev->allocator, sz, 16);
+ if(!buf->accums) {
+ log_err(dev, "%s: could not allocate a buffer accumulations of %lux%lu.\n",
+ FUNC_NAME, (unsigned long)width, (unsigned long)height);
+ res = RES_MEM_ERR;
+ goto error;
+ }
+
+exit:
+ if(out_buf) *out_buf = buf;
+ return res;
+error:
+ if(buf) {
+ SDIS(accum_buffer_ref_put(buf));
+ buf = NULL;
+ }
+ goto exit;
+}
+
+res_T
+sdis_accum_buffer_ref_get(struct sdis_accum_buffer* buf)
+{
+ if(!buf) return RES_BAD_ARG;
+ ref_get(&buf->ref);
+ return RES_OK;
+}
+
+res_T
+sdis_accum_buffer_ref_put(struct sdis_accum_buffer* buf)
+{
+ if(!buf) return RES_BAD_ARG;
+ ref_put(&buf->ref, accum_buffer_release);
+ return RES_OK;
+}
+
+res_T
+sdis_accum_buffer_get_layout
+ (const struct sdis_accum_buffer* buf,
+ struct sdis_accum_buffer_layout* layout)
+{
+ if(!buf || !layout) return RES_BAD_ARG;
+ layout->width = buf->width;
+ layout->height = buf->height;
+ return RES_OK;
+}
+
+res_T
+sdis_accum_buffer_map
+ (const struct sdis_accum_buffer* buf, const struct sdis_accum** accums)
+{
+ if(!buf || !accums) return RES_BAD_ARG;
+ *accums = buf->accums;
+ return RES_OK;
+}
+
+res_T
+sdis_accum_buffer_unmap(const struct sdis_accum_buffer* buf)
+{
+ if(!buf) return RES_BAD_ARG;
+ /* Do nothing */
+ return RES_OK;
+}
+
+res_T
+sdis_accum_buffer_write
+ (void* buffer,
+ const size_t origin[2],
+ const size_t size[2],
+ const struct sdis_accum* accums)
+{
+ struct sdis_accum_buffer* buf = buffer;
+ const struct sdis_accum* src_row = accums;
+ size_t dst_iy;
+
+ if(UNLIKELY(!buffer || !origin || !size || !accums)) return RES_BAD_ARG;
+ if(UNLIKELY((origin[0] + size[0]) > buf->width)) return RES_BAD_ARG;
+ if(UNLIKELY((origin[1] + size[1]) > buf->height)) return RES_BAD_ARG;
+
+ FOR_EACH(dst_iy, origin[1], origin[1] + size[1]) {
+ const size_t dst_irow = dst_iy*buf->width + origin[0];
+ memcpy(buf->accums + dst_irow, src_row, size[0] * sizeof(struct sdis_accum));
+ src_row += size[0];
+ }
+ return RES_OK;
+}
+
diff --git a/src/sdis_camera.c b/src/sdis_camera.c
@@ -0,0 +1,138 @@
+/* Copyright (C) |Meso|Star> 2016-2018 (contact@meso-star.com)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. */
+
+#include "sdis.h"
+#include "sdis_camera.h"
+#include "sdis_device_c.h"
+
+#include <rsys/mem_allocator.h>
+
+/*******************************************************************************
+ * Helper functions
+ ******************************************************************************/
+static void
+camera_release(ref_T* ref)
+{
+ struct sdis_camera* cam = CONTAINER_OF(ref, struct sdis_camera, ref);
+ struct sdis_device* dev;
+ ASSERT(ref);
+ dev = cam->dev;
+ MEM_RM(dev->allocator, cam);
+ SDIS(device_ref_put(dev));
+}
+
+/*******************************************************************************
+ * Exported functions
+ ******************************************************************************/
+res_T
+sdis_camera_create(struct sdis_device* dev, struct sdis_camera** out_cam)
+{
+ const double pos[3] = {0, 0, 0};
+ const double tgt[3] = {0, 0,-1};
+ const double up[3] = {0, 1, 0};
+ struct sdis_camera* cam = NULL;
+ res_T res = RES_OK;
+
+ if(!dev || !out_cam) {
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ cam = MEM_CALLOC(dev->allocator, 1, sizeof(struct sdis_camera));
+ if(!cam) {
+ res = RES_MEM_ERR;
+ goto error;
+ }
+ SDIS(device_ref_get(dev));
+ cam->dev = dev;
+ ref_init(&cam->ref);
+ cam->rcp_proj_ratio = 1.0;
+ cam->fov_x = PI/2.0;
+ SDIS(camera_look_at(cam, pos, tgt, up));
+
+exit:
+ if(out_cam) *out_cam = cam;
+ return res;
+error:
+ if(cam) {
+ SDIS(camera_ref_put(cam));
+ cam = NULL;
+ }
+ goto exit;
+}
+
+res_T
+sdis_camera_ref_get(struct sdis_camera* cam)
+{
+ if(!cam) return RES_BAD_ARG;
+ ref_get(&cam->ref);
+ return RES_OK;
+}
+
+res_T
+sdis_camera_ref_put(struct sdis_camera* cam)
+{
+ if(!cam) return RES_BAD_ARG;
+ ref_put(&cam->ref, camera_release);
+ return RES_OK;
+}
+
+res_T
+sdis_camera_set_proj_ratio(struct sdis_camera* cam, const double ratio)
+{
+ double y[3] = {0};
+ if(!cam || ratio <= 0) return RES_BAD_ARG;
+ if(d3_normalize(y, cam->axis_y) <= 0) return RES_BAD_ARG;
+ cam->rcp_proj_ratio = 1.0 / ratio;
+ d3_muld(cam->axis_y, y, cam->rcp_proj_ratio);
+ return RES_OK;
+}
+
+res_T
+sdis_camera_set_fov(struct sdis_camera* cam, const double fov_x)
+{
+ double z[3] = {0};
+ double img_plane_depth;
+ if(!cam || (float)fov_x <= 0) return RES_BAD_ARG;
+ if(d3_normalize(z, cam->axis_z) <= 0) return RES_BAD_ARG;
+ img_plane_depth = 1.0/tan(fov_x*0.5);
+ d3_muld(cam->axis_z, z, img_plane_depth);
+ cam->fov_x = fov_x;
+ return RES_OK;
+}
+
+res_T
+sdis_camera_look_at
+ (struct sdis_camera* cam,
+ const double pos[3],
+ const double tgt[3],
+ const double up[3])
+{
+ double x[3], y[3], z[3];
+ double img_plane_depth;
+ if(!cam || !pos || !tgt || !up) return RES_BAD_ARG;
+
+ if(d3_normalize(z, d3_sub(z, tgt, pos)) <= 0) return RES_BAD_ARG;
+ if(d3_normalize(x, d3_cross(x, z, up)) <= 0) return RES_BAD_ARG;
+ if(d3_normalize(y, d3_cross(y, z, x)) <= 0) return RES_BAD_ARG;
+ img_plane_depth = 1.0/tan(cam->fov_x*0.5);
+
+ d3_set(cam->axis_x, x);
+ d3_muld(cam->axis_y, y, cam->rcp_proj_ratio);
+ d3_muld(cam->axis_z, z, img_plane_depth);
+ d3_set(cam->position, pos);
+ return RES_OK;
+}
+
diff --git a/src/sdis_camera.h b/src/sdis_camera.h
@@ -0,0 +1,60 @@
+/* Copyright (C) |Meso|Star> 2016-2018 (contact@meso-star.com)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. */
+
+#ifndef SDIS_CAMERA_H
+#define SDIS_CAMERA_H
+
+#include <rsys/double3.h>
+#include <rsys/ref_count.h>
+
+struct sdis_device;
+
+struct sdis_camera {
+ /* Orthogonal camera frame */
+ double axis_x[3];
+ double axis_y[3];
+ double axis_z[3];
+
+ double position[3];
+ double fov_x; /* Field of view in radians */
+ double rcp_proj_ratio; /* height / width */
+
+ ref_T ref;
+ struct sdis_device* dev;
+};
+
+static FINLINE void
+camera_ray
+ (const struct sdis_camera* cam,
+ const double sample[2], /* In [0, 1[ */
+ double org[3],
+ double dir[3])
+{
+ double x[3], y[3], len;
+ (void)len; /* Avoid warning in debug */
+ ASSERT(cam && sample && org && dir);
+ ASSERT(sample[0] >= 0.0 || sample[0] < 1.0);
+ ASSERT(sample[1] >= 0.0 || sample[1] < 1.0);
+
+ d3_muld(x, cam->axis_x, sample[0]*2.0 - 1.0);
+ d3_muld(y, cam->axis_y, sample[1]*2.0 - 1.0);
+ d3_add(dir, d3_add(dir, x, y), cam->axis_z);
+ len = d3_normalize(dir, dir);
+ ASSERT(len >= 1.e-6);
+ d3_set(org, cam->position);
+}
+
+#endif /* SDIS_CAMERA_H */
+
diff --git a/src/sdis_device.c b/src/sdis_device.c
@@ -52,6 +52,7 @@ device_release(ref_T* ref)
if(dev->s3d) S3D(device_ref_put(dev->s3d));
ASSERT(flist_name_is_empty(&dev->names));
flist_name_release(&dev->names);
+ darray_tile_release(&dev->tiles);
MEM_RM(dev->allocator, dev);
}
@@ -94,11 +95,19 @@ sdis_device_create
dev->nthreads = MMIN(nthreads_hint, (unsigned)omp_get_num_procs());
ref_init(&dev->ref);
flist_name_init(allocator, &dev->names);
+ darray_tile_init(allocator, &dev->tiles);
+
+ res = darray_tile_resize(&dev->tiles, dev->nthreads);
+ if(res != RES_OK) {
+ log_err(dev,
+ "%s: could not allocate the per thread buffer of estimations.\n", FUNC_NAME);
+ goto error;
+ }
res = s2d_device_create(log, allocator, 0, &dev->s2d);
if(res != RES_OK) {
log_err(dev,
- "%s, could not create the Star-2D device on Stardis.\n", FUNC_NAME);
+ "%s: could not create the Star-2D device on Stardis.\n", FUNC_NAME);
}
res = s3d_device_create(log, allocator, 0, &dev->s3d);
diff --git a/src/sdis_device_c.h b/src/sdis_device_c.h
@@ -16,6 +16,7 @@
#ifndef SDIS_DEVICE_C_H
#define SDIS_DEVICE_C_H
+#include <rsys/dynamic_array.h>
#include <rsys/free_list.h>
#include <rsys/ref_count.h>
@@ -23,6 +24,18 @@ struct name { FITEM; };
#define FITEM_TYPE name
#include <rsys/free_list.h>
+#define DARRAY_NAME accum
+#define DARRAY_DATA struct sdis_accum
+#include <rsys/dynamic_array.h>
+
+#define DARRAY_NAME tile
+#define DARRAY_DATA struct darray_accum
+#define DARRAY_FUNCTOR_INIT darray_accum_init
+#define DARRAY_FUNCTOR_RELEASE darray_accum_release
+#define DARRAY_FUNCTOR_COPY darray_accum_copy
+#define DARRAY_FUNCTOR_COPY_AND_RELEASE darray_accum_copy_and_release
+#include <rsys/dynamic_array.h>
+
struct sdis_device {
struct logger* logger;
struct mem_allocator* allocator;
@@ -30,6 +43,7 @@ struct sdis_device {
int verbose;
struct flist_name names;
+ struct darray_tile tiles;
struct s2d_device* s2d;
struct s3d_device* s3d;
diff --git a/src/sdis_interface.c b/src/sdis_interface.c
@@ -41,15 +41,12 @@ check_interface_shader
type1 = sdis_medium_get_type(back);
/* Fluid<->solid interface */
- if(type0 != type1 && shader->convection_coef == NULL) {
- return 0;
- }
-
- /* Solid<->solid interface */
- if(type0 == SDIS_MEDIUM_SOLID
- && type1 == SDIS_MEDIUM_SOLID
- && shader->convection_coef) {
- return 0;
+ if(type0 != type1) {
+ if(shader->convection_coef == NULL
+ || shader->emissivity == NULL
+ || shader->specular_fraction == NULL) {
+ return 0;
+ }
}
return 1;
diff --git a/src/sdis_interface_c.h b/src/sdis_interface_c.h
@@ -76,5 +76,23 @@ interface_get_convection_coef
return interf->shader.convection_coef(frag, interf->data);
}
+static INLINE double
+interface_get_emissivity
+ (const struct sdis_interface* interf,
+ const struct sdis_interface_fragment* frag)
+{
+ ASSERT(interf && frag);
+ return interf->shader.emissivity(frag, interf->data);
+}
+
+static INLINE double
+interface_get_specular_fraction
+ (const struct sdis_interface* interf,
+ const struct sdis_interface_fragment* frag)
+{
+ ASSERT(interf && frag);
+ return interf->shader.specular_fraction(frag, interf->data);
+}
+
#endif /* SDIS_INTERFACE_C_H */
diff --git a/src/sdis_scene.c b/src/sdis_scene.c
@@ -361,6 +361,7 @@ scene_create
ref_init(&scn->ref);
SDIS(device_ref_get(dev));
scn->dev = dev;
+ scn->ambient_radiative_temperature = -1;
darray_interf_init(dev->allocator, &scn->interfaces);
darray_interf_init(dev->allocator, &scn->prim_interfaces);
diff --git a/src/sdis_scene_c.h b/src/sdis_scene_c.h
@@ -40,6 +40,8 @@ struct sdis_scene {
struct s2d_scene_view* s2d_view;
struct s3d_scene_view* s3d_view;
+ double ambient_radiative_temperature; /* In Kelvin */
+
ref_T ref;
struct sdis_device* dev;
};
diff --git a/src/sdis_solve.c b/src/sdis_solve.c
@@ -0,0 +1,405 @@
+/* Copyright (C) |Meso|Star> 2016-2018 (contact@meso-star.com)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. */
+
+#include "sdis.h"
+#include "sdis_camera.h"
+#include "sdis_device_c.h"
+#include "sdis_estimator_c.h"
+#include "sdis_solve_Xd.h"
+
+/* Generate the 2D solver */
+#define SDIS_SOLVE_DIMENSION 2
+#include "sdis_solve_Xd.h"
+
+/* Generate the 3D solver */
+#define SDIS_SOLVE_DIMENSION 3
+#include "sdis_solve_Xd.h"
+
+#include <star/ssp.h>
+#include <omp.h>
+
+/*******************************************************************************
+ * Helper functions
+ ******************************************************************************/
+static FINLINE uint16_t
+morton2D_decode(const uint32_t u32)
+{
+ uint32_t x = u32 & 0x55555555;
+ x = (x | (x >> 1)) & 0x33333333;
+ x = (x | (x >> 2)) & 0x0F0F0F0F;
+ x = (x | (x >> 4)) & 0x00FF00FF;
+ x = (x | (x >> 8)) & 0x0000FFFF;
+ return (uint16_t)x;
+}
+
+static res_T
+solve_pixel
+ (struct sdis_scene* scn,
+ struct ssp_rng* rng,
+ const struct sdis_medium* mdm,
+ const struct sdis_camera* cam,
+ const double time, /* Observation time */
+ const double fp_to_meter, /* Scale from floating point units to meters */
+ const double Tarad, /* In Kelvin */
+ const double Tref, /* In Kelvin */
+ const size_t ipix[2], /* Pixel coordinate in the image plane */
+ const size_t nrealisations,
+ const double pix_sz[2], /* Pixel size in the normalized image plane */
+ struct sdis_accum* accum)
+{
+ double sum_weights = 0;
+ double sum_weights_sqr = 0;
+ size_t N = 0; /* #realisations that do not fail */
+ size_t irealisation;
+ res_T res = RES_OK;
+ ASSERT(scn && mdm && rng && cam && ipix && nrealisations);
+ ASSERT(pix_sz && pix_sz[0] > 0 && pix_sz[1] > 0);
+
+ FOR_EACH(irealisation, 0, nrealisations) {
+ double samp[2]; /* Pixel sample */
+ double ray_pos[3];
+ double ray_dir[3];
+ double w = 0;
+
+ /* Generate a sample into the pixel to estimate */
+ samp[0] = ((double)ipix[0] + ssp_rng_canonical(rng)) * pix_sz[0];
+ samp[1] = ((double)ipix[1] + ssp_rng_canonical(rng)) * pix_sz[1];
+
+ /* Generate a ray starting from the camera position and passing through
+ * pixel sample */
+ camera_ray(cam, samp, ray_pos, ray_dir);
+
+ /* Launch the realisation */
+ res = ray_realisation_3d(scn, rng, mdm, ray_pos, ray_dir,
+ time, fp_to_meter, Tarad, Tref, &w);
+ if(res == RES_OK) {
+ sum_weights += w;
+ sum_weights_sqr += w*w;
+ ++N;
+ } else if(res != RES_BAD_OP) {
+ goto error;
+ }
+ }
+
+ accum->sum_weights = sum_weights;
+ accum->sum_weights_sqr = sum_weights_sqr;
+ accum->nweights = N;
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
+static res_T
+solve_tile
+ (struct sdis_scene* scn,
+ struct ssp_rng* rng,
+ const struct sdis_medium* mdm,
+ const struct sdis_camera* cam,
+ const double time,
+ const double fp_to_meter,
+ const double Tarad,
+ const double Tref,
+ const size_t origin[2], /* Tile origin in image plane */
+ const size_t size[2], /* #pixels in X and Y */
+ const size_t spp, /* #samples per pixel */
+ const double pix_sz[2], /* Pixel size in the normalized image plane */
+ struct sdis_accum* accums)
+{
+ size_t mcode; /* Morton code of the tile pixel */
+ size_t npixels;
+ res_T res = RES_OK;
+ ASSERT(scn && rng && mdm && cam && spp && origin && accums);
+ ASSERT(size &&size[0] && size[1]);
+ ASSERT(pix_sz && pix_sz[0] > 0 && pix_sz[1] > 0);
+
+ /* Adjust the #pixels to process them wrt a morton order */
+ npixels = round_up_pow2(MMAX(size[0], size[1]));
+ npixels *= npixels;
+
+ FOR_EACH(mcode, 0, npixels) {
+ size_t ipix[2];
+ struct sdis_accum* accum;
+
+ ipix[0] = morton2D_decode((uint32_t)(mcode>>0));
+ if(ipix[0] >= size[0]) continue;
+ ipix[1] = morton2D_decode((uint32_t)(mcode>>1));
+ if(ipix[1] >= size[1]) continue;
+
+ accum = accums + ipix[1]*size[0] + ipix[0];
+ ipix[0] = ipix[0] + origin[0];
+ ipix[1] = ipix[1] + origin[1];
+
+ res = solve_pixel(scn, rng, mdm, cam, time, fp_to_meter, Tarad, Tref, ipix,
+ spp, pix_sz, accum);
+ if(res != RES_OK) goto error;
+ }
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
+/*******************************************************************************
+ * Exported functions
+ ******************************************************************************/
+res_T
+sdis_solve_probe
+ (struct sdis_scene* scn,
+ const size_t nrealisations,
+ const double position[3],
+ const double time,
+ const double fp_to_meter,/* Scale factor from floating point unit to meter */
+ const double Tarad, /* Ambient radiative temperature */
+ const double Tref, /* Reference temperature */
+ struct sdis_estimator** out_estimator)
+{
+ const struct sdis_medium* medium = NULL;
+ struct sdis_estimator* estimator = NULL;
+ struct ssp_rng_proxy* rng_proxy = NULL;
+ struct ssp_rng** rngs = NULL;
+ double weight = 0;
+ double sqr_weight = 0;
+ size_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) {
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ /* Create the proxy RNG */
+ res = ssp_rng_proxy_create(scn->dev->allocator, &ssp_rng_mt19937_64,
+ scn->dev->nthreads, &rng_proxy);
+ if(res != RES_OK) goto error;
+
+ /* Create the per thread RNG */
+ rngs = MEM_CALLOC
+ (scn->dev->allocator, scn->dev->nthreads, sizeof(struct ssp_rng*));
+ if(!rngs) {
+ res = RES_MEM_ERR;
+ goto error;
+ }
+ FOR_EACH(i, 0, scn->dev->nthreads) {
+ res = ssp_rng_proxy_create_rng(rng_proxy, i, rngs+i);
+ if(res != RES_OK) goto error;
+ }
+
+ /* Create the estimator */
+ res = estimator_create(scn->dev, &estimator);
+ if(res != RES_OK) goto error;
+
+ /* Retrieve the medium in which the submitted position lies */
+ res = scene_get_medium(scn, position, &medium);
+ if(res != RES_OK) goto error;
+
+ /* Here we go! Launch the Monte Carlo estimation */
+ #pragma omp parallel for schedule(static) reduction(+:weight,sqr_weight,N)
+ for(irealisation = 0; irealisation < nrealisations; ++irealisation) {
+ res_T res_local;
+ double w;
+ const int ithread = omp_get_thread_num();
+ struct ssp_rng* rng = rngs[ithread];
+
+ if(ATOMIC_GET(&res) != RES_OK) continue; /* An error occured */
+
+ if(scene_is_2d(scn)) {
+ res_local = probe_realisation_2d
+ (scn, rng, medium, position, time, fp_to_meter, Tarad, Tref, &w);
+ } else {
+ res_local = probe_realisation_3d
+ (scn, rng, medium, position, time, fp_to_meter, Tarad, Tref, &w);
+ }
+ if(res_local != RES_OK) {
+ if(res_local != RES_BAD_OP) {
+ ATOMIC_SET(&res, res_local);
+ continue;
+ }
+ } else {
+ weight += w;
+ sqr_weight += w*w;
+ ++N;
+ }
+ }
+
+ estimator->nrealisations = N;
+ estimator->nfailures = nrealisations - N;
+ estimator->temperature.E = weight / (double)N;
+ estimator->temperature.V =
+ sqr_weight / (double)N
+ - estimator->temperature.E * estimator->temperature.E;
+ estimator->temperature.SE = sqrt(estimator->temperature.V / (double)N);
+
+exit:
+ if(rngs) {
+ FOR_EACH(i, 0, scn->dev->nthreads) {
+ if(rngs[i]) SSP(rng_ref_put(rngs[i]));
+ }
+ MEM_RM(scn->dev->allocator, rngs);
+ }
+ if(rng_proxy) SSP(rng_proxy_ref_put(rng_proxy));
+ if(out_estimator) *out_estimator = estimator;
+ return (res_T)res;
+error:
+ if(estimator) {
+ SDIS(estimator_ref_put(estimator));
+ estimator = NULL;
+ }
+ goto exit;
+}
+
+res_T
+sdis_solve_camera
+ (struct sdis_scene* scn,
+ const struct sdis_camera* cam,
+ const double time,
+ const double fp_to_meter, /* Scale from floating point units to meters */
+ const double Tarad, /* In Kelvin */
+ const double Tref, /* In Kelvin */
+ const size_t width, /* #pixels in X */
+ const size_t height, /* #pixels in Y */
+ const size_t spp, /* #samples per pixel */
+ sdis_write_accums_T writer,
+ void* writer_data)
+{
+ #define TILE_SIZE 32 /* definition in X & Y of a tile */
+ STATIC_ASSERT(IS_POW2(TILE_SIZE), TILE_SIZE_must_be_a_power_of_2);
+
+ const struct sdis_medium* medium = NULL;
+ struct darray_accum* tiles = NULL;
+ struct ssp_rng_proxy* rng_proxy = NULL;
+ struct ssp_rng** rngs = NULL;
+ size_t ntiles_x, ntiles_y, ntiles;
+ double pix_sz[2]; /* Size of a pixel in the normalized image plane */
+ int64_t mcode; /* Morton code of a tile */
+ size_t i;
+ ATOMIC res = RES_OK;
+
+ if(!scn || !cam || time < 0 || fp_to_meter <= 0 || Tref < 0 || !width
+ || !height || !spp || !writer) {
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ if(scene_is_2d(scn)) {
+ log_err(scn->dev, "%s: 2D scene are not supported.\n", FUNC_NAME);
+ goto error;
+ }
+
+ /* Retrieve the medium in which the submitted position lies */
+ res = scene_get_medium(scn, cam->position, &medium);
+ if(res != RES_OK) goto error;
+
+ if(medium->type != SDIS_MEDIUM_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;
+ goto error;
+ }
+
+ /* Create the proxy RNG */
+ res = ssp_rng_proxy_create(scn->dev->allocator, &ssp_rng_mt19937_64,
+ scn->dev->nthreads, &rng_proxy);
+ if(res != RES_OK) goto error;
+
+ /* Create the per thread RNG */
+ rngs = MEM_CALLOC
+ (scn->dev->allocator, scn->dev->nthreads, sizeof(struct ssp_rng*));
+ if(!rngs) {
+ res = RES_MEM_ERR;
+ goto error;
+ }
+ FOR_EACH(i, 0, scn->dev->nthreads) {
+ res = ssp_rng_proxy_create_rng(rng_proxy, i, rngs+i);
+ if(res != RES_OK) goto error;
+ }
+
+ /* Allocate per thread buffer of accumulations */
+ tiles = darray_tile_data_get(&scn->dev->tiles);
+ ASSERT(darray_tile_size_get(&scn->dev->tiles) == scn->dev->nthreads);
+ FOR_EACH(i, 0, scn->dev->nthreads) {
+ const size_t naccums = TILE_SIZE * TILE_SIZE;
+ res = darray_accum_resize(tiles+i, naccums);
+ if(res != RES_OK) goto error;
+ }
+
+ ntiles_x = (width + (TILE_SIZE-1)/*ceil*/)/TILE_SIZE;
+ ntiles_y = (height + (TILE_SIZE-1)/*ceil*/)/TILE_SIZE;
+ ntiles = round_up_pow2(MMAX(ntiles_x, ntiles_y));
+ ntiles *= ntiles;
+
+ pix_sz[0] = 1.0 / (double)width;
+ pix_sz[1] = 1.0 / (double)height;
+
+ omp_set_num_threads((int)scn->dev->nthreads);
+ #pragma omp parallel for schedule(static, 1/*chunki size*/)
+ for(mcode = 0; mcode < (int64_t)ntiles; ++mcode) {
+ size_t tile_org[2] = {0, 0};
+ size_t tile_sz[2] = {0, 0};
+ const int ithread = omp_get_thread_num();
+ struct sdis_accum* accums = NULL;
+ struct ssp_rng* rng = rngs[ithread];
+ res_T res_local = RES_OK;
+
+ if(ATOMIC_GET(&res) != RES_OK) continue;
+
+ tile_org[0] = morton2D_decode((uint32_t)(mcode>>0));
+ if(tile_org[0] >= ntiles_x) continue; /* Discard tile */
+ tile_org[1] = morton2D_decode((uint32_t)(mcode>>1));
+ if(tile_org[1] >= ntiles_y) continue; /* Disaard tile */
+
+ /* Setup the tile coordinates in the image plane */
+ 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]);
+
+ /* Fetch the accumulations buffer */
+ accums = darray_accum_data_get(tiles+ithread);
+
+ /* Draw the tile */
+ res_local = solve_tile(scn, rng, medium, cam, time, fp_to_meter, Tarad,
+ Tref, tile_org, tile_sz, spp, pix_sz, accums);
+ if(res_local != RES_OK) {
+ ATOMIC_SET(&res, res_local);
+ continue;
+ }
+
+ /* Write the accumulations */
+ res_local = writer(writer_data, tile_org, tile_sz, accums);
+ if(res_local != RES_OK) {
+ ATOMIC_SET(&res, res_local);
+ continue;
+ }
+ }
+
+exit:
+ if(rngs) {
+ FOR_EACH(i, 0, scn->dev->nthreads) {
+ if(rngs[i]) SSP(rng_ref_put(rngs[i]));
+ }
+ MEM_RM(scn->dev->allocator, rngs);
+ }
+ if(rng_proxy) SSP(rng_proxy_ref_put(rng_proxy));
+ return (res_T)res;
+error:
+ goto exit;
+}
+
diff --git a/src/sdis_solve_Xd.h b/src/sdis_solve_Xd.h
@@ -0,0 +1,841 @@
+/* Copyright (C) |Meso|Star> 2016-2018 (contact@meso-star.com)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. */
+
+#ifndef SDIS_SOLVE_DIMENSION
+#ifndef SDIS_SOLVE_XD_H
+#define SDIS_SOLVE_XD_H
+
+#include "sdis_device_c.h"
+#include "sdis_interface_c.h"
+#include "sdis_medium_c.h"
+#include "sdis_scene_c.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
+ * to handle numerical imprecisions */
+#define RAY_RANGE_MAX_SCALE 1.0001f
+
+#define BOLTZMANN_CONSTANT 5.6696e-8 /* W/m^2/K^4 */
+
+struct rwalk_context {
+ double Tarad; /* Ambient radiative temperature */
+ double Tref3; /* Reference temperature ^ 3 */
+};
+
+/* 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])
+{
+ float tmp[3];
+ float cos_V_N;
+ ASSERT(res && V && N);
+ ASSERT(f3_is_normalized(V) && f3_is_normalized(N));
+ cos_V_N = f3_dot(V, N);
+ f3_mulf(tmp, N, 2*cos_V_N);
+ f3_sub(res, tmp, V);
+ return res;
+}
+
+#endif /* SDIS_SOLVE_XD_H */
+#else
+
+#if (SDIS_SOLVE_DIMENSION == 2)
+ #include <rsys/double2.h>
+ #include <rsys/float2.h>
+ #include <star/s2d.h>
+#elif (SDIS_SOLVE_DIMENSION == 3)
+ #include <rsys/double2.h>
+ #include <rsys/double3.h>
+ #include <rsys/float3.h>
+ #include <star/s3d.h>
+#else
+ #error "Invalid SDIS_SOLVE_DIMENSION value."
+#endif
+
+/* Syntactic sugar */
+#define DIM SDIS_SOLVE_DIMENSION
+
+/* Star-XD macros generic to SDIS_SOLVE_DIMENSION */
+#define sXd(Name) CONCAT(CONCAT(CONCAT(s, DIM), d_), Name)
+#define SXD_HIT_NONE CONCAT(CONCAT(S,DIM), D_HIT_NONE)
+#define SXD_HIT_NULL CONCAT(CONCAT(S,DIM), D_HIT_NULL)
+#define SXD_HIT_NULL__ CONCAT(CONCAT(S, DIM), D_HIT_NULL__)
+#define SXD CONCAT(CONCAT(S, DIM), D)
+
+/* Vector macros generic to SDIS_SOLVE_DIMENSION */
+#define dX(Func) CONCAT(CONCAT(CONCAT(d, DIM), _), Func)
+#define fX(Func) CONCAT(CONCAT(CONCAT(f, DIM), _), Func)
+#define fX_set_dX CONCAT(CONCAT(CONCAT(f, DIM), _set_d), DIM)
+#define dX_set_fX CONCAT(CONCAT(CONCAT(d, DIM), _set_f), DIM)
+
+/* Macro making generic its subimitted name to SDIS_SOLVE_DIMENSION */
+#define XD(Name) CONCAT(CONCAT(CONCAT(Name, _), DIM), d)
+
+/* Current state of the random walk */
+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 */
+};
+static const struct XD(rwalk) XD(RWALK_NULL) = {
+ SDIS_RWALK_VERTEX_NULL__, NULL, SXD_HIT_NULL__
+};
+
+struct XD(temperature) {
+ res_T (*func)/* Next function to invoke in order to compute the temperature */
+ (const struct sdis_scene* scn,
+ const double fp_to_meter,
+ const struct rwalk_context* ctx,
+ struct XD(rwalk)* rwalk,
+ struct ssp_rng* rng,
+ struct XD(temperature)* temp);
+ double value; /* Current value of the temperature */
+ int done;
+};
+static const struct XD(temperature) XD(TEMPERATURE_NULL) = { NULL, 0, 0 };
+
+static res_T
+XD(boundary_temperature)
+ (const struct sdis_scene* scn,
+ const double fp_to_meter,
+ const struct rwalk_context* ctx,
+ struct XD(rwalk)* rwalk,
+ struct ssp_rng* rng,
+ struct XD(temperature)* T);
+
+static res_T
+XD(solid_temperature)
+ (const struct sdis_scene* scn,
+ const double fp_to_meter,
+ const struct rwalk_context* ctx,
+ struct XD(rwalk)* rwalk,
+ struct ssp_rng* rng,
+ struct XD(temperature)* T);
+
+static res_T
+XD(fluid_temperature)
+ (const struct sdis_scene* scn,
+ const double fp_to_meter,
+ const struct rwalk_context* ctx,
+ struct XD(rwalk)* rwalk,
+ struct ssp_rng* rng,
+ struct XD(temperature)* T);
+
+static res_T
+XD(radiative_temperature)
+ (const struct sdis_scene* scn,
+ const double fp_to_meter,
+ const struct rwalk_context* ctx,
+ struct XD(rwalk)* rwalk,
+ struct ssp_rng* rng,
+ struct XD(temperature)* T);
+
+/*******************************************************************************
+ * Helper functions
+ ******************************************************************************/
+
+static FINLINE void
+XD(move_pos)(double pos[DIM], const float dir[DIM], const float delta)
+{
+ ASSERT(pos && dir);
+ pos[0] += dir[0] * delta;
+ pos[1] += dir[1] * delta;
+#if(SDIS_SOLVE_DIMENSION == 3)
+ pos[2] += dir[2] * delta;
+#endif
+}
+
+/* Check that the interface fragment is consistent with the current state of
+ * the random walk */
+static INLINE int
+XD(check_rwalk_fragment_consistency)
+ (const struct XD(rwalk)* rwalk,
+ const struct sdis_interface_fragment* frag)
+{
+ double N[DIM];
+ double uv[2] = {0, 0};
+ ASSERT(rwalk && frag);
+ dX(normalize)(N, dX_set_fX(N, rwalk->hit.normal));
+ if( SXD_HIT_NONE(&rwalk->hit)
+ || !dX(eq_eps)(rwalk->vtx.P, frag->P, 1.e-6)
+ || !dX(eq_eps)(N, frag->Ng, 1.e-6)
+ || !( (IS_INF(rwalk->vtx.time) && IS_INF(frag->time))
+ || eq_eps(rwalk->vtx.time, frag->time, 1.e-6))) {
+ return 0;
+ }
+#if (SDIS_SOLVE_DIMENSION == 2)
+ uv[0] = rwalk->hit.u;
+#else
+ d2_set_f2(uv, rwalk->hit.uv);
+#endif
+ return d2_eq_eps(uv, frag->uv, 1.e-6);
+}
+
+static res_T
+XD(trace_radiative_path)
+ (const struct sdis_scene* scn,
+ const float ray_dir[3],
+ const double fp_to_meter,
+ const struct rwalk_context* ctx,
+ struct XD(rwalk)* rwalk,
+ struct ssp_rng* rng,
+ struct XD(temperature)* T)
+{
+ /* 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};
+ float dir[3] = {0, 0, 0};
+ res_T res = RES_OK;
+
+ ASSERT(scn && ray_dir && fp_to_meter > 0 && ctx && rwalk && rng && T);
+ (void)fp_to_meter;
+
+ f3_set(dir, ray_dir);
+
+ /* Launch the radiative random walk */
+ for(;;) {
+ const struct sdis_interface* interf = NULL;
+ struct sdis_interface_fragment frag = SDIS_INTERFACE_FRAGMENT_NULL;
+ const struct sdis_medium* chk_mdm = NULL;
+ double alpha;
+ double epsilon;
+ double r;
+ float pos[DIM];
+ const float range[2] = { 0, FLT_MAX };
+
+ fX_set_dX(pos, rwalk->vtx.P);
+
+ /* Trace the radiative ray */
+#if (SDIS_SOLVE_DIMENSION == 2)
+ SXD(scene_view_trace_ray_3d
+ (scn->sXd(view), pos, dir, range, &rwalk->hit, &rwalk->hit));
+#else
+ SXD(scene_view_trace_ray
+ (scn->sXd(view), pos, dir, range, &rwalk->hit, &rwalk->hit));
+#endif
+ if(SXD_HIT_NONE(&rwalk->hit)) { /* Fetch the ambient radiative temperature */
+ 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",
+ FUNC_NAME,
+ ctx->Tarad,
+ SPLIT3(rwalk->vtx.P));
+ res = RES_BAD_ARG;
+ goto error;
+ }
+ }
+
+ /* 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);
+
+ /* Fetch the interface emissivity */
+ epsilon = interface_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));
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ /* Switch in boundary temperature ? */
+ r = ssp_rng_canonical(rng);
+ if(r < epsilon) {
+ T->func = XD(boundary_temperature);
+ break;
+ }
+
+ /* 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) {
+ chk_mdm = interf->medium_back;
+ fX(minus)(N, N);
+ } else {
+ chk_mdm = interf->medium_front;
+ }
+
+ 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;
+ goto error;
+ }
+ alpha = interface_get_specular_fraction(interf, &frag);
+ r = ssp_rng_canonical(rng);
+ if(r < alpha) { /* Sample specular part */
+ reflect(dir, f3_minus(dir, dir), N);
+ } else { /* Sample diffuse part */
+ ssp_ran_hemisphere_cos_float(rng, N, dir, NULL);
+ }
+ }
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
+res_T
+XD(radiative_temperature)
+ (const struct sdis_scene* scn,
+ const double fp_to_meter,
+ const struct rwalk_context* ctx,
+ struct XD(rwalk)* rwalk,
+ 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};
+ float dir[3] = {0, 0, 0};
+ res_T res = RES_OK;
+
+ ASSERT(scn && fp_to_meter > 0 && ctx && rwalk && rng && T);
+ 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) {
+ fX(minus(N, N));
+ }
+
+ /* Cosine weighted sampling of a direction around the surface normal */
+ ssp_ran_hemisphere_cos_float(rng, N, dir, NULL);
+
+ /* Launch the radiative random walk */
+ res = XD(trace_radiative_path)(scn, dir, fp_to_meter, ctx, rwalk, rng, T);
+ if(res != RES_OK) goto error;
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
+res_T
+XD(fluid_temperature)
+ (const struct sdis_scene* scn,
+ const double fp_to_meter,
+ const struct rwalk_context* ctx,
+ struct XD(rwalk)* rwalk,
+ struct ssp_rng* rng,
+ struct XD(temperature)* T)
+{
+ double tmp;
+ (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);
+
+ tmp = fluid_get_temperature(rwalk->mdm, &rwalk->vtx);
+ if(tmp < 0) {
+ log_err(scn->dev, "%s: unknown fluid temperature is not supported.\n",
+ FUNC_NAME);
+ return RES_BAD_OP;
+ }
+ T->value += tmp;
+ T->done = 1;
+ return RES_OK;
+}
+
+static void
+XD(solid_solid_boundary_temperature)
+ (const struct sdis_scene* scn,
+ const double fp_to_meter,
+ const struct rwalk_context* ctx,
+ const struct sdis_interface_fragment* frag,
+ struct XD(rwalk)* rwalk,
+ struct ssp_rng* rng,
+ struct XD(temperature)* T)
+{
+ const struct sdis_interface* interf = NULL;
+ const struct sdis_medium* solid_front = NULL;
+ const struct sdis_medium* solid_back = NULL;
+ double lambda_front, lambda_back;
+ double delta_front_boundary, delta_back_boundary;
+ double delta_front_boundary_meter, delta_back_boundary_meter;
+ double delta_boundary;
+ double proba;
+ double tmp;
+ double r;
+ float pos[DIM], dir[DIM], range[2];
+ ASSERT(scn && fp_to_meter > 0 && ctx && frag && rwalk && rng && T);
+ ASSERT(XD(check_rwalk_fragment_consistency)(rwalk, frag));
+ (void)frag, (void)ctx;
+
+ /* Retrieve the current boundary media */
+ 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);
+
+ /* 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);
+
+ r = ssp_rng_canonical(rng);
+ fX(normalize)(dir, rwalk->hit.normal);
+ if(r < proba) { /* Reinject in front */
+ delta_boundary = delta_front_boundary;
+ rwalk->mdm = solid_front;
+ } else { /* Reinject in back */
+ delta_boundary = delta_back_boundary;
+ fX(minus)(dir, dir);
+ rwalk->mdm = solid_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);
+}
+
+static void
+XD(solid_fluid_boundary_temperature)
+ (const struct sdis_scene* scn,
+ const double fp_to_meter,
+ const struct rwalk_context* ctx,
+ const struct sdis_interface_fragment* frag,
+ struct XD(rwalk)* rwalk,
+ struct ssp_rng* rng,
+ struct XD(temperature)* T)
+{
+ const struct sdis_interface* interf = NULL;
+ const struct sdis_medium* mdm_front = NULL;
+ const struct sdis_medium* mdm_back = NULL;
+ const struct sdis_medium* solid = NULL;
+ const struct sdis_medium* fluid = NULL;
+ double hc;
+ double hr;
+ double epsilon; /* Interface emissivity */
+ double lambda;
+ double fluid_proba;
+ double radia_proba;
+ double delta_boundary;
+ double r;
+ double tmp;
+ float dir[DIM], pos[DIM], range[2];
+
+ 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 */
+ 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) {
+ solid = mdm_front;
+ fluid = mdm_back;
+ } else {
+ solid = mdm_back;
+ fluid = mdm_front;
+ }
+
+ /* Fetch the solid properties */
+ lambda = solid_get_thermal_conductivity(solid, &rwalk->vtx);
+ delta_boundary = solid_get_delta_boundary(solid, &rwalk->vtx);
+
+ /* Fetch the boundary properties */
+ epsilon = interface_get_emissivity(interf, frag);
+ 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);
+ 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;
+ T->func = XD(fluid_temperature);
+ } else { /* Solid random walk */
+ rwalk->mdm = solid;
+ fX(normalize)(dir, rwalk->hit.normal);
+ if(solid == mdm_back) fX(minus)(dir, dir);
+
+ /* "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);
+
+ /* Switch in solid random walk */
+ T->func = XD(solid_temperature);
+ }
+}
+
+res_T
+XD(boundary_temperature)
+ (const struct sdis_scene* scn,
+ const double fp_to_meter,
+ const struct rwalk_context* ctx,
+ struct XD(rwalk)* rwalk,
+ struct ssp_rng* rng,
+ struct XD(temperature)* T)
+{
+ struct sdis_interface_fragment frag = SDIS_INTERFACE_FRAGMENT_NULL;
+ const struct sdis_interface* interf = NULL;
+ const struct sdis_medium* mdm_front = NULL;
+ const struct sdis_medium* mdm_back = NULL;
+ double tmp;
+ ASSERT(scn && fp_to_meter > 0 && ctx && rwalk && rng && T);
+ ASSERT(!SXD_HIT_NONE(&rwalk->hit));
+
+ XD(setup_interface_fragment)(&frag, &rwalk->vtx, &rwalk->hit);
+
+ /* 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);
+ if(tmp >= 0) {
+ T->value += tmp;
+ T->done = 1;
+ return RES_OK;
+ }
+
+ mdm_front = interface_get_medium(interf, SDIS_FRONT);
+ mdm_back = interface_get_medium(interf, SDIS_BACK);
+
+ if(mdm_front->type == mdm_back->type) {
+ XD(solid_solid_boundary_temperature)
+ (scn, fp_to_meter, ctx, &frag, rwalk, rng, T);
+ } else {
+ XD(solid_fluid_boundary_temperature)
+ (scn, fp_to_meter, ctx, &frag, rwalk, rng, T);
+ }
+ return RES_OK;
+}
+
+res_T
+XD(solid_temperature)
+ (const struct sdis_scene* scn,
+ const double fp_to_meter,
+ const struct rwalk_context* ctx,
+ struct XD(rwalk)* rwalk,
+ struct ssp_rng* rng,
+ struct XD(temperature)* T)
+{
+ 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);
+ (void)ctx;
+
+ /* Check the random walk consistency */
+ CHK(scene_get_medium(scn, rwalk->vtx.P, &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;
+ }
+ /* Save the submitted position */
+ dX(set)(position_start, rwalk->vtx.P);
+
+ do { /* Solid random walk */
+ struct sXd(hit) hit0, hit1;
+ double lambda; /* Thermal conductivity */
+ double rho; /* Volumic mass */
+ double cp; /* Calorific capacity */
+ double tau, mu;
+ double tmp;
+ float delta, delta_solid; /* Random walk numerical parameter */
+ float range[2];
+ float dir0[DIM], dir1[DIM];
+ float org[DIM];
+
+ /* Check the limit condition */
+ tmp = solid_get_temperature(mdm, &rwalk->vtx);
+ if(tmp >= 0) {
+ T->value += tmp;
+ T->done = 1;
+ return RES_OK;
+ }
+
+ /* Fetch solid properties */
+ delta_solid = (float)solid_get_delta(mdm, &rwalk->vtx);
+ lambda = solid_get_thermal_conductivity(mdm, &rwalk->vtx);
+ rho = solid_get_volumic_mass(mdm, &rwalk->vtx);
+ cp = solid_get_calorific_capacity(mdm, &rwalk->vtx);
+
+#if (SDIS_SOLVE_DIMENSION == 2)
+ /* Sample a direction around 2PI */
+ ssp_ran_circle_uniform_float(rng, dir0, NULL);
+#else
+ /* Sample a direction around 4PI */
+ ssp_ran_sphere_uniform_float(rng, dir0, NULL);
+#endif
+
+ /* Trace a ray along the sampled direction and its opposite to check if a
+ * surface is hit in [0, delta_solid]. */
+ fX_set_dX(org, rwalk->vtx.P);
+ fX(minus)(dir1, dir0);
+ hit0 = hit1 = SXD_HIT_NULL;
+ range[0] = 0.f, range[1] = delta_solid*RAY_RANGE_MAX_SCALE;
+ SXD(scene_view_trace_ray(scn->sXd(view), org, dir0, range, NULL, &hit0));
+ SXD(scene_view_trace_ray(scn->sXd(view), org, dir1, range, NULL, &hit1));
+
+ if(SXD_HIT_NONE(&hit0) && SXD_HIT_NONE(&hit1)) {
+ /* Hit nothing: move along dir0 of the original delta */
+ delta = delta_solid;
+ } else {
+ /* Hit something: move along dir0 of the minimum hit distance */
+ delta = MMIN(hit0.distance, hit1.distance);
+ }
+
+ /* Sample the time */
+ mu = (2*DIM*lambda) / (rho*cp*delta*fp_to_meter*delta*fp_to_meter);
+ tau = ssp_ran_exp(rng, mu);
+ rwalk->vtx.time -= tau;
+
+ /* Check the initial condition */
+ tmp = solid_get_temperature(mdm, &rwalk->vtx);
+ if(tmp >= 0) {
+ T->value += tmp;
+ T->done = 1;
+ return RES_OK;
+ }
+
+ /* Define if the random walk hits something along dir0 */
+ rwalk->hit = hit0.distance > delta ? SXD_HIT_NULL : hit0;
+
+ /* 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);
+ } 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);
+ }
+
+ /* 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));
+ }
+ return RES_BAD_OP;
+ }
+
+ /* Keep going while the solid random walk does not hit an interface */
+ } while(SXD_HIT_NONE(&rwalk->hit));
+
+ T->func = XD(boundary_temperature);
+ return RES_OK;
+}
+
+static res_T
+XD(compute_temperature)
+ (struct sdis_scene* scn,
+ const double fp_to_meter,
+ const struct rwalk_context* ctx,
+ struct XD(rwalk)* rwalk,
+ struct ssp_rng* rng,
+ struct XD(temperature)* T)
+{
+#ifndef NDEBUG
+ struct XD(temperature)* stack = NULL;
+ size_t istack = 0;
+#endif
+ 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;
+
+#ifndef NDEBUG
+ sa_push(stack, *T);
+ ++istack;
+#endif
+ } while(!T->done);
+
+exit:
+#ifndef NDEBUG
+ sa_release(stack);
+#endif
+ return res;
+error:
+ goto exit;
+}
+
+static res_T
+XD(probe_realisation)
+ (struct sdis_scene* scn,
+ struct ssp_rng* rng,
+ const struct sdis_medium* medium,
+ const double position[],
+ const double time,
+ const double fp_to_meter,/* Scale factor from floating point unit to meter */
+ const double ambient_radiative_temperature,
+ const double reference_temperature,
+ double* weight)
+{
+ struct rwalk_context ctx;
+ struct XD(rwalk) rwalk = XD(RWALK_NULL);
+ struct XD(temperature) T = XD(TEMPERATURE_NULL);
+ res_T res = RES_OK;
+ 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;
+ default: FATAL("Unreachable code\n"); break;
+ }
+
+ dX(set)(rwalk.vtx.P, position);
+ rwalk.vtx.time = time;
+ rwalk.hit = SXD_HIT_NULL;
+ rwalk.mdm = medium;
+
+ ctx.Tarad = ambient_radiative_temperature;
+ ctx.Tref3 =
+ reference_temperature
+ * reference_temperature
+ * reference_temperature;
+
+ res = XD(compute_temperature)(scn, fp_to_meter, &ctx, &rwalk, rng, &T);
+ if(res != RES_OK) return res;
+
+ *weight = T.value;
+ return RES_OK;
+}
+
+#if SDIS_SOLVE_DIMENSION == 3
+static res_T
+XD(ray_realisation)
+ (struct sdis_scene* scn,
+ struct ssp_rng* rng,
+ const struct sdis_medium* medium,
+ const double position[],
+ const double direction[],
+ const double time,
+ const double fp_to_meter,
+ const double Tarad,
+ const double Tref,
+ double* weight)
+{
+ struct rwalk_context ctx;
+ struct XD(rwalk) rwalk = XD(RWALK_NULL);
+ struct XD(temperature) T = XD(TEMPERATURE_NULL);
+ 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);
+
+ dX(set)(rwalk.vtx.P, position);
+ rwalk.vtx.time = time;
+ rwalk.hit = SXD_HIT_NULL;
+ rwalk.mdm = medium;
+
+ ctx.Tarad = Tarad;
+ ctx.Tref3 = Tref*Tref*Tref;
+
+ f3_set_d3(dir, direction);
+
+ res = XD(trace_radiative_path)(scn, dir, fp_to_meter, &ctx, &rwalk, rng, &T);
+ if(res != RES_OK) goto error;
+
+ if(!T.done) {
+ res = XD(compute_temperature)(scn, fp_to_meter, &ctx, &rwalk, rng, &T);
+ if(res != RES_OK) goto error;
+ }
+
+ *weight = T.value;
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+#endif /* SDIS_SOLVE_DIMENSION == 3 */
+
+#undef SDIS_SOLVE_DIMENSION
+#undef DIM
+#undef sXd
+#undef SXD_HIT_NONE
+#undef SXD_HIT_NULL
+#undef SXD_HIT_NULL__
+#undef SXD
+#undef dX
+#undef fX
+#undef fX_set_dX
+#undef XD
+
+#endif /* !SDIS_SOLVE_DIMENSION */
+
diff --git a/src/sdis_solve_probe.c b/src/sdis_solve_probe.c
@@ -1,138 +0,0 @@
-/* Copyright (C) |Meso|Star> 2016-2018 (contact@meso-star.com)
- *
- * This program is free software: you can redistribute it and/or modify
- * it under the terms of the GNU General Public License as published by
- * the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * This program is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU General Public License
- * along with this program. If not, see <http://www.gnu.org/licenses/>. */
-
-#include "sdis.h"
-#include "sdis_device_c.h"
-#include "sdis_estimator_c.h"
-#include "sdis_solve_probe_Xd.h"
-
-/* Generate the 2D solver */
-#define SDIS_SOLVE_PROBE_DIMENSION 2
-#include "sdis_solve_probe_Xd.h"
-
-/* Generate the 3D solver */
-#define SDIS_SOLVE_PROBE_DIMENSION 3
-#include "sdis_solve_probe_Xd.h"
-
-#include <star/ssp.h>
-#include <omp.h>
-
-res_T
-sdis_solve_probe
- (struct sdis_scene* scn,
- const size_t nrealisations,
- const double position[3],
- const double time,
- const double fp_to_meter,/* Scale factor from floating point unit to meter */
- struct sdis_estimator** out_estimator)
-{
- const struct sdis_medium* medium = NULL;
- struct sdis_estimator* estimator = NULL;
- struct ssp_rng_proxy* rng_proxy = NULL;
- struct ssp_rng** rngs = NULL;
- double weight = 0;
- double sqr_weight = 0;
- size_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
- || !out_estimator) {
- res = RES_BAD_ARG;
- goto error;
- }
-
- /* Create the proxy RNG */
- res = ssp_rng_proxy_create(scn->dev->allocator, &ssp_rng_mt19937_64,
- scn->dev->nthreads, &rng_proxy);
- if(res != RES_OK) goto error;
-
- /* Create the per thread RNG */
- rngs = MEM_CALLOC
- (scn->dev->allocator, scn->dev->nthreads, sizeof(struct ssp_rng*));
- if(!rngs) {
- res = RES_MEM_ERR;
- goto error;
- }
- FOR_EACH(i, 0, scn->dev->nthreads) {
- res = ssp_rng_proxy_create_rng(rng_proxy, i, rngs+i);
- if(res != RES_OK) goto error;
- }
-
- /* Create the estimator */
- res = estimator_create(scn->dev, &estimator);
- if(res != RES_OK) goto error;
-
- /* Retrieve the medium in which the submitted position lies */
- res = scene_get_medium(scn, position, &medium);
- if(res != RES_OK) goto error;
-
- /* Here we go! Launch the Monte Carlo estimation */
- #pragma omp parallel for schedule(static) reduction(+:weight,sqr_weight,N)
- for(irealisation = 0; irealisation < nrealisations; ++irealisation) {
- res_T res_local;
- double w;
- const int ithread = omp_get_thread_num();
- struct ssp_rng* rng = rngs[ithread];
-
- if(ATOMIC_GET(&res) != RES_OK) continue; /* An error occured */
-
- if(scene_is_2d(scn)) {
- res_local = probe_realisation_2d
- (scn, rng, medium, position, time, fp_to_meter, &w);
- } else {
- res_local = probe_realisation_3d
- (scn, rng, medium, position, time, fp_to_meter, &w);
- }
- if(res_local != RES_OK) {
- if(res_local == RES_BAD_OP) {
- ++estimator->nfailures;
- } else {
- ATOMIC_SET(&res, res_local);
- continue;
- }
- } else {
- weight += w;
- sqr_weight += w*w;
- ++N;
- }
- }
-
- estimator->nrealisations = N;
- estimator->temperature.E = weight / (double)N;
- estimator->temperature.V =
- sqr_weight / (double)N
- - estimator->temperature.E * estimator->temperature.E;
- estimator->temperature.SE = sqrt(estimator->temperature.V / (double)N);
-
-exit:
- if(rngs) {
- FOR_EACH(i, 0, scn->dev->nthreads) {
- if(rngs[i]) SSP(rng_ref_put(rngs[i]));
- }
- MEM_RM(scn->dev->allocator, rngs);
- }
- if(rng_proxy) SSP(rng_proxy_ref_put(rng_proxy));
- if(out_estimator) *out_estimator = estimator;
- return (res_T)res;
-error:
- if(estimator) {
- SDIS(estimator_ref_put(estimator));
- estimator = NULL;
- }
- goto exit;
-}
-
diff --git a/src/sdis_solve_probe_Xd.h b/src/sdis_solve_probe_Xd.h
@@ -1,565 +0,0 @@
-/* Copyright (C) |Meso|Star> 2016-2018 (contact@meso-star.com)
- *
- * This program is free software: you can redistribute it and/or modify
- * it under the terms of the GNU General Public License as published by
- * the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * This program is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU General Public License
- * along with this program. If not, see <http://www.gnu.org/licenses/>. */
-
-#ifndef SDIS_SOLVE_PROBE_DIMENSION
-#ifndef SDIS_SOLVE_PROBE_XD_H
-#define SDIS_SOLVE_PROBE_XD_H
-
-#include "sdis_device_c.h"
-#include "sdis_interface_c.h"
-#include "sdis_medium_c.h"
-#include "sdis_scene_c.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
- * to handle numerical imprecisions */
-#define RAY_RANGE_MAX_SCALE 1.0001f
-
-#endif /* SDIS_SOLVE_PROBE_XD_H */
-#else
-
-#if (SDIS_SOLVE_PROBE_DIMENSION == 2)
- #include <rsys/double2.h>
- #include <rsys/float2.h>
- #include <star/s2d.h>
-#elif (SDIS_SOLVE_PROBE_DIMENSION == 3)
- #include <rsys/double2.h>
- #include <rsys/double3.h>
- #include <rsys/float3.h>
- #include <star/s3d.h>
-#else
- #error "Invalid SDIS_SOLVE_PROBE_DIMENSION value."
-#endif
-
-/* Syntactic sugar */
-#define DIM SDIS_SOLVE_PROBE_DIMENSION
-
-/* Star-XD macros generic to SDIS_SOLVE_PROBE_DIMENSION */
-#define sXd(Name) CONCAT(CONCAT(CONCAT(s, DIM), d_), Name)
-#define SXD_HIT_NONE CONCAT(CONCAT(S,DIM), D_HIT_NONE)
-#define SXD_HIT_NULL CONCAT(CONCAT(S,DIM), D_HIT_NULL)
-#define SXD_HIT_NULL__ CONCAT(CONCAT(S, DIM), D_HIT_NULL__)
-#define SXD CONCAT(CONCAT(S, DIM), D)
-
-/* Vector macros generic to SDIS_SOLVE_PROBE_DIMENSION */
-#define dX(Func) CONCAT(CONCAT(CONCAT(d, DIM), _), Func)
-#define fX(Func) CONCAT(CONCAT(CONCAT(f, DIM), _), Func)
-#define fX_set_dX CONCAT(CONCAT(CONCAT(f, DIM), _set_d), DIM)
-#define dX_set_fX CONCAT(CONCAT(CONCAT(d, DIM), _set_f), DIM)
-
-/* Macro making generic its subimitted name to SDIS_SOLVE_PROBE_DIMENSION */
-#define XD(Name) CONCAT(CONCAT(CONCAT(Name, _), DIM), d)
-
-/* Current state of the random walk */
-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 */
-};
-static const struct XD(rwalk) XD(RWALK_NULL) = {
- SDIS_RWALK_VERTEX_NULL__, NULL, SXD_HIT_NULL__
-};
-
-struct XD(temperature) {
- res_T (*func)/* Next function to invoke in order to compute the temperature */
- (const struct sdis_scene* scn,
- const double fp_to_meter,
- struct XD(rwalk)* rwalk,
- struct ssp_rng* rng,
- struct XD(temperature)* temp);
- double value; /* Current value of the temperature */
- int done;
-};
-static const struct XD(temperature) XD(TEMPERATURE_NULL) = { NULL, 0, 0 };
-
-static res_T
-XD(boundary_temperature)
- (const struct sdis_scene* scn,
- const double fp_to_meter,
- struct XD(rwalk)* rwalk,
- struct ssp_rng* rng,
- struct XD(temperature)* T);
-
-static res_T
-XD(solid_temperature)
- (const struct sdis_scene* scn,
- const double fp_to_meter,
- struct XD(rwalk)* rwalk,
- struct ssp_rng* rng,
- struct XD(temperature)* T);
-
-static res_T
-XD(fluid_temperature)
- (const struct sdis_scene* scn,
- const double fp_to_meter,
- struct XD(rwalk)* rwalk,
- struct ssp_rng* rng,
- struct XD(temperature)* T);
-
-/*******************************************************************************
- * Helper functions
- ******************************************************************************/
-static FINLINE void
-XD(move_pos)(double pos[DIM], const float dir[DIM], const float delta)
-{
- ASSERT(pos && dir);
- pos[0] += dir[0] * delta;
- pos[1] += dir[1] * delta;
-#if (SDIS_SOLVE_PROBE_DIMENSION == 3)
- pos[2] += dir[2] * delta;
-#endif
-
-}
-
-/* Check that the interface fragment is consistent with the current state of
- * the random walk */
-static INLINE int
-XD(check_rwalk_fragment_consistency)
- (const struct XD(rwalk)* rwalk,
- const struct sdis_interface_fragment* frag)
-{
- double N[DIM];
- double uv[2] = {0, 0};
- ASSERT(rwalk && frag);
- dX(normalize)(N, dX_set_fX(N, rwalk->hit.normal));
- if( SXD_HIT_NONE(&rwalk->hit)
- || !dX(eq_eps)(rwalk->vtx.P, frag->P, 1.e-6)
- || !dX(eq_eps)(N, frag->Ng, 1.e-6)
- || !( (IS_INF(rwalk->vtx.time) && IS_INF(frag->time))
- || eq_eps(rwalk->vtx.time, frag->time, 1.e-6))) {
- return 0;
- }
-#if (SDIS_SOLVE_PROBE_DIMENSION == 2)
- uv[0] = rwalk->hit.u;
-#else
- d2_set_f2(uv, rwalk->hit.uv);
-#endif
- return d2_eq_eps(uv, frag->uv, 1.e-6);
-}
-
-res_T
-XD(fluid_temperature)
- (const struct sdis_scene* scn,
- const double fp_to_meter,
- struct XD(rwalk)* rwalk,
- struct ssp_rng* rng,
- struct XD(temperature)* T)
-{
- double tmp;
- (void)rng, (void)fp_to_meter;
- ASSERT(scn && fp_to_meter > 0 && rwalk && rng && T);
- ASSERT(rwalk->mdm->type == SDIS_MEDIUM_FLUID);
-
- tmp = fluid_get_temperature(rwalk->mdm, &rwalk->vtx);
- if(tmp < 0) {
- log_err(scn->dev, "%s: unknown fluid temperature is not supported.\n",
- FUNC_NAME);
- return RES_BAD_OP;
- }
- T->value += tmp;
- T->done = 1;
- return RES_OK;
-}
-
-static void
-XD(solid_solid_boundary_temperature)
- (const struct sdis_scene* scn,
- const double fp_to_meter,
- const struct sdis_interface_fragment* frag,
- struct XD(rwalk)* rwalk,
- struct ssp_rng* rng,
- struct XD(temperature)* T)
-{
- const struct sdis_interface* interf = NULL;
- const struct sdis_medium* solid_front = NULL;
- const struct sdis_medium* solid_back = NULL;
- double lambda_front, lambda_back;
- double delta_front_boundary, delta_back_boundary;
- double delta_front_boundary_meter, delta_back_boundary_meter;
- double delta_boundary;
- double proba;
- double tmp;
- double r;
- float pos[DIM], dir[DIM], range[2];
- ASSERT(scn && fp_to_meter > 0 && frag && rwalk && rng && T);
- ASSERT(XD(check_rwalk_fragment_consistency)(rwalk, frag));
- (void)frag;
-
- /* Retrieve the current boundary media */
- 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);
-
- /* 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);
-
- r = ssp_rng_canonical(rng);
- fX(normalize)(dir, rwalk->hit.normal);
- if(r < proba) { /* Reinject in front */
- delta_boundary = delta_front_boundary;
- rwalk->mdm = solid_front;
- } else { /* Reinject in back */
- delta_boundary = delta_back_boundary;
- fX(minus)(dir, dir);
- rwalk->mdm = solid_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);
-}
-
-static void
-XD(solid_fluid_boundary_temperature)
- (const struct sdis_scene* scn,
- const double fp_to_meter,
- const struct sdis_interface_fragment* frag,
- struct XD(rwalk)* rwalk,
- struct ssp_rng* rng,
- struct XD(temperature)* T)
-{
- const struct sdis_interface* interf = NULL;
- const struct sdis_medium* mdm_front = NULL;
- const struct sdis_medium* mdm_back = NULL;
- const struct sdis_medium* solid = NULL;
- const struct sdis_medium* fluid = NULL;
- double hc;
- double lambda;
- double fluid_proba;
- double delta_boundary;
- double r;
- double tmp;
- float dir[DIM], pos[DIM], range[2];
-
- ASSERT(scn && fp_to_meter > 0 && rwalk && rng && T);
- ASSERT(XD(check_rwalk_fragment_consistency)(rwalk, frag));
-
- /* 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) {
- solid = mdm_front;
- fluid = mdm_back;
- } else {
- solid = mdm_back;
- fluid = mdm_front;
- }
-
- /* Fetch the solid properties */
- lambda = solid_get_thermal_conductivity(solid, &rwalk->vtx);
- delta_boundary = solid_get_delta_boundary(solid, &rwalk->vtx);
- hc = interface_get_convection_coef(interf, frag);
-
- /* Compute the probas to switch in solid or fluid random walk */
- tmp = lambda / (delta_boundary*fp_to_meter);
- fluid_proba = hc / (tmp + hc);
-
- r = ssp_rng_canonical(rng);
- if(r < fluid_proba) { /* Switch to fluid random walk */
- rwalk->mdm = fluid;
- T->func = XD(fluid_temperature);
- } else { /* Solid random walk */
- rwalk->mdm = solid;
- fX(normalize)(dir, rwalk->hit.normal);
- if(solid == mdm_back) fX(minus)(dir, dir);
-
- /* "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);
-
- /* Switch in solid random walk */
- T->func = XD(solid_temperature);
- }
-}
-
-res_T
-XD(boundary_temperature)
- (const struct sdis_scene* scn,
- const double fp_to_meter,
- struct XD(rwalk)* rwalk,
- struct ssp_rng* rng,
- struct XD(temperature)* T)
-{
- struct sdis_interface_fragment frag = SDIS_INTERFACE_FRAGMENT_NULL;
- const struct sdis_interface* interf = NULL;
- const struct sdis_medium* mdm_front = NULL;
- const struct sdis_medium* mdm_back = NULL;
- double tmp;
- ASSERT(scn && fp_to_meter > 0 && rwalk && rng && T);
- ASSERT(!SXD_HIT_NONE(&rwalk->hit));
-
- XD(setup_interface_fragment)(&frag, &rwalk->vtx, &rwalk->hit);
-
- /* 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);
- if(tmp >= 0) {
- T->value += tmp;
- T->done = 1;
- return RES_OK;
- }
-
- mdm_front = interface_get_medium(interf, SDIS_FRONT);
- mdm_back = interface_get_medium(interf, SDIS_BACK);
-
- if(mdm_front->type == mdm_back->type) {
- XD(solid_solid_boundary_temperature)(scn, fp_to_meter, &frag, rwalk, rng, T);
- } else {
- XD(solid_fluid_boundary_temperature)(scn, fp_to_meter, &frag, rwalk, rng, T);
- }
- return RES_OK;
-}
-
-res_T
-XD(solid_temperature)
- (const struct sdis_scene* scn,
- const double fp_to_meter,
- struct XD(rwalk)* rwalk,
- struct ssp_rng* rng,
- struct XD(temperature)* T)
-{
- 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);
-
- /* Check the random walk consistency */
- CHK(scene_get_medium(scn, rwalk->vtx.P, &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_ARG;
- }
- /* Save the submitted position */
- dX(set)(position_start, rwalk->vtx.P);
-
- do { /* Solid random walk */
- struct sXd(hit) hit0, hit1;
- double lambda; /* Thermal conductivity */
- double rho; /* Volumic mass */
- double cp; /* Calorific capacity */
- double tau, mu;
- double tmp;
- float delta, delta_solid; /* Random walk numerical parameter */
- float range[2];
- float dir0[DIM], dir1[DIM];
- float org[DIM];
-
- /* Check the limit condition */
- tmp = solid_get_temperature(mdm, &rwalk->vtx);
- if(tmp >= 0) {
- T->value += tmp;
- T->done = 1;
- return RES_OK;
- }
-
- /* Fetch solid properties */
- delta_solid = (float)solid_get_delta(mdm, &rwalk->vtx);
- lambda = solid_get_thermal_conductivity(mdm, &rwalk->vtx);
- rho = solid_get_volumic_mass(mdm, &rwalk->vtx);
- cp = solid_get_calorific_capacity(mdm, &rwalk->vtx);
-
-#if (SDIS_SOLVE_PROBE_DIMENSION == 2)
- /* Sample a direction around 2PI */
- ssp_ran_circle_uniform_float(rng, dir0, NULL);
-#else
- /* Sample a direction around 4PI */
- ssp_ran_sphere_uniform_float(rng, dir0, NULL);
-#endif
-
- /* Trace a ray along the sampled direction and its opposite to check if a
- * surface is hit in [0, delta_solid]. */
- fX_set_dX(org, rwalk->vtx.P);
- fX(minus)(dir1, dir0);
- hit0 = hit1 = SXD_HIT_NULL;
- range[0] = 0.f, range[1] = delta_solid*RAY_RANGE_MAX_SCALE;
- SXD(scene_view_trace_ray(scn->sXd(view), org, dir0, range, NULL, &hit0));
- SXD(scene_view_trace_ray(scn->sXd(view), org, dir1, range, NULL, &hit1));
-
- if(SXD_HIT_NONE(&hit0) && SXD_HIT_NONE(&hit1)) {
- /* Hit nothing: move along dir0 of the original delta */
- delta = delta_solid;
- } else {
- /* Hit something: move along dir0 of the minimum hit distance */
- delta = MMIN(hit0.distance, hit1.distance);
- }
-
- /* Sample the time */
- mu = (2*DIM*lambda) / (rho*cp*delta*fp_to_meter*delta*fp_to_meter);
- tau = ssp_ran_exp(rng, mu);
- rwalk->vtx.time -= tau;
-
- /* Check the initial condition */
- tmp = solid_get_temperature(mdm, &rwalk->vtx);
- if(tmp >= 0) {
- T->value += tmp;
- T->done = 1;
- return RES_OK;
- }
-
- /* Define if the random walk hits something along dir0 */
- rwalk->hit = hit0.distance > delta ? SXD_HIT_NULL : hit0;
-
- /* 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);
- } 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);
- }
-
- /* 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));
- }
- return RES_BAD_OP;
- }
-
- /* Keep going while the solid random walk does not hit an interface */
- } while(SXD_HIT_NONE(&rwalk->hit));
-
- T->func = XD(boundary_temperature);
- return RES_OK;
-}
-
-static res_T
-XD(compute_temperature)
- (struct sdis_scene* scn,
- const double fp_to_meter,
- struct XD(rwalk)* rwalk,
- struct ssp_rng* rng,
- struct XD(temperature)* T)
-{
-#ifndef NDEBUG
- struct XD(temperature)* stack = NULL;
- size_t istack = 0;
-#endif
- res_T res = RES_OK;
- ASSERT(scn && fp_to_meter && rwalk && rng && T);
-
- do {
- res = T->func(scn, fp_to_meter, rwalk, rng, T);
- if(res != RES_OK) goto error;
-
-#ifndef NDEBUG
- sa_push(stack, *T);
- ++istack;
-#endif
- } while(!T->done);
-
-exit:
-#ifndef NDEBUG
- sa_release(stack);
-#endif
- return res;
-error:
- goto exit;
-}
-
-static res_T
-XD(probe_realisation)
- (struct sdis_scene* scn,
- struct ssp_rng* rng,
- const struct sdis_medium* medium,
- const double position[],
- const double time,
- const double fp_to_meter,/* Scale factor from floating point unit to meter */
- double* weight)
-{
- struct XD(rwalk) rwalk = XD(RWALK_NULL);
- struct XD(temperature) T = XD(TEMPERATURE_NULL);
- res_T res = RES_OK;
- 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;
- default: FATAL("Unreachable code\n"); break;
- }
-
- dX(set)(rwalk.vtx.P, position);
- rwalk.vtx.time = time;
- rwalk.hit = SXD_HIT_NULL;
- rwalk.mdm = medium;
-
- res = XD(compute_temperature)(scn, fp_to_meter, &rwalk, rng, &T);
- if(res != RES_OK) return res;
-
- *weight = T.value;
- return RES_OK;
-}
-
-
-#undef SDIS_SOLVE_PROBE_DIMENSION
-#undef DIM
-#undef sXd
-#undef SXD_HIT_NONE
-#undef SXD_HIT_NULL
-#undef SXD_HIT_NULL__
-#undef SXD
-#undef dX
-#undef fX
-#undef fX_set_dX
-#undef XD
-
-#endif /* !SDIS_SOLVE_PROBE_DIMENSION */
-
diff --git a/src/test_sdis_accum_buffer.c b/src/test_sdis_accum_buffer.c
@@ -0,0 +1,77 @@
+/* Copyright (C) |Meso|Star> 2016-2018 (contact@meso-star.com)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. */
+
+#include "test_sdis_utils.h"
+#include <string.h>
+
+int
+main(int argc, char** argv)
+{
+ struct mem_allocator allocator;
+ struct sdis_accum_buffer* buf = NULL;
+ struct sdis_device* dev = NULL;
+ struct sdis_accum_buffer_layout layout;
+ const struct sdis_accum* accums = NULL;
+ struct sdis_accum* accums_tmp = NULL;
+ (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);
+
+ CHK(sdis_accum_buffer_create(NULL, 4 ,4, &buf) == RES_BAD_ARG);
+ CHK(sdis_accum_buffer_create(dev, 0 ,4, &buf) == RES_BAD_ARG);
+ CHK(sdis_accum_buffer_create(dev, 4 ,0, &buf) == RES_BAD_ARG);
+ CHK(sdis_accum_buffer_create(dev, 4 ,0, NULL) == RES_BAD_ARG);
+ CHK(sdis_accum_buffer_create(dev, 4 ,4, &buf) == RES_OK);
+
+ CHK(sdis_accum_buffer_ref_get(NULL) == RES_BAD_ARG);
+ CHK(sdis_accum_buffer_ref_get(buf) == RES_OK);
+ CHK(sdis_accum_buffer_ref_put(NULL) == RES_BAD_ARG);
+ CHK(sdis_accum_buffer_ref_put(buf) == RES_OK);
+ CHK(sdis_accum_buffer_ref_put(buf) == RES_OK);
+
+ CHK(sdis_accum_buffer_create(dev, 16, 8, &buf) == RES_OK);
+
+ CHK(sdis_accum_buffer_get_layout(NULL, &layout) == RES_BAD_ARG);
+ CHK(sdis_accum_buffer_get_layout(buf, NULL) == RES_BAD_ARG);
+ CHK(sdis_accum_buffer_get_layout(buf, &layout) == RES_OK);
+
+ CHK(layout.width == 16);
+ CHK(layout.height == 8);
+
+ CHK(sdis_accum_buffer_map(NULL, &accums) == RES_BAD_ARG);
+ CHK(sdis_accum_buffer_map(buf, NULL) == RES_BAD_ARG);
+ CHK(sdis_accum_buffer_map(buf, &accums) == RES_OK);
+
+ /* Check the accessibility to the mapped data */
+ accums_tmp = MEM_CALLOC
+ (&allocator, layout.width*layout.height, sizeof(struct sdis_accum));
+ CHK(accums_tmp != NULL);
+ memcpy(accums_tmp, accums_tmp,
+ layout.width*layout.height*sizeof(struct sdis_accum));
+ MEM_RM(&allocator, accums_tmp);
+
+ CHK(sdis_accum_buffer_unmap(NULL) == RES_BAD_ARG);
+ CHK(sdis_accum_buffer_unmap(buf) == RES_OK);
+
+ CHK(sdis_accum_buffer_ref_put(buf) == 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_camera.c b/src/test_sdis_camera.c
@@ -0,0 +1,94 @@
+/* Copyright (C) |Meso|Star> 2016-2018 (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>
+
+int
+main(int argc, char** argv)
+{
+ struct sdis_device* dev;
+ struct sdis_camera* cam;
+ struct mem_allocator allocator;
+ double pos[3] = {0};
+ double tgt[3] = {0};
+ double up[3] = {0};
+ (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);
+
+ CHK(sdis_camera_create(NULL, NULL) == RES_BAD_ARG);
+ CHK(sdis_camera_create(dev, NULL) == RES_BAD_ARG);
+ CHK(sdis_camera_create(NULL, &cam) == RES_BAD_ARG);
+ CHK(sdis_camera_create(dev, &cam) == RES_OK);
+
+ CHK(sdis_camera_ref_get(NULL) == RES_BAD_ARG);
+ CHK(sdis_camera_ref_get(cam) == RES_OK);
+ CHK(sdis_camera_ref_put(NULL) == RES_BAD_ARG);
+ CHK(sdis_camera_ref_put(cam) == RES_OK);
+ CHK(sdis_camera_ref_put(cam) == RES_OK);
+
+ CHK(sdis_camera_create(dev, &cam) == RES_OK);
+ CHK(sdis_camera_set_proj_ratio(NULL, 0) == RES_BAD_ARG);
+ CHK(sdis_camera_set_proj_ratio(cam, 0) == RES_BAD_ARG);
+ CHK(sdis_camera_set_proj_ratio(NULL, 4.0/3.0) == RES_BAD_ARG);
+ CHK(sdis_camera_set_proj_ratio(cam, 4.0/3.0) == RES_OK);
+ CHK(sdis_camera_set_proj_ratio(cam, -4.0/3.0) == RES_BAD_ARG);
+
+ CHK(sdis_camera_set_fov(NULL, 0) == RES_BAD_ARG);
+ CHK(sdis_camera_set_fov(cam, 0) == RES_BAD_ARG);
+ CHK(sdis_camera_set_fov(NULL, PI/4.0) == RES_BAD_ARG);
+ CHK(sdis_camera_set_fov(cam, PI/4.0) == RES_OK);
+ CHK(sdis_camera_set_fov(cam, -PI/4.0) == RES_BAD_ARG);
+
+ pos[0] = 0, pos[1] = 0, pos[2] = 0;
+ tgt[0] = 0, tgt[1] = 0, tgt[2] = -1;
+ up[0] = 0, up[1] = 1, up[2] = 0;
+ CHK(sdis_camera_look_at(NULL, NULL, NULL, NULL) == RES_BAD_ARG);
+ CHK(sdis_camera_look_at(cam, NULL, NULL, NULL) == RES_BAD_ARG);
+ CHK(sdis_camera_look_at(NULL, pos, NULL, NULL) == RES_BAD_ARG);
+ CHK(sdis_camera_look_at(cam, pos, NULL, NULL) == RES_BAD_ARG);
+ CHK(sdis_camera_look_at(NULL, NULL, tgt, NULL) == RES_BAD_ARG);
+ CHK(sdis_camera_look_at(cam, NULL, tgt, NULL) == RES_BAD_ARG);
+ CHK(sdis_camera_look_at(NULL, pos, tgt, NULL) == RES_BAD_ARG);
+ CHK(sdis_camera_look_at(cam, pos, tgt, NULL) == RES_BAD_ARG);
+ CHK(sdis_camera_look_at(NULL, NULL, NULL, up) == RES_BAD_ARG);
+ CHK(sdis_camera_look_at(cam, NULL, NULL, up) == RES_BAD_ARG);
+ CHK(sdis_camera_look_at(NULL, pos, NULL, up) == RES_BAD_ARG);
+ CHK(sdis_camera_look_at(cam, pos, NULL, up) == RES_BAD_ARG);
+ CHK(sdis_camera_look_at(NULL, NULL, tgt, up) == RES_BAD_ARG);
+ CHK(sdis_camera_look_at(cam, NULL, tgt, up) == RES_BAD_ARG);
+ CHK(sdis_camera_look_at(NULL, pos, tgt, up) == RES_BAD_ARG);
+ CHK(sdis_camera_look_at(cam, pos, tgt, up) == RES_OK);
+ tgt[0] = 0, tgt[1] = 0, tgt[2] = 0;
+ CHK(sdis_camera_look_at(cam, pos, tgt, up) == RES_BAD_ARG);
+ tgt[0] = 0, tgt[1] = 0, tgt[2] = -1;
+ up[0] = 0, up[1] = 0, up[2] = 0;
+ CHK(sdis_camera_look_at(cam, pos, tgt, up) == RES_BAD_ARG);
+
+ CHK(sdis_device_ref_put(dev) == RES_OK);
+ CHK(sdis_camera_ref_put(cam) == RES_OK);
+
+ check_memory_allocator(&allocator);
+ mem_shutdown_proxy_allocator(&allocator);
+ CHK(mem_allocated_size() == 0);
+ return 0;
+}
+
diff --git a/src/test_sdis_conducto_radiative.c b/src/test_sdis_conducto_radiative.c
@@ -0,0 +1,430 @@
+/* Copyright (C) |Meso|Star> 2016-2018 (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>
+#include <star/ssp.h>
+
+#define UNKNOWN_TEMPERATURE -1
+
+/* The scene is composed of a solid cube whose temperature is unknown. The cube
+ * faces on +/-X are in contact with a fluid and their convection coefficient
+ * is null while their emissivity is 1. The left and right fluids are enclosed
+ * by surfaces whose emissivity are null excepted for the faces orthogonal to
+ * the X axis that are fully emissive and whose temperature is known. The
+ * medium that surrounds the solid cube and the 2 fluids is a solid with a null
+ * conductivity.
+ *
+ * Y (1, 1, 1)
+ * | +------+----------+------+ (1.5,1,1)
+ * o--- X /' /##########/' /|
+ * / +------+----------+------+ |
+ * Z | ' |##########|*' | | 310K
+ * | ' |##########|*' | |
+ * 300K | ' E=1|##########|*'E=1 | |
+ * | +....|##########|*+....|.+
+ * |/ |##########|/ |/
+ * (-1.5,-1,-1) +------+----------+------+
+ * (-1,-1,-1)
+ */
+
+/*******************************************************************************
+ * Geometry
+ ******************************************************************************/
+struct geometry {
+ const double* positions;
+ const size_t* indices;
+ struct sdis_interface** interfaces;
+};
+
+static const double vertices[16/*#vertices*/*3/*#coords per vertex*/] = {
+ -1.0,-1.0,-1.0,
+ 1.0,-1.0,-1.0,
+ -1.0, 1.0,-1.0,
+ 1.0, 1.0,-1.0,
+ -1.0,-1.0, 1.0,
+ 1.0,-1.0, 1.0,
+ -1.0, 1.0, 1.0,
+ 1.0, 1.0, 1.0,
+ -1.5,-1.0,-1.0,
+ 1.5,-1.0,-1.0,
+ -1.5, 1.0,-1.0,
+ 1.5, 1.0,-1.0,
+ -1.5,-1.0, 1.0,
+ 1.5,-1.0, 1.0,
+ -1.5, 1.0, 1.0,
+ 1.5, 1.0, 1.0,
+};
+static const size_t nvertices = sizeof(vertices) / sizeof(double[3]);
+
+static const size_t indices[32/*#triangles*/*3/*#indices per triangle*/] = {
+ 0, 2, 1, 1, 2, 3, /* Solid back face */
+ 0, 4, 2, 2, 4, 6, /* Solid left face*/
+ 4, 5, 6, 6, 5, 7, /* Solid front face */
+ 3, 7, 1, 1, 7, 5, /* Solid right face */
+ 2, 6, 3, 3, 6, 7, /* Solid top face */
+ 0, 1, 4, 4, 1, 5, /* Solid bottom face */
+
+ 8, 10, 0, 0, 10, 2, /* Left fluid back face */
+ 8, 12, 10, 10, 12, 14, /* Left fluid left face */
+ 12, 4, 14, 14, 4, 6, /* Left fluid front face */
+ 10, 14, 2, 2, 14, 6, /* Left fluid top face */
+ 8, 0, 12, 12, 0, 4, /* Left fluid bottom face */
+
+ 1, 3, 9, 9, 3, 11, /* Right fluid back face */
+ 5, 13, 7, 7, 13, 15, /* Right fluid front face */
+ 11, 15, 9, 9, 15, 13, /* Right fluid right face */
+ 3, 7, 11, 11, 7, 15, /* Right fluid top face */
+ 1, 9, 5, 5, 9, 13 /* Right fluid bottom face */
+};
+static const size_t ntriangles = sizeof(indices) / sizeof(size_t[3]);
+
+static void
+get_indices(const size_t itri, size_t ids[3], void* ctx)
+{
+ struct geometry* geom = ctx;
+ CHK(ctx != NULL);
+ ids[0] = geom->indices[itri*3+0];
+ ids[1] = geom->indices[itri*3+1];
+ ids[2] = geom->indices[itri*3+2];
+}
+
+static void
+get_position(const size_t ivert, double pos[3], void* ctx)
+{
+ struct geometry* geom = ctx;
+ CHK(ctx != NULL);
+ pos[0] = geom->positions[ivert*3+0];
+ pos[1] = geom->positions[ivert*3+1];
+ pos[2] = geom->positions[ivert*3+2];
+}
+
+static void
+get_interface(const size_t itri, struct sdis_interface** bound, void* ctx)
+{
+ struct geometry* geom = ctx;
+ CHK(ctx != NULL);
+ *bound = geom->interfaces[itri];
+}
+
+/*******************************************************************************
+ * Media
+ ******************************************************************************/
+struct solid {
+ double lambda;
+};
+
+static double
+temperature_unknown(const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ CHK(vtx != NULL); (void)data;
+ return -1;
+}
+
+static double
+solid_get_calorific_capacity
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ CHK(vtx != NULL); (void)data;
+ return 1;
+}
+
+static double
+solid_get_thermal_conductivity
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ CHK(vtx != NULL);
+ CHK(data != 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(vtx != NULL); (void)data;
+ return 1;
+}
+
+static double
+solid_get_delta
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ CHK(vtx != NULL); (void)data;
+ 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;
+ double convection_coef;
+ double emissivity;
+ double specular_fraction;
+};
+
+static double
+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;
+}
+
+static double
+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;
+}
+
+static double
+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;
+}
+
+static double
+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;
+}
+
+/*******************************************************************************
+ * Helper functions
+ ******************************************************************************/
+static void
+create_interface
+ (struct sdis_device* dev,
+ struct sdis_medium* front,
+ struct sdis_medium* back,
+ const struct interface* interf,
+ struct sdis_interface** out_interf)
+{
+ struct sdis_interface_shader shader = DUMMY_INTERFACE_SHADER;
+ 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;
+
+ CHK(sdis_data_create(dev, sizeof(struct interface), ALIGNOF(struct interface),
+ NULL, &data) == RES_OK);
+ *((struct interface*)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
+ ******************************************************************************/
+int
+main(int argc, char** argv)
+{
+ struct mem_allocator allocator;
+ struct interface interf;
+ struct geometry geom;
+ struct sdis_data* data = NULL;
+ struct sdis_device* dev = NULL;
+ struct sdis_medium* fluid = NULL;
+ struct sdis_medium* solid = NULL;
+ struct sdis_medium* solid2 = NULL;
+ struct sdis_interface* interfaces[5] = {NULL};
+ struct sdis_interface* prim_interfaces[32/*#triangles*/];
+ struct sdis_fluid_shader fluid_shader = DUMMY_FLUID_SHADER;
+ struct sdis_solid_shader solid_shader = DUMMY_SOLID_SHADER;
+ struct sdis_scene* scn = NULL;
+ struct ssp_rng* rng = NULL;
+ const size_t nsimuls = 4;
+ size_t isimul;
+ const double emissivity = 1;/* Emissivity of the side +/-X of the solid */
+ const double lambda = 0.1; /* Conductivity of the solid */
+ const double Tref = 300; /* Reference temperature */
+ const double T0 = 300; /* Fixed temperature on the left side of the system */
+ const double T1 = 310; /* Fixed temperature on the right side of the system */
+ const double thickness = 2.0; /* Thickness of the solid along X */
+ double Ts0, Ts1, hr, tmp;
+ (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 = temperature_unknown;
+ CHK(sdis_fluid_create(dev, &fluid_shader, NULL, &fluid) == RES_OK);
+
+ /* Create the solid medium */
+ CHK(sdis_data_create(dev, sizeof(struct solid), ALIGNOF(struct solid),
+ NULL, &data) == RES_OK);
+ ((struct solid*)sdis_data_get(data))->lambda = lambda;
+ 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 = temperature_unknown;
+ CHK(sdis_solid_create(dev, &solid_shader, data, &solid) == RES_OK);
+ CHK(sdis_data_ref_put(data) == RES_OK);
+
+ /* Create the surrounding solid medium */
+ CHK(sdis_data_create(dev, sizeof(struct solid), ALIGNOF(struct solid),
+ NULL, &data) == RES_OK);
+ ((struct solid*)sdis_data_get(data))->lambda = 0;
+ 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 = temperature_unknown;
+ 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;
+ 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;
+ 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;
+ 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;
+ 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;
+ create_interface(dev, fluid, solid2, &interf, interfaces+4);
+
+ /* Setup the per primitive interface of the solid medium */
+ prim_interfaces[0] = prim_interfaces[1] = interfaces[0];
+ prim_interfaces[2] = prim_interfaces[3] = interfaces[1];
+ prim_interfaces[4] = prim_interfaces[5] = interfaces[0];
+ prim_interfaces[6] = prim_interfaces[7] = interfaces[1];
+ prim_interfaces[8] = prim_interfaces[9] = interfaces[0];
+ prim_interfaces[10] = prim_interfaces[11] = interfaces[0];
+
+ /* Setup the per primitive interface of the fluid on the left of the medium */
+ prim_interfaces[12] = prim_interfaces[13] = interfaces[2];
+ prim_interfaces[14] = prim_interfaces[15] = interfaces[3];
+ prim_interfaces[16] = prim_interfaces[17] = interfaces[2];
+ prim_interfaces[18] = prim_interfaces[19] = interfaces[2];
+ prim_interfaces[20] = prim_interfaces[21] = interfaces[2];
+
+ /* Setup the per primitive interface of the fluid on the right of the medium */
+ prim_interfaces[22] = prim_interfaces[23] = interfaces[2];
+ prim_interfaces[24] = prim_interfaces[25] = interfaces[2];
+ prim_interfaces[26] = prim_interfaces[27] = interfaces[4];
+ prim_interfaces[28] = prim_interfaces[29] = interfaces[2];
+ prim_interfaces[30] = prim_interfaces[31] = interfaces[2];
+
+ /* Create the scene */
+ geom.positions = vertices;
+ geom.indices = indices;
+ geom.interfaces = prim_interfaces;
+ CHK(sdis_scene_create(dev, ntriangles, get_indices, get_interface, nvertices,
+ get_position, &geom, &scn) == RES_OK);
+
+ hr = 4.0 * BOLTZMANN_CONSTANT * Tref*Tref*Tref * emissivity;
+ tmp = lambda/(2*lambda + thickness*hr) * (T1 - T0);
+ Ts0 = T0 + tmp;
+ Ts1 = T1 - tmp;
+
+ /* Run the simulations */
+ CHK(ssp_rng_create(&allocator, &ssp_rng_kiss, &rng) == RES_OK);
+ FOR_EACH(isimul, 0, nsimuls) {
+ struct sdis_mc T = SDIS_MC_NULL;
+ struct sdis_estimator* estimator;
+ double pos[3];
+ double ref, u;
+ size_t nreals = 0;
+ size_t nfails = 0;
+
+ 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_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);
+
+ u = (pos[0] + 1) / thickness;
+ ref = u * Ts1 + (1-u) * Ts0;
+ printf("Temperature at (%g, %g, %g) = %g ~ %g +/- %g\n",
+ SPLIT3(pos), ref, T.E, T.SE);
+
+ CHK(eq_eps(T.E, ref, 2*T.SE) == 1);
+
+ CHK(sdis_estimator_ref_put(estimator) == RES_OK);
+ }
+
+ /* Release memory */
+ CHK(sdis_scene_ref_put(scn) == RES_OK);
+ CHK(sdis_interface_ref_put(interfaces[0]) == RES_OK);
+ CHK(sdis_interface_ref_put(interfaces[1]) == RES_OK);
+ CHK(sdis_interface_ref_put(interfaces[2]) == RES_OK);
+ CHK(sdis_interface_ref_put(interfaces[3]) == RES_OK);
+ CHK(sdis_interface_ref_put(interfaces[4]) == RES_OK);
+ CHK(sdis_medium_ref_put(fluid) == RES_OK);
+ CHK(sdis_medium_ref_put(solid) == RES_OK);
+ CHK(sdis_medium_ref_put(solid2) == RES_OK);
+ CHK(sdis_device_ref_put(dev) == RES_OK);
+ CHK(ssp_rng_ref_put(rng) == RES_OK);
+
+ check_memory_allocator(&allocator);
+ mem_shutdown_proxy_allocator(&allocator);
+ CHK(mem_allocated_size() == 0);
+ return 0;
+}
diff --git a/src/test_sdis_conducto_radiative_2d.c b/src/test_sdis_conducto_radiative_2d.c
@@ -0,0 +1,405 @@
+/* Copyright (C) |Meso|Star> 2016-2018 (contact@meso-star.com)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. */
+
+#include "sdis.h"
+#include "test_sdis_utils.h"
+
+#include <star/ssp.h>
+
+#define UNKNOWN_TEMPERATURE -1
+
+/* The scene is composed of a solid square whose temperature is unknown. The
+ * square segments on +/-X are in contact with a fluid and their convection
+ * coefficient is null while their emissivity is 1. The left and right fluids
+ * are enclosed by segments whose emissivity are null excepted for the segments
+ * orthogonal to the X axis that are fully emissive and whose temperature is
+ * known. The medium that surrounds the solid square and the 2 fluids is a
+ * solid with a null conductivity.
+ *
+ * (1, 1)
+ * +-----+----------+-----+ (1.5,1,1)
+ * | |##########| |
+ * | |##########| |
+ * 300K | E=1 |##########| E=1 | 310K
+ * | |##########| |
+ * | |##########| |
+ * (-1.5,-1) +-----+----------+-----+
+ * (-1,-1)
+ */
+
+/*******************************************************************************
+ * Geometry
+ ******************************************************************************/
+struct geometry {
+ const double* positions;
+ const size_t* indices;
+ struct sdis_interface** interfaces;
+};
+
+static const double vertices[8/*#vertices*/*2/*#coords par vertex*/] = {
+ 1.0, -1.0,
+ -1.0, -1.0,
+ -1.0, 1.0,
+ 1.0, 1.0,
+ 1.5, -1.0,
+ -1.5, -1.0,
+ -1.5, 1.0,
+ 1.5, 1.0
+};
+static const size_t nvertices = sizeof(vertices) / sizeof(double[2]);
+
+static const size_t indices[10/*#segments*/*2/*#indices per segment*/] = {
+ 0, 1, /* Solid bottom segment */
+ 1, 2, /* Solid left segment */
+ 2, 3, /* Solid top segment */
+ 3, 0, /* Solid right segment */
+
+ 1, 5, /* Left fluid bottom segment */
+ 5, 6, /* Left fluid left segment */
+ 6, 2, /* Left fluid top segment */
+
+ 4, 0, /* Right fluid bottom segment */
+ 3, 7, /* Right fluid top segment */
+ 7, 4 /* Right fluid right segment */
+};
+static const size_t nsegments = sizeof(indices) / sizeof(size_t[2]);
+
+static void
+get_indices(const size_t iseg, size_t ids[2], void* ctx)
+{
+ struct geometry* geom = ctx;
+ CHK(ctx != NULL);
+ ids[0] = geom->indices[iseg*2+0];
+ ids[1] = geom->indices[iseg*2+1];
+}
+
+static void
+get_position(const size_t ivert, double pos[2], void* ctx)
+{
+ struct geometry* geom = ctx;
+ CHK(ctx != NULL);
+ pos[0] = geom->positions[ivert*2+0];
+ pos[1] = geom->positions[ivert*2+1];
+}
+
+static void
+get_interface(const size_t iseg, struct sdis_interface** bound, void* ctx)
+{
+ struct geometry* geom = ctx;
+ CHK(ctx != NULL);
+ *bound = geom->interfaces[iseg];
+}
+
+/*******************************************************************************
+ * Media
+ ******************************************************************************/
+struct solid {
+ double lambda;
+};
+
+static double
+temperature_unknown(const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ CHK(vtx != NULL); (void)data;
+ return -1;
+}
+
+static double
+solid_get_calorific_capacity
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ CHK(vtx != NULL); (void)data;
+ return 1;
+}
+
+static double
+solid_get_thermal_conductivity
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ CHK(vtx != NULL);
+ CHK(data != 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(vtx != NULL); (void)data;
+ return 1;
+}
+
+static double
+solid_get_delta
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ CHK(vtx != NULL); (void)data;
+ 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;
+ double convection_coef;
+ double emissivity;
+ double specular_fraction;
+};
+
+static double
+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;
+}
+
+static double
+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;
+}
+
+static double
+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;
+}
+
+static double
+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;
+}
+
+/*******************************************************************************
+ * Helper functions
+ ******************************************************************************/
+static void
+create_interface
+ (struct sdis_device* dev,
+ struct sdis_medium* front,
+ struct sdis_medium* back,
+ const struct interface* interf,
+ struct sdis_interface** out_interf)
+{
+ struct sdis_interface_shader shader = DUMMY_INTERFACE_SHADER;
+ 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;
+
+ CHK(sdis_data_create(dev, sizeof(struct interface), ALIGNOF(struct interface),
+ NULL, &data) == RES_OK);
+ *((struct interface*)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
+ ******************************************************************************/
+int
+main(int argc, char** argv)
+{
+ struct mem_allocator allocator;
+ struct interface interf;
+ struct geometry geom;
+ struct ssp_rng* rng = NULL;
+ struct sdis_scene* scn = NULL;
+ struct sdis_data* data = NULL;
+ struct sdis_device* dev = NULL;
+ struct sdis_medium* fluid = NULL;
+ struct sdis_medium* solid = NULL;
+ struct sdis_medium* solid2 = NULL;
+ struct sdis_interface* interfaces[5] = {NULL};
+ struct sdis_interface* prim_interfaces[10/*#segment*/];
+ struct sdis_fluid_shader fluid_shader = DUMMY_FLUID_SHADER;
+ struct sdis_solid_shader solid_shader = DUMMY_SOLID_SHADER;
+ const size_t nsimuls = 4;
+ size_t isimul;
+ const double emissivity = 1;/* Emissivity of the side +/-X of the solid */
+ const double lambda = 0.1; /* Conductivity of the solid */
+ const double Tref = 300; /* Reference temperature */
+ const double T0 = 300; /* Fixed temperature on the left side of the system */
+ const double T1 = 310; /* Fixed temperature on the right side of the system */
+ const double thickness = 2.0; /* Thickness of the solid along X */
+ double Ts0, Ts1, hr, tmp;
+ (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 = temperature_unknown;
+ CHK(sdis_fluid_create(dev, &fluid_shader, NULL, &fluid) == RES_OK);
+
+ /* Create the solid medium */
+ CHK(sdis_data_create
+ (dev, sizeof(struct solid), ALIGNOF(struct solid), NULL, &data) == RES_OK);
+ ((struct solid*)sdis_data_get(data))->lambda = lambda;
+ 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 = temperature_unknown;
+ CHK(sdis_solid_create(dev, &solid_shader, data, &solid) == RES_OK);
+ CHK(sdis_data_ref_put(data) == RES_OK);
+
+ /* Create the surrounding solid medium */
+ CHK(sdis_data_create(dev, sizeof(struct solid), ALIGNOF(struct solid),
+ NULL, &data) == RES_OK);
+ ((struct solid*)sdis_data_get(data))->lambda = 0;
+ solid_shader.calorific_capacity = solid_get_thermal_conductivity;
+ 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;
+ 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;
+ 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;
+ 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;
+ 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;
+ create_interface(dev, fluid, solid2, &interf, interfaces+4);
+
+ /* Setup the per primitive interface of the solid medium */
+ prim_interfaces[0] = interfaces[0];
+ prim_interfaces[1] = interfaces[1];
+ prim_interfaces[2] = interfaces[0];
+ prim_interfaces[3] = interfaces[1];
+
+ /* Setup the per primitive interface of the fluid on the left of the medium */
+ prim_interfaces[4] = interfaces[2];
+ prim_interfaces[5] = interfaces[3];
+ prim_interfaces[6] = interfaces[2];
+
+ /* Setup the per primitive interface of the fluid on the right of the medium */
+ prim_interfaces[7] = interfaces[2];
+ prim_interfaces[8] = interfaces[2];
+ prim_interfaces[9] = interfaces[4];
+
+ /* Create the scene */
+ geom.positions = vertices;
+ geom.indices = indices;
+ geom.interfaces = prim_interfaces;
+ CHK(sdis_scene_2d_create(dev, nsegments, get_indices, get_interface, nvertices,
+ get_position, &geom, &scn) == RES_OK);
+
+ hr = 4*BOLTZMANN_CONSTANT * Tref*Tref*Tref * emissivity;
+ tmp = lambda/(2*lambda + thickness*hr) * (T1 - T0);
+ Ts0 = T0 + tmp;
+ Ts1 = T1 - tmp;
+
+ /* Run the simulations */
+ CHK(ssp_rng_create(&allocator, &ssp_rng_kiss, &rng) == RES_OK);
+ FOR_EACH(isimul, 0, nsimuls) {
+ struct sdis_mc T = SDIS_MC_NULL;
+ struct sdis_estimator* estimator;
+ double pos[2];
+ double ref, u;
+ size_t nreals = 0;
+ size_t nfails = 0;
+
+ pos[0] = ssp_rng_uniform_double(rng, -0.9, 0.9);
+ pos[1] = 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_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);
+
+ u = (pos[0] + 1) / thickness;
+ ref = u * Ts1 + (1-u) * Ts0;
+ printf("Temperature at (%g, %g) = %g ~ %g +/- %g\n",
+ SPLIT2(pos), ref, T.E, T.SE);
+
+ CHK(eq_eps(T.E, ref, 2*T.SE) == 1);
+
+ CHK(sdis_estimator_ref_put(estimator) == RES_OK);
+ }
+
+ /* Release memory */
+ CHK(sdis_scene_ref_put(scn) == RES_OK);
+ CHK(sdis_interface_ref_put(interfaces[0]) == RES_OK);
+ CHK(sdis_interface_ref_put(interfaces[1]) == RES_OK);
+ CHK(sdis_interface_ref_put(interfaces[2]) == RES_OK);
+ CHK(sdis_interface_ref_put(interfaces[3]) == RES_OK);
+ CHK(sdis_interface_ref_put(interfaces[4]) == RES_OK);
+ CHK(sdis_medium_ref_put(fluid) == RES_OK);
+ CHK(sdis_medium_ref_put(solid) == RES_OK);
+ CHK(sdis_medium_ref_put(solid2) == RES_OK);
+ CHK(ssp_rng_ref_put(rng) == 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
@@ -76,8 +76,11 @@ main(int argc, char** argv)
CHK(sdis_interface_ref_put(interf) == RES_OK);
CHK(sdis_interface_ref_put(interf) == RES_OK);
- CHK(CREATE(dev, solid, solid, &shader, NULL, &interf) == RES_BAD_ARG);
+ 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;
CHK(CREATE(dev, solid, solid, &shader, NULL, &interf) == RES_OK);
CHK(sdis_interface_ref_put(interf) == RES_OK);
@@ -87,6 +90,10 @@ main(int argc, char** argv)
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;
CHK(CREATE(dev, solid, fluid, &shader, NULL, &interf) == RES_OK);
CHK(sdis_interface_ref_put(interf) == RES_OK);
#undef CREATE
diff --git a/src/test_sdis_solve_camera.c b/src/test_sdis_solve_camera.c
@@ -0,0 +1,623 @@
+/* Copyright (C) |Meso|Star> 2016-2018 (contact@meso-star.com)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. */
+
+#include "sdis.h"
+#include "test_sdis_utils.h"
+
+#include <star/s3dut.h>
+
+#include <rsys/algorithm.h>
+#include <rsys/image.h>
+#include <rsys/double3.h>
+#include <rsys/double33.h>
+#include <rsys/math.h>
+#include <rsys/stretchy_array.h>
+
+#include <string.h>
+
+#define UNKOWN_TEMPERATURE -1
+#define IMG_WIDTH 640
+#define IMG_HEIGHT 480
+#define SPP 4 /* #Samples per pixel, i.e. #realisations per pixel */
+
+/*
+ * The scene is composed of a solid cube whose temperature is unknown. The
+ * emissivity of the cube is 1 and its convection coefficient with the
+ * surrounding fluid is null. At the center of the cube there is a spherical
+ * fluid cavity whose temperature is 350K. The convection coefficient between
+ * the solid and the cavity is 1 and the emissivity of this interface is null.
+ * The ambient radiative temperature of the system is 300K.
+ *
+ * In this test, we compute the radiative temperature that reaches a camera
+ * that looks the cube and we `dump' a normalized image of the result.
+ *
+ * (1,1,1)
+ * +----------------+
+ * /' # # /|
+ * +----*--------*--+ |
+ * | ' # # | |
+ * | ' # 350K # | |
+ * | ' # # | |
+ * | +.....#..#.....|.+
+ * |/ |/
+ * +----------------+
+ * (-1,-1,-1)
+ */
+
+/*******************************************************************************
+ * Geometry
+ ******************************************************************************/
+struct map_interf {
+ size_t key;
+ struct sdis_interface* interf;
+};
+
+static int
+cmp_map_inter(const void* key, const void* elmt)
+{
+ const size_t* iprim = key;
+ const struct map_interf* interf = elmt;
+ if(*iprim < interf->key) return -1;
+ else if(*iprim > interf->key) return 1;
+ else return 0;
+}
+
+struct geometry {
+ double* positions;
+ size_t* indices;
+ struct map_interf* interfaces;
+};
+static const struct geometry GEOMETRY_NULL = {NULL, NULL, NULL};
+
+static void
+geometry_release(struct geometry* geom)
+{
+ CHK(geom != NULL);
+ sa_release(geom->positions);
+ sa_release(geom->indices);
+ sa_release(geom->interfaces);
+ *geom = GEOMETRY_NULL;
+}
+
+static void
+geometry_get_indices(const size_t itri, size_t ids[3], void* ctx)
+{
+ struct geometry* geom = ctx;
+
+ CHK(ids != NULL);
+ CHK(geom != NULL);
+ CHK(itri < sa_size(geom->indices)/3/*#indices per triangle*/);
+
+ /* Fetch the indices */
+ ids[0] = geom->indices[itri*3+0];
+ ids[1] = geom->indices[itri*3+1];
+ ids[2] = geom->indices[itri*3+2];
+}
+
+static void
+geometry_get_position(const size_t ivert, double pos[3], void* ctx)
+{
+ struct geometry* geom = ctx;
+
+ CHK(pos != NULL);
+ CHK(geom != NULL);
+ CHK(ivert < sa_size(geom->positions)/3/*#coords per triangle*/);
+
+ /* Fetch the vertices */
+ pos[0] = geom->positions[ivert*3+0];
+ pos[1] = geom->positions[ivert*3+1];
+ pos[2] = geom->positions[ivert*3+2];
+}
+
+static void
+geometry_get_interface
+ (const size_t itri,
+ struct sdis_interface** bound,
+ void* ctx)
+{
+ struct geometry* geom = ctx;
+ struct map_interf* interf = NULL;
+
+ CHK(bound != NULL);
+ CHK(geom != NULL);
+ CHK(itri < sa_size(geom->indices)/3/*#indices per triangle*/);
+
+ /* Find the interface of the triangle */
+ interf = search_lower_bound(&itri, geom->interfaces,
+ sa_size(geom->interfaces), sizeof(struct map_interf), cmp_map_inter);
+
+ CHK(interf != NULL);
+ CHK(interf->key >= itri);
+ *bound = interf->interf;
+}
+
+/*******************************************************************************
+ * Fluid medium
+ ******************************************************************************/
+struct fluid {
+ double cp;
+ double rho;
+ double temperature;
+};
+static const struct fluid FLUID_NULL = {0, 0, UNKOWN_TEMPERATURE};
+
+static double
+fluid_get_calorific_capacity
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ CHK(data != NULL && vtx != NULL);
+ return ((const struct fluid*)sdis_data_cget(data))->cp;
+}
+
+static double
+fluid_get_volumic_mass
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ CHK(data != NULL && vtx != NULL);
+ return ((const struct fluid*)sdis_data_cget(data))->rho;
+}
+
+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;
+}
+
+/*******************************************************************************
+ * Solid medium
+ ******************************************************************************/
+struct solid {
+ double cp;
+ double lambda;
+ double rho;
+ double delta;
+ double temperature;
+};
+static const struct solid SOLID_NULL = {0, 0, 0, 0, UNKOWN_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_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)
+{
+ CHK(data != NULL && vtx != NULL);
+ return ((const struct solid*)sdis_data_cget(data))->temperature;
+}
+
+/*******************************************************************************
+ * Interface
+ ******************************************************************************/
+struct interf {
+ double hc;
+ double epsilon;
+ double specular_fraction;
+ double temperature;
+};
+static const struct interf INTERF_NULL = {0, 0, 0, UNKOWN_TEMPERATURE};
+
+static double
+interface_get_convection_coef
+ (const struct sdis_interface_fragment* frag, struct sdis_data* data)
+{
+ CHK(data != NULL && frag != NULL);
+ return ((const struct interf*)sdis_data_cget(data))->hc;
+}
+
+static double
+interface_get_emissivity
+ (const struct sdis_interface_fragment* frag, struct sdis_data* data)
+{
+ CHK(data != NULL && frag != NULL);
+ return ((const struct interf*)sdis_data_cget(data))->epsilon;
+}
+
+static double
+interface_get_specular_fraction
+ (const struct sdis_interface_fragment* frag, struct sdis_data* data)
+{
+ CHK(data != NULL && frag != NULL);
+ return ((const struct interf*)sdis_data_cget(data))->specular_fraction;
+}
+
+static double
+interface_get_temperature
+ (const struct sdis_interface_fragment* frag, struct sdis_data* data)
+{
+ CHK(data != NULL && frag != NULL);
+ return ((const struct interf*)sdis_data_cget(data))->temperature;
+}
+
+/*******************************************************************************
+ * Helper functions
+ ******************************************************************************/
+static void
+create_solid
+ (struct sdis_device* dev,
+ const struct solid* param,
+ struct sdis_medium** solid)
+{
+ struct sdis_data* data = NULL;
+ struct solid* solid_param = NULL;
+ struct sdis_solid_shader solid_shader = DUMMY_SOLID_SHADER;
+
+ CHK(param != NULL);
+ CHK(solid != NULL);
+
+ /* Copy the solid parameters into the Stardis memory space */
+ CHK(sdis_data_create
+ (dev, sizeof(struct solid), ALIGNOF(struct solid), NULL, &data) == RES_OK);
+ solid_param = sdis_data_get(data);
+ memcpy(solid_param, param, sizeof(struct solid));
+
+ /* Setup the solid shader */
+ solid_shader.calorific_capacity = solid_get_calorific_capacity;
+ solid_shader.calorific_capacity = solid_get_calorific_capacity;
+ solid_shader.thermal_conductivity = solid_get_thermal_conductivity;
+ solid_shader.volumic_mass = solid_get_volumic_mass;
+ solid_shader.delta_solid = solid_get_delta;
+ solid_shader.delta_boundary = solid_get_delta_boundary;
+ solid_shader.temperature = solid_get_temperature;
+
+ /* Create the solid medium */
+ CHK(sdis_solid_create(dev, &solid_shader, data, solid) == RES_OK);
+
+ /* Release the ownership onto the Stardis memory space storing the solid
+ * parameters */
+ CHK(sdis_data_ref_put(data) == RES_OK);
+}
+
+static void
+create_fluid
+ (struct sdis_device* dev,
+ const struct fluid* param,
+ struct sdis_medium** fluid)
+{
+ struct sdis_data* data = NULL;
+ struct fluid* fluid_param = NULL;
+ struct sdis_fluid_shader fluid_shader = DUMMY_FLUID_SHADER;
+
+ CHK(param != NULL);
+ CHK(fluid != NULL);
+
+ /* Copy the fluid parameters into the Stardis memory space */
+ CHK(sdis_data_create
+ (dev, sizeof(struct fluid), ALIGNOF(struct fluid), NULL, &data) == RES_OK);
+ fluid_param = sdis_data_get(data);
+ memcpy(fluid_param, param, sizeof(struct fluid));
+
+ /* Setup the fluid shader */
+ fluid_shader.calorific_capacity = fluid_get_calorific_capacity;
+ fluid_shader.volumic_mass = fluid_get_volumic_mass;
+ fluid_shader.temperature = fluid_get_temperature;
+
+ /* Create the fluid medium */
+ CHK(sdis_fluid_create(dev, &fluid_shader, data, fluid) == RES_OK);
+
+ /* Release the ownership onto the Stardis memory space storing the fluid
+ * parameters */
+ CHK(sdis_data_ref_put(data) == RES_OK);
+}
+
+static void
+create_interface
+ (struct sdis_device* dev,
+ struct sdis_medium* mdm_front,
+ struct sdis_medium* mdm_back,
+ const struct interf* param,
+ struct sdis_interface** interf)
+{
+ struct sdis_data* data = NULL;
+ struct interf* interface_param = NULL;
+ struct sdis_interface_shader interface_shader = DUMMY_INTERFACE_SHADER;
+
+ CHK(mdm_front != NULL);
+ CHK(mdm_back != NULL);
+ CHK(param != NULL);
+ CHK(interf != NULL);
+
+ /* Copy the interface parameters into the Stardis memory space */
+ CHK(sdis_data_create
+ (dev, sizeof(struct interf), ALIGNOF(struct interf), NULL, &data) == RES_OK);
+ interface_param = sdis_data_get(data);
+ memcpy(interface_param, param, sizeof(struct interf));
+
+ /* 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;
+
+ /* Create the interface */
+ CHK(sdis_interface_create
+ (dev, mdm_front, mdm_back, &interface_shader, data, interf) == RES_OK);
+
+ /* Release the ownership onto the Stardis memory space storing the interface
+ * parameters */
+ CHK(sdis_data_ref_put(data) == RES_OK);
+}
+
+static void
+geometry_add_shape
+ (struct geometry* geom,
+ const double* positions,
+ const size_t nverts,
+ const size_t* indices,
+ const size_t nprims,
+ const double transform[12], /* May be NULL <=> no transformation */
+ struct sdis_interface* interf)
+{
+ struct map_interf* geom_interf = NULL;
+ size_t nverts_prev = 0;
+ size_t i;
+
+ CHK(geom != NULL);
+ CHK(positions != NULL);
+ CHK(indices != NULL);
+ CHK(nverts != 0);
+ CHK(nprims != 0);
+ CHK(interf != NULL);
+
+ /* Save the previous number of vertices/primitives of the geometry */
+ nverts_prev = sa_size(geom->positions) / 3;
+
+ /* Add the vertices */
+ FOR_EACH(i, 0, nverts) {
+ double* pos = sa_add(geom->positions, 3);
+ d3_set(pos, positions + i*3);
+ if(transform) {
+ d33_muld3(pos, transform, pos);
+ d3_add(pos, transform+9, pos);
+ }
+ }
+
+ /* Add the indices */
+ FOR_EACH(i, 0, nprims) {
+ sa_push(geom->indices, indices[i*3+0] + nverts_prev);
+ sa_push(geom->indices, indices[i*3+1] + nverts_prev);
+ sa_push(geom->indices, indices[i*3+2] + nverts_prev);
+ }
+
+ geom_interf = sa_add(geom->interfaces, 1);
+ geom_interf->key = sa_size(geom->indices) / 3 - 1;
+ geom_interf->interf = interf;
+}
+
+static void
+dump_image(const struct sdis_accum_buffer* buf)
+{
+ struct sdis_accum_buffer_layout layout = SDIS_ACCUM_BUFFER_LAYOUT_NULL;
+ struct image img;
+ double* temps = NULL;
+ const struct sdis_accum* accums = NULL;
+ double Tmax = -DBL_MAX;
+ double Tmin = DBL_MAX;
+ double norm;
+ size_t ix, iy;
+
+ CHK(buf != NULL);
+ CHK(sdis_accum_buffer_get_layout(buf, &layout) == RES_OK);
+
+ 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);
+ FOR_EACH(iy, 0, layout.height) {
+ const struct sdis_accum* row_accums = accums + iy * layout.width;
+ double* row = temps + iy * layout.width;
+ FOR_EACH(ix, 0, layout.width) {
+ row[ix] = row_accums[ix].sum_weights / (double)row_accums[ix].nweights;
+ Tmax = MMAX(row[ix], Tmax);
+ Tmin = MMIN(row[ix], Tmin);
+ }
+ }
+ if(Tmax != Tmin) {
+ norm = Tmax - Tmin;
+ } else {
+ Tmin = 0;
+ norm = 1;
+ }
+
+ /* Allocate the image memory space */
+ CHK(image_init(NULL, &img) == RES_OK);
+ CHK(image_setup(&img, IMG_WIDTH, IMG_HEIGHT, IMG_WIDTH*3, IMAGE_RGB8, NULL)
+ == RES_OK);
+
+ FOR_EACH(iy, 0, layout.height) {
+ const double* src_row = temps + iy*layout.width;
+ char* dst_row = img.pixels + iy*img.pitch;
+ FOR_EACH(ix, 0, layout.width) {
+ unsigned char* pixels = (unsigned char*)
+ (dst_row + ix * sizeof_image_format(img.format));
+ const unsigned char T = (unsigned char)
+ ((src_row[ix] - Tmin)/ norm * 255.0);
+ pixels[0] = T;
+ pixels[1] = T;
+ pixels[2] = T;
+ }
+ }
+ CHK(image_write_ppm_stream(&img, 0/*binary?*/, stdout) == RES_OK);
+ image_release(&img);
+ mem_rm(temps);
+}
+
+/*******************************************************************************
+ * Test
+ ******************************************************************************/
+int
+main(int argc, char** argv)
+{
+ struct mem_allocator allocator;
+ struct geometry geom = GEOMETRY_NULL;
+ struct s3dut_mesh* msh = NULL;
+ struct s3dut_mesh_data msh_data;
+ struct sdis_accum_buffer* buf = NULL;
+ struct sdis_camera* cam = NULL;
+ struct sdis_device* dev = NULL;
+ struct sdis_medium* solid = NULL;
+ struct sdis_medium* fluid0 = NULL;
+ struct sdis_medium* fluid1 = NULL;
+ struct sdis_interface* interf0 = NULL;
+ struct sdis_interface* interf1 = NULL;
+ struct sdis_scene* scn = NULL;
+ struct fluid fluid_param = FLUID_NULL;
+ struct solid solid_param = SOLID_NULL;
+ struct interf interface_param = INTERF_NULL;
+ size_t ntris, npos;
+ double pos[3];
+ double tgt[3];
+ double up[3];
+ (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 fluid0 */
+ fluid_param.temperature = 350;
+ fluid_param.rho = 0;
+ fluid_param.cp = 0;
+ create_fluid(dev, &fluid_param, &fluid0);
+
+ /* Create the fluid1 */
+ fluid_param.temperature = UNKOWN_TEMPERATURE;
+ fluid_param.rho = 0;
+ fluid_param.cp = 0;
+ create_fluid(dev, &fluid_param, &fluid1);
+
+ /* Create the solid medium */
+ solid_param.cp = 1.0;
+ solid_param.lambda = 0.1;
+ solid_param.rho = 1.0;
+ solid_param.delta = 1.0/20.0;
+ solid_param.temperature = UNKOWN_TEMPERATURE;
+ create_solid(dev, &solid_param, &solid);
+
+ /* Create the fluid0/solid interface */
+ interface_param.hc = 1;
+ interface_param.epsilon = 0;
+ interface_param.specular_fraction = 0;
+ interface_param.temperature = UNKOWN_TEMPERATURE;
+ create_interface(dev, solid, fluid0, &interface_param, &interf0);
+
+ /* Create the fluid1/solid interface */
+ interface_param.hc = 0;
+ interface_param.epsilon = 1;
+ interface_param.specular_fraction = 1;
+ interface_param.temperature = UNKOWN_TEMPERATURE;
+ create_interface(dev, fluid1, solid, &interface_param, &interf1);
+
+ /* Release the ownership onto the media */
+ CHK(sdis_medium_ref_put(solid) == RES_OK);
+ CHK(sdis_medium_ref_put(fluid0) == RES_OK);
+ CHK(sdis_medium_ref_put(fluid1) == RES_OK);
+
+ /* Setup the cube geometry */
+ CHK(s3dut_create_cuboid(&allocator, 2, 2, 2, &msh) == RES_OK);
+ CHK(s3dut_mesh_get_data(msh, &msh_data) == RES_OK);
+ geometry_add_shape(&geom, msh_data.positions, msh_data.nvertices,
+ msh_data.indices, msh_data.nprimitives, NULL, interf1);
+ CHK(s3dut_mesh_ref_put(msh) == RES_OK);
+
+ /* Setup the sphere geometry */
+ CHK(s3dut_create_sphere(&allocator, 0.5, 32, 16, &msh) == RES_OK);
+ CHK(s3dut_mesh_get_data(msh, &msh_data) == RES_OK);
+ geometry_add_shape(&geom, msh_data.positions, msh_data.nvertices,
+ msh_data.indices, msh_data.nprimitives, NULL, interf0);
+ CHK(s3dut_mesh_ref_put(msh) == RES_OK);
+
+ /* Setup the scene */
+ ntris = sa_size(geom.indices) / 3; /* #primitives */
+ npos = sa_size(geom.positions) / 3; /* #positions */
+ CHK(sdis_scene_create(dev, ntris, geometry_get_indices,
+ geometry_get_interface, npos, geometry_get_position,
+ &geom, &scn) == RES_OK);
+
+ /* Setup the camera */
+ d3(pos, 3, 3, 3);
+ d3(tgt, 0, 0, 0);
+ d3(up, 0, 0, 1);
+ CHK(sdis_camera_create(dev, &cam) == RES_OK);
+ CHK(sdis_camera_set_proj_ratio
+ (cam, (double)IMG_WIDTH/(double)IMG_HEIGHT) == RES_OK);
+ CHK(sdis_camera_set_fov(cam, MDEG2RAD(70)) == RES_OK);
+ CHK(sdis_camera_look_at(cam, pos, tgt, up) == RES_OK);
+
+ /* Create the accum buffer */
+ CHK(sdis_accum_buffer_create(dev, IMG_WIDTH, IMG_HEIGHT, &buf) == RES_OK);
+
+ /* 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);
+
+ /* Write the image */
+ dump_image(buf);
+
+ /* Release memory */
+ CHK(sdis_scene_ref_put(scn) == RES_OK);
+ CHK(sdis_camera_ref_put(cam) == RES_OK);
+ CHK(sdis_interface_ref_put(interf0) == RES_OK);
+ CHK(sdis_interface_ref_put(interf1) == RES_OK);
+ CHK(sdis_device_ref_put(dev) == RES_OK);
+ CHK(sdis_accum_buffer_ref_put(buf) == RES_OK);
+ geometry_release(&geom);
+
+ check_memory_allocator(&allocator);
+ mem_shutdown_proxy_allocator(&allocator);
+ CHK(mem_allocated_size() == 0);
+ return 0;
+}
+
diff --git a/src/test_sdis_solve_probe.c b/src/test_sdis_solve_probe.c
@@ -146,6 +146,8 @@ solid_get_temperature
******************************************************************************/
struct interf {
double hc;
+ double epsilon;
+ double specular_fraction;
};
static double
@@ -156,6 +158,22 @@ interface_get_convection_coef
return ((const struct interf*)sdis_data_cget(data))->hc;
}
+static double
+interface_get_emissivity
+ (const struct sdis_interface_fragment* frag, struct sdis_data* data)
+{
+ CHK(data != NULL && frag != NULL);
+ return ((const struct interf*)sdis_data_cget(data))->epsilon;
+}
+
+static double
+interface_get_specular_fraction
+ (const struct sdis_interface_fragment* frag, struct sdis_data* data)
+{
+ CHK(data != NULL && frag != NULL);
+ return ((const struct interf*)sdis_data_cget(data))->specular_fraction;
+}
+
/*******************************************************************************
* Test
******************************************************************************/
@@ -222,8 +240,12 @@ main(int argc, char** argv)
ALIGNOF(struct interf), NULL, &data) == RES_OK);
interface_param = sdis_data_get(data);
interface_param->hc = 0.5;
+ 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;
CHK(sdis_interface_create
(dev, solid, fluid, &interface_shader, data, &interf) == RES_OK);
CHK(sdis_data_ref_put(data) == RES_OK);
@@ -246,12 +268,13 @@ main(int argc, char** argv)
pos[1] = 0.5;
pos[2] = 0.5;
time = INF;
- CHK(sdis_solve_probe(NULL, N, pos, time, 1.0, &estimator) == RES_BAD_ARG);
- CHK(sdis_solve_probe(scn, 0, pos, time, 1.0, &estimator) == RES_BAD_ARG);
- CHK(sdis_solve_probe(scn, N, NULL, time, 1.0, &estimator) == RES_BAD_ARG);
- CHK(sdis_solve_probe(scn, N, pos, time, 0, &estimator) == RES_BAD_ARG);
- CHK(sdis_solve_probe(scn, N, pos, time, 1.0, NULL) == RES_BAD_ARG);
- CHK(sdis_solve_probe(scn, N, pos, time, 1.0, &estimator) == RES_OK);
+ CHK(sdis_solve_probe(NULL, N, pos, time, 1.0, 0, 0, &estimator) == RES_BAD_ARG);
+ CHK(sdis_solve_probe(scn, 0, pos, time, 1.0, 0, 0, &estimator) == RES_BAD_ARG);
+ CHK(sdis_solve_probe(scn, N, NULL, time, 1.0, 0, 0, &estimator) == RES_BAD_ARG);
+ CHK(sdis_solve_probe(scn, N, pos, time, 0, 0, 0, &estimator) == RES_BAD_ARG);
+ CHK(sdis_solve_probe(scn, N, pos, time, 0, 0, -1, &estimator) == RES_BAD_ARG);
+ CHK(sdis_solve_probe(scn, N, pos, time, 1.0, 0, 0, NULL) == RES_BAD_ARG);
+ CHK(sdis_solve_probe(scn, N, pos, time, 1.0, 0, 0, &estimator) == RES_OK);
CHK(sdis_estimator_get_realisation_count(estimator, NULL) == RES_BAD_ARG);
CHK(sdis_estimator_get_realisation_count(NULL, &nreals) == RES_BAD_ARG);
diff --git a/src/test_sdis_solve_probe2.c b/src/test_sdis_solve_probe2.c
@@ -131,7 +131,7 @@ struct interf {
};
static double
-null_convection_coef
+null_interface_value
(const struct sdis_interface_fragment* frag, struct sdis_data* data)
{
CHK(frag != NULL);
@@ -196,8 +196,10 @@ main(int argc, char** argv)
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_convection_coef;
+ interface_shader.convection_coef = null_interface_value;
interface_shader.temperature = NULL;
+ interface_shader.emissivity = null_interface_value;
+ interface_shader.specular_fraction = null_interface_value;
CHK(sdis_interface_create
(dev, solid, fluid, &interface_shader, NULL, &Tnone) == RES_OK);
@@ -206,8 +208,10 @@ 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_convection_coef;
+ 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);
@@ -217,8 +221,10 @@ 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_convection_coef;
+ 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);
@@ -253,7 +259,7 @@ main(int argc, char** argv)
pos[1] = 0.5;
pos[2] = 0.5;
time = INF;
- CHK(sdis_solve_probe( scn, N, pos, time, 1.0, &estimator) == RES_OK);
+ CHK(sdis_solve_probe( scn, N, pos, time, 1.0, -1, 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);
diff --git a/src/test_sdis_solve_probe2_2d.c b/src/test_sdis_solve_probe2_2d.c
@@ -128,7 +128,7 @@ struct interf {
};
static double
-null_convection_coef
+null_interface_value
(const struct sdis_interface_fragment* frag, struct sdis_data* data)
{
CHK(frag != NULL);
@@ -193,8 +193,10 @@ main(int argc, char** argv)
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_convection_coef;
+ interface_shader.convection_coef = null_interface_value;
interface_shader.temperature = NULL;
+ interface_shader.emissivity = null_interface_value;
+ interface_shader.specular_fraction = null_interface_value;
CHK(sdis_interface_create
(dev, solid, fluid, &interface_shader, NULL, &Tnone) == RES_OK);
@@ -203,8 +205,10 @@ 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_convection_coef;
+ 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);
@@ -214,8 +218,10 @@ 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_convection_coef;
+ 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);
@@ -247,7 +253,7 @@ main(int argc, char** argv)
pos[0] = 0.5;
pos[1] = 0.5;
time = INF;
- CHK(sdis_solve_probe( scn, N, pos, time, 1.0, &estimator) == RES_OK);
+ CHK(sdis_solve_probe( scn, N, pos, time, 1.0, -1, 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);
diff --git a/src/test_sdis_solve_probe3.c b/src/test_sdis_solve_probe3.c
@@ -52,7 +52,7 @@ struct context {
struct sdis_interface* solid_fluid_T350;
struct sdis_interface* solid_solid;
};
-static const struct context CONTEXT_NULL = { NULL };
+static const struct context CONTEXT_NULL = {NULL, NULL, NULL, NULL, NULL, NULL};
static void
get_indices(const size_t itri, size_t ids[3], void* context)
@@ -153,7 +153,7 @@ struct interf {
};
static double
-null_convection_coef
+null_interface_value
(const struct sdis_interface_fragment* frag, struct sdis_data* data)
{
CHK(frag != NULL);
@@ -223,8 +223,10 @@ main(int argc, char** argv)
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_convection_coef;
+ interface_shader.convection_coef = null_interface_value;
interface_shader.temperature = NULL;
+ interface_shader.emissivity = null_interface_value;
+ interface_shader.specular_fraction = null_interface_value;
CHK(sdis_interface_create
(dev, solid, fluid, &interface_shader, NULL, &Tnone) == RES_OK);
@@ -233,8 +235,10 @@ 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_convection_coef;
+ 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);
@@ -244,8 +248,10 @@ 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_convection_coef;
+ 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);
@@ -253,6 +259,8 @@ main(int argc, char** argv)
/* Create the solid/solid interface */
interface_shader.convection_coef = NULL;
interface_shader.temperature = NULL;
+ interface_shader.specular_fraction = NULL;
+ interface_shader.emissivity = NULL;
CHK(sdis_interface_create
(dev, solid, solid, &interface_shader, NULL, &solid_solid) == RES_OK);
@@ -310,7 +318,7 @@ main(int argc, char** argv)
pos[1] = 0.5;
pos[2] = 0.5;
time = INF;
- CHK(sdis_solve_probe( scn, N, pos, time, 1.0, &estimator) == RES_OK);
+ CHK(sdis_solve_probe( scn, N, pos, time, 1.0, -1, 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);
diff --git a/src/test_sdis_solve_probe3_2d.c b/src/test_sdis_solve_probe3_2d.c
@@ -150,7 +150,7 @@ struct interf {
};
static double
-null_convection_coef
+null_interface_value
(const struct sdis_interface_fragment* frag, struct sdis_data* data)
{
CHK(frag != NULL);
@@ -218,8 +218,10 @@ main(int argc, char** argv)
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_convection_coef;
+ interface_shader.convection_coef = null_interface_value;
interface_shader.temperature = NULL;
+ interface_shader.emissivity = null_interface_value;
+ interface_shader.specular_fraction = null_interface_value;
CHK(sdis_interface_create
(dev, solid, fluid, &interface_shader, NULL, &Tnone) == RES_OK);
@@ -228,8 +230,10 @@ 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_convection_coef;
+ 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);
@@ -239,8 +243,10 @@ 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_convection_coef;
+ 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);
@@ -248,6 +254,8 @@ main(int argc, char** argv)
/* Create the solid/solid interface */
interface_shader.convection_coef = NULL;
interface_shader.temperature = NULL;
+ interface_shader.specular_fraction = NULL;
+ interface_shader.emissivity = NULL;
CHK(sdis_interface_create
(dev, solid, solid, &interface_shader, NULL, &solid_solid) == RES_OK);
@@ -300,7 +308,7 @@ main(int argc, char** argv)
pos[0] = 0.5;
pos[1] = 0.5;
time = INF;
- CHK(sdis_solve_probe( scn, N, pos, time, 1.0, &estimator) == RES_OK);
+ CHK(sdis_solve_probe( scn, N, pos, time, 1.0, -1, 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);
@@ -325,5 +333,4 @@ main(int argc, char** argv)
mem_shutdown_proxy_allocator(&allocator);
CHK(mem_allocated_size() == 0);
return 0;
-
}
diff --git a/src/test_sdis_solve_probe_2d.c b/src/test_sdis_solve_probe_2d.c
@@ -131,6 +131,14 @@ 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
******************************************************************************/
@@ -177,6 +185,8 @@ main(int argc, char** argv)
/* 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;
CHK(sdis_interface_create
(dev, solid, fluid, &interface_shader, NULL, &interf) == RES_OK);
@@ -197,7 +207,7 @@ main(int argc, char** argv)
pos[0] = 0.5;
pos[1] = 0.5;
time = INF;
- CHK(sdis_solve_probe(scn, N, pos, time, 1.0, &estimator) == RES_OK);
+ CHK(sdis_solve_probe(scn, N, pos, time, 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);
diff --git a/src/test_sdis_utils.h b/src/test_sdis_utils.h
@@ -16,9 +16,13 @@
#ifndef TEST_SDIS_UTILS_H
#define TEST_SDIS_UTILS_H
+#include "sdis.h"
+
#include <rsys/mem_allocator.h>
#include <stdio.h>
+#define BOLTZMANN_CONSTANT 5.6696e-8 /* W/m^2/K^4 */
+
/*******************************************************************************
* Box geometry
******************************************************************************/
@@ -82,7 +86,7 @@ dummy_medium_getter
{
(void)data;
CHK(vert != NULL);
- return 1;
+ return 0;
}
static INLINE double
@@ -91,7 +95,7 @@ dummy_interface_getter
{
(void)data;
CHK(frag != NULL);
- return 1;
+ return 0;
}
static const struct sdis_solid_shader DUMMY_SOLID_SHADER = {
@@ -111,6 +115,8 @@ static const struct sdis_fluid_shader DUMMY_FLUID_SHADER = {
static const struct sdis_interface_shader DUMMY_INTERFACE_SHADER = {
dummy_interface_getter,
+ dummy_interface_getter,
+ dummy_interface_getter,
dummy_interface_getter
};
