commit db52f0b369426b18d4335b8a5c80c34ab7fb0704
parent 03a18686244416c338e1eb2aae8babcc7f133a9c
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Tue, 13 Feb 2018 11:17:30 +0100
Merge branch 'release_0.0'
Diffstat:
20 files changed, 2082 insertions(+), 755 deletions(-)
diff --git a/README.md b/README.md
@@ -8,7 +8,8 @@ thermal problems.
Stardis relies on the [CMake](http://www.cmake.org) and the
[RCMake](https://gitlab.com/vaplv/rcmake/) package to build.
It also depends on the
-[RSys](https://gitlab.com/vaplv/rsys/) and
+[RSys](https://gitlab.com/vaplv/rsys/),
+[Star-2D](https://gitlab.com/meso-star/star-2d/),
[Star-3D](https://gitlab.com/meso-star/star-3d/) and
[Star-SP](https://gitlab.com/meso-star/star-sp/) libraries as well as on the
[OpenMP](http://www.openmp.org) 1.2 specification to parallelize its
@@ -20,6 +21,17 @@ well as all the aforementioned prerequisites. Finally generate the project from
the `cmake/CMakeLists.txt` file by appending to the `CMAKE_PREFIX_PATH`
variable the install directories of its dependencies.
+## Release notes
+
+### Version 0.0
+
+First version and implementation of the Stardis solver API.
+
+- Support fluid/solid and solid/solid interfaces.
+- Only conduction is currently fully supported: convection and radiative
+ temperature are not computed yet. Fluid media can be added to the system but
+ currently, Stardis assumes that their temperature are known.
+
## License
Stardis is Copyright (C) |Meso|Star> 2016-2018 (<contact@meso-star.com>). It is
diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt
@@ -24,8 +24,9 @@ option(NO_TEST "Do not build tests" OFF)
# Check dependencies
################################################################################
find_package(RCMake 0.3 REQUIRED)
+find_package(Star2D 0.1 REQUIRED)
find_package(Star3D 0.4 REQUIRED)
-find_package(StarSP 0.5 REQUIRED)
+find_package(StarSP 0.7 REQUIRED)
find_package(RSys 0.6 REQUIRED)
find_package(OpenMP 1.2 REQUIRED)
@@ -38,7 +39,7 @@ rcmake_append_runtime_dirs(_runtime_dirs RSys Star3D StarSP)
################################################################################
# Configure and define targets
################################################################################
-set(VERSION_MAJOR 1)
+set(VERSION_MAJOR 0)
set(VERSION_MINOR 0)
set(VERSION_PATCH 0)
set(VERSION ${VERSION_MAJOR}.${VERSION_MINOR}.${VERSION_PATCH})
@@ -74,7 +75,7 @@ add_library(sdis SHARED
${SDIS_FILES_SRC}
${SDIS_FILES_INC}
${SDIS_FILES_INC_API})
-target_link_libraries(sdis RSys Star3D StarSP)
+target_link_libraries(sdis RSys Star2D Star3D StarSP)
set_target_properties(sdis PROPERTIES
DEFINE_SYMBOL SDIS_SHARED_BUILD
@@ -118,8 +119,14 @@ if(NOT NO_TEST)
new_test(test_sdis_solve_probe)
new_test(test_sdis_solve_probe2)
new_test(test_sdis_solve_probe3)
+ new_test(test_sdis_solve_probe_2d)
+ new_test(test_sdis_solve_probe2_2d)
+ new_test(test_sdis_solve_probe3_2d)
target_link_libraries(test_sdis_solve_probe3 Star3DUT)
+ if(CMAKE_COMPILER_IS_GNUCC)
+ target_link_libraries(test_sdis_solve_probe3_2d m)
+ endif()
endif()
################################################################################
diff --git a/src/sdis.h b/src/sdis.h
@@ -70,20 +70,21 @@ enum sdis_medium_type {
SDIS_MEDIUM_TYPES_COUNT__
};
-/* Random walk vertex */
+/* Random walk vertex, i.e. a spatiotemporal position at a given step of the
+ * random walk. */
struct sdis_rwalk_vertex {
double P[3]; /* World space position */
- double time; /* Time of the vertex */
+ double time; /* "Time" of the vertex */
};
#define SDIS_RWALK_VERTEX_NULL__ {{0}, -1}
static const struct sdis_rwalk_vertex SDIS_RWALK_VERTEX_NULL =
SDIS_RWALK_VERTEX_NULL__;
-/* Point onto an interface */
+/* Spatiotemporal position onto an interface */
struct sdis_interface_fragment {
double P[3]; /* World space position */
- double Ng[3]; /* Normalized world space geometry normal */
- double uv[2]; /* Texture coordinates */
+ double Ng[3]; /* Normalized world space geometry normal at the interface */
+ double uv[2]; /* Parametric coordinates of the interface */
double time; /* Current time */
};
#define SDIS_INTERFACE_FRAGMENT_NULL__ {{0}, {0}, {0}, -1}
@@ -105,7 +106,6 @@ typedef double
(const struct sdis_rwalk_vertex* vert,
struct sdis_data* data);
-
/* Functor type to retrieve the interface properties. */
typedef double
(*sdis_interface_getter_T)
@@ -160,7 +160,7 @@ sdis_device_create
(struct logger* logger, /* May be NULL <=> use default logger */
struct mem_allocator* allocator, /* May be NULL <=> use default allocator */
const unsigned nthreads_hint, /* Hint on the number of threads to use */
- const int verbose,
+ const int verbose, /* Verbosity level */
struct sdis_device** dev);
SDIS_API res_T
@@ -172,12 +172,13 @@ sdis_device_ref_put
(struct sdis_device* dev);
/*******************************************************************************
- * A data stores in the Stardis memory space a set of user defined data.
+ * A data stores in the Stardis memory space a set of user defined data. It can
+ * be seen as a ref counted memory space allocated by Stardis.
******************************************************************************/
SDIS_API res_T
sdis_data_create
(struct sdis_device* dev,
- const size_t size, /* Size in bytes */
+ const size_t size, /* Size in bytes of user defined data */
const size_t align, /* Data alignment. Must be a power of 2 */
void (*release)(void*),/* Invoked priorly to data destruction. May be NULL */
struct sdis_data** data);
@@ -237,29 +238,46 @@ sdis_interface_create
struct sdis_medium* back,
const struct sdis_interface_shader* shader,
struct sdis_data* data, /* Data sent to the shader. May be NULL */
- struct sdis_interface** interface);
+ struct sdis_interface** interf);
SDIS_API res_T
sdis_interface_ref_get
- (struct sdis_interface* interface);
+ (struct sdis_interface* interf);
SDIS_API res_T
sdis_interface_ref_put
- (struct sdis_interface* interface);
+ (struct sdis_interface* interf);
/*******************************************************************************
- * A scene is a collection of triangles. Each triangle is the support of the
- * interface between 2 mediums.
+ * A scene is a collection of primitives. Each primitive is the geometric
+ * support of the interface between 2 mediums.
******************************************************************************/
SDIS_API res_T
sdis_scene_create
(struct sdis_device* dev,
const size_t ntris, /* #triangles */
- void (*indices)(const size_t itri, size_t ids[3], void*),
- void (*interface)(const size_t itri, struct sdis_interface** bound, void*),
+ void (*indices) /* Retrieve the indices toward the vertices of `itri' */
+ (const size_t itri, size_t ids[3], void*),
+ void (*interf) /* Get the interface of the triangle `itri' */
+ (const size_t itri, struct sdis_interface** bound, void*),
+ const size_t nverts, /* #vertices */
+ void (*position) /* Retrieve the position of the vertex `ivert' */
+ (const size_t ivert, double pos[3], void* ctx),
+ void* ctx, /* Client side data sent as input of the previous callbacks */
+ struct sdis_scene** scn);
+
+SDIS_API res_T
+sdis_scene_2d_create
+ (struct sdis_device* dev,
+ const size_t nsegs, /* #segments */
+ void (*indices) /* Retrieve the indices toward the vertices of `iseg' */
+ (const size_t iseg, size_t ids[2], void*),
+ void (*interf) /* Get the interface of the segment `iseg' */
+ (const size_t iseg, struct sdis_interface** bound, void*),
const size_t nverts, /* #vertices */
- void (*position)(const size_t ivert, double pos[3], void* ctx),
- void* ctx,
+ void (*position) /* Retrieve the position of the vertex `ivert' */
+ (const size_t ivert, double pos[2], void* ctx),
+ void* ctx, /* Client side data sent as input of the previous callbacks */
struct sdis_scene** scn);
SDIS_API res_T
@@ -270,6 +288,7 @@ SDIS_API res_T
sdis_scene_ref_put
(struct sdis_scene* scn);
+/* Retrieve the Axis Aligned Bounding Box of the scene */
SDIS_API res_T
sdis_scene_get_aabb
(const struct sdis_scene* scn,
diff --git a/src/sdis_device.c b/src/sdis_device.c
@@ -19,6 +19,7 @@
#include <rsys/logger.h>
#include <rsys/mem_allocator.h>
+#include <star/s2d.h>
#include <star/s3d.h>
#include <omp.h>
@@ -47,6 +48,7 @@ device_release(ref_T* ref)
struct sdis_device* dev;
ASSERT(ref);
dev = CONTAINER_OF(ref, struct sdis_device, ref);
+ if(dev->s2d) S2D(device_ref_put(dev->s2d));
if(dev->s3d) S3D(device_ref_put(dev->s3d));
ASSERT(flist_name_is_empty(&dev->names));
flist_name_release(&dev->names);
@@ -93,9 +95,15 @@ sdis_device_create
ref_init(&dev->ref);
flist_name_init(allocator, &dev->names);
+ 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);
+ }
+
res = s3d_device_create(log, allocator, 0, &dev->s3d);
if(res != RES_OK) {
- log_err(dev,
+ log_err(dev,
"%s: could not create the Star-3D device on Stardis.\n", FUNC_NAME);
goto error;
}
diff --git a/src/sdis_device_c.h b/src/sdis_device_c.h
@@ -31,6 +31,7 @@ struct sdis_device {
struct flist_name names;
+ struct s2d_device* s2d;
struct s3d_device* s3d;
ref_T ref;
diff --git a/src/sdis_interface.c b/src/sdis_interface.c
@@ -21,6 +21,7 @@
#include <rsys/double3.h>
#include <rsys/mem_allocator.h>
+#include <star/s2d.h>
#include <star/s3d.h>
/*******************************************************************************
@@ -57,16 +58,16 @@ check_interface_shader
static void
interface_release(ref_T* ref)
{
- struct sdis_interface* interface = NULL;
+ struct sdis_interface* interf = NULL;
struct sdis_device* dev = NULL;
ASSERT(ref);
- interface = CONTAINER_OF(ref, struct sdis_interface, ref);
- dev = interface->dev;
- if(interface->medium_front) SDIS(medium_ref_put(interface->medium_front));
- if(interface->medium_back) SDIS(medium_ref_put(interface->medium_back));
- if(interface->data) SDIS(data_ref_put(interface->data));
- flist_name_del(&dev->names, interface->id);
- MEM_RM(dev->allocator, interface);
+ interf = CONTAINER_OF(ref, struct sdis_interface, ref);
+ dev = interf->dev;
+ if(interf->medium_front) SDIS(medium_ref_put(interf->medium_front));
+ if(interf->medium_back) SDIS(medium_ref_put(interf->medium_back));
+ if(interf->data) SDIS(data_ref_put(interf->data));
+ flist_name_del(&dev->names, interf->id);
+ MEM_RM(dev->allocator, interf);
SDIS(device_ref_put(dev));
}
@@ -82,7 +83,7 @@ sdis_interface_create
struct sdis_data* data,
struct sdis_interface** out_interface)
{
- struct sdis_interface* interface = NULL;
+ struct sdis_interface* interf = NULL;
res_T res = RES_OK;
if(!dev || !front || !back || !shader || !out_interface) {
@@ -103,51 +104,51 @@ sdis_interface_create
goto error;
}
- interface = MEM_CALLOC(dev->allocator, 1, sizeof(struct sdis_interface));
- if(!interface) {
+ interf = MEM_CALLOC(dev->allocator, 1, sizeof(struct sdis_interface));
+ if(!interf) {
log_err(dev, "%s: could not create the interface.\n", FUNC_NAME);
res = RES_MEM_ERR;
goto error;
}
- ref_init(&interface->ref);
+ ref_init(&interf->ref);
SDIS(medium_ref_get(front));
SDIS(medium_ref_get(back));
SDIS(device_ref_get(dev));
- interface->medium_front = front;
- interface->medium_back = back;
- interface->dev = dev;
- interface->shader = *shader;
- interface->id = flist_name_add(&dev->names);
+ interf->medium_front = front;
+ interf->medium_back = back;
+ interf->dev = dev;
+ interf->shader = *shader;
+ interf->id = flist_name_add(&dev->names);
if(data) {
SDIS(data_ref_get(data));
- interface->data = data;
+ interf->data = data;
}
exit:
- if(out_interface) *out_interface = interface;
+ if(out_interface) *out_interface = interf;
return res;
error:
- if(interface) {
- SDIS(interface_ref_put(interface));
- interface = NULL;
+ if(interf) {
+ SDIS(interface_ref_put(interf));
+ interf = NULL;
}
goto exit;
}
res_T
-sdis_interface_ref_get(struct sdis_interface* interface)
+sdis_interface_ref_get(struct sdis_interface* interf)
{
- if(!interface) return RES_BAD_ARG;
- ref_get(&interface->ref);
+ if(!interf) return RES_BAD_ARG;
+ ref_get(&interf->ref);
return RES_OK;
}
res_T
-sdis_interface_ref_put(struct sdis_interface* interface)
+sdis_interface_ref_put(struct sdis_interface* interf)
{
- if(!interface) return RES_BAD_ARG;
- ref_put(&interface->ref, interface_release);
+ if(!interf) return RES_BAD_ARG;
+ ref_put(&interf->ref, interface_release);
return RES_OK;
}
@@ -156,27 +157,42 @@ sdis_interface_ref_put(struct sdis_interface* interface)
******************************************************************************/
const struct sdis_medium*
interface_get_medium
- (const struct sdis_interface* interface, const enum sdis_side_flag side)
+ (const struct sdis_interface* interf, const enum sdis_side_flag side)
{
struct sdis_medium* mdm = NULL;
- ASSERT(interface);
+ ASSERT(interf);
switch(side) {
- case SDIS_FRONT: mdm = interface->medium_front; break;
- case SDIS_BACK: mdm = interface->medium_back; break;
+ case SDIS_FRONT: mdm = interf->medium_front; break;
+ case SDIS_BACK: mdm = interf->medium_back; break;
default: FATAL("Unreachable code.\n"); break;
}
return mdm;
}
unsigned
-interface_get_id(const struct sdis_interface* interface)
+interface_get_id(const struct sdis_interface* interf)
{
- ASSERT(interface);
- return interface->id.index;
+ ASSERT(interf);
+ return interf->id.index;
}
void
-setup_interface_fragment
+setup_interface_fragment_2d
+ (struct sdis_interface_fragment* frag,
+ const struct sdis_rwalk_vertex* vertex,
+ const struct s2d_hit* hit)
+{
+ ASSERT(frag && vertex && hit && !S2D_HIT_NONE(hit));
+ d2_set(frag->P, vertex->P);
+ frag->P[2] = 0;
+ d2_normalize(frag->Ng, d2_set_f2(frag->Ng, hit->normal));
+ frag->Ng[2] = 0;
+ frag->uv[0] = hit->u;
+ frag->time = vertex->time;
+}
+
+void
+setup_interface_fragment_3d
(struct sdis_interface_fragment* frag,
const struct sdis_rwalk_vertex* vertex,
const struct s3d_hit* hit)
diff --git a/src/sdis_interface_c.h b/src/sdis_interface_c.h
@@ -22,6 +22,7 @@
#include <float.h>
/* Forward declaration of external type */
+struct s2d_hit;
struct s3d_hit;
struct sdis_interface {
@@ -37,36 +38,42 @@ struct sdis_interface {
extern LOCAL_SYM const struct sdis_medium*
interface_get_medium
- (const struct sdis_interface* interface,
+ (const struct sdis_interface* interf,
const enum sdis_side_flag side);
extern LOCAL_SYM unsigned
interface_get_id
- (const struct sdis_interface* interface);
+ (const struct sdis_interface* interf);
extern LOCAL_SYM void
-setup_interface_fragment
+setup_interface_fragment_2d
+ (struct sdis_interface_fragment* frag,
+ const struct sdis_rwalk_vertex* vertex,
+ const struct s2d_hit* hit);
+
+extern LOCAL_SYM void
+setup_interface_fragment_3d
(struct sdis_interface_fragment* frag,
const struct sdis_rwalk_vertex* vertex,
const struct s3d_hit* hit);
static INLINE double
interface_get_temperature
- (const struct sdis_interface* interface,
+ (const struct sdis_interface* interf,
const struct sdis_interface_fragment* frag)
{
- ASSERT(interface && frag);
- if(!interface->shader.temperature) return -DBL_MAX;
- return interface->shader.temperature(frag, interface->data);
+ ASSERT(interf && frag);
+ if(!interf->shader.temperature) return -DBL_MAX;
+ return interf->shader.temperature(frag, interf->data);
}
static INLINE double
interface_get_convection_coef
- (const struct sdis_interface* interface,
+ (const struct sdis_interface* interf,
const struct sdis_interface_fragment* frag)
{
- ASSERT(interface && frag);
- return interface->shader.convection_coef(frag, interface->data);
+ ASSERT(interf && frag);
+ return interf->shader.convection_coef(frag, interf->data);
}
#endif /* SDIS_INTERFACE_C_H */
diff --git a/src/sdis_scene.c b/src/sdis_scene.c
@@ -18,19 +18,122 @@
#include "sdis_interface_c.h"
#include "sdis_scene_c.h"
+#include <rsys/float2.h>
#include <rsys/float3.h>
+#include <rsys/double2.h>
#include <rsys/double3.h>
#include <rsys/mem_allocator.h>
+#include <star/s2d.h>
#include <star/s3d.h>
#include <limits.h>
+/* Context used to wrap the user geometry to Star-XD. */
+struct geometry_context {
+ void (*indices)(const size_t iprim, size_t ids[], void*);
+ void (*position)(const size_t ivert, double pos[], void*);
+ void* data;
+};
+
/*******************************************************************************
* Helper function
******************************************************************************/
+/* Check that `hit' roughly lies on a vertex. For segments, a simple but
+ * approximative way is to test that its position have at least one barycentric
+ * coordinate roughly equal to 0 or 1. */
+static FINLINE int
+hit_on_vertex(const struct s2d_hit* hit)
+{
+ const float on_vertex_eps = 1.e-4f;
+ float v;
+ ASSERT(hit && !S2D_HIT_NONE(hit));
+ v = 1.f - hit->u;
+ return eq_epsf(hit->u, 0.f, on_vertex_eps)
+ || eq_epsf(hit->u, 1.f, on_vertex_eps)
+ || eq_epsf(v, 0.f, on_vertex_eps)
+ || eq_epsf(v, 1.f, on_vertex_eps);
+}
+
+/* Check that `hit' roughly lies on an edge. For triangular primitives, a
+ * simple but approximative way is to test that its position have at least one
+ * barycentric coordinate roughly equal to 0 or 1. */
+static FINLINE int
+hit_on_edge(const struct s3d_hit* hit)
+{
+ const float on_edge_eps = 1.e-4f;
+ float w;
+ ASSERT(hit && !S3D_HIT_NONE(hit));
+ w = 1.f - hit->uv[0] - hit->uv[1];
+ return eq_epsf(hit->uv[0], 0.f, on_edge_eps)
+ || eq_epsf(hit->uv[0], 1.f, on_edge_eps)
+ || eq_epsf(hit->uv[1], 0.f, on_edge_eps)
+ || eq_epsf(hit->uv[1], 1.f, on_edge_eps)
+ || eq_epsf(w, 0.f, on_edge_eps)
+ || eq_epsf(w, 1.f, on_edge_eps);
+}
+
+static int
+hit_filter_function_3d
+ (const struct s3d_hit* hit,
+ const float org[3],
+ const float dir[3],
+ void* ray_data,
+ void* filter_data)
+{
+ const struct s3d_hit* hit_from = ray_data;
+ (void)org, (void)dir, (void)filter_data;
+
+ if(!hit_from || S3D_HIT_NONE(hit_from)) return 0; /* No filtering */
+
+ if(S3D_PRIMITIVE_EQ(&hit_from->prim, &hit->prim)) return 1;
+
+ if(eq_epsf(hit->distance, 0, 1.e-6f)) {
+ /* If the targeted point is near of the origin, check that it lies on an
+ * edge shared by the 2 primitives. */
+ return hit_on_edge(hit_from) && hit_on_edge(hit);
+ }
+
+ return 0;
+}
+
+static int
+hit_filter_function_2d
+ (const struct s2d_hit* hit,
+ const float org[2],
+ const float dir[2],
+ void* ray_data,
+ void* filter_data)
+{
+ const struct s2d_hit* hit_from = ray_data;
+ (void)org, (void)dir, (void)filter_data;
+
+ if(!hit_from || S2D_HIT_NONE(hit_from)) return 0; /* No filtering */
+
+ if(S2D_PRIMITIVE_EQ(&hit_from->prim, &hit->prim)) return 1;
+
+ if(eq_epsf(hit->distance, 0, 1.e-6f)) {
+ /* If the targeted point is near of the origin, check that it lies on a
+ * vertex shared by the 2 segments. */
+ return hit_on_vertex(hit_from) && hit_on_vertex(hit);
+ }
+
+ return 0;
+}
+
+static void
+get_indices_2d(const unsigned iseg, unsigned out_ids[2], void* data)
+{
+ struct geometry_context* ctx = data;
+ size_t ids[2];
+ ASSERT(ctx);
+ ctx->indices(iseg, ids, ctx->data);
+ out_ids[0] = (unsigned)ids[0];
+ out_ids[1] = (unsigned)ids[1];
+}
+
static void
-get_indices(const unsigned itri, unsigned out_ids[3], void* data)
+get_indices_3d(const unsigned itri, unsigned out_ids[3], void* data)
{
struct geometry_context* ctx = data;
size_t ids[3];
@@ -42,7 +145,18 @@ get_indices(const unsigned itri, unsigned out_ids[3], void* data)
}
static void
-get_position(const unsigned ivert, float out_pos[3], void* data)
+get_position_2d(const unsigned ivert, float out_pos[2], void* data)
+{
+ struct geometry_context* ctx = data;
+ double pos[2];
+ ASSERT(ctx);
+ ctx->position(ivert, pos, ctx->data);
+ out_pos[0] = (float)pos[0];
+ out_pos[1] = (float)pos[1];
+}
+
+static void
+get_position_3d(const unsigned ivert, float out_pos[3], void* data)
{
struct geometry_context* ctx = data;
double pos[3];
@@ -58,25 +172,25 @@ clear_interfaces(struct sdis_scene* scn)
{
size_t i;
ASSERT(scn);
- FOR_EACH(i, 0, darray_interface_size_get(&scn->interfaces)) {
- if(darray_interface_cdata_get(&scn->interfaces)[i]) {
- SDIS(interface_ref_put(darray_interface_data_get(&scn->interfaces)[i]));
+ FOR_EACH(i, 0, darray_interf_size_get(&scn->interfaces)) {
+ if(darray_interf_cdata_get(&scn->interfaces)[i]) {
+ SDIS(interface_ref_put(darray_interf_data_get(&scn->interfaces)[i]));
}
}
- darray_interface_clear(&scn->interfaces);
- darray_interface_clear(&scn->prim_interfaces);
+ darray_interf_clear(&scn->interfaces);
+ darray_interf_clear(&scn->prim_interfaces);
}
static res_T
setup_interfaces
(struct sdis_scene* scn,
const size_t ntris, /* #triangles */
- void (*interface)(const size_t itri, struct sdis_interface**, void*),
+ void (*interf)(const size_t itri, struct sdis_interface**, void*),
void* ctx)
{
size_t itri;
res_T res = RES_OK;
- ASSERT(ntris && interface);
+ ASSERT(ntris && interf);
clear_interfaces(scn);
@@ -86,24 +200,24 @@ setup_interfaces
unsigned id;
/* Retrieve the interface of the primitive */
- interface(itri, &itface, ctx);
+ interf(itri, &itface, ctx);
id = interface_get_id(itface);
/* Check that the interface is already registered against the scene */
- ninterfaces = darray_interface_size_get(&scn->interfaces);
+ ninterfaces = darray_interf_size_get(&scn->interfaces);
if(id >= ninterfaces) {
- res = darray_interface_resize(&scn->interfaces, id + 1);
+ res = darray_interf_resize(&scn->interfaces, id + 1);
if(res != RES_OK) goto error;
}
- if(darray_interface_cdata_get(&scn->interfaces)[id]) {
- ASSERT(darray_interface_cdata_get(&scn->interfaces)[id] == itface);
+ if(darray_interf_cdata_get(&scn->interfaces)[id]) {
+ ASSERT(darray_interf_cdata_get(&scn->interfaces)[id] == itface);
} else {
SDIS(interface_ref_get(itface));
- darray_interface_data_get(&scn->interfaces)[id] = itface;
+ darray_interf_data_get(&scn->interfaces)[id] = itface;
}
/* Register the primitive interface */
- res = darray_interface_push_back(&scn->prim_interfaces, &itface);
+ res = darray_interf_push_back(&scn->prim_interfaces, &itface);
if(res != RES_OK) goto error;
}
@@ -115,7 +229,59 @@ error:
}
static res_T
-setup_geometry
+setup_geometry_2d
+ (struct sdis_scene* scn,
+ const size_t nsegs, /* #segments */
+ void (*indices)(const size_t itri, size_t ids[2], void*),
+ const size_t nverts, /* #vertices */
+ void (*position)(const size_t ivert, double pos[2], void* ctx),
+ void* ctx)
+{
+ struct geometry_context context;
+ struct s2d_shape* s2d_msh = NULL;
+ struct s2d_scene* s2d_scn = NULL;
+ struct s2d_vertex_data vdata = S2D_VERTEX_DATA_NULL;
+ res_T res = RES_OK;
+ ASSERT(scn && nsegs && indices && nverts && position);
+
+ /* Setup the intermediary geometry context */
+ context.indices = indices;
+ context.position = position;
+ context.data = ctx;
+
+ /* Setup the vertex data */
+ vdata.usage = S2D_POSITION;
+ vdata.type = S2D_FLOAT2;
+ vdata.get = get_position_2d;
+
+ /* Create the Star-2D geometry */
+ res = s2d_scene_create(scn->dev->s2d, &s2d_scn);
+ if(res != RES_OK) goto error;
+ res = s2d_shape_create_line_segments(scn->dev->s2d, &s2d_msh);
+ if(res != RES_OK) goto error;
+ res = s2d_line_segments_set_hit_filter_function
+ (s2d_msh, hit_filter_function_2d, NULL);
+ if(res != RES_OK) goto error;
+ res = s2d_scene_attach_shape(s2d_scn, s2d_msh);
+ if(res != RES_OK) goto error;
+ res = s2d_line_segments_setup_indexed_vertices(s2d_msh, (unsigned)nsegs,
+ get_indices_2d, (unsigned)nverts, &vdata, 1, &context);
+ if(res != RES_OK) goto error;
+ res = s2d_scene_view_create(s2d_scn, S2D_SAMPLE|S2D_TRACE|S2D_GET_PRIMITIVE,
+ &scn->s2d_view);
+ if(res != RES_OK) goto error;
+
+exit:
+ if(s2d_msh) S2D(shape_ref_put(s2d_msh));
+ if(s2d_scn) S2D(scene_ref_put(s2d_scn));
+ return res;
+error:
+ if(scn->s2d_view) S2D(scene_view_ref_put(scn->s2d_view));
+ goto exit;
+}
+
+static res_T
+setup_geometry_3d
(struct sdis_scene* scn,
const size_t ntris, /* #triangles */
void (*indices)(const size_t itri, size_t ids[3], void*),
@@ -138,19 +304,19 @@ setup_geometry
/* Setup the vertex data */
vdata.usage = S3D_POSITION;
vdata.type = S3D_FLOAT3;
- vdata.get = get_position;
+ vdata.get = get_position_3d;
/* Create the Star-3D geometry */
res = s3d_scene_create(scn->dev->s3d, &s3d_scn);
if(res != RES_OK) goto error;
res = s3d_shape_create_mesh(scn->dev->s3d, &s3d_msh);
if(res != RES_OK) goto error;
- res = s3d_mesh_set_hit_filter_function(s3d_msh, hit_filter_function, NULL);
+ res = s3d_mesh_set_hit_filter_function(s3d_msh, hit_filter_function_3d, NULL);
if(res != RES_OK) goto error;
res = s3d_scene_attach_shape(s3d_scn, s3d_msh);
if(res != RES_OK) goto error;
- res = s3d_mesh_setup_indexed_vertices(s3d_msh, (unsigned)ntris, get_indices,
- (unsigned)nverts, &vdata, 1, &context);
+ res = s3d_mesh_setup_indexed_vertices(s3d_msh, (unsigned)ntris,
+ get_indices_3d, (unsigned)nverts, &vdata, 1, &context);
if(res != RES_OK) goto error;
res = s3d_scene_view_create(s3d_scn, S3D_SAMPLE|S3D_TRACE|S3D_GET_PRIMITIVE,
&scn->s3d_view);
@@ -165,59 +331,23 @@ error:
goto exit;
}
-/* Check that `hit' roughly lies on an edge. For triangular primitives, a
- * simple but approximative way is to test that its position have at least one
- * barycentric coordinate roughly equal to 0 or 1. */
-static FINLINE int
-hit_on_edge(const struct s3d_hit* hit)
-{
- const float on_edge_eps = 1.e-4f;
- float w;
- ASSERT(hit && !S3D_HIT_NONE(hit));
- w = 1.f - hit->uv[0] - hit->uv[1];
- return eq_epsf(hit->uv[0], 0.f, on_edge_eps)
- || eq_epsf(hit->uv[0], 1.f, on_edge_eps)
- || eq_epsf(hit->uv[1], 0.f, on_edge_eps)
- || eq_epsf(hit->uv[1], 1.f, on_edge_eps)
- || eq_epsf(w, 0.f, on_edge_eps)
- || eq_epsf(w, 1.f, on_edge_eps);
-}
-
-static void
-scene_release(ref_T * ref)
-{
- struct sdis_device* dev = NULL;
- struct sdis_scene* scn = NULL;
- ASSERT(ref);
- scn = CONTAINER_OF(ref, struct sdis_scene, ref);
- dev = scn->dev;
- clear_interfaces(scn);
- darray_interface_release(&scn->interfaces);
- darray_interface_release(&scn->prim_interfaces);
- if(scn->s3d_view) S3D(scene_view_ref_put(scn->s3d_view));
- MEM_RM(dev->allocator, scn);
- SDIS(device_ref_put(dev));
-}
-
-/*******************************************************************************
- * Exported functions
- ******************************************************************************/
-res_T
-sdis_scene_create
+static res_T
+scene_create
(struct sdis_device* dev,
- const size_t ntris, /* #triangles */
- void (*indices)(const size_t itri, size_t ids[3], void*),
- void (*interface)(const size_t itri, struct sdis_interface** bound, void*),
+ const int is_2d,
+ const size_t nprims, /* #primitives */
+ void (*indices)(const size_t iprim, size_t ids[], void*),
+ void (*interf)(const size_t iprim, struct sdis_interface** bound, void*),
const size_t nverts, /* #vertices */
- void (*position)(const size_t ivert, double pos[3], void* ctx),
+ void (*position)(const size_t ivert, double pos[], void* ctx),
void* ctx,
struct sdis_scene** out_scn)
{
struct sdis_scene* scn = NULL;
res_T res = RES_OK;
- if(!dev || !out_scn || !ntris || !indices || !interface || !nverts
- || !position || ntris > UINT_MAX || nverts > UINT_MAX) {
+ if(!dev || !out_scn || !nprims || !indices || !interf || !nverts
+ || !position || nprims > UINT_MAX || nverts > UINT_MAX) {
res = RES_BAD_ARG;
goto error;
}
@@ -231,16 +361,20 @@ sdis_scene_create
ref_init(&scn->ref);
SDIS(device_ref_get(dev));
scn->dev = dev;
- darray_interface_init(dev->allocator, &scn->interfaces);
- darray_interface_init(dev->allocator, &scn->prim_interfaces);
+ darray_interf_init(dev->allocator, &scn->interfaces);
+ darray_interf_init(dev->allocator, &scn->prim_interfaces);
- res = setup_interfaces(scn, ntris, interface, ctx);
+ res = setup_interfaces(scn, nprims, interf, ctx);
if(res != RES_OK) {
log_err(dev, "%s: could not setup the scene interfaces.\n", FUNC_NAME);
goto error;
}
- res = setup_geometry(scn, ntris, indices, nverts, position, ctx);
+ if(is_2d) {
+ res = setup_geometry_2d(scn, nprims, indices, nverts, position, ctx);
+ } else {
+ res = setup_geometry_3d(scn, nprims, indices, nverts, position, ctx);
+ }
if(res != RES_OK) {
log_err(dev, "%s: could not setup the scene geometry.\n", FUNC_NAME);
goto error;
@@ -257,48 +391,72 @@ error:
goto exit;
}
-res_T
-sdis_scene_ref_get(struct sdis_scene* scn)
+static INLINE res_T
+scene_get_medium_2d
+ (const struct sdis_scene* scn,
+ const double pos[2],
+ const struct sdis_medium** out_medium)
{
- if(!scn) return RES_BAD_ARG;
- ref_get(&scn->ref);
- return RES_OK;
-}
+ const struct sdis_medium* medium = NULL;
+ size_t iprim, nprims;
+ size_t nfailures = 0;
+ const size_t max_failures = 10;
+ res_T res = RES_OK;
+ ASSERT(scn && pos);
-res_T
-sdis_scene_ref_put(struct sdis_scene* scn)
-{
- if(!scn) return RES_BAD_ARG;
- ref_put(&scn->ref, scene_release);
- return RES_OK;
-}
+ S2D(scene_view_primitives_count(scn->s2d_view, &nprims));
+ FOR_EACH(iprim, 0, nprims) {
+ struct s2d_hit hit;
+ struct s2d_attrib attr;
+ struct s2d_primitive prim;
+ float s;
+ const float range[2] = {0.f, FLT_MAX};
+ float N[2], P[2], dir[2], cos_N_dir;
+ s = 1.f / 3.f;
-res_T
-sdis_scene_get_aabb
- (const struct sdis_scene* scn, double lower[3], double upper[3])
-{
- float low[3], upp[3];
- res_T res = RES_OK;
- if(!scn || !lower || !upper) return RES_BAD_ARG;
- res = s3d_scene_view_get_aabb(scn->s3d_view, low, upp);
- if(res != RES_OK) return res;
- d3_set_f3(lower, low);
- d3_set_f3(upper, upp);
- return RES_OK;
-}
+ /* Retrieve a position onto the primitive */
+ S2D(scene_view_get_primitive(scn->s2d_view, (unsigned)iprim, &prim));
+ S2D(primitive_get_attrib(&prim, S2D_POSITION, s, &attr));
-/*******************************************************************************
- * Local miscellaneous function
- ******************************************************************************/
-const struct sdis_interface*
-scene_get_interface(const struct sdis_scene* scn, const unsigned iprim)
-{
- ASSERT(scn && iprim < darray_interface_size_get(&scn->prim_interfaces));
- return darray_interface_cdata_get(&scn->prim_interfaces)[iprim];
+ /* Trace a ray from the randomw walk vertex toward the retrieved primitive
+ * position */
+ f2_normalize(dir, f2_sub(dir, attr.value, f2_set_d2(P, pos)));
+ S2D(scene_view_trace_ray(scn->s2d_view, P, dir, range, NULL, &hit));
+
+ f2_normalize(N, hit.normal);
+ cos_N_dir = f2_dot(N, dir);
+
+ /* Unforeseen error. One has to intersect a primitive ! */
+ if(S2D_HIT_NONE(&hit)) {
+ ++nfailures;
+ if(nfailures < max_failures) {
+ continue;
+ } else {
+ res = RES_BAD_ARG;
+ goto error;
+ }
+ }
+
+ if(absf(cos_N_dir) > 1.e-1f) { /* Not roughly orthognonal */
+ const struct sdis_interface* interf;
+ interf = scene_get_interface(scn, hit.prim.prim_id);
+ medium = interface_get_medium
+ (interf, cos_N_dir < 0 ? SDIS_FRONT : SDIS_BACK);
+ break;
+ }
+ }
+
+exit:
+ *out_medium = medium;
+ return res;
+error:
+ log_err(scn->dev, "%s: could not retrieve the medium at {%g, %g}.\n",
+ FUNC_NAME, SPLIT2(pos));
+ goto exit;
}
-res_T
-scene_get_medium
+static INLINE res_T
+scene_get_medium_3d
(const struct sdis_scene* scn,
const double pos[3],
const struct sdis_medium** out_medium)
@@ -315,9 +473,10 @@ scene_get_medium
struct s3d_hit hit;
struct s3d_attrib attr;
struct s3d_primitive prim;
- const float st[2] = { 1.f/3.f, 1.f/3.f };
+ float st[2];
const float range[2] = {0.f, FLT_MAX};
float N[3], P[3], dir[3], cos_N_dir;
+ st[0] = st[1] = 1.f / 3.f; /* Or MSVC will issue a warning */
/* Retrieve a position onto the primitive */
S3D(scene_view_get_primitive(scn->s3d_view, (unsigned)iprim, &prim));
@@ -343,10 +502,10 @@ scene_get_medium
}
if(absf(cos_N_dir) > 1.e-1f) { /* Not roughly orthognonal */
- const struct sdis_interface* interface;
- interface = scene_get_interface(scn, hit.prim.prim_id);
+ const struct sdis_interface* interf;
+ interf = scene_get_interface(scn, hit.prim.prim_id);
medium = interface_get_medium
- (interface, cos_N_dir < 0 ? SDIS_FRONT : SDIS_BACK);
+ (interf, cos_N_dir < 0 ? SDIS_FRONT : SDIS_BACK);
break;
}
}
@@ -360,27 +519,112 @@ error:
goto exit;
}
-int
-hit_filter_function
- (const struct s3d_hit* hit,
- const float org[3],
- const float dir[3],
- void* ray_data,
- void* filter_data)
+static void
+scene_release(ref_T * ref)
{
- const struct s3d_hit* hit_from = ray_data;
- (void)org, (void)dir, (void)filter_data;
+ struct sdis_device* dev = NULL;
+ struct sdis_scene* scn = NULL;
+ ASSERT(ref);
+ scn = CONTAINER_OF(ref, struct sdis_scene, ref);
+ dev = scn->dev;
+ clear_interfaces(scn);
+ darray_interf_release(&scn->interfaces);
+ darray_interf_release(&scn->prim_interfaces);
+ if(scn->s2d_view) S2D(scene_view_ref_put(scn->s2d_view));
+ if(scn->s3d_view) S3D(scene_view_ref_put(scn->s3d_view));
+ MEM_RM(dev->allocator, scn);
+ SDIS(device_ref_put(dev));
+}
- if(!hit_from || S3D_HIT_NONE(hit_from)) return 0; /* No filtering */
+/*******************************************************************************
+ * Exported functions
+ ******************************************************************************/
+res_T
+sdis_scene_create
+ (struct sdis_device* dev,
+ const size_t ntris, /* #triangles */
+ void (*indices)(const size_t itri, size_t ids[3], void*),
+ void (*interf)(const size_t itri, struct sdis_interface** bound, void*),
+ const size_t nverts, /* #vertices */
+ void (*position)(const size_t ivert, double pos[3], void* ctx),
+ void* ctx,
+ struct sdis_scene** out_scn)
+{
+ return scene_create
+ (dev, 0/*is_2D*/, ntris, indices, interf, nverts, position, ctx, out_scn);
+}
- if(S3D_PRIMITIVE_EQ(&hit_from->prim, &hit->prim)) return 1;
+res_T
+sdis_scene_2d_create
+ (struct sdis_device* dev,
+ const size_t nsegs, /* #segments */
+ void (*indices)(const size_t itri, size_t ids[2], void*),
+ void (*interf)(const size_t itri, struct sdis_interface** bound, void*),
+ const size_t nverts, /* #vertices */
+ void (*position)(const size_t ivert, double pos[2], void* ctx),
+ void* ctx,
+ struct sdis_scene** out_scn)
+{
+ return scene_create
+ (dev, 1/*is_2D*/, nsegs, indices, interf, nverts, position, ctx, out_scn);
+}
- if(eq_epsf(hit->distance, 0, 1.e-6f)) {
- /* If the targeted point is near of the origin, check that it lies on an
- * edge shared by the 2 primitives. */
- return hit_on_edge(hit_from) && hit_on_edge(hit);
+res_T
+sdis_scene_ref_get(struct sdis_scene* scn)
+{
+ if(!scn) return RES_BAD_ARG;
+ ref_get(&scn->ref);
+ return RES_OK;
+}
+
+res_T
+sdis_scene_ref_put(struct sdis_scene* scn)
+{
+ if(!scn) return RES_BAD_ARG;
+ ref_put(&scn->ref, scene_release);
+ return RES_OK;
+}
+
+res_T
+sdis_scene_get_aabb
+ (const struct sdis_scene* scn, double lower[], double upper[])
+{
+ float low[3], upp[3];
+ res_T res = RES_OK;
+ if(!scn || !lower || !upper) return RES_BAD_ARG;
+
+ if(scene_is_2d(scn)) {
+ res = s2d_scene_view_get_aabb(scn->s2d_view, low, upp);
+ if(res != RES_OK) return res;
+ d2_set_f2(lower, low);
+ d2_set_f2(upper, upp);
+ } else {
+ res = s3d_scene_view_get_aabb(scn->s3d_view, low, upp);
+ if(res != RES_OK) return res;
+ d3_set_f3(lower, low);
+ d3_set_f3(upper, upp);
}
+ return RES_OK;
+}
- return 0;
+/*******************************************************************************
+ * Local miscellaneous function
+ ******************************************************************************/
+const struct sdis_interface*
+scene_get_interface(const struct sdis_scene* scn, const unsigned iprim)
+{
+ ASSERT(scn && iprim < darray_interf_size_get(&scn->prim_interfaces));
+ return darray_interf_cdata_get(&scn->prim_interfaces)[iprim];
+}
+
+res_T
+scene_get_medium
+ (const struct sdis_scene* scn,
+ const double pos[],
+ const struct sdis_medium** out_medium)
+{
+ return scene_is_2d(scn)
+ ? scene_get_medium_2d(scn, pos, out_medium)
+ : scene_get_medium_3d(scn, pos, out_medium);
}
diff --git a/src/sdis_scene_c.h b/src/sdis_scene_c.h
@@ -19,34 +19,25 @@
#include <rsys/dynamic_array.h>
#include <rsys/ref_count.h>
-/* Forward declaration of external types */
-struct s3d_hit;
-
-/* Context used to wrap the user geometry to Star-3D. */
-struct geometry_context {
- void (*indices)(const size_t itri, size_t ids[3], void*);
- void (*position)(const size_t ivert, double pos[3], void*);
- void* data;
-};
-
static INLINE void
interface_init
(struct mem_allocator* allocator,
- struct sdis_interface** interface)
+ struct sdis_interface** interf)
{
(void)allocator;
- *interface = NULL;
+ *interf = NULL;
}
/* Declare the array of interfaces */
-#define DARRAY_NAME interface
+#define DARRAY_NAME interf
#define DARRAY_DATA struct sdis_interface*
#define DARRAY_FUNCTOR_INIT interface_init
#include <rsys/dynamic_array.h>
struct sdis_scene {
- struct darray_interface interfaces; /* List of interfaces own by the scene */
- struct darray_interface prim_interfaces; /* Per primitive interface */
+ struct darray_interf interfaces; /* List of interfaces own by the scene */
+ struct darray_interf prim_interfaces; /* Per primitive interface */
+ struct s2d_scene_view* s2d_view;
struct s3d_scene_view* s3d_view;
ref_T ref;
@@ -61,16 +52,15 @@ scene_get_interface
extern LOCAL_SYM res_T
scene_get_medium
(const struct sdis_scene* scene,
- const double position[3],
+ const double position[],
const struct sdis_medium** medium);
-extern LOCAL_SYM int
-hit_filter_function
- (const struct s3d_hit* hit,
- const float ray_org[3],
- const float ray_dir[3],
- void* ray_data, /* struct s3d_hit* */
- void* filter_data); /* NULL */
+static FINLINE int
+scene_is_2d(const struct sdis_scene* scn)
+{
+ ASSERT(scn && (scn->s2d_view || scn->s3d_view));
+ return scn->s2d_view != NULL;
+}
#endif /* SDIS_SCENE_C_H */
diff --git a/src/sdis_solve_probe.c b/src/sdis_solve_probe.c
@@ -16,454 +16,19 @@
#include "sdis.h"
#include "sdis_device_c.h"
#include "sdis_estimator_c.h"
-#include "sdis_interface_c.h"
-#include "sdis_medium_c.h"
-#include "sdis_scene_c.h"
+#include "sdis_solve_probe_Xd.h"
-#include <rsys/double2.h>
-#include <rsys/double3.h>
-#include <rsys/float3.h>
-#include <rsys/stretchy_array.h>
+/* Generate the 2D solver */
+#define SDIS_SOLVE_PROBE_DIMENSION 2
+#include "sdis_solve_probe_Xd.h"
-#include <star/s3d.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
-
-/* Current state of the random walk */
-struct 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 s3d_hit hit; /* Hit of the random walk */
-};
-static const struct rwalk RWALK_NULL = {
- SDIS_RWALK_VERTEX_NULL__, NULL, S3D_HIT_NULL__
-};
-
-struct 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 rwalk* rwalk,
- struct ssp_rng* rng,
- struct temperature* temp);
- double value; /* Current value of the temperature */
- int done;
-};
-static const struct temperature TEMPERATURE_NULL = { NULL, 0, 0 };
-
-static res_T
-boundary_temperature
- (const struct sdis_scene* scn,
- const double fp_to_meter,
- struct rwalk* rwalk,
- struct ssp_rng* rng,
- struct temperature* T);
-
-static res_T
-solid_temperature
- (const struct sdis_scene* scn,
- const double fp_to_meter,
- struct rwalk* rwalk,
- struct ssp_rng* rng,
- struct temperature* T);
-
-static res_T
-fluid_temperature
- (const struct sdis_scene* scn,
- const double fp_to_meter,
- struct rwalk* rwalk,
- struct ssp_rng* rng,
- struct temperature* T);
-
-/*******************************************************************************
- * Helper functions
- ******************************************************************************/
-static FINLINE void
-move_pos(double pos[3], const float dir[3], const float delta)
-{
- ASSERT(pos && dir);
- pos[0] += dir[0] * delta;
- pos[1] += dir[1] * delta;
- pos[2] += dir[2] * delta;
-}
-
-/* Check that the interface fragment is consistent with the current state of
- * the random walk */
-static INLINE int
-check_rwalk_fragment_consistency
- (const struct rwalk* rwalk,
- const struct sdis_interface_fragment* frag)
-{
- double N[3];
- double uv[2];
- ASSERT(rwalk && frag);
- d3_normalize(N, d3_set_f3(N, rwalk->hit.normal));
- d2_set_f2(uv, rwalk->hit.uv);
- return !S3D_HIT_NONE(&rwalk->hit)
- && d3_eq_eps(rwalk->vtx.P, frag->P, 1.e-6)
- && d3_eq_eps(N, frag->Ng, 1.e-6)
- && d2_eq_eps(uv, frag->uv, 1.e-6)
- && ( (IS_INF(rwalk->vtx.time) && IS_INF(frag->time))
- || eq_eps(rwalk->vtx.time, frag->time, 1.e-6));
-}
-
-res_T
-fluid_temperature
- (const struct sdis_scene* scn,
- const double fp_to_meter,
- struct rwalk* rwalk,
- struct ssp_rng* rng,
- struct 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
-solid_solid_boundary_temperature
- (const struct sdis_scene* scn,
- const double fp_to_meter,
- const struct sdis_interface_fragment* frag,
- struct rwalk* rwalk,
- struct ssp_rng* rng,
- struct temperature* T)
-{
- const struct sdis_interface* interface = 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[3], dir[3], range[2];
- ASSERT(scn && fp_to_meter > 0 && frag && rwalk && rng && T);
- ASSERT(check_rwalk_fragment_consistency(rwalk, frag));
- (void)frag;
-
- /* Retrieve the current boundary media */
- interface = scene_get_interface(scn, rwalk->hit.prim.prim_id);
- solid_front = interface_get_medium(interface, SDIS_FRONT);
- solid_back = interface_get_medium(interface, 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);
- f3_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;
- f3_minus(dir, dir);
- rwalk->mdm = solid_back;
- }
-
- /* "Reinject" the path into the solid along the surface normal. */
- f3_set_d3(pos, rwalk->vtx.P);
- range[0] = 0, range[1] = (float)delta_boundary*RAY_RANGE_MAX_SCALE;
- S3D(scene_view_trace_ray
- (scn->s3d_view, pos, dir, range, &rwalk->hit, &rwalk->hit));
- if(!S3D_HIT_NONE(&rwalk->hit)) delta_boundary = rwalk->hit.distance * 0.5;
- move_pos(rwalk->vtx.P, dir, (float)delta_boundary);
-
- /* Switch in solid random walk */
- T->func = solid_temperature;
-}
-
-static void
-solid_fluid_boundary_temperature
- (const struct sdis_scene* scn,
- const double fp_to_meter,
- const struct sdis_interface_fragment* frag,
- struct rwalk* rwalk,
- struct ssp_rng* rng,
- struct temperature* T)
-{
- const struct sdis_interface* interface = 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[3], pos[3], range[2];
-
- ASSERT(scn && fp_to_meter > 0 && rwalk && rng && T);
- ASSERT(check_rwalk_fragment_consistency(rwalk, frag));
-
- /* Retrieve the solid and the fluid split by the boundary */
- interface = scene_get_interface(scn, rwalk->hit.prim.prim_id);
- mdm_front = interface_get_medium(interface, SDIS_FRONT);
- mdm_back = interface_get_medium(interface, 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(interface, 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 = fluid_temperature;
- } else { /* Solid random walk */
- rwalk->mdm = solid;
- f3_normalize(dir, rwalk->hit.normal);
- if(solid == mdm_back) f3_minus(dir, dir);
-
- /* "Reinject" the random walk into the solid */
- f3_set_d3(pos, rwalk->vtx.P);
- range[0] = 0, range[1] = (float)delta_boundary*RAY_RANGE_MAX_SCALE;
- S3D(scene_view_trace_ray
- (scn->s3d_view, pos, dir, range, &rwalk->hit, &rwalk->hit));
- if(!S3D_HIT_NONE(&rwalk->hit)) delta_boundary = rwalk->hit.distance * 0.5;
- move_pos(rwalk->vtx.P, dir, (float)delta_boundary);
-
- /* Switch in solid random walk */
- T->func = solid_temperature;
- }
-}
-
-res_T
-boundary_temperature
- (const struct sdis_scene* scn,
- const double fp_to_meter,
- struct rwalk* rwalk,
- struct ssp_rng* rng,
- struct temperature* T)
-{
- struct sdis_interface_fragment frag = SDIS_INTERFACE_FRAGMENT_NULL;
- const struct sdis_interface* interface = 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(!S3D_HIT_NONE(&rwalk->hit));
-
- setup_interface_fragment(&frag, &rwalk->vtx, &rwalk->hit);
-
- /* Retrieve the current interface */
- interface = scene_get_interface(scn, rwalk->hit.prim.prim_id);
-
- /* Check if the boundary condition is known */
- tmp = interface_get_temperature(interface, &frag);
- if(tmp >= 0) {
- T->value += tmp;
- T->done = 1;
- return RES_OK;
- }
-
- mdm_front = interface_get_medium(interface, SDIS_FRONT);
- mdm_back = interface_get_medium(interface, SDIS_BACK);
-
- if(mdm_front->type == mdm_back->type) {
- solid_solid_boundary_temperature(scn, fp_to_meter, &frag, rwalk, rng, T);
- } else {
- solid_fluid_boundary_temperature(scn, fp_to_meter, &frag, rwalk, rng, T);
- }
- return RES_OK;
-}
+/* Generate the 3D solver */
+#define SDIS_SOLVE_PROBE_DIMENSION 3
+#include "sdis_solve_probe_Xd.h"
-res_T
-solid_temperature
- (const struct sdis_scene* scn,
- const double fp_to_meter,
- struct rwalk* rwalk,
- struct ssp_rng* rng,
- struct temperature* T)
-{
- double position_start[3];
- 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 */
- d3_set(position_start, rwalk->vtx.P);
-
- do { /* Solid random walk */
- struct s3d_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[3], dir1[3];
- float org[3];
-
- /* 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);
-
- /* Sample a direction around 4PI */
- ssp_ran_sphere_uniform_float(rng, dir0, NULL);
-
- /* Trace a ray along the sampled direction and its opposite to check if a
- * surface is hit in [0, delta_solid]. */
- f3_set_d3(org, rwalk->vtx.P);
- f3_minus(dir1, dir0);
- hit0 = hit1 = S3D_HIT_NULL;
- range[0] = 0.f, range[1] = delta_solid*RAY_RANGE_MAX_SCALE;
- S3D(scene_view_trace_ray(scn->s3d_view, org, dir0, range, NULL, &hit0));
- S3D(scene_view_trace_ray(scn->s3d_view, org, dir1, range, NULL, &hit1));
-
- if(S3D_HIT_NONE(&hit0) && S3D_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 = (6 * 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 ? S3D_HIT_NULL : hit0;
-
- /* Update the random walk position */
- move_pos(rwalk->vtx.P, dir0, delta);
-
- /* Fetch the current medium */
- if(S3D_HIT_NONE(&rwalk->hit)) {
- CHK(scene_get_medium(scn, rwalk->vtx.P, &mdm) == RES_OK);
- } else {
- const struct sdis_interface* interface;
- interface = scene_get_interface(scn, rwalk->hit.prim.prim_id);
- mdm = interface_get_medium
- (interface,
- f3_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"
- " start position: %g %g %g; current position: %g %g %g\n",
- FUNC_NAME, SPLIT3(position_start), SPLIT3(rwalk->vtx.P));
- return RES_BAD_OP;
- }
-
- /* Keep going while the solid random walk does not hit an interface */
- } while(S3D_HIT_NONE(&rwalk->hit));
-
- T->func = boundary_temperature;
- return RES_OK;
-}
-
-static res_T
-compute_temperature
- (struct sdis_scene* scn,
- const double fp_to_meter,
- struct rwalk* rwalk,
- struct ssp_rng* rng,
- struct temperature* T)
-{
-#ifndef NDEBUG
- struct temperature* stack = NULL;
- size_t istack;
-#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;
-}
+#include <star/ssp.h>
+#include <omp.h>
-/*******************************************************************************
- * Exported functions
- ******************************************************************************/
res_T
sdis_solve_probe
(struct sdis_scene* scn,
@@ -475,11 +40,14 @@ sdis_solve_probe
{
const struct sdis_medium* medium = NULL;
struct sdis_estimator* estimator = NULL;
- struct ssp_rng* rng = NULL;
+ struct ssp_rng_proxy* rng_proxy = NULL;
+ struct ssp_rng** rngs = NULL;
double weight = 0;
double sqr_weight = 0;
size_t irealisation = 0;
- res_T res = RES_OK;
+ 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) {
@@ -487,56 +55,79 @@ sdis_solve_probe
goto error;
}
- res = scene_get_medium(scn, position, &medium);
+ /* 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;
- res = ssp_rng_create(scn->dev->allocator, &ssp_rng_mt19937_64, &rng);
- 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;
- FOR_EACH(irealisation, 0, nrealisations) {
- struct rwalk rwalk = RWALK_NULL;
- struct temperature T = TEMPERATURE_NULL;
+ /* Retrieve the medium in which the submitted position lies */
+ res = scene_get_medium(scn, position, &medium);
+ if(res != RES_OK) goto error;
- switch(medium->type) {
- case SDIS_MEDIUM_FLUID: T.func = fluid_temperature; break;
- case SDIS_MEDIUM_SOLID: T.func = solid_temperature; break;
- default: FATAL("Unreachable code\n"); break;
- }
+ /* 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];
- rwalk.vtx.P[0] = position[0];
- rwalk.vtx.P[1] = position[1];
- rwalk.vtx.P[2] = position[2];
- rwalk.vtx.time = time;
- rwalk.hit = S3D_HIT_NULL;
- rwalk.mdm = medium;
+ if(ATOMIC_GET(&res) != RES_OK) continue; /* An error occured */
- res = compute_temperature(scn, fp_to_meter, &rwalk, rng, &T);
- if(res != RES_OK) {
- if(res == RES_BAD_OP) {
+ 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 {
- goto error;
+ ATOMIC_SET(&res, res_local);
+ continue;
}
} else {
- weight += T.value;
- sqr_weight += T.value*T.value;
- ++estimator->nrealisations;
+ weight += w;
+ sqr_weight += w*w;
+ ++N;
}
}
- estimator->temperature.E = weight / (double)estimator->nrealisations;
+ estimator->nrealisations = N;
+ estimator->temperature.E = weight / (double)N;
estimator->temperature.V =
- sqr_weight / (double)estimator->nrealisations
+ sqr_weight / (double)N
- estimator->temperature.E * estimator->temperature.E;
- estimator->temperature.SE =
- sqrt(estimator->temperature.V / (double)estimator->nrealisations);
+ estimator->temperature.SE = sqrt(estimator->temperature.V / (double)N);
exit:
- if(rng) SSP(rng_ref_put(rng));
+ 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;
+ return (res_T)res;
error:
if(estimator) {
SDIS(estimator_ref_put(estimator));
@@ -544,3 +135,4 @@ error:
}
goto exit;
}
+
diff --git a/src/sdis_solve_probe_Xd.h b/src/sdis_solve_probe_Xd.h
@@ -0,0 +1,565 @@
+/* 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_interface.c b/src/test_sdis_interface.c
@@ -23,7 +23,7 @@ main(int argc, char** argv)
struct sdis_device* dev = NULL;
struct sdis_medium* fluid = NULL;
struct sdis_medium* solid = NULL;
- struct sdis_interface* interface = NULL;
+ struct sdis_interface* interf = NULL;
struct sdis_fluid_shader fluid_shader = DUMMY_FLUID_SHADER;
struct sdis_solid_shader solid_shader = DUMMY_SOLID_SHADER;
struct sdis_interface_shader shader = DUMMY_INTERFACE_SHADER;
@@ -53,42 +53,42 @@ main(int argc, char** argv)
CHK(CREATE(dev, NULL, fluid, &shader, NULL, NULL) == RES_BAD_ARG);
CHK(CREATE(NULL, solid, fluid, &shader, NULL, NULL) == RES_BAD_ARG);
CHK(CREATE(dev, solid, fluid, &shader, NULL, NULL) == RES_BAD_ARG);
- CHK(CREATE(NULL, NULL, NULL, NULL, NULL, &interface) == RES_BAD_ARG);
- CHK(CREATE(dev, NULL, NULL, NULL, NULL, &interface) == RES_BAD_ARG);
- CHK(CREATE(NULL, solid, NULL, NULL, NULL, &interface) == RES_BAD_ARG);
- CHK(CREATE(dev, solid, NULL, NULL, NULL, &interface) == RES_BAD_ARG);
- CHK(CREATE(NULL, NULL, fluid, NULL, NULL, &interface) == RES_BAD_ARG);
- CHK(CREATE(dev, NULL, fluid, NULL, NULL, &interface) == RES_BAD_ARG);
- CHK(CREATE(NULL, solid, fluid, NULL, NULL, &interface) == RES_BAD_ARG);
- CHK(CREATE(dev, solid, fluid, NULL, NULL, &interface) == RES_BAD_ARG);
- CHK(CREATE(NULL, NULL, NULL, &shader, NULL, &interface) == RES_BAD_ARG);
- CHK(CREATE(dev, NULL, NULL, &shader, NULL, &interface) == RES_BAD_ARG);
- CHK(CREATE(NULL, solid, NULL, &shader, NULL, &interface) == RES_BAD_ARG);
- CHK(CREATE(dev, solid, NULL, &shader, NULL, &interface) == RES_BAD_ARG);
- CHK(CREATE(NULL, NULL, fluid, &shader, NULL, &interface) == RES_BAD_ARG);
- CHK(CREATE(dev, NULL, fluid, &shader, NULL, &interface) == RES_BAD_ARG);
- CHK(CREATE(NULL, solid, fluid, &shader, NULL, &interface) == RES_BAD_ARG);
- CHK(CREATE(dev, solid, fluid, &shader, NULL, &interface) == RES_OK);
+ CHK(CREATE(NULL, NULL, NULL, NULL, NULL, &interf) == RES_BAD_ARG);
+ CHK(CREATE(dev, NULL, NULL, NULL, NULL, &interf) == RES_BAD_ARG);
+ CHK(CREATE(NULL, solid, NULL, NULL, NULL, &interf) == RES_BAD_ARG);
+ CHK(CREATE(dev, solid, NULL, NULL, NULL, &interf) == RES_BAD_ARG);
+ CHK(CREATE(NULL, NULL, fluid, NULL, NULL, &interf) == RES_BAD_ARG);
+ CHK(CREATE(dev, NULL, fluid, NULL, NULL, &interf) == RES_BAD_ARG);
+ CHK(CREATE(NULL, solid, fluid, NULL, NULL, &interf) == RES_BAD_ARG);
+ CHK(CREATE(dev, solid, fluid, NULL, NULL, &interf) == RES_BAD_ARG);
+ CHK(CREATE(NULL, NULL, NULL, &shader, NULL, &interf) == RES_BAD_ARG);
+ CHK(CREATE(dev, NULL, NULL, &shader, NULL, &interf) == RES_BAD_ARG);
+ CHK(CREATE(NULL, solid, NULL, &shader, NULL, &interf) == RES_BAD_ARG);
+ CHK(CREATE(dev, solid, NULL, &shader, NULL, &interf) == RES_BAD_ARG);
+ CHK(CREATE(NULL, NULL, fluid, &shader, NULL, &interf) == RES_BAD_ARG);
+ CHK(CREATE(dev, NULL, fluid, &shader, NULL, &interf) == RES_BAD_ARG);
+ CHK(CREATE(NULL, solid, fluid, &shader, NULL, &interf) == RES_BAD_ARG);
+ CHK(CREATE(dev, solid, fluid, &shader, NULL, &interf) == RES_OK);
CHK(sdis_interface_ref_get(NULL) == RES_BAD_ARG);
- CHK(sdis_interface_ref_get(interface) == RES_OK);
+ CHK(sdis_interface_ref_get(interf) == RES_OK);
CHK(sdis_interface_ref_put(NULL) == RES_BAD_ARG);
- CHK(sdis_interface_ref_put(interface) == RES_OK);
- CHK(sdis_interface_ref_put(interface) == RES_OK);
+ CHK(sdis_interface_ref_put(interf) == RES_OK);
+ CHK(sdis_interface_ref_put(interf) == RES_OK);
- CHK(CREATE(dev, solid, solid, &shader, NULL, &interface) == RES_BAD_ARG);
+ CHK(CREATE(dev, solid, solid, &shader, NULL, &interf) == RES_BAD_ARG);
shader.convection_coef = NULL;
- CHK(CREATE(dev, solid, solid, &shader, NULL, &interface) == RES_OK);
- CHK(sdis_interface_ref_put(interface) == RES_OK);
+ CHK(CREATE(dev, solid, solid, &shader, NULL, &interf) == RES_OK);
+ CHK(sdis_interface_ref_put(interf) == RES_OK);
shader.temperature = NULL;
- CHK(CREATE(dev, solid, solid, &shader, NULL, &interface) == RES_OK);
- CHK(sdis_interface_ref_put(interface) == RES_OK);
+ CHK(CREATE(dev, solid, solid, &shader, NULL, &interf) == RES_OK);
+ CHK(sdis_interface_ref_put(interf) == RES_OK);
- CHK(CREATE(dev, solid, fluid, &shader, NULL, &interface) == RES_BAD_ARG);
+ 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, &interface) == RES_OK);
- CHK(sdis_interface_ref_put(interface) == RES_OK);
+ CHK(CREATE(dev, solid, fluid, &shader, NULL, &interf) == RES_OK);
+ CHK(sdis_interface_ref_put(interf) == RES_OK);
#undef CREATE
CHK(sdis_device_ref_put(dev) == RES_OK);
diff --git a/src/test_sdis_scene.c b/src/test_sdis_scene.c
@@ -20,7 +20,7 @@
struct context {
const double* positions;
const size_t* indices;
- struct sdis_interface* interface;
+ struct sdis_interface* interf;
};
static INLINE void
@@ -52,7 +52,7 @@ get_interface(const size_t itri, struct sdis_interface** bound, void* context)
CHK(ctx != NULL);
CHK(itri < box_ntriangles);
CHK(bound != NULL);
- *bound = ctx->interface;
+ *bound = ctx->interf;
}
int
@@ -62,7 +62,7 @@ main(int argc, char** argv)
struct sdis_device* dev = NULL;
struct sdis_medium* solid = NULL;
struct sdis_medium* fluid = NULL;
- struct sdis_interface* interface = NULL;
+ struct sdis_interface* interf = NULL;
struct sdis_scene* scn = NULL;
struct sdis_fluid_shader fluid_shader = DUMMY_FLUID_SHADER;
struct sdis_solid_shader solid_shader = DUMMY_SOLID_SHADER;
@@ -78,11 +78,11 @@ main(int argc, char** argv)
CHK(sdis_fluid_create(dev, &fluid_shader, NULL, &fluid) == RES_OK);
CHK(sdis_solid_create(dev, &solid_shader, NULL, &solid) == RES_OK);
CHK(sdis_interface_create
- (dev, solid, fluid, &interface_shader, NULL, &interface) == RES_OK);
+ (dev, solid, fluid, &interface_shader, NULL, &interf) == RES_OK);
ctx.positions = box_vertices;
ctx.indices = box_indices;
- ctx.interface = interface;
+ ctx.interf = interf;
ntris = box_ntriangles;
npos = box_nvertices;
@@ -122,7 +122,7 @@ main(int argc, char** argv)
CHK(sdis_scene_ref_put(scn) == RES_OK);
CHK(sdis_device_ref_put(dev) == RES_OK);
- CHK(sdis_interface_ref_put(interface) == RES_OK);
+ CHK(sdis_interface_ref_put(interf) == RES_OK);
CHK(sdis_medium_ref_put(solid) == RES_OK);
CHK(sdis_medium_ref_put(fluid) == RES_OK);
diff --git a/src/test_sdis_solve_probe.c b/src/test_sdis_solve_probe.c
@@ -38,7 +38,7 @@
struct context {
const double* positions;
const size_t* indices;
- struct sdis_interface* interface;
+ struct sdis_interface* interf;
};
static void
@@ -64,7 +64,7 @@ get_interface(const size_t itri, struct sdis_interface** bound, void* context)
{
struct context* ctx = context;
(void)itri;
- *bound = ctx->interface;
+ *bound = ctx->interf;
}
/*******************************************************************************
@@ -144,7 +144,7 @@ solid_get_temperature
/*******************************************************************************
* Interface
******************************************************************************/
-struct interface {
+struct interf {
double hc;
};
@@ -153,7 +153,7 @@ interface_get_convection_coef
(const struct sdis_interface_fragment* frag, struct sdis_data* data)
{
CHK(data != NULL && frag != NULL);
- return ((const struct interface*)sdis_data_cget(data))->hc;
+ return ((const struct interf*)sdis_data_cget(data))->hc;
}
/*******************************************************************************
@@ -167,7 +167,7 @@ main(int argc, char** argv)
struct sdis_device* dev = NULL;
struct sdis_medium* solid = NULL;
struct sdis_medium* fluid = NULL;
- struct sdis_interface* interface = NULL;
+ struct sdis_interface* interf = NULL;
struct sdis_scene* scn = NULL;
struct sdis_data* data = NULL;
struct sdis_estimator* estimator = NULL;
@@ -177,7 +177,7 @@ main(int argc, char** argv)
struct context ctx;
struct fluid* fluid_param;
struct solid* solid_param;
- struct interface* interface_param;
+ struct interf* interface_param;
double pos[3];
double time;
double ref;
@@ -218,14 +218,14 @@ main(int argc, char** argv)
CHK(sdis_data_ref_put(data) == RES_OK);
/* Create the solid/fluid interface */
- CHK(sdis_data_create(dev, sizeof(struct interface),
- ALIGNOF(struct interface), NULL, &data) == RES_OK);
+ CHK(sdis_data_create(dev, sizeof(struct interf),
+ ALIGNOF(struct interf), NULL, &data) == RES_OK);
interface_param = sdis_data_get(data);
interface_param->hc = 0.5;
interface_shader.convection_coef = interface_get_convection_coef;
interface_shader.temperature = NULL;
CHK(sdis_interface_create
- (dev, solid, fluid, &interface_shader, data, &interface) == RES_OK);
+ (dev, solid, fluid, &interface_shader, data, &interf) == RES_OK);
CHK(sdis_data_ref_put(data) == RES_OK);
/* Release the media */
@@ -235,11 +235,11 @@ main(int argc, char** argv)
/* Create the scene */
ctx.positions = box_vertices;
ctx.indices = box_indices;
- ctx.interface = interface;
+ ctx.interf = interf;
CHK(sdis_scene_create(dev, box_ntriangles, get_indices, get_interface,
box_nvertices, get_position, &ctx, &scn) == RES_OK);
- CHK(sdis_interface_ref_put(interface) == RES_OK);
+ CHK(sdis_interface_ref_put(interf) == RES_OK);
/* Test the solver */
pos[0] = 0.5;
diff --git a/src/test_sdis_solve_probe2.c b/src/test_sdis_solve_probe2.c
@@ -126,7 +126,7 @@ solid_get_delta_boundary
/*******************************************************************************
* Interface
******************************************************************************/
-struct interface {
+struct interf {
double temperature;
};
@@ -144,7 +144,7 @@ interface_get_temperature
(const struct sdis_interface_fragment* frag, struct sdis_data* data)
{
CHK(data != NULL && frag != NULL);
- return ((const struct interface*)sdis_data_cget(data))->temperature;
+ return ((const struct interf*)sdis_data_cget(data))->temperature;
}
/*******************************************************************************
@@ -169,7 +169,7 @@ main(int argc, char** argv)
struct sdis_interface_shader interface_shader = DUMMY_INTERFACE_SHADER;
struct sdis_interface* interfaces[12];
struct context ctx;
- struct interface* interface_param = NULL;
+ struct interf* interface_param = NULL;
double pos[3];
double time;
double ref;
@@ -202,8 +202,8 @@ main(int argc, char** argv)
(dev, solid, fluid, &interface_shader, NULL, &Tnone) == RES_OK);
/* Create the fluid/solid interface with a fixed temperature of 300K */
- CHK(sdis_data_create(dev, sizeof(struct interface),
- ALIGNOF(struct interface), NULL, &data) == RES_OK);
+ CHK(sdis_data_create(dev, sizeof(struct interf),
+ ALIGNOF(struct interf), NULL, &data) == RES_OK);
interface_param = sdis_data_get(data);
interface_param->temperature = 300;
interface_shader.convection_coef = null_convection_coef;
@@ -213,8 +213,8 @@ main(int argc, char** argv)
CHK(sdis_data_ref_put(data) == RES_OK);
/* Create the fluid/solid interface with a fixed temperature of 350K */
- CHK(sdis_data_create(dev, sizeof(struct interface),
- ALIGNOF(struct interface), NULL, &data) == RES_OK);
+ CHK(sdis_data_create(dev, sizeof(struct interf),
+ ALIGNOF(struct interf), NULL, &data) == RES_OK);
interface_param = sdis_data_get(data);
interface_param->temperature = 350;
interface_shader.convection_coef = null_convection_coef;
diff --git a/src/test_sdis_solve_probe2_2d.c b/src/test_sdis_solve_probe2_2d.c
@@ -0,0 +1,275 @@
+/* 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>
+
+/*
+ * The scene is composed of a solid square whose temperature is unknown. The
+ * convection coefficient with the surrounding fluid is null. The temperature
+ * is fixed at the left and right segment.
+ *
+ * (1,1)
+ * +-------+
+ * | |350K
+ * | |
+ * 300K| |
+ * +-------+
+ * (0,0)
+ */
+
+/*******************************************************************************
+ * Geometry
+ ******************************************************************************/
+struct context {
+ const double* positions;
+ const size_t* indices;
+ struct sdis_interface** interfaces; /* Per primitive interfaces */
+};
+
+static void
+get_indices(const size_t iseg, size_t ids[2], void* context)
+{
+ struct context* ctx = context;
+ ids[0] = ctx->indices[iseg*2+0];
+ ids[1] = ctx->indices[iseg*2+1];
+}
+
+static void
+get_position(const size_t ivert, double pos[2], void* context)
+{
+ struct context* ctx = context;
+ pos[0] = ctx->positions[ivert*2+0];
+ pos[1] = ctx->positions[ivert*2+1];
+}
+
+static void
+get_interface(const size_t iseg, struct sdis_interface** bound, void* context)
+{
+ struct context* ctx = context;
+ *bound = ctx->interfaces[iseg];
+}
+
+/*******************************************************************************
+ * Medium data
+ ******************************************************************************/
+static double
+temperature_unknown(const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ (void)data;
+ CHK(vtx != NULL);
+ return -1;
+}
+
+static double
+solid_get_calorific_capacity
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ (void)data;
+ CHK(vtx != NULL && data == NULL);
+ return 2.0;
+}
+
+static double
+solid_get_thermal_conductivity
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ (void)data;
+ CHK(vtx != NULL);
+ return 50.0;
+}
+
+static double
+solid_get_volumic_mass
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ (void)data;
+ CHK(vtx != NULL);
+ return 25.0;
+}
+
+static double
+solid_get_delta
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ (void)data;
+ CHK(vtx != NULL);
+ return 1.0/20.0;
+}
+
+static double
+solid_get_delta_boundary
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ (void)data;
+ CHK(vtx != NULL);
+ return 2.1/20.0;
+}
+
+/*******************************************************************************
+ * Interface
+ ******************************************************************************/
+struct interf {
+ double temperature;
+};
+
+static double
+null_convection_coef
+ (const struct sdis_interface_fragment* frag, struct sdis_data* data)
+{
+ CHK(frag != NULL);
+ (void)data;
+ return 0;
+}
+
+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;
+}
+
+/*******************************************************************************
+ * Test
+ ******************************************************************************/
+int
+main(int argc, char** argv)
+{
+ struct mem_allocator allocator;
+ struct sdis_mc T = SDIS_MC_NULL;
+ struct sdis_device* dev = NULL;
+ struct sdis_data* data = NULL;
+ struct sdis_estimator* estimator = NULL;
+ struct sdis_medium* solid = NULL;
+ struct sdis_medium* fluid = NULL;
+ struct sdis_interface* Tnone = NULL;
+ struct sdis_interface* T300 = NULL;
+ struct sdis_interface* T350 = NULL;
+ struct sdis_scene* scn = NULL;
+ struct sdis_fluid_shader fluid_shader = DUMMY_FLUID_SHADER;
+ struct sdis_solid_shader solid_shader = DUMMY_SOLID_SHADER;
+ struct sdis_interface_shader interface_shader = DUMMY_INTERFACE_SHADER;
+ struct sdis_interface* interfaces[4];
+ struct context ctx;
+ struct interf* interface_param = NULL;
+ double pos[2];
+ double time;
+ double ref;
+ const size_t N = 10000;
+ size_t nreals;
+ size_t nfails;
+ (void)argc, (void)argv;
+
+ CHK(mem_init_proxy_allocator(&allocator, &mem_default_allocator) == RES_OK);
+ CHK(sdis_device_create
+ (NULL, &allocator, SDIS_NTHREADS_DEFAULT, 0, &dev) == RES_OK);
+
+ /* Create the fluid medium */
+ fluid_shader.temperature = temperature_unknown;
+ CHK(sdis_fluid_create(dev, &fluid_shader, NULL, &fluid) == RES_OK);
+
+ /* Create the solid medium */
+ solid_shader.calorific_capacity = solid_get_calorific_capacity;
+ solid_shader.thermal_conductivity = solid_get_thermal_conductivity;
+ solid_shader.volumic_mass = solid_get_volumic_mass;
+ solid_shader.delta_solid = solid_get_delta;
+ solid_shader.delta_boundary = solid_get_delta_boundary;
+ solid_shader.temperature = temperature_unknown;
+ CHK(sdis_solid_create(dev, &solid_shader, NULL, &solid) == RES_OK);
+
+ /* Create the fluid/solid interface with no limit conidition */
+ interface_shader.convection_coef = null_convection_coef;
+ interface_shader.temperature = NULL;
+ CHK(sdis_interface_create
+ (dev, solid, fluid, &interface_shader, NULL, &Tnone) == RES_OK);
+
+ /* Create the fluid/solid interface with a fixed temperature of 300K */
+ CHK(sdis_data_create(dev, sizeof(struct interf),
+ ALIGNOF(struct interf), NULL, &data) == RES_OK);
+ interface_param = sdis_data_get(data);
+ interface_param->temperature = 300;
+ interface_shader.convection_coef = null_convection_coef;
+ interface_shader.temperature = interface_get_temperature;
+ CHK(sdis_interface_create
+ (dev, solid, fluid, &interface_shader, data, &T300) == RES_OK);
+ CHK(sdis_data_ref_put(data) == RES_OK);
+
+ /* Create the fluid/solid interface with a fixed temperature of 350K */
+ CHK(sdis_data_create(dev, sizeof(struct interf),
+ 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.temperature = interface_get_temperature;
+ CHK(sdis_interface_create
+ (dev, solid, fluid, &interface_shader, data, &T350) == RES_OK);
+ CHK(sdis_data_ref_put(data) == RES_OK);
+
+ /* Release the media */
+ CHK(sdis_medium_ref_put(solid) == RES_OK);
+ CHK(sdis_medium_ref_put(fluid) == RES_OK);
+
+ /* Setup the per primitive scene interfaces */
+ CHK(sizeof(interfaces)/sizeof(struct sdis_interface*) == square_nsegments);
+ interfaces[0] = Tnone; /* Bottom segment */
+ interfaces[1] = T300; /* Left segment */
+ interfaces[2] = Tnone; /* Top segment */
+ interfaces[3] = T350; /* Right segment */
+
+ /* Create the scene */
+ ctx.positions = square_vertices;
+ ctx.indices = square_indices;
+ ctx.interfaces = interfaces;
+ CHK(sdis_scene_2d_create(dev, square_nsegments, get_indices, get_interface,
+ square_nvertices, get_position, &ctx, &scn) == RES_OK);
+
+ /* Release the interfaces */
+ CHK(sdis_interface_ref_put(Tnone) == RES_OK);
+ CHK(sdis_interface_ref_put(T300) == RES_OK);
+ CHK(sdis_interface_ref_put(T350) == RES_OK);
+
+ /* Launch the solver */
+ 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_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);
+
+ /* Print the estimation results */
+ ref = 350 * pos[0] + (1-pos[0]) * 300;
+ printf("Temperature at (%g, %g) = %g ~ %g +/- %g\n",
+ SPLIT2(pos), ref, T.E, T.SE);
+ printf("#realisations: %lu; #failures: %lu\n",
+ (unsigned long)nreals, (unsigned long)nfails);
+
+ /* Check the results */
+ CHK(nfails + nreals == N);
+ CHK(eq_eps(T.E, ref, T.SE));
+
+ /* Release data */
+ CHK(sdis_estimator_ref_put(estimator) == RES_OK);
+ CHK(sdis_scene_ref_put(scn) == RES_OK);
+ CHK(sdis_device_ref_put(dev) == RES_OK);
+
+ check_memory_allocator(&allocator);
+ mem_shutdown_proxy_allocator(&allocator);
+ CHK(mem_allocated_size() == 0);
+ return 0;
+}
diff --git a/src/test_sdis_solve_probe3.c b/src/test_sdis_solve_probe3.c
@@ -148,7 +148,7 @@ solid_get_delta_boundary
/*******************************************************************************
* Interface
******************************************************************************/
-struct interface {
+struct interf {
double temperature;
};
@@ -166,7 +166,7 @@ interface_get_temperature
(const struct sdis_interface_fragment* frag, struct sdis_data* data)
{
CHK(data != NULL && frag != NULL);
- return ((const struct interface*)sdis_data_cget(data))->temperature;
+ return ((const struct interf*)sdis_data_cget(data))->temperature;
}
/*******************************************************************************
@@ -193,7 +193,7 @@ main(int argc, char** argv)
struct s3dut_mesh* msh = NULL;
struct s3dut_mesh_data msh_data;
struct context ctx = CONTEXT_NULL;
- struct interface* interface_param = NULL;
+ struct interf* interface_param = NULL;
double pos[3];
double time;
double ref;
@@ -229,8 +229,8 @@ main(int argc, char** argv)
(dev, solid, fluid, &interface_shader, NULL, &Tnone) == RES_OK);
/* Create the fluid/solid interface with a fixed temperature of 300K */
- CHK(sdis_data_create(dev, sizeof(struct interface),
- ALIGNOF(struct interface), NULL, &data) == RES_OK);
+ CHK(sdis_data_create(dev, sizeof(struct interf),
+ ALIGNOF(struct interf), NULL, &data) == RES_OK);
interface_param = sdis_data_get(data);
interface_param->temperature = 300;
interface_shader.convection_coef = null_convection_coef;
@@ -240,8 +240,8 @@ main(int argc, char** argv)
CHK(sdis_data_ref_put(data) == RES_OK);
/* Create the fluid/solid interface with a fixed temperature of 350K */
- CHK(sdis_data_create(dev, sizeof(struct interface),
- ALIGNOF(struct interface), NULL, &data) == RES_OK);
+ CHK(sdis_data_create(dev, sizeof(struct interf),
+ ALIGNOF(struct interf), NULL, &data) == RES_OK);
interface_param = sdis_data_get(data);
interface_param->temperature = 350;
interface_shader.convection_coef = null_convection_coef;
@@ -324,7 +324,7 @@ main(int argc, char** argv)
/* Check the results */
CHK(nfails + nreals == N);
- CHK(eq_eps(T.E, ref, T.SE));
+ CHK(eq_eps(T.E, ref, 2*T.SE));
/* Release data */
CHK(sdis_estimator_ref_put(estimator) == RES_OK);
diff --git a/src/test_sdis_solve_probe3_2d.c b/src/test_sdis_solve_probe3_2d.c
@@ -0,0 +1,329 @@
+/* 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>
+
+#include <rsys/stretchy_array.h>
+#include <rsys/math.h>
+
+/*
+ * The scene is composed of a solid square whose temperature is unknown. The
+ * convection coefficient with the surrounding fluid is null. The temperature
+ * is fixed at the left and right face. At the center of the cube there is a
+ * solid disk whose physical properties are the same of the solid square; i.e.
+ * the disk influences the random walks but not the result.
+ *
+ * (1,1)
+ * +--------------+
+ * | # # |
+ * | # # |350K
+ * | # # |
+ * | # # |
+ * 300K| # # |
+ * | # # |
+ * +--------------+
+ * (0,0)
+ */
+
+/*******************************************************************************
+ * Geometry
+ ******************************************************************************/
+struct context {
+ double* positions;
+ size_t* indices;
+ struct sdis_interface* solid_fluid_Tnone;
+ struct sdis_interface* solid_fluid_T300;
+ struct sdis_interface* solid_fluid_T350;
+ struct sdis_interface* solid_solid;
+};
+static const struct context CONTEXT_NULL = { NULL };
+
+static void
+get_indices(const size_t iseg, size_t ids[2], void* context)
+{
+ struct context* ctx = context;
+ ids[0] = ctx->indices[iseg*2+0];
+ ids[1] = ctx->indices[iseg*2+1];
+}
+
+static void
+get_position(const size_t ivert, double pos[2], void* context)
+{
+ struct context* ctx = context;
+ pos[0] = ctx->positions[ivert*2+0];
+ pos[1] = ctx->positions[ivert*2+1];
+}
+
+static void
+get_interface(const size_t iseg, struct sdis_interface** bound, void* context)
+{
+ struct context* ctx = context;
+ CHK(bound != NULL && context != NULL);
+
+ if(iseg == 1) { /* Square left segment */
+ *bound = ctx->solid_fluid_T300;
+ } else if(iseg == 3 ) { /* Square right segment */
+ *bound = ctx->solid_fluid_T350;
+ } else if(iseg < square_nsegments) { /* Square remaining segments */
+ *bound = ctx->solid_fluid_Tnone;
+ } else { /* Faces of the internal geometry */
+ *bound = ctx->solid_solid;
+ }
+}
+
+/*******************************************************************************
+ * Medium data
+ ******************************************************************************/
+static double
+temperature_unknown(const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ (void)data;
+ CHK(vtx != NULL);
+ return -1;
+}
+
+static double
+solid_get_calorific_capacity
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ (void)data;
+ CHK(vtx != NULL && data == NULL);
+ return 2.0;
+}
+
+static double
+solid_get_thermal_conductivity
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ (void)data;
+ CHK(vtx != NULL);
+ return 50.0;
+}
+
+static double
+solid_get_volumic_mass
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ (void)data;
+ CHK(vtx != NULL);
+ return 25.0;
+}
+
+static double
+solid_get_delta
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ (void)data;
+ CHK(vtx != NULL);
+ return 1.0/20.0;
+}
+
+static double
+solid_get_delta_boundary
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ (void)data;
+ CHK(vtx != NULL);
+ return 2.1/20.0;
+}
+
+/*******************************************************************************
+ * Interface
+ ******************************************************************************/
+struct interf {
+ double temperature;
+};
+
+static double
+null_convection_coef
+ (const struct sdis_interface_fragment* frag, struct sdis_data* data)
+{
+ CHK(frag != NULL);
+ (void)data;
+ return 0;
+}
+
+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;
+}
+
+/*******************************************************************************
+ * Test
+ ******************************************************************************/
+int
+main(int argc, char** argv)
+{
+ struct mem_allocator allocator;
+ struct sdis_mc T = SDIS_MC_NULL;
+ struct sdis_device* dev = NULL;
+ struct sdis_data* data = NULL;
+ struct sdis_estimator* estimator = NULL;
+ struct sdis_medium* solid = NULL;
+ struct sdis_medium* fluid = NULL;
+ struct sdis_interface* Tnone = NULL;
+ struct sdis_interface* T300 = NULL;
+ struct sdis_interface* T350 = NULL;
+ struct sdis_interface* solid_solid = NULL;
+ struct sdis_scene* scn = NULL;
+ struct sdis_fluid_shader fluid_shader = DUMMY_FLUID_SHADER;
+ struct sdis_solid_shader solid_shader = DUMMY_SOLID_SHADER;
+ struct sdis_interface_shader interface_shader = DUMMY_INTERFACE_SHADER;
+ struct context ctx = CONTEXT_NULL;
+ struct interf* interface_param = NULL;
+ double pos[2];
+ double time;
+ double ref;
+ const size_t N = 10000;
+ size_t nsegs;
+ size_t nverts;
+ size_t nreals;
+ size_t nfails;
+ size_t i;
+ (void)argc, (void)argv;
+
+ CHK(mem_init_proxy_allocator(&allocator, &mem_default_allocator) == RES_OK);
+ CHK(sdis_device_create
+ (NULL, &allocator, SDIS_NTHREADS_DEFAULT, 0, &dev) == RES_OK);
+
+ /* Create the fluid medium */
+ fluid_shader.temperature = temperature_unknown;
+ CHK(sdis_fluid_create(dev, &fluid_shader, NULL, &fluid) == RES_OK);
+
+ /* Create the solid medium */
+ solid_shader.calorific_capacity = solid_get_calorific_capacity;
+ solid_shader.thermal_conductivity = solid_get_thermal_conductivity;
+ solid_shader.volumic_mass = solid_get_volumic_mass;
+ solid_shader.delta_solid = solid_get_delta;
+ solid_shader.delta_boundary = solid_get_delta_boundary;
+ solid_shader.temperature = temperature_unknown;
+ CHK(sdis_solid_create(dev, &solid_shader, NULL, &solid) == RES_OK);
+
+ /* Create the fluid/solid interface with no limit conidition */
+ interface_shader.convection_coef = null_convection_coef;
+ interface_shader.temperature = NULL;
+ CHK(sdis_interface_create
+ (dev, solid, fluid, &interface_shader, NULL, &Tnone) == RES_OK);
+
+ /* Create the fluid/solid interface with a fixed temperature of 300K */
+ CHK(sdis_data_create(dev, sizeof(struct interf),
+ ALIGNOF(struct interf), NULL, &data) == RES_OK);
+ interface_param = sdis_data_get(data);
+ interface_param->temperature = 300;
+ interface_shader.convection_coef = null_convection_coef;
+ interface_shader.temperature = interface_get_temperature;
+ CHK(sdis_interface_create
+ (dev, solid, fluid, &interface_shader, data, &T300) == RES_OK);
+ CHK(sdis_data_ref_put(data) == RES_OK);
+
+ /* Create the fluid/solid interface with a fixed temperature of 350K */
+ CHK(sdis_data_create(dev, sizeof(struct interf),
+ 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.temperature = interface_get_temperature;
+ CHK(sdis_interface_create
+ (dev, solid, fluid, &interface_shader, data, &T350) == RES_OK);
+ CHK(sdis_data_ref_put(data) == RES_OK);
+
+ /* Create the solid/solid interface */
+ interface_shader.convection_coef = NULL;
+ interface_shader.temperature = NULL;
+ CHK(sdis_interface_create
+ (dev, solid, solid, &interface_shader, NULL, &solid_solid) == RES_OK);
+
+ /* Release the media */
+ CHK(sdis_medium_ref_put(solid) == RES_OK);
+ CHK(sdis_medium_ref_put(fluid) == RES_OK);
+
+ /* Register the square geometry */
+ FOR_EACH(i, 0, square_nvertices) {
+ sa_push(ctx.positions, square_vertices[i*2+0]);
+ sa_push(ctx.positions, square_vertices[i*2+1]);
+ }
+ FOR_EACH(i, 0, square_nsegments) {
+ sa_push(ctx.indices, square_indices[i*2+0]);
+ sa_push(ctx.indices, square_indices[i*2+1]);
+ }
+
+ /* Setup a disk at the center of the square */
+ nverts = 64;
+ FOR_EACH(i, 0, nverts) {
+ const float theta = (float)i * (2*(float)PI)/(float)nverts;
+ const float r = 0.25f; /* radius */
+ sa_push(ctx.positions, (float)cos(theta) * r + 0.5f);
+ sa_push(ctx.positions, (float)sin(theta) * r + 0.5f);
+ }
+ FOR_EACH(i, 0, nverts) {
+ sa_push(ctx.indices, i + square_nvertices);
+ sa_push(ctx.indices, (i+1)%nverts + square_nvertices);
+ }
+
+ /* Create the scene */
+ ctx.solid_fluid_Tnone = Tnone;
+ ctx.solid_fluid_T300 = T300;
+ ctx.solid_fluid_T350 = T350;
+ ctx.solid_solid = solid_solid;
+ nverts = sa_size(ctx.positions) / 2;
+ nsegs = sa_size(ctx.indices) / 2;
+ CHK(sdis_scene_2d_create(dev, nsegs, get_indices, get_interface, nverts,
+ get_position, &ctx, &scn) == RES_OK);
+
+ /* Release the scene data */
+ CHK(sdis_interface_ref_put(Tnone) == RES_OK);
+ CHK(sdis_interface_ref_put(T300) == RES_OK);
+ CHK(sdis_interface_ref_put(T350) == RES_OK);
+ CHK(sdis_interface_ref_put(solid_solid) == RES_OK);
+ sa_release(ctx.positions);
+ sa_release(ctx.indices);
+
+ /* Launch the solver */
+ 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_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);
+
+ /* Print the estimation results */
+ ref = 350 * pos[0] + (1-pos[0]) * 300;
+ printf("Temperature at (%g, %g) = %g ~ %g +/- %g\n",
+ SPLIT2(pos), ref, T.E, T.SE);
+ printf("#realisations: %lu; #failures: %lu\n",
+ (unsigned long)nreals, (unsigned long)nfails);
+
+ /* Check the results */
+ CHK(nfails + nreals == N);
+ CHK(eq_eps(T.E, ref, 3*T.SE));
+
+ /* Release data */
+ CHK(sdis_estimator_ref_put(estimator) == RES_OK);
+ CHK(sdis_scene_ref_put(scn) == RES_OK);
+ CHK(sdis_device_ref_put(dev) == RES_OK);
+
+ check_memory_allocator(&allocator);
+ mem_shutdown_proxy_allocator(&allocator);
+ CHK(mem_allocated_size() == 0);
+ return 0;
+
+}
diff --git a/src/test_sdis_solve_probe_2d.c b/src/test_sdis_solve_probe_2d.c
@@ -0,0 +1,222 @@
+/* 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>
+
+/*
+ * The scene is composed of a solid square with unknown temperature. The
+ * surrounding fluid has a fixed constant temperature.
+ *
+ * (1,1)
+ * +-------+ _\
+ * | | / /
+ * | | \__/ 300K
+ * | |
+ * +-------+
+ * (0,0)
+ */
+
+/*******************************************************************************
+ * Geometry
+ ******************************************************************************/
+struct context {
+ const double* positions;
+ const size_t* indices;
+ struct sdis_interface* interf;
+};
+
+static void
+get_indices(const size_t iseg, size_t ids[2], void* context)
+{
+ struct context* ctx = context;
+ ids[0] = ctx->indices[iseg*2+0];
+ ids[1] = ctx->indices[iseg*2+1];
+}
+
+static void
+get_position(const size_t ivert, double pos[2], void* context)
+{
+ struct context* ctx = context;
+ pos[0] = ctx->positions[ivert*2+0];
+ pos[1] = ctx->positions[ivert*2+1];
+}
+
+static void
+get_interface(const size_t iseg, struct sdis_interface** bound, void* context)
+{
+ struct context* ctx = context;
+ (void)iseg;
+ *bound = ctx->interf;
+}
+
+/*******************************************************************************
+ * Media & interface
+ ******************************************************************************/
+static double
+fluid_get_temperature
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ (void)vtx, (void)data;
+ return 300.0;
+}
+
+static double
+solid_get_calorific_capacity
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ (void)vtx, (void)data;
+ return 1.0;
+}
+
+static double
+solid_get_thermal_conductivity
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ (void)vtx, (void)data;
+ return 0.1;
+}
+
+static double
+solid_get_volumic_mass
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ (void)vtx, (void)data;
+ return 1.0;
+}
+
+static double
+solid_get_delta
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ (void)vtx, (void)data;
+ return 1.0/20.0;
+}
+
+static double
+solid_get_delta_boundary
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ (void)vtx, (void)data;
+ return 2.1/20.0;
+}
+
+static double
+solid_get_temperature
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ (void)vtx, (void)data;
+ return -1;
+}
+
+static double
+interface_get_convection_coef
+ (const struct sdis_interface_fragment* frag, struct sdis_data* data)
+{
+ (void)frag, (void)data;
+ return 0.5;
+}
+
+/*******************************************************************************
+ * Main test
+ ******************************************************************************/
+int
+main(int argc, char** argv)
+{
+ struct mem_allocator allocator;
+ struct sdis_mc T = SDIS_MC_NULL;
+ struct sdis_device* dev = NULL;
+ struct sdis_medium* solid = NULL;
+ struct sdis_medium* fluid = NULL;
+ struct sdis_interface* interf = NULL;
+ struct sdis_scene* scn = NULL;
+ struct sdis_estimator* estimator = NULL;
+ struct sdis_fluid_shader fluid_shader = DUMMY_FLUID_SHADER;
+ struct sdis_solid_shader solid_shader = DUMMY_SOLID_SHADER;
+ struct sdis_interface_shader interface_shader = DUMMY_INTERFACE_SHADER;
+ struct context ctx;
+ double pos[2];
+ double time;
+ double ref;
+ const size_t N = 1000;
+ size_t nreals;
+ size_t nfails;
+ (void)argc, (void)argv;
+
+ CHK(mem_init_proxy_allocator(&allocator, &mem_default_allocator) == RES_OK);
+ CHK(sdis_device_create
+ (NULL, &allocator, SDIS_NTHREADS_DEFAULT, 0, &dev) == RES_OK);
+
+ /* Create the fluid medium */
+ fluid_shader.temperature = fluid_get_temperature;
+ CHK(sdis_fluid_create(dev, &fluid_shader, NULL, &fluid) == RES_OK);
+
+ /* Create the solid medium */
+ solid_shader.calorific_capacity = solid_get_calorific_capacity;
+ solid_shader.thermal_conductivity = solid_get_thermal_conductivity;
+ solid_shader.volumic_mass = solid_get_volumic_mass;
+ solid_shader.delta_solid = solid_get_delta;
+ solid_shader.delta_boundary = solid_get_delta_boundary;
+ solid_shader.temperature = solid_get_temperature;
+ CHK(sdis_solid_create(dev, &solid_shader, NULL, &solid) == RES_OK);
+
+ /* Create the solid/fluid interface */
+ interface_shader.convection_coef = interface_get_convection_coef;
+ interface_shader.temperature = NULL;
+ CHK(sdis_interface_create
+ (dev, solid, fluid, &interface_shader, NULL, &interf) == RES_OK);
+
+ /* Release the media */
+ CHK(sdis_medium_ref_put(solid) == RES_OK);
+ CHK(sdis_medium_ref_put(fluid) == RES_OK);
+
+ /* Create the scene */
+ ctx.positions = square_vertices;
+ ctx.indices = square_indices;
+ ctx.interf = interf;
+ CHK(sdis_scene_2d_create(dev, square_nsegments, get_indices, get_interface,
+ square_nvertices, get_position, &ctx, &scn) == RES_OK);
+
+ CHK(sdis_interface_ref_put(interf) == RES_OK);
+
+ /* Test the solver */
+ 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_estimator_get_realisation_count(estimator, &nreals) == RES_OK);
+ CHK(sdis_estimator_get_failure_count(estimator, &nfails) == RES_OK);
+
+ CHK(nfails + nreals == N);
+
+ CHK(sdis_estimator_get_temperature(estimator, &T) == RES_OK);
+
+ ref = 300;
+ printf("Temperature at (%g, %g) = %g ~ %g +/- %g\n",
+ SPLIT2(pos), ref, T.E, T.SE);
+ CHK(eq_eps(T.E, ref, T.SE));
+
+ CHK(sdis_estimator_ref_put(estimator) == RES_OK);
+
+ CHK(sdis_scene_ref_put(scn) == RES_OK);
+ CHK(sdis_device_ref_put(dev) == RES_OK);
+
+ check_memory_allocator(&allocator);
+ mem_shutdown_proxy_allocator(&allocator);
+ CHK(mem_allocated_size() == 0);
+ return 0;
+}
diff --git a/src/test_sdis_utils.h b/src/test_sdis_utils.h
@@ -20,7 +20,7 @@
#include <stdio.h>
/*******************************************************************************
- * Geometry
+ * Box geometry
******************************************************************************/
static const double box_vertices[8/*#vertices*/*3/*#coords per vertex*/] = {
0.0, 0.0, 0.0,
@@ -55,6 +55,25 @@ static const size_t box_indices[12/*#triangles*/*3/*#indices per triangle*/] = {
static const size_t box_ntriangles = sizeof(box_indices) / sizeof(size_t[3]);
/*******************************************************************************
+ * Square geometry
+ ******************************************************************************/
+static const double square_vertices[4/*#vertices*/*2/*#coords per vertex*/] = {
+ 1.0, 0.0,
+ 0.0, 0.0,
+ 0.0, 1.0,
+ 1.0, 1.0
+};
+static const size_t square_nvertices = sizeof(square_vertices)/sizeof(double[2]);
+
+static const size_t square_indices[4/*#triangles*/*2/*#indices per segment*/]= {
+ 0, 1, /* Bottom */
+ 1, 2, /* Left */
+ 2, 3, /* Top */
+ 3, 0 /* Right */
+};
+static const size_t square_nsegments = sizeof(square_indices)/sizeof(size_t[2]);
+
+/*******************************************************************************
* Medium & interface
******************************************************************************/
static INLINE double
@@ -121,6 +140,27 @@ dump_mesh
}
static INLINE void
+dump_segments
+ (FILE* stream,
+ const double* pos,
+ const size_t npos,
+ const size_t* ids,
+ const size_t nids)
+{
+ size_t i;
+ CHK(pos != NULL && npos != 0);
+ CHK(ids != NULL && nids != 0);
+ FOR_EACH(i, 0, npos) {
+ fprintf(stream, "v %g %g 0\n", SPLIT2(pos+i*2));
+ }
+ FOR_EACH(i, 0, nids) {
+ fprintf(stream, "l %lu %lu\n",
+ (unsigned long)(ids[i*2+0] + 1),
+ (unsigned long)(ids[i*2+1] + 1));
+ }
+}
+
+static INLINE void
check_memory_allocator(struct mem_allocator* allocator)
{
if(MEM_ALLOCATED_SIZE(allocator)) {
