commit 7afcf6b3a1a3649c1a4d3bbed51d2d4204aa4411
parent 06b2fb21996751b6c9d7d0cda8543b6229b26d8e
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Fri, 5 Oct 2018 10:12:20 +0200
Merge branch 'feature_convection' into develop
Diffstat:
21 files changed, 2000 insertions(+), 771 deletions(-)
diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt
@@ -28,18 +28,20 @@ CMAKE_DEPENDENT_OPTION(ALL_TESTS
################################################################################
# Check dependencies
################################################################################
-find_package(RCMake 0.3 REQUIRED)
-find_package(Star2D 0.1 REQUIRED)
-find_package(Star3D 0.4 REQUIRED)
+find_package(RCMake 0.4 REQUIRED)
+find_package(Star2D 0.1.1 REQUIRED)
+find_package(Star3D 0.5 REQUIRED)
find_package(StarSP 0.7 REQUIRED)
-find_package(RSys 0.6 REQUIRED)
-find_package(OpenMP 1.2 REQUIRED)
+find_package(StarEnc 0.2.0 REQUIRED)
+find_package(StarEnc2D 0.2.0 REQUIRED)
+find_package(RSys 0.6.1 REQUIRED)
+find_package(OpenMP 2.0 REQUIRED)
set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} ${RCMAKE_SOURCE_DIR})
include(rcmake)
include(rcmake_runtime)
-rcmake_append_runtime_dirs(_runtime_dirs RSys Star3D StarSP)
+rcmake_append_runtime_dirs(_runtime_dirs RSys Star3D StarSP StarEnc StarEnc2D)
################################################################################
# Configure and define targets
@@ -70,6 +72,7 @@ set(SDIS_FILES_INC
sdis_interface_c.h
sdis_medium_c.h
sdis_scene_c.h
+ sdis_scene_Xd.h
sdis_solve_Xd.h)
set(SDIS_FILES_DOC COPYING README.md)
@@ -97,7 +100,7 @@ add_library(sdis SHARED
${SDIS_FILES_SRC}
${SDIS_FILES_INC}
${SDIS_FILES_INC_API})
-target_link_libraries(sdis RSys Star2D Star3D StarSP ${MATH_LIB})
+target_link_libraries(sdis RSys Star2D Star3D StarSP StarEnc StarEnc2D ${MATH_LIB})
set_target_properties(sdis PROPERTIES
DEFINE_SYMBOL SDIS_SHARED_BUILD
@@ -137,6 +140,8 @@ if(NOT NO_TEST)
new_test(test_sdis_camera)
new_test(test_sdis_conducto_radiative)
new_test(test_sdis_conducto_radiative_2d)
+ new_test(test_sdis_convection)
+ new_test(test_sdis_convection_non_uniform)
new_test(test_sdis_data)
new_test(test_sdis_device)
new_test(test_sdis_flux)
diff --git a/src/sdis.h b/src/sdis.h
@@ -188,12 +188,15 @@ struct sdis_interface_shader {
/* May be NULL for solid/solid or if the convection coefficient is 0 onto
* the whole interface. */
sdis_interface_getter_T convection_coef; /* In W.K^-1.m^-2 */
+ /* Under no circumstance can convection_coef() return outside of
+ * [0 convection_coef_upper_bound] */
+ double convection_coef_upper_bound;
struct sdis_interface_side_shader front;
struct sdis_interface_side_shader back;
};
#define SDIS_INTERFACE_SHADER_NULL__ \
- {NULL, SDIS_INTERFACE_SIDE_SHADER_NULL__, SDIS_INTERFACE_SIDE_SHADER_NULL__}
+ {NULL, 0, SDIS_INTERFACE_SIDE_SHADER_NULL__, SDIS_INTERFACE_SIDE_SHADER_NULL__}
static const struct sdis_interface_shader SDIS_INTERFACE_SHADER_NULL =
SDIS_INTERFACE_SHADER_NULL__;
diff --git a/src/sdis_device.c b/src/sdis_device.c
@@ -50,8 +50,10 @@ device_release(ref_T* 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);
+ ASSERT(flist_name_is_empty(&dev->interfaces_names));
+ ASSERT(flist_name_is_empty(&dev->media_names));
+ flist_name_release(&dev->interfaces_names);
+ flist_name_release(&dev->media_names);
darray_tile_release(&dev->tiles);
MEM_RM(dev->allocator, dev);
}
@@ -94,7 +96,8 @@ sdis_device_create
dev->verbose = verbose;
dev->nthreads = MMIN(nthreads_hint, (unsigned)omp_get_num_procs());
ref_init(&dev->ref);
- flist_name_init(allocator, &dev->names);
+ flist_name_init(allocator, &dev->interfaces_names);
+ flist_name_init(allocator, &dev->media_names);
darray_tile_init(allocator, &dev->tiles);
res = darray_tile_resize(&dev->tiles, dev->nthreads);
diff --git a/src/sdis_device_c.h b/src/sdis_device_c.h
@@ -16,6 +16,8 @@
#ifndef SDIS_DEVICE_C_H
#define SDIS_DEVICE_C_H
+#include "sdis.h"
+
#include <rsys/dynamic_array.h>
#include <rsys/free_list.h>
#include <rsys/ref_count.h>
@@ -42,7 +44,8 @@ struct sdis_device {
unsigned nthreads;
int verbose;
- struct flist_name names;
+ struct flist_name interfaces_names;
+ struct flist_name media_names;
struct darray_tile tiles;
struct s2d_device* s2d;
diff --git a/src/sdis_interface.c b/src/sdis_interface.c
@@ -54,6 +54,13 @@ check_interface_shader
"function of the interface shader should be NULL.\n", caller_name);
}
+ if(shader->convection_coef_upper_bound < 0) {
+ log_warn(dev,
+ "%s: Invalid upper bound for convection coefficient (%g).\n",
+ caller_name, shader->convection_coef_upper_bound);
+ if(type[0] == SDIS_FLUID || type[1] == SDIS_FLUID) return 0;
+ }
+
FOR_EACH(i, 0, 2) {
switch(type[i]) {
case SDIS_SOLID:
@@ -89,7 +96,7 @@ interface_release(ref_T* ref)
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);
+ flist_name_del(&dev->interfaces_names, interf->id);
MEM_RM(dev->allocator, interf);
SDIS(device_ref_put(dev));
}
@@ -141,7 +148,7 @@ sdis_interface_create
interf->medium_back = back;
interf->dev = dev;
interf->shader = *shader;
- interf->id = flist_name_add(&dev->names);
+ interf->id = flist_name_add(&dev->interfaces_names);
if(data) {
SDIS(data_ref_get(data));
diff --git a/src/sdis_interface_c.h b/src/sdis_interface_c.h
@@ -70,6 +70,14 @@ interface_get_convection_coef
}
static INLINE double
+interface_get_convection_coef_upper_bound
+ (const struct sdis_interface* interf)
+{
+ ASSERT(interf);
+ return interf->shader.convection_coef_upper_bound;
+}
+
+static INLINE double
interface_side_get_temperature
(const struct sdis_interface* interf,
const struct sdis_interface_fragment* frag)
diff --git a/src/sdis_medium.c b/src/sdis_medium.c
@@ -64,6 +64,7 @@ medium_create
SDIS(device_ref_get(dev));
medium->dev = dev;
medium->type = type;
+ medium->id = flist_name_add(&dev->media_names);
exit:
if(out_medium) *out_medium = medium;
@@ -85,6 +86,7 @@ medium_release(ref_T* ref)
medium = CONTAINER_OF(ref, struct sdis_medium, ref);
dev = medium->dev;
if(medium->data) SDIS(data_ref_put(medium->data));
+ flist_name_del(&dev->media_names, medium->id);
MEM_RM(dev->allocator, medium);
SDIS(device_ref_put(dev));
}
diff --git a/src/sdis_medium_c.h b/src/sdis_medium_c.h
@@ -26,11 +26,19 @@ struct sdis_medium {
} shader;
struct sdis_data* data;
+ struct fid id; /* Unique identifier of the medium */
ref_T ref;
struct sdis_device* dev;
};
+static FINLINE unsigned
+medium_get_id(const struct sdis_medium* mdm)
+{
+ ASSERT(mdm);
+ return mdm->id.index;
+}
+
/*******************************************************************************
* Fluid local functions
******************************************************************************/
diff --git a/src/sdis_scene.c b/src/sdis_scene.c
@@ -13,26 +13,21 @@
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>. */
+#include "sdis_scene_Xd.h"
+
+/* Generate the 2D functions of the scene */
+#define SDIS_SCENE_DIMENSION 2
+#include "sdis_scene_Xd.h"
+
+/* Generate the 3D functions of the scene */
+#define SDIS_SCENE_DIMENSION 3
+#include "sdis_scene_Xd.h"
+
#include "sdis.h"
-#include "sdis_device_c.h"
-#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 <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
******************************************************************************/
@@ -95,532 +90,6 @@ project_position
}
}
-/* 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_3d(const unsigned itri, unsigned out_ids[3], void* data)
-{
- struct geometry_context* ctx = data;
- size_t ids[3];
- ASSERT(ctx);
- ctx->indices(itri, ids, ctx->data);
- out_ids[0] = (unsigned)ids[0];
- out_ids[1] = (unsigned)ids[1];
- out_ids[2] = (unsigned)ids[2];
-}
-
-static void
-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];
- ASSERT(ctx);
- ctx->position(ivert, pos, ctx->data);
- out_pos[0] = (float)pos[0];
- out_pos[1] = (float)pos[1];
- out_pos[2] = (float)pos[2];
-}
-
-static void
-clear_interfaces(struct sdis_scene* scn)
-{
- size_t i;
- ASSERT(scn);
- 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_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 (*interf)(const size_t itri, struct sdis_interface**, void*),
- void* ctx)
-{
- size_t itri;
- res_T res = RES_OK;
- ASSERT(ntris && interf);
-
- clear_interfaces(scn);
-
- FOR_EACH(itri, 0, ntris) {
- struct sdis_interface* itface;
- size_t ninterfaces;
- unsigned id;
-
- /* Retrieve the interface of the primitive */
- interf(itri, &itface, ctx);
- id = interface_get_id(itface);
-
- /* Check that the interface is already registered against the scene */
- ninterfaces = darray_interf_size_get(&scn->interfaces);
- if(id >= ninterfaces) {
- res = darray_interf_resize(&scn->interfaces, id + 1);
- if(res != RES_OK) goto error;
- }
- if(darray_interf_cdata_get(&scn->interfaces)[id]) {
- ASSERT(darray_interf_cdata_get(&scn->interfaces)[id] == itface);
- } else {
- SDIS(interface_ref_get(itface));
- darray_interf_data_get(&scn->interfaces)[id] = itface;
- }
-
- /* Register the primitive interface */
- res = darray_interf_push_back(&scn->prim_interfaces, &itface);
- if(res != RES_OK) goto error;
- }
-
-exit:
- return res;
-error:
- clear_interfaces(scn);
- goto exit;
-}
-
-static res_T
-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*),
- const size_t nverts, /* #vertices */
- void (*position)(const size_t ivert, double pos[3], void* ctx),
- void* ctx)
-{
- struct geometry_context context;
- struct s3d_shape* s3d_msh = NULL;
- struct s3d_scene* s3d_scn = NULL;
- struct s3d_vertex_data vdata = S3D_VERTEX_DATA_NULL;
- res_T res = RES_OK;
- ASSERT(scn && ntris && indices && nverts && position);
-
- /* Setup the intermediary geometry context */
- context.indices = indices;
- context.position = position;
- context.data = ctx;
-
- /* Setup the vertex data */
- vdata.usage = S3D_POSITION;
- vdata.type = S3D_FLOAT3;
- 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_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_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);
- if(res != RES_OK) goto error;
-
-exit:
- if(s3d_msh) S3D(shape_ref_put(s3d_msh));
- if(s3d_scn) S3D(scene_ref_put(s3d_scn));
- return res;
-error:
- if(scn->s3d_view) S3D(scene_view_ref_put(scn->s3d_view));
- goto exit;
-}
-
-static res_T
-scene_create
- (struct sdis_device* dev,
- 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[], void* ctx),
- void* ctx,
- struct sdis_scene** out_scn)
-{
- struct sdis_scene* scn = NULL;
- res_T res = RES_OK;
-
- if(!dev || !out_scn || !nprims || !indices || !interf || !nverts
- || !position || nprims > UINT_MAX || nverts > UINT_MAX) {
- res = RES_BAD_ARG;
- goto error;
- }
-
- scn = MEM_CALLOC(dev->allocator, 1, sizeof(struct sdis_scene));
- if(!scn) {
- log_err(dev, "%s: could not allocate the Stardis scene.\n", FUNC_NAME);
- res = RES_MEM_ERR;
- goto error;
- }
- ref_init(&scn->ref);
- SDIS(device_ref_get(dev));
- scn->dev = dev;
- scn->ambient_radiative_temperature = -1;
- darray_interf_init(dev->allocator, &scn->interfaces);
- darray_interf_init(dev->allocator, &scn->prim_interfaces);
-
- 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;
- }
-
- 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;
- }
-
-exit:
- if(out_scn) *out_scn = scn;
- return res;
-error:
- if(scn) {
- SDIS(scene_ref_put(scn));
- scn = NULL;
- }
- goto exit;
-}
-
-static INLINE res_T
-scene_get_medium_2d
- (const struct sdis_scene* scn,
- const double pos[2],
- struct get_medium_info* info, /* May be NULL */
- const struct sdis_medium** out_medium)
-{
- const struct sdis_medium* medium = NULL;
- size_t iprim, nprims;
- size_t nfailures = 0;
- const size_t max_failures = 10;
- /* Range of the parametric coordinate into which positions are challenged */
- const float s_range[2] = {0.25, 0.75};
- const size_t s_nsteps = 3; /* #challenges per primitive into the range */
- float s_step;
- res_T res = RES_OK;
- ASSERT(scn && pos);
-
- s_step = (s_range[1] - s_range[0]) / (float)(s_nsteps-1);
-
- 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;
- const float range[2] = {0.f, FLT_MAX};
- float N[2], P[2], dir[2], cos_N_dir;
- float s = s_range[0];
-
- do {
- /* 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));
-
- /* Trace a ray from the random 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));
-
- /* 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;
- }
- }
- /* Discard the hit if it is on a vertex, i.e. between 2 segments, and
- * target a new position onto the current primitive */
- } while((S2D_HIT_NONE(&hit) || hit_on_vertex(&hit))
- && (s+=s_step) <= s_range[1]);
-
- /* The hits of all targeted positions on the current primitive are on
- * vertices. Challenge positions on another primitive. */
- if(s > s_range[1]) continue;
-
- f2_normalize(N, hit.normal);
- cos_N_dir = f2_dot(N, dir);
-
- if(absf(cos_N_dir) > 1.e-1f) { /* Not roughly orthognonal */
- const struct sdis_interface* interf;
- interf = scene_get_interface(scn, hit.prim.prim_id);
- medium = interface_get_medium
- (interf, cos_N_dir < 0 ? SDIS_FRONT : SDIS_BACK);
-
- /* Register the get_medium_info */
- if(info) {
- f2_set(info->pos_tgt, attr.value);
- f2_set(info->ray_org, P);
- f2_set(info->ray_dir, dir);
- info->hit_2d = hit;
- }
- 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;
-}
-
-static INLINE res_T
-scene_get_medium_3d
- (const struct sdis_scene* scn,
- const double pos[3],
- struct get_medium_info* info,
- const struct sdis_medium** out_medium)
-{
- const struct sdis_medium* medium = NULL;
- size_t iprim, nprims;
- size_t nfailures = 0;
- const size_t max_failures = 10;
- float st[3][2]; /* Position to challenge onto the primitive */
- const size_t nsteps = sizeof(st)/sizeof(float[2]);
- res_T res = RES_OK;
- ASSERT(scn && pos);
-
- /* Setup the position to challenge onto the primitives */
- f2(st[0], 1.f/6.f, 5.f/12.f);
- f2(st[1], 5.f/12.f, 1.f/6.f);
- f2(st[2], 5.f/12.f, 5.f/12.f);
-
- S3D(scene_view_primitives_count(scn->s3d_view, &nprims));
- FOR_EACH(iprim, 0, nprims) {
- struct s3d_hit hit;
- struct s3d_attrib attr;
- struct s3d_primitive prim;
- const float range[2] = {0.f, FLT_MAX};
- float N[3], P[3], dir[3], cos_N_dir;
- size_t istep = 0;
-
- do {
- /* Retrieve a position onto the primitive */
- S3D(scene_view_get_primitive(scn->s3d_view, (unsigned)iprim, &prim));
- S3D(primitive_get_attrib(&prim, S3D_POSITION, st[istep], &attr));
-
- /* Trace a ray from the randomw walk vertex toward the retrieved primitive
- * position */
- f3_normalize(dir, f3_sub(dir, attr.value, f3_set_d3(P, pos)));
- S3D(scene_view_trace_ray(scn->s3d_view, P, dir, range, NULL, &hit));
-
- /* Unforeseen error. One has to intersect a primitive ! */
- if(S3D_HIT_NONE(&hit)) {
- ++nfailures;
- if(nfailures < max_failures) {
- continue;
- } else {
- res = RES_BAD_ARG;
- goto error;
- }
- }
- /* Discard the hit if it is on an edge, i.e. between 2 triangles, and
- * target a new position onto the current primitive */
- } while((S3D_HIT_NONE(&hit) || hit_on_edge(&hit)) && ++istep < nsteps);
-
- /* The hits of all targeted positions on the current primitive are on
- * edges. Challenge positions on another primitive. */
- if(istep >= nsteps) continue;
-
- f3_normalize(N, hit.normal);
- cos_N_dir = f3_dot(N, dir);
-
- 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);
-
- if(info) {
- f3_set(info->pos_tgt, attr.value);
- f3_set(info->ray_org, P);
- f3_set(info->ray_dir, dir);
- info->hit_3d = hit;
- }
- break;
- }
- }
-
-exit:
- *out_medium = medium;
- return res;
-error:
- log_err(scn->dev, "%s: could not retrieve the medium at {%g, %g, %g}.\n",
- FUNC_NAME, SPLIT3(pos));
- goto exit;
-}
-
static void
scene_release(ref_T * ref)
{
@@ -629,9 +98,12 @@ scene_release(ref_T * ref)
ASSERT(ref);
scn = CONTAINER_OF(ref, struct sdis_scene, ref);
dev = scn->dev;
- clear_interfaces(scn);
+ clear_properties(scn);
darray_interf_release(&scn->interfaces);
- darray_interf_release(&scn->prim_interfaces);
+ darray_medium_release(&scn->media);
+ darray_prim_prop_release(&scn->prim_props);
+ htable_enclosure_release(&scn->enclosures);
+ htable_d_release(&scn->tmp_hc_ub);
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);
@@ -652,8 +124,8 @@ sdis_scene_create
void* ctx,
struct sdis_scene** out_scn)
{
- return scene_create
- (dev, 0/*is_2D*/, ntris, indices, interf, nverts, position, ctx, out_scn);
+ return scene_create_3d
+ (dev, ntris, indices, interf, nverts, position, ctx, out_scn);
}
res_T
@@ -667,8 +139,8 @@ sdis_scene_2d_create
void* ctx,
struct sdis_scene** out_scn)
{
- return scene_create
- (dev, 1/*is_2D*/, nsegs, indices, interf, nverts, position, ctx, out_scn);
+ return scene_create_2d
+ (dev, nsegs, indices, interf, nverts, position, ctx, out_scn);
}
res_T
@@ -809,8 +281,8 @@ sdis_scene_boundary_project_position
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];
+ ASSERT(scn && iprim < darray_prim_prop_size_get(&scn->prim_props));
+ return darray_prim_prop_cdata_get(&scn->prim_props)[iprim].interf;
}
res_T
diff --git a/src/sdis_scene_Xd.h b/src/sdis_scene_Xd.h
@@ -0,0 +1,878 @@
+/* Copyright (C) 2016-2018 |Meso|Star> (contact@meso-star.com)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. */
+
+#ifndef SDIS_SCENE_DIMENSION
+
+#ifndef SDIS_SCENE_XD_H
+#define SDIS_SCENE_XD_H
+
+#include "sdis_interface_c.h"
+#include "sdis_medium_c.h"
+#include "sdis_scene_c.h"
+
+#include <rsys/rsys.h>
+
+/*******************************************************************************
+ * Define the helper functions and the data types used by the scene
+ * independently of its dimension, i.e. 2D or 3D.
+ ******************************************************************************/
+
+/* Context used to wrap the user geometry and interfaces to Star-Enc */
+struct geometry {
+ void (*indices)(const size_t iprim, size_t ids[], void*);
+ void (*interf)(const size_t iprim, struct sdis_interface**, void*);
+ void (*position)(const size_t ivert, double pos[], void*);
+ void* data;
+};
+
+/* Fetch the media split by the primitive `iprim'. This first and second media
+ * are the media from the front face side and back face side of the primitive,
+ * respectively. */
+static void
+geometry_media(const unsigned iprim, unsigned media[2], void* data)
+{
+ struct geometry* ctx = data;
+ struct sdis_interface* interf;
+ ASSERT(ctx && media);
+ ctx->interf(iprim, &interf, ctx->data);
+ media[0] = medium_get_id(interf->medium_front);
+ media[1] = medium_get_id(interf->medium_back);
+}
+
+/* Register the submitted medium against the scene if it is not already
+ * registered. On registration, no reference is taken onto the medium; the
+ * scene references its media through its interfaces and it is thus useless to
+ * take another reference onto them. */
+static res_T
+register_medium(struct sdis_scene* scn, struct sdis_medium* mdm)
+{
+ unsigned id;
+ size_t nmedia;
+ res_T res = RES_OK;
+ ASSERT(scn && mdm);
+
+ /* Check that the front medium is already registered against the scene */
+ id = medium_get_id(mdm);
+ nmedia = darray_medium_size_get(&scn->media);
+ if(id >= nmedia) {
+ res = darray_medium_resize(&scn->media, id + 1);
+ if(res != RES_OK) return res;
+ }
+ if(darray_medium_cdata_get(&scn->media)[id]) {
+ ASSERT(darray_medium_cdata_get(&scn->media)[id] == mdm);
+ } else {
+ /* Do not take a reference onto the medium since we already take a
+ * reference onto at least one interface that uses it, and thus that has a
+ * reference onto it */
+ darray_medium_data_get(&scn->media)[id] = mdm;
+ }
+ return RES_OK;
+}
+
+/* Release the reference toward the interfaces and thus clear the list of scene
+ * interfaces, the list of scene media, and the list of per-primitive
+ * interface. */
+static void
+clear_properties(struct sdis_scene* scn)
+{
+ size_t i;
+ ASSERT(scn);
+ 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_interf_clear(&scn->interfaces);
+ darray_medium_clear(&scn->media);
+ darray_prim_prop_clear(&scn->prim_props);
+}
+
+#endif /* SDIS_SCENE_XD_H */
+#else /* !SDIS_SCENE_DIMENSION */
+
+#include "sdis_device_c.h"
+
+#include <rsys/float2.h>
+#include <rsys/float3.h>
+#include <rsys/double2.h>
+#include <rsys/double3.h>
+#include <rsys/mem_allocator.h>
+
+#include <limits.h>
+
+/* Check the submitted dimension and include its specific headers */
+#if (SDIS_SCENE_DIMENSION == 2)
+ #include <star/senc2d.h>
+ #include <star/s2d.h>
+#elif (SDIS_SCENE_DIMENSION == 3)
+ #include <star/senc.h>
+ #include <star/s3d.h>
+#else
+ #error "Invalid SDIS_SCENE_DIMENSION value."
+#endif
+
+/* Syntactic sugar */
+#define DIM SDIS_SCENE_DIMENSION
+
+/* Star-Enc macros generic to the SDIS_SCENE_DIMENSION */
+#if DIM == 2
+ #define sencXd(Name) CONCAT(senc2d_, Name)
+ #define SENCXD SENC2D
+#else
+ #define sencXd(Name) CONCAT(senc_, Name)
+ #define SENCXD SENC
+#endif
+
+/* Star-XD macros generic to SDIS_SCENE_DIMENSION */
+#define sXd(Name) CONCAT(CONCAT(CONCAT(s, DIM), d_), Name)
+#define SXD CONCAT(CONCAT(S, DIM), D)
+#define SXD_VERTEX_DATA_NULL CONCAT(CONCAT(S,DIM),D_VERTEX_DATA_NULL)
+#define SXD_POSITION CONCAT(CONCAT(S, DIM), D_POSITION)
+#define SXD_TRACE CONCAT(CONCAT(S,DIM), D_TRACE)
+#define SXD_SAMPLE CONCAT(CONCAT(S,DIM), D_SAMPLE)
+#define SXD_GET_PRIMITIVE CONCAT(CONCAT(S,DIM), D_GET_PRIMITIVE)
+#define SXD_HIT_NONE CONCAT(CONCAT(S,DIM), D_HIT_NONE)
+#define SXD_PRIMITIVE_EQ CONCAT(CONCAT(S,DIM), D_PRIMITIVE_EQ)
+
+/* Vector macros generic to SDIS_SCENE_DIMENSION */
+#define fX(Func) CONCAT(CONCAT(CONCAT(f, DIM), _), Func)
+#define fX_set_dX CONCAT(CONCAT(CONCAT(f, DIM), _set_d), DIM)
+
+/* Macro making generic its subimitted name to SDIS_SCENE_DIMENSION */
+#define XD(Name) CONCAT(CONCAT(CONCAT(Name, _), DIM), d)
+
+#if DIM == 2
+ #define HIT_ON_BOUNDARY hit_on_vertex
+#else
+ #define HIT_ON_BOUNDARY hit_on_edge
+#endif
+
+/*******************************************************************************
+ * Helper functions
+ ******************************************************************************/
+#if DIM == 2
+/* 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);
+}
+
+#else /* DIM == 3 */
+/* 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);
+}
+#endif /* DIM == 2 */
+
+/* Avoid self-intersection for a ray starting from a planar primitive, i.e. a
+ * triangle or a line segment */
+static int
+XD(hit_filter_function)
+ (const struct sXd(hit)* hit,
+ const float org[DIM],
+ const float dir[DIM],
+ void* ray_data,
+ void* filter_data)
+{
+ const struct sXd(hit)* hit_from = ray_data;
+ (void)org, (void)dir, (void)filter_data;
+
+ if(!hit_from || SXD_HIT_NONE(hit_from)) return 0; /* No filtering */
+
+ if(SXD_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/vertex shared by the 2 primitives. */
+#if DIM == 2
+ return hit_on_vertex(hit_from) && hit_on_vertex(hit);
+#else
+ return hit_on_edge(hit_from) && hit_on_edge(hit);
+#endif
+ }
+ return 0;
+}
+
+/* Retrieve the indices of `struct geometry' primitive */
+static void
+XD(geometry_indices)(const unsigned iprim, unsigned out_ids[DIM], void* data)
+{
+ struct geometry* ctx = data;
+ size_t ids[DIM];
+ int i;
+ ASSERT(ctx && out_ids);
+ ctx->indices(iprim, ids, ctx->data);
+ FOR_EACH(i, 0, DIM) out_ids[i] = (unsigned)ids[i];
+}
+
+/* Retrieve the coordinates of `struct geometry' vertex */
+static void
+XD(geometry_position)(const unsigned ivert, double out_pos[DIM], void* data)
+{
+ struct geometry* ctx = data;
+ double pos[DIM];
+ int i;
+ ASSERT(ctx && out_pos);
+ ctx->position(ivert, pos, ctx->data);
+ FOR_EACH(i, 0, DIM) out_pos[i] = pos[i];
+}
+
+/* Retrieve the indices of a primitive of a Star-EncXD descriptor */
+static void
+XD(descriptor_indices)(const unsigned iprim, unsigned ids[DIM], void* data)
+{
+ struct sencXd(descriptor)* desc = data;
+#if DIM == 2
+ SENCXD(descriptor_get_global_segment(desc, iprim, ids));
+#else
+ SENCXD(descriptor_get_global_triangle(desc, iprim, ids));
+#endif
+}
+
+/* Retrieve the coordinates of a vertex of a Star-EncXD descriptor */
+static void
+XD(descriptor_position)(const unsigned ivert, float out_pos[DIM], void* data)
+{
+ struct sencXd(descriptor)* desc = data;
+ double pos[3];
+ int i;
+ SENCXD(descriptor_get_global_vertex(desc, ivert, pos));
+ FOR_EACH(i, 0, DIM) out_pos[i] = (float)pos[i];
+}
+
+/* Retrieve the indices of a primitive of a Star-EncXD enclosure */
+static void
+XD(enclosure_indices)(const unsigned iprim, unsigned ids[DIM], void* data)
+{
+ struct sencXd(enclosure)* enc = data;
+#if DIM == 2
+ SENCXD(enclosure_get_segment(enc, iprim, ids));
+#else
+ SENCXD(enclosure_get_triangle(enc, iprim, ids));
+#endif
+}
+
+/* Retrieve the coordinates of a vertex of a Star-EncXD encolsure */
+static void
+XD(enclosure_position)(const unsigned ivert, float out_pos[DIM], void* data)
+{
+ struct sencXd(enclosure)* enc = data;
+ double pos[DIM];
+ int i;
+ ASSERT(out_pos);
+ SENCXD(enclosure_get_vertex(enc, ivert, pos));
+ FOR_EACH(i, 0, DIM) out_pos[i] = (float)pos[i];
+}
+
+/* Use Star-EncXD to analyze the user defined data. It essentially cleans-up
+ * the geometry and extracts the enclosures wrt to the submitted media. Note
+ * that data inconsistencies are also detected by this function. In this case
+ * an error is returned. */
+static res_T
+XD(run_analyze)
+ (struct sdis_scene* scn,
+ const size_t nprims, /* #primitives */
+ void (*indices)(const size_t iprim, size_t ids[], void*),
+ void (interf)(const size_t iprim, struct sdis_interface**, void*),
+ const size_t nverts, /* #vertices */
+ void (*position)(const size_t ivert, double pos[], void*),
+ void* ctx,
+ struct sencXd(descriptor)** out_desc)
+{
+ struct geometry geom;
+ struct sencXd(device)* senc = NULL;
+ struct sencXd(scene)* senc_scn = NULL;
+ struct sencXd(descriptor)* desc = NULL;
+ res_T res = RES_OK;
+ ASSERT(scn && nprims && indices && interf && nverts && position && out_desc);
+
+ res = sencXd(device_create)(scn->dev->logger, scn->dev->allocator,
+ scn->dev->nthreads, scn->dev->verbose, &senc);
+ if(res != RES_OK) goto error;
+
+ res = sencXd(scene_create)(senc,
+#if DIM == 2
+ SENC2D_CONVENTION_NORMAL_BACK | SENC2D_CONVENTION_NORMAL_OUTSIDE,
+#else
+ SENC_CONVENTION_NORMAL_BACK | SENC_CONVENTION_NORMAL_OUTSIDE,
+#endif
+ &senc_scn);
+ if(res != RES_OK) goto error;
+
+ /* Setup the geometry data */
+ geom.indices = indices;
+ geom.interf = interf;
+ geom.position = position;
+ geom.data = ctx;
+ res = sencXd(scene_add_geometry)
+ (senc_scn, (unsigned)nprims, XD(geometry_indices), geometry_media, NULL,
+ (unsigned)nverts, XD(geometry_position), &geom);
+ if(res != RES_OK) goto error;
+
+ /* Launch the scene analyze */
+ res = sencXd(scene_analyze)(senc_scn, &desc);
+ if(res != RES_OK) goto error;
+
+exit:
+ if(senc) SENCXD(device_ref_put(senc));
+ if(senc_scn) SENCXD(scene_ref_put(senc_scn));
+ if(out_desc) *out_desc = desc;
+ return res;
+error:
+ if(desc) {
+ SENCXD(descriptor_ref_put(desc));
+ desc = NULL;
+ }
+ goto exit;
+
+}
+
+/* Register the media and the interfaces, map each primitive to its interface
+ * and associated to each primitive the identifier of the enclosures that it
+ * splits. */
+static res_T
+XD(setup_properties)
+ (struct sdis_scene* scn,
+ struct sencXd(descriptor)* desc,
+ void (*interf)(const size_t itri, struct sdis_interface**, void*),
+ void* ctx)
+{
+ unsigned iprim, nprims;
+ res_T res = RES_OK;
+ ASSERT(scn && interf);
+
+ clear_properties(scn);
+
+#if DIM == 2
+ SENCXD(descriptor_get_global_segments_count(desc, &nprims));
+#else
+ SENCXD(descriptor_get_global_triangles_count(desc, &nprims));
+#endif
+ FOR_EACH(iprim, 0, nprims) {
+ struct prim_prop* prim_prop;
+ struct sdis_interface* itface;
+ unsigned enclosures[2];
+ unsigned iprim_adjusted; /* Primitive id in user space */
+ unsigned id;
+ int i;
+ double* enc_upper_bound;
+ size_t ninterfaces;
+
+#if DIM == 2
+ /* Retrieve the triangle id in user space */
+ SENCXD(descriptor_get_global_segment_global_id(desc, iprim, &iprim_adjusted));
+ /* Fetch the enclosures that the segment splits */
+ SENCXD(descriptor_get_global_segment_enclosures(desc, iprim, enclosures));
+#else
+ /* Retrieve the triangle id in user space */
+ SENCXD(descriptor_get_global_triangle_global_id(desc, iprim, &iprim_adjusted));
+ /* Fetch the enclosures that the triangle splits */
+ SENCXD(descriptor_get_global_triangle_enclosures(desc, iprim, enclosures));
+#endif
+
+ /* Fetch the interface of the primitive */
+ interf(iprim, &itface, ctx);
+
+ /* Check that the interface is already registered against the scene */
+ id = interface_get_id(itface);
+ ninterfaces = darray_interf_size_get(&scn->interfaces);
+ if(id >= ninterfaces) {
+ res = darray_interf_resize(&scn->interfaces, id + 1);
+ if(res != RES_OK) goto error;
+ }
+ if(darray_interf_cdata_get(&scn->interfaces)[id]) {
+ ASSERT(darray_interf_cdata_get(&scn->interfaces)[id] == itface);
+ } else {
+ SDIS(interface_ref_get(itface));
+ darray_interf_data_get(&scn->interfaces)[id] = itface;
+ }
+
+ /* Register the interface media against the scene */
+ res = register_medium(scn, itface->medium_front);
+ if(res != RES_OK) goto error;
+ res = register_medium(scn, itface->medium_back);
+ if(res != RES_OK) goto error;
+
+ /* Allocate primitive properties */
+ res = darray_prim_prop_resize(&scn->prim_props, iprim+1);
+ if(res != RES_OK) goto error;
+
+ /* Setup primitive properties */
+ prim_prop = darray_prim_prop_data_get(&scn->prim_props) + iprim;
+ prim_prop->interf = itface;
+ prim_prop->front_enclosure = enclosures[0];
+ prim_prop->back_enclosure = enclosures[1];
+
+ /* Build per-interface hc upper bounds in a tmp table */
+ FOR_EACH(i, 0, 2) {
+ double hc_ub = interface_get_convection_coef_upper_bound(itface);
+ enc_upper_bound = htable_d_find(&scn->tmp_hc_ub, enclosures+i);
+ if(!enc_upper_bound) {
+ res = htable_d_set(&scn->tmp_hc_ub, enclosures+i, &hc_ub);
+ } else {
+ *enc_upper_bound = MMAX(*enc_upper_bound, hc_ub);
+ }
+ }
+ }
+
+exit:
+ return res;
+error:
+ clear_properties(scn);
+ goto exit;
+}
+
+/* Build the Star-XD scene view of the whole scene */
+static res_T
+XD(setup_scene_geometry)(struct sdis_scene* scn, struct sencXd(descriptor)* desc)
+{
+ struct sXd(device)* sXd_dev = NULL;
+ struct sXd(shape)* sXd_shape = NULL;
+ struct sXd(scene)* sXd_scn = NULL;
+ struct sXd(vertex_data) vdata = SXD_VERTEX_DATA_NULL;
+ unsigned nprims, nverts;
+ res_T res = RES_OK;
+ ASSERT(scn && desc);
+
+ SENCXD(descriptor_get_global_vertices_count(desc, &nverts));
+
+ /* Setup the vertex data */
+ vdata.usage = SXD_POSITION;
+#if DIM == 2
+ vdata.type = S2D_FLOAT2;
+#else
+ vdata.type = S3D_FLOAT3;
+#endif
+ vdata.get = XD(descriptor_position);
+
+ /* Create the Star-XD geometry of the whole scene */
+ #define CALL(Func) { if(RES_OK != (res = Func)) goto error; } (void)0
+#if DIM == 2
+ sXd_dev = scn->dev->s2d;
+ SENCXD(descriptor_get_global_segments_count(desc, &nprims));
+ CALL(sXd(shape_create_line_segments)(sXd_dev, &sXd_shape));
+ CALL(sXd(line_segments_set_hit_filter_function)(sXd_shape,
+ XD(hit_filter_function), NULL));
+ CALL(sXd(line_segments_setup_indexed_vertices)(sXd_shape, nprims,
+ XD(descriptor_indices), nverts, &vdata, 1, desc));
+#else
+ sXd_dev = scn->dev->s3d;
+ SENCXD(descriptor_get_global_triangles_count(desc, &nprims));
+ CALL(sXd(shape_create_mesh)(sXd_dev, &sXd_shape));
+ CALL(sXd(mesh_set_hit_filter_function)(sXd_shape, XD(hit_filter_function), NULL));
+ CALL(sXd(mesh_setup_indexed_vertices)(sXd_shape, nprims, XD(descriptor_indices),
+ nverts, &vdata, 1, desc));
+#endif
+ CALL(sXd(scene_create)(sXd_dev, &sXd_scn));
+ CALL(sXd(scene_attach_shape)(sXd_scn, sXd_shape));
+ CALL(sXd(scene_view_create)(sXd_scn, SXD_TRACE|SXD_GET_PRIMITIVE,
+ &scn->sXd(view)));
+ #undef CALL
+
+exit:
+ if(sXd_shape) SXD(shape_ref_put(sXd_shape));
+ if(sXd_scn) SXD(scene_ref_put(sXd_scn));
+ return res;
+error:
+ if(scn->sXd(view)) SXD(scene_view_ref_put(scn->sXd(view)));
+ goto exit;
+}
+
+/* Build the Star-XD scene view of a specific enclosure and map their local
+ * primitive id to their primitive id in the whole scene */
+static res_T
+XD(setup_enclosure_geometry)(struct sdis_scene* scn, struct sencXd(enclosure)* enc)
+{
+ struct sXd(device)* sXd_dev = NULL;
+ struct sXd(scene)* sXd_scn = NULL;
+ struct sXd(shape)* sXd_shape = NULL;
+ struct sXd(vertex_data) vdata = SXD_VERTEX_DATA_NULL;
+ struct enclosure enc_dummy;
+ struct enclosure* enc_data;
+ float S, V;
+ double* p_ub;
+ unsigned iprim, nprims, nverts;
+#if DIM == 2
+ struct senc2d_enclosure_header header;
+#else
+ struct senc_enclosure_header header;
+#endif
+ res_T res = RES_OK;
+ ASSERT(scn && enc);
+
+ enclosure_init(scn->dev->allocator, &enc_dummy);
+
+ SENCXD(enclosure_get_header(enc, &header));
+#if DIM == 2
+ sXd_dev = scn->dev->s2d;
+ nprims = header.segment_count;
+#else
+ sXd_dev = scn->dev->s3d;
+ nprims = header.triangle_count;
+#endif
+ nverts = header.vertices_count;
+
+ /* Register the enclosure into the scene. Use a dummy data on their
+ * registration; in order to avoid a costly copy, we are going to setup the
+ * registered data rather than a local data that would be then registered. In
+ * other words, the following hash table registration can be seen as an
+ * allocation of the enclosure data to setup. */
+ res = htable_enclosure_set(&scn->enclosures, &header.enclosure_id, &enc_dummy);
+ if(res != RES_OK) goto error;
+
+ /* Fetch the data of the registered enclosure */
+ enc_data = htable_enclosure_find(&scn->enclosures, &header.enclosure_id);
+ ASSERT(enc_data != NULL);
+
+ /* Setup the vertex data */
+ vdata.usage = SXD_POSITION;
+#if DIM == 2
+ vdata.type = S2D_FLOAT2;
+#else
+ vdata.type = S2D_FLOAT3;
+#endif
+ vdata.get = XD(enclosure_position);
+
+ /* Create the Star-XD geometry */
+ #define CALL(Func) { if(RES_OK != (res = Func)) goto error; } (void)0
+#if DIM == 2
+ CALL(sXd(shape_create_line_segments)(sXd_dev, &sXd_shape));
+ CALL(sXd(line_segments_setup_indexed_vertices)(sXd_shape, nprims,
+ XD(enclosure_indices), nverts, &vdata, 1, enc));
+#else
+ CALL(sXd(shape_create_mesh)(sXd_dev, &sXd_shape));
+ CALL(sXd(mesh_setup_indexed_vertices)(sXd_shape, nprims, XD(enclosure_indices),
+ nverts, &vdata, 1, enc));
+#endif
+ CALL(sXd(scene_create)(sXd_dev, &sXd_scn));
+ CALL(sXd(scene_attach_shape)(sXd_scn, sXd_shape));
+ CALL(sXd(scene_view_create)(sXd_scn, SXD_SAMPLE, &enc_data->sXd(view)));
+
+ /* Compute the S/V ratio */
+#if DIM == 2
+ CALL(sXd(scene_view_compute_contour_length)(enc_data->sXd(view), &S));
+ CALL(sXd(scene_view_compute_area)(enc_data->sXd(view), &V));
+#else
+ CALL(sXd(scene_view_compute_area)(enc_data->sXd(view), &S));
+ CALL(sXd(scene_view_compute_volume)(enc_data->sXd(view), &V));
+#endif
+ /* The volume of the enclosure is actually negative since Star-Enc ensures
+ * that the normal of its primitives point outward the enclosure. Take its
+ * absolute value in order to ensure a postive value. */
+ enc_data->S_over_V = S/absf(V);
+ #undef CALL
+
+ /* Set enclosure hc upper bound regardless of its media being a fluid */
+ p_ub = htable_d_find(&scn->tmp_hc_ub, &header.enclosure_id);
+ ASSERT(p_ub);
+ enc_data->hc_upper_bound = *p_ub;
+
+ /* Define the identifier of the enclosure primitives in the whole scene */
+ res = darray_uint_resize(&enc_data->local2global, nprims);
+ if(res != RES_OK) goto error;
+ FOR_EACH(iprim, 0, nprims) {
+#if DIM == 2
+ SENCXD(enclosure_get_segment_global_id
+ (enc, iprim, darray_uint_data_get(&enc_data->local2global)+iprim));
+#else
+ SENCXD(enclosure_get_triangle_global_id
+ (enc, iprim, darray_uint_data_get(&enc_data->local2global)+iprim));
+#endif
+ }
+
+exit:
+ enclosure_release(&enc_dummy);
+ if(sXd_shape) SXD(shape_ref_put(sXd_shape));
+ if(sXd_scn) SXD(scene_ref_put(sXd_scn));
+ return res;
+error:
+ htable_enclosure_erase(&scn->enclosures, &header.enclosure_id);
+ goto exit;
+}
+
+/* Build the Star-XD scene view and define its associated data of the finite
+ * fluid enclosures */
+static res_T
+XD(setup_enclosures)(struct sdis_scene* scn, struct sencXd(descriptor)* desc)
+{
+ struct sencXd(enclosure)* enc = NULL;
+ unsigned ienc, nencs;
+ unsigned enclosed_medium;
+ res_T res = RES_OK;
+ ASSERT(scn && desc);
+
+ SENCXD(descriptor_get_enclosure_count(desc, &nencs));
+ FOR_EACH(ienc, 0, nencs) {
+#if DIM == 2
+ struct senc2d_enclosure_header header;
+#else
+ struct senc_enclosure_header header;
+#endif
+ const struct sdis_medium* mdm;
+
+ SENCXD(descriptor_get_enclosure(desc, ienc, &enc));
+ SENCXD(enclosure_get_header(enc, &header));
+
+ /* As paths don't go in infinite enclosures we can accept models are broken
+ * there. But nowhere else. */
+ if(header.enclosed_media_count != 1 && !header.is_infinite) {
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ SENCXD(enclosure_get_medium(enc, 0, &enclosed_medium));
+ if(res != RES_OK) goto error;
+ ASSERT(enclosed_medium < darray_medium_size_get(&scn->media));
+ mdm = darray_medium_cdata_get(&scn->media)[enclosed_medium];
+ ASSERT(mdm);
+
+ /* Silently discard the solid and infinite enclosures */
+ if(mdm->type == SDIS_FLUID && !header.is_infinite) {
+ res = XD(setup_enclosure_geometry)(scn, enc);
+ if(res != RES_OK) goto error;
+ }
+ SENCXD(enclosure_ref_put(enc));
+ enc = NULL;
+ }
+
+ /* tmp table no more useful */
+ htable_d_purge(&scn->tmp_hc_ub);
+exit:
+ if(enc) SENCXD(enclosure_ref_put(enc));
+ return res;
+error:
+ goto exit;
+}
+
+/* Create a Stardis scene */
+static res_T
+XD(scene_create)
+ (struct sdis_device* dev,
+ 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[], void* ctx),
+ void* ctx,
+ struct sdis_scene** out_scn)
+{
+ struct sencXd(descriptor)* desc = NULL;
+ struct sdis_scene* scn = NULL;
+ res_T res = RES_OK;
+
+ if(!dev || !out_scn || !nprims || !indices || !interf || !nverts
+ || !position || nprims > UINT_MAX || nverts > UINT_MAX) {
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ scn = MEM_CALLOC(dev->allocator, 1, sizeof(struct sdis_scene));
+ if(!scn) {
+ log_err(dev, "%s: could not allocate the Stardis scene.\n", FUNC_NAME);
+ res = RES_MEM_ERR;
+ goto error;
+ }
+ ref_init(&scn->ref);
+ SDIS(device_ref_get(dev));
+ scn->dev = dev;
+ scn->ambient_radiative_temperature = -1;
+ darray_interf_init(dev->allocator, &scn->interfaces);
+ darray_medium_init(dev->allocator, &scn->media);
+ darray_prim_prop_init(dev->allocator, &scn->prim_props);
+ htable_enclosure_init(dev->allocator, &scn->enclosures);
+ htable_d_init(dev->allocator, &scn->tmp_hc_ub);
+
+ res = XD(run_analyze)(scn, nprims, indices, interf, nverts, position, ctx, &desc);
+ if(res != RES_OK) {
+ log_err(dev, "%s: error during the scene analysis.\n", FUNC_NAME);
+ goto error;
+ }
+ res = XD(setup_properties)(scn, desc, interf, ctx);
+ if(res != RES_OK) {
+ log_err(dev, "%s: could not setup the scene interfaces and their media.\n",
+ FUNC_NAME);
+ goto error;
+ }
+ res = XD(setup_scene_geometry)(scn, desc);
+ if(res != RES_OK) {
+ log_err(dev, "%s: could not setup the scene geometry.\n", FUNC_NAME);
+ goto error;
+ }
+ res = XD(setup_enclosures)(scn, desc);
+ if(res != RES_OK) {
+ log_err(dev, "%s: could not setup the enclosures.\n", FUNC_NAME);
+ goto error;
+ }
+
+exit:
+ if(out_scn) *out_scn = scn;
+ if(desc) SENCXD(descriptor_ref_put(desc));
+ return res;
+error:
+ if(scn) {
+ SDIS(scene_ref_put(scn));
+ scn = NULL;
+ }
+ goto exit;
+}
+
+/*******************************************************************************
+ * Local functions
+ ******************************************************************************/
+static INLINE res_T
+XD(scene_get_medium)
+ (const struct sdis_scene* scn,
+ const double pos[2],
+ struct get_medium_info* info, /* May be NULL */
+ const struct sdis_medium** out_medium)
+{
+ const struct sdis_medium* medium = NULL;
+ size_t iprim, nprims;
+ size_t nfailures = 0;
+ const size_t max_failures = 10;
+ /* Range of the parametric coordinate into which positions are challenged */
+#if DIM == 2
+ float st[3];
+#else
+ float st[3][2];
+#endif
+ size_t nsteps = 3;
+ res_T res = RES_OK;
+ ASSERT(scn && pos);
+
+#if DIM == 2
+ st[0] = 0.25f;
+ st[1] = 0.50f;
+ st[2] = 0.75f;
+#else
+ f2(st[0], 1.f/6.f, 5.f/12.f);
+ f2(st[1], 5.f/12.f, 1.f/6.f);
+ f2(st[2], 5.f/12.f, 5.f/12.f);
+#endif
+
+ SXD(scene_view_primitives_count(scn->sXd(view), &nprims));
+ FOR_EACH(iprim, 0, nprims) {
+ struct sXd(hit) hit;
+ struct sXd(attrib) attr;
+ struct sXd(primitive) prim;
+ const float range[2] = {0.f, FLT_MAX};
+ float N[DIM], P[DIM], dir[DIM], cos_N_dir;
+ size_t istep = 0;
+
+ do {
+ /* Retrieve a position onto the primitive */
+ SXD(scene_view_get_primitive(scn->sXd(view), (unsigned)iprim, &prim));
+ SXD(primitive_get_attrib(&prim, SXD_POSITION, st[istep], &attr));
+
+ /* Trace a ray from the random walk vertex toward the retrieved primitive
+ * position */
+ fX(normalize)(dir, fX(sub)(dir, attr.value, fX_set_dX(P, pos)));
+ SXD(scene_view_trace_ray(scn->sXd(view), P, dir, range, NULL, &hit));
+
+ /* Unforeseen error. One has to intersect a primitive ! */
+ if(SXD_HIT_NONE(&hit)) {
+ ++nfailures;
+ if(nfailures < max_failures) {
+ continue;
+ } else {
+ res = RES_BAD_ARG;
+ goto error;
+ }
+ }
+ /* Discard the hit if it is on a vertex, i.e. between 2 segments, and
+ * target a new position onto the current primitive */
+ } while((SXD_HIT_NONE(&hit) || HIT_ON_BOUNDARY(&hit))
+ && ++istep < nsteps);
+
+ /* The hits of all targeted positions on the current primitive are on
+ * vertices. Challenge positions on another primitive. */
+ if(istep > nsteps) continue;
+
+ fX(normalize)(N, hit.normal);
+ cos_N_dir = fX(dot)(N, dir);
+
+ if(absf(cos_N_dir) > 1.e-1f) { /* Not roughly orthognonal */
+ const struct sdis_interface* interf;
+ interf = scene_get_interface(scn, hit.prim.prim_id);
+ medium = interface_get_medium
+ (interf, cos_N_dir < 0 ? SDIS_FRONT : SDIS_BACK);
+
+ /* Register the get_medium_info */
+ if(info) {
+ fX(set)(info->pos_tgt, attr.value);
+ fX(set)(info->ray_org, P);
+ fX(set)(info->ray_dir, dir);
+ info->XD(hit) = hit;
+ }
+ break;
+ }
+ }
+
+exit:
+ *out_medium = medium;
+ return res;
+error:
+#if DIM == 2
+ log_err(scn->dev, "%s: could not retrieve the medium at {%g, %g}.\n",
+ FUNC_NAME, SPLIT2(pos));
+#else
+ log_err(scn->dev, "%s: could not retrieve the medium at {%g, %g, %g}.\n",
+ FUNC_NAME, SPLIT3(pos));
+#endif
+ goto exit;
+}
+#undef SDIS_SCENE_DIMENSION
+#undef DIM
+#undef sencXd
+#undef SENCXD
+#undef sXd
+#undef SXD
+#undef SXD_VERTEX_DATA_NULL
+#undef SXD_POSITION
+#undef SXD_FLOAT3
+#undef SXD_TRACE
+#undef SXD_TRACE
+#undef SXD_GET_PRIMITIVE
+#undef SXD_HIT_NONE
+#undef SXD_PRIMITIVE_EQ
+#undef fX
+#undef fX_set_dX
+#undef XD
+#undef HIT_ON_BOUNDARY
+
+#endif /* !SDIS_SCENE_DIMENSION */
diff --git a/src/sdis_scene_c.h b/src/sdis_scene_c.h
@@ -16,12 +16,21 @@
#ifndef SDIS_SCENE_C_H
#define SDIS_SCENE_C_H
-#include <rsys/dynamic_array.h>
-#include <rsys/ref_count.h>
-
#include <star/s2d.h>
#include <star/s3d.h>
+#include <rsys/dynamic_array_uint.h>
+#include <rsys/hash_table.h>
+#include <rsys/ref_count.h>
+
+#include <limits.h>
+
+struct prim_prop {
+ struct sdis_interface* interf;
+ unsigned front_enclosure; /* Id of the front facing enclosure */
+ unsigned back_enclosure; /* Id of the back facing enclosure */
+};
+
struct get_medium_info {
/* Targeted position */
float pos_tgt[3];
@@ -34,26 +43,139 @@ struct get_medium_info {
};
static INLINE void
-interface_init
- (struct mem_allocator* allocator,
- struct sdis_interface** interf)
+prim_prop_init(struct mem_allocator* allocator, struct prim_prop* prim)
+{
+ (void)allocator;
+ prim->interf = NULL;
+ prim->front_enclosure = UINT_MAX;
+ prim->back_enclosure = UINT_MAX;
+}
+
+static INLINE void
+interface_init(struct mem_allocator* allocator, struct sdis_interface** interf)
{
(void)allocator;
*interf = NULL;
}
+static INLINE void
+medium_init(struct mem_allocator* allocator, struct sdis_medium** medium)
+{
+ (void)allocator;
+ *medium = NULL;
+}
+
+struct enclosure {
+ struct s2d_scene_view* s2d_view;
+ struct s3d_scene_view* s3d_view;
+ /* Map the id of the enclosure primitives to their primitive id into the
+ * whole scene */
+ struct darray_uint local2global;
+
+ double hc_upper_bound;
+ double S_over_V; /* in 3D = surface/volume; in 2D = perimeter/area */
+};
+
+static INLINE void
+enclosure_init(struct mem_allocator* allocator, struct enclosure* enc)
+{
+ ASSERT(allocator && enc);
+ enc->s2d_view = NULL;
+ enc->s3d_view = NULL;
+ darray_uint_init(allocator, &enc->local2global);
+ enc->S_over_V = 0;
+ enc->hc_upper_bound = 0;
+}
+
+static INLINE void
+enclosure_release(struct enclosure* enc)
+{
+ if(enc->s2d_view) S2D(scene_view_ref_put(enc->s2d_view));
+ if(enc->s3d_view) S3D(scene_view_ref_put(enc->s3d_view));
+ darray_uint_release(&enc->local2global);
+}
+
+static INLINE res_T
+enclosure_copy(struct enclosure* dst, const struct enclosure* src)
+{
+ if(src->s3d_view) {
+ S3D(scene_view_ref_get(src->s3d_view));
+ dst->s3d_view = src->s3d_view;
+ }
+ if(src->s2d_view) {
+ S2D(scene_view_ref_get(src->s2d_view));
+ dst->s2d_view = src->s2d_view;
+ }
+ dst->S_over_V = src->S_over_V;
+ dst->hc_upper_bound = src->hc_upper_bound;
+ return darray_uint_copy(&dst->local2global, &src->local2global);
+}
+
+static INLINE res_T
+enclosure_copy_and_release(struct enclosure* dst, struct enclosure* src)
+{
+ res_T res = RES_OK;
+ res = darray_uint_copy_and_release(&dst->local2global, &src->local2global);
+ if(res != RES_OK) return res;
+ if(src->s3d_view) {
+ /* Only transfer ownership */
+ dst->s3d_view = src->s3d_view;
+ src->s3d_view = NULL;
+ }
+ if(src->s2d_view) {
+ /* Only transfer ownership */
+ dst->s2d_view = src->s2d_view;
+ src->s2d_view = NULL;
+ }
+ dst->S_over_V = src->S_over_V;
+ dst->hc_upper_bound = src->hc_upper_bound;
+ return RES_OK;
+}
+
/* Declare the array of interfaces */
#define DARRAY_NAME interf
#define DARRAY_DATA struct sdis_interface*
#define DARRAY_FUNCTOR_INIT interface_init
#include <rsys/dynamic_array.h>
+/* Declare the array of medium */
+#define DARRAY_NAME medium
+#define DARRAY_DATA struct sdis_medium*
+#define DARRAY_FUNCTOR_INIT medium_init
+#include <rsys/dynamic_array.h>
+
+/* Declare the array of primitive */
+#define DARRAY_NAME prim_prop
+#define DARRAY_DATA struct prim_prop
+#define DARRAY_FUNCTOR_INIT prim_prop_init
+#include <rsys/dynamic_array.h>
+
+/* Declare the hash table that maps an enclosure id to its data */
+#define HTABLE_NAME enclosure
+#define HTABLE_KEY unsigned
+#define HTABLE_DATA struct enclosure
+#define HTABLE_DATA_FUNCTOR_INIT enclosure_init
+#define HTABLE_DATA_FUNCTOR_RELEASE enclosure_release
+#define HTABLE_DATA_FUNCTOR_COPY enclosure_copy
+#define HTABLE_DATA_FUNCTOR_COPY_AND_RELEASE enclosure_copy_and_release
+#include <rsys/hash_table.h>
+
+/* Declare the hash table that maps an enclosure id to its data */
+#define HTABLE_NAME d
+#define HTABLE_KEY unsigned
+#define HTABLE_DATA double
+#include <rsys/hash_table.h>
+
struct sdis_scene {
struct darray_interf interfaces; /* List of interfaces own by the scene */
- struct darray_interf prim_interfaces; /* Per primitive interface */
+ struct darray_medium media; /* List of media own by the scene */
+ struct darray_prim_prop prim_props; /* Per primitive properties */
struct s2d_scene_view* s2d_view;
struct s3d_scene_view* s3d_view;
+ struct htable_d tmp_hc_ub; /* Map an enclosure id to its hc upper bound */
+ struct htable_enclosure enclosures; /* Map an enclosure id to its data */
+
double ambient_radiative_temperature; /* In Kelvin */
ref_T ref;
@@ -64,7 +186,7 @@ static FINLINE size_t
scene_get_primitives_count(const struct sdis_scene* scn)
{
ASSERT(scn);
- return darray_interf_size_get(&scn->prim_interfaces);
+ return darray_prim_prop_size_get(&scn->prim_props);
}
extern LOCAL_SYM const struct sdis_interface*
@@ -79,6 +201,27 @@ scene_get_medium
struct get_medium_info* info, /* May be NULL */
const struct sdis_medium** medium);
+static INLINE void
+scene_get_enclosure_ids
+ (const struct sdis_scene* scn,
+ const unsigned iprim,
+ unsigned encs[2]) /* Front and Back enclosure identifiers */
+{
+ ASSERT(scn && iprim < darray_prim_prop_size_get(&scn->prim_props));
+ ASSERT(encs);
+ encs[0] = darray_prim_prop_cdata_get(&scn->prim_props)[iprim].front_enclosure;
+ encs[1] = darray_prim_prop_cdata_get(&scn->prim_props)[iprim].back_enclosure;
+}
+
+static INLINE const struct enclosure*
+scene_get_enclosure(struct sdis_scene* scn, const unsigned ienc)
+{
+ const struct enclosure* enc = NULL;
+ ASSERT(scn);
+ enc = htable_enclosure_find(&scn->enclosures, &ienc);
+ return enc;
+}
+
static FINLINE int
scene_is_2d(const struct sdis_scene* scn)
{
diff --git a/src/sdis_solve_Xd.h b/src/sdis_solve_Xd.h
@@ -130,7 +130,7 @@ static const struct XD(rwalk) XD(RWALK_NULL) = {
struct XD(temperature) {
res_T (*func)/* Next function to invoke in order to compute the temperature */
- (const struct sdis_scene* scn,
+ (struct sdis_scene* scn,
const double fp_to_meter,
const struct rwalk_context* ctx,
struct XD(rwalk)* rwalk,
@@ -143,7 +143,7 @@ static const struct XD(temperature) XD(TEMPERATURE_NULL) = { NULL, 0, 0 };
static res_T
XD(boundary_temperature)
- (const struct sdis_scene* scn,
+ (struct sdis_scene* scn,
const double fp_to_meter,
const struct rwalk_context* ctx,
struct XD(rwalk)* rwalk,
@@ -152,7 +152,7 @@ XD(boundary_temperature)
static res_T
XD(solid_temperature)
- (const struct sdis_scene* scn,
+ (struct sdis_scene* scn,
const double fp_to_meter,
const struct rwalk_context* ctx,
struct XD(rwalk)* rwalk,
@@ -161,7 +161,7 @@ XD(solid_temperature)
static res_T
XD(fluid_temperature)
- (const struct sdis_scene* scn,
+ (struct sdis_scene* scn,
const double fp_to_meter,
const struct rwalk_context* ctx,
struct XD(rwalk)* rwalk,
@@ -170,7 +170,7 @@ XD(fluid_temperature)
static res_T
XD(radiative_temperature)
- (const struct sdis_scene* scn,
+ (struct sdis_scene* scn,
const double fp_to_meter,
const struct rwalk_context* ctx,
struct XD(rwalk)* rwalk,
@@ -261,7 +261,7 @@ XD(check_rwalk_fragment_consistency)
static res_T
XD(trace_radiative_path)
- (const struct sdis_scene* scn,
+ (struct sdis_scene* scn,
const float ray_dir[3],
const double fp_to_meter,
const struct rwalk_context* ctx,
@@ -385,7 +385,7 @@ error:
res_T
XD(radiative_temperature)
- (const struct sdis_scene* scn,
+ (struct sdis_scene* scn,
const double fp_to_meter,
const struct rwalk_context* ctx,
struct XD(rwalk)* rwalk,
@@ -424,26 +424,192 @@ error:
res_T
XD(fluid_temperature)
- (const struct sdis_scene* scn,
+ (struct sdis_scene* scn,
const double fp_to_meter,
const struct rwalk_context* ctx,
struct XD(rwalk)* rwalk,
struct ssp_rng* rng,
struct XD(temperature)* T)
{
+ struct sXd(attrib) attr_P, attr_N;
+ struct sdis_interface_fragment frag;
+ const struct sdis_interface* interf;
+ const struct enclosure* enc;
+ unsigned enc_ids[2];
+ unsigned enc_id;
+ double rho; /* Volumic mass */
+ double hc; /* Convection coef */
+ double cp; /* Calorific capacity */
+ double mu;
+ double tau;
double tmp;
+ double r;
+#if DIM == 2
+ float st;
+#else
+ float st[2];
+#endif
(void)rng, (void)fp_to_meter, (void)ctx;
ASSERT(scn && fp_to_meter > 0 && ctx && rwalk && rng && T);
ASSERT(rwalk->mdm->type == SDIS_FLUID);
tmp = fluid_get_temperature(rwalk->mdm, &rwalk->vtx);
- if(tmp < 0) {
- log_err(scn->dev, "%s: unknown fluid temperature is not supported.\n",
+ if(tmp >= 0) { /* T is known. */
+ T->value += tmp;
+ T->done = 1;
+ return RES_OK;
+ }
+
+ if(SXD_HIT_NONE(&rwalk->hit)) { /* The path begins in the fluid */
+ const float range[2] = {0, FLT_MAX};
+ float dir[DIM] = {0};
+ float org[DIM];
+
+ dir[DIM-1] = 1;
+ fX_set_dX(org, rwalk->vtx.P);
+
+ /* Init the path hit field required to define the current enclosure and
+ * fetch the interface data */
+ SXD(scene_view_trace_ray(scn->sXd(view), org, dir, range, NULL, &rwalk->hit));
+ rwalk->hit_side = fX(dot)(rwalk->hit.normal, dir) < 0 ? SDIS_FRONT : SDIS_BACK;
+
+ if(SXD_HIT_NONE(&rwalk->hit)) {
+ log_err(scn->dev,
+"%s: the position %g %g %g lies in the surrounding fluid whose temperature must \n"
+"be known.\n",
+ FUNC_NAME, SPLIT3(rwalk->vtx.P));
+ return RES_BAD_OP;
+ }
+ }
+
+ /* Fetch the current interface and its associated enclosures */
+ interf = scene_get_interface(scn, rwalk->hit.prim.prim_id);
+ scene_get_enclosure_ids(scn, rwalk->hit.prim.prim_id, enc_ids);
+
+ /* Define the enclosure identifier of the current medium */
+ ASSERT(interf->medium_front != interf->medium_back);
+ if(rwalk->mdm == interf->medium_front) {
+ enc_id = enc_ids[0];
+ ASSERT(rwalk->hit_side == SDIS_FRONT);
+ } else {
+ ASSERT(rwalk->mdm == interf->medium_back);
+ enc_id = enc_ids[1];
+ ASSERT(rwalk->hit_side == SDIS_BACK);
+ }
+
+ /* Fetch the enclosure data */
+ enc = scene_get_enclosure(scn, enc_id);
+ if(!enc) {
+ log_err(scn->dev,
+"%s: invalid enclosure. The position %g %g %g may lie in the surrounding fluid.\n",
+ FUNC_NAME, SPLIT3(rwalk->vtx.P));
+ return RES_BAD_OP;
+ }
+
+ /* The hc upper bound can be 0 is h is uniformly 0. In that case the result
+ * is the initial condition. */
+ if(enc->hc_upper_bound == 0) {
+ /* Cannot be in the fluid without starting there. */
+ ASSERT(SXD_HIT_NONE(&rwalk->hit));
+ rwalk->vtx.time = 0;
+ tmp = fluid_get_temperature(rwalk->mdm, &rwalk->vtx);
+ if(tmp >= 0) {
+ T->value += tmp;
+ T->done = 1;
+ return RES_OK;
+ }
+
+ /* At t=0, the initial condition should have been reached. */
+ log_err(scn->dev,
+"%s: undefined initial condition. "
+"Time is 0 but the temperature remains unknown.\n",
FUNC_NAME);
return RES_BAD_OP;
}
- T->value += tmp;
- T->done = 1;
+
+ /* A trick to force first r test result. */
+ r = 1;
+
+ /* Sample time until init condition is reached or a true convection occurs. */
+ for(;;) {
+ /* Setup the fragment of the interface. */
+ XD(setup_interface_fragment)(&frag, &rwalk->vtx, &rwalk->hit, rwalk->hit_side);
+
+ /* Fetch hc. */
+ hc = interface_get_convection_coef(interf, &frag);
+ if(hc > enc->hc_upper_bound) {
+ log_err(scn->dev,
+ "%s: hc (%g) exceeds its provided upper bound (%g) at %g %g %g.\n",
+ FUNC_NAME, hc, enc->hc_upper_bound, SPLIT3(rwalk->vtx.P));
+ return RES_BAD_OP;
+ }
+
+ if(r < hc / enc->hc_upper_bound) {
+ /* True convection. Always true if hc == bound. */
+ break;
+ }
+
+ /* Fetch other physical properties. */
+ cp = fluid_get_calorific_capacity(rwalk->mdm, &rwalk->vtx);
+ rho = fluid_get_volumic_mass(rwalk->mdm, &rwalk->vtx);
+
+ /* Sample the time using the upper bound. */
+ mu = enc->hc_upper_bound / (rho * cp) * enc->S_over_V;
+ tau = ssp_ran_exp(rng, mu);
+ rwalk->vtx.time = MMAX(rwalk->vtx.time - tau, 0);
+
+ /* Check the initial condition. */
+ tmp = fluid_get_temperature(rwalk->mdm, &rwalk->vtx);
+ if(tmp >= 0) {
+ T->value += tmp;
+ T->done = 1;
+ return RES_OK;
+ }
+
+ if(rwalk->vtx.time <= 0) {
+ /* The initial condition should have been reached. */
+ log_err(scn->dev,
+"%s: undefined initial condition. "
+"Time is 0 but the temperature remains unknown.\n",
+ FUNC_NAME);
+ return RES_BAD_OP;
+ }
+
+ /* Uniformly sample the enclosure. */
+#if DIM == 2
+ SXD(scene_view_sample
+ (enc->sXd(view),
+ ssp_rng_canonical_float(rng),
+ ssp_rng_canonical_float(rng),
+ &rwalk->hit.prim,
+ &rwalk->hit.u));
+ st = rwalk->hit.u;
+#else
+ SXD(scene_view_sample
+ (enc->sXd(view),
+ ssp_rng_canonical_float(rng),
+ ssp_rng_canonical_float(rng),
+ ssp_rng_canonical_float(rng),
+ &rwalk->hit.prim,
+ rwalk->hit.uv));
+ f2_set(st, rwalk->hit.uv);
+#endif
+
+ SXD(primitive_get_attrib(&rwalk->hit.prim, SXD_POSITION, st, &attr_P));
+ SXD(primitive_get_attrib(&rwalk->hit.prim, SXD_GEOMETRY_NORMAL, st, &attr_N));
+ dX_set_fX(rwalk->vtx.P, attr_P.value);
+ fX(set)(rwalk->hit.normal, attr_N.value);
+
+ /* Fetch the interface of the sampled point. */
+ interf = scene_get_interface(scn, rwalk->hit.prim.prim_id);
+
+ /* Renew r for next loop. */
+ r = ssp_rng_canonical_float(rng);
+ }
+
+ rwalk->hit.distance = 0;
+ T->func = XD(boundary_temperature);
+ rwalk->mdm = NULL; /* The random walk is at an interface between 2 media */
return RES_OK;
}
@@ -861,7 +1027,7 @@ XD(solid_boundary_with_flux_temperature)
res_T
XD(boundary_temperature)
- (const struct sdis_scene* scn,
+ (struct sdis_scene* scn,
const double fp_to_meter,
const struct rwalk_context* ctx,
struct XD(rwalk)* rwalk,
@@ -920,7 +1086,7 @@ XD(boundary_temperature)
res_T
XD(solid_temperature)
- (const struct sdis_scene* scn,
+ (struct sdis_scene* scn,
const double fp_to_meter,
const struct rwalk_context* ctx,
struct XD(rwalk)* rwalk,
@@ -1049,7 +1215,7 @@ XD(solid_temperature)
/* 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;
+ rwalk->vtx.time = MMAX(rwalk->vtx.time - tau, 0);
/* Check the initial condition */
tmp = solid_get_temperature(mdm, &rwalk->vtx);
@@ -1059,6 +1225,15 @@ XD(solid_temperature)
return RES_OK;
}
+ /* The initial condition should be reached */
+ if(rwalk->vtx.time <=0) {
+ log_err(scn->dev,
+ "%s: undefined initial condition. "
+ "The time is null but the temperature remains unknown.\n",
+ FUNC_NAME);
+ return RES_BAD_OP;
+ }
+
/* Define if the random walk hits something along dir0 */
if(hit0.distance > delta) {
rwalk->hit = SXD_HIT_NULL;
diff --git a/src/test_sdis_conducto_radiative.c b/src/test_sdis_conducto_radiative.c
@@ -239,6 +239,7 @@ create_interface
shader.back.emissivity = interface_get_emissivity;
shader.back.specular_fraction = interface_get_specular_fraction;
}
+ shader.convection_coef_upper_bound = MMAX(0, interf->convection_coef);
CHK(sdis_data_create(dev, sizeof(struct interfac), ALIGNOF(struct interfac),
NULL, &data) == RES_OK);
diff --git a/src/test_sdis_conducto_radiative_2d.c b/src/test_sdis_conducto_radiative_2d.c
@@ -255,6 +255,8 @@ create_interface
shader.back.emissivity = interface_get_emissivity;
shader.back.specular_fraction = interface_get_specular_fraction;
}
+ shader.convection_coef_upper_bound = MMAX(0, interf->convection_coef);
+
CHK(sdis_data_create(dev, sizeof(struct interfac), ALIGNOF(struct interfac),
NULL, &data) == RES_OK);
*((struct interfac*)sdis_data_get(data)) = *interf;
diff --git a/src/test_sdis_convection.c b/src/test_sdis_convection.c
@@ -0,0 +1,310 @@
+/* Copyright (C) 2016-2018 |Meso|Star> (contact@meso-star.com)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. */
+
+#include "sdis.h"
+#include "test_sdis_utils.h"
+
+#include <rsys/double3.h>
+#include <rsys/math.h>
+
+/*
+ * The scene is composed of an unit fluid cube/square whose temperature is
+ * unknown. The convection coefficient with the surrounding solid is H
+ * everywhere the temperature of the -/+X, -/+Y and -/+Z faces are fixed to T0
+ * and T1, T2, T3, T4 and T5, respectively. This test computes the temperature
+ * of the fluid Tf at an observation time t. This temperature is equal to:
+ *
+ * Tf(t) = Tf(0) * e^(-nu*t) + Tinf*(1-e^(-nu*t))
+ *
+ * nu = (Sum_{i=0..5}(H*Si)) / (RHO*CP*V)
+ * Tinf = (Sum_{i=0..5}(H*Si*Ti) / (Sum_{i=0..5}(H*Si));
+ *
+ * with Si surface of the faces (i.e. one), V the volume of the cube (i.e.
+ * one), RHO the volumic mass of the fluid and CP its calorific capacity.
+ *
+ * 3D 2D
+ *
+ * (1,1,1) (1,1)
+ * +---------+ +-----T3----+
+ * /' T3 /|T4 | |
+ * +---------+ | | H _\ |
+ * | ' H _\ |T1 T0 / / T1
+ * |T0 / / | | | \__/ |
+ * | +..\__/.|.+ | |
+ * T5|, T2 |/ +-----------+
+ * +---------+ (0,0)
+ * (0,0,0)
+ */
+
+#define UNKNOWN_TEMPERATURE -1
+#define N 100000 /* #realisations */
+
+#define Tf_0 280.0
+
+#define T0 300.0
+#define T1 310.0
+#define T2 320.0
+#define T3 330.0
+#define T4 340.0
+#define T5 350.0
+
+#define H 10.0
+#define RHO 25.0
+#define CP 2.0
+
+/*******************************************************************************
+ * Media
+ ******************************************************************************/
+static double
+fluid_get_temperature
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ (void)data;
+ CHK(vtx != NULL);
+ return vtx->time <= 0 ? Tf_0 : UNKNOWN_TEMPERATURE;
+}
+
+static double
+fluid_get_volumic_mass
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ (void)data;
+ CHK(vtx != NULL);
+ return RHO;
+}
+
+static double
+fluid_get_calorific_capacity
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ (void)data;
+ CHK(vtx != NULL);
+ return CP;
+}
+
+/*******************************************************************************
+ * Interface
+ ******************************************************************************/
+struct interf {
+ double temperature;
+};
+
+static double
+interface_get_temperature
+ (const struct sdis_interface_fragment* frag, struct sdis_data* data)
+{
+ const struct interf* interf = sdis_data_cget(data);
+ CHK(frag && data);
+ return interf->temperature;
+}
+
+static double
+interface_get_convection_coef
+ (const struct sdis_interface_fragment* frag, struct sdis_data* data)
+{
+ CHK(frag && data);
+ return H;
+}
+
+static double
+interface_get_emissivity
+ (const struct sdis_interface_fragment* frag, struct sdis_data* data)
+{
+ CHK(frag && data);
+ return 0;
+}
+
+static double
+interface_get_specular_fraction
+ (const struct sdis_interface_fragment* frag, struct sdis_data* data)
+{
+ CHK(frag && data);
+ return 0;
+}
+
+static struct sdis_interface*
+create_interface
+ (struct sdis_device* dev,
+ struct sdis_medium* front,
+ struct sdis_medium* back,
+ const struct sdis_interface_shader* interf_shader,
+ const double temperature)
+{
+ struct sdis_data* data;
+ struct sdis_interface* interf;
+ struct interf* interf_props;
+
+ CHK(sdis_data_create
+ (dev, sizeof(struct interf), ALIGNOF(struct interf), NULL, &data) == RES_OK);
+ interf_props = sdis_data_get(data);
+ interf_props->temperature = temperature;
+ CHK(sdis_interface_create
+ (dev, front, back, interf_shader, data, &interf) == RES_OK);
+ CHK(sdis_data_ref_put(data) == RES_OK);
+ return interf;
+}
+
+/*******************************************************************************
+ * 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* fluid = NULL;
+ struct sdis_medium* solid = NULL;
+ struct sdis_interface* interf_T0 = NULL;
+ struct sdis_interface* interf_T1 = NULL;
+ struct sdis_interface* interf_T2 = NULL;
+ struct sdis_interface* interf_T3 = NULL;
+ struct sdis_interface* interf_T4 = NULL;
+ struct sdis_interface* interf_T5 = NULL;
+ struct sdis_scene* box_scn = NULL;
+ struct sdis_scene* square_scn = NULL;
+ struct sdis_estimator* estimator = NULL;
+ struct sdis_fluid_shader fluid_shader = DUMMY_FLUID_SHADER;
+ struct sdis_solid_shader solid_shader = DUMMY_SOLID_SHADER;
+ struct sdis_interface_shader interf_shader = DUMMY_INTERFACE_SHADER;
+ struct sdis_interface* box_interfaces[12/*#triangles*/];
+ struct sdis_interface* square_interfaces[4/*#segments*/];
+ double pos[3];
+ double ref;
+ double Tinf;
+ double nu;
+ size_t nreals;
+ size_t nfails;
+ int 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 = fluid_get_temperature;
+ fluid_shader.calorific_capacity = fluid_get_calorific_capacity;
+ fluid_shader.volumic_mass = fluid_get_volumic_mass;
+ CHK(sdis_fluid_create(dev, &fluid_shader, NULL, &fluid) == RES_OK);
+
+ /* Create the solid_medium */
+ CHK(sdis_solid_create(dev, &solid_shader, NULL, &solid) == RES_OK);
+
+ /* Setup the interface shader */
+ interf_shader.convection_coef = interface_get_convection_coef;
+ interf_shader.front.temperature = interface_get_temperature;
+ interf_shader.front.emissivity = interface_get_emissivity;
+ interf_shader.front.specular_fraction = interface_get_specular_fraction;
+ interf_shader.convection_coef_upper_bound = H;
+
+ /* Create the interfaces */
+ interf_T0 = create_interface(dev, fluid, solid, &interf_shader, T0);
+ interf_T1 = create_interface(dev, fluid, solid, &interf_shader, T1);
+ interf_T2 = create_interface(dev, fluid, solid, &interf_shader, T2);
+ interf_T3 = create_interface(dev, fluid, solid, &interf_shader, T3);
+ interf_T4 = create_interface(dev, fluid, solid, &interf_shader, T4);
+ interf_T5 = create_interface(dev, fluid, solid, &interf_shader, T5);
+
+ /* Release the media */
+ CHK(sdis_medium_ref_put(solid) == RES_OK);
+ CHK(sdis_medium_ref_put(fluid) == RES_OK);
+
+ /* Map the interfaces to their box triangles */
+ box_interfaces[0] = box_interfaces[1] = interf_T5; /* Front */
+ box_interfaces[2] = box_interfaces[3] = interf_T0; /* Left */
+ box_interfaces[4] = box_interfaces[5] = interf_T4; /* Back */
+ box_interfaces[6] = box_interfaces[7] = interf_T1; /* Right */
+ box_interfaces[8] = box_interfaces[9] = interf_T3; /* Top */
+ box_interfaces[10]= box_interfaces[11]= interf_T2; /* Bottom */
+
+ /* Map the interfaces to their square segments */
+ square_interfaces[0] = interf_T2; /* Bottom */
+ square_interfaces[1] = interf_T0; /* Left */
+ square_interfaces[2] = interf_T3; /* Top */
+ square_interfaces[3] = interf_T1; /* Right */
+
+ /* Create the box scene */
+ CHK(sdis_scene_create(dev, box_ntriangles, box_get_indices,
+ box_get_interface, box_nvertices, box_get_position, box_interfaces,
+ &box_scn) == RES_OK);
+
+ /* Create the square scene */
+ CHK(sdis_scene_2d_create(dev, square_nsegments, square_get_indices,
+ square_get_interface, square_nvertices, square_get_position,
+ square_interfaces, &square_scn) == RES_OK);
+
+ /* Release the interfaces */
+ CHK(sdis_interface_ref_put(interf_T0) == RES_OK);
+ CHK(sdis_interface_ref_put(interf_T1) == RES_OK);
+ CHK(sdis_interface_ref_put(interf_T2) == RES_OK);
+ CHK(sdis_interface_ref_put(interf_T3) == RES_OK);
+ CHK(sdis_interface_ref_put(interf_T4) == RES_OK);
+ CHK(sdis_interface_ref_put(interf_T5) == RES_OK);
+
+ d3_splat(pos, 0.25);
+
+ /* Test in 3D for various time values. */
+ nu = (6 * H) / (RHO*CP);
+ Tinf = (H*(T0 + T1 + T2 + T3 + T4 + T5)) / (6 * H);
+ printf("Temperature of the box at (%g %g %g)\n", SPLIT3(pos));
+ FOR_EACH(i, 0, 5) {
+ double time = i ? (double) i / nu : INF;
+ ref = Tf_0 * exp(-nu * time) + Tinf * (1 - exp(-nu * time));
+
+ /* Solve in 3D */
+ CHK(sdis_solve_probe(box_scn, N, pos, time, 1.0, 0, 0, &estimator) == RES_OK);
+ CHK(sdis_estimator_get_realisation_count(estimator, &nreals) == RES_OK);
+ CHK(sdis_estimator_get_failure_count(estimator, &nfails) == RES_OK);
+ CHK(nfails + nreals == N);
+ CHK(sdis_estimator_get_temperature(estimator, &T) == RES_OK);
+ CHK(sdis_estimator_ref_put(estimator) == RES_OK);
+ printf(" t=%g : %g ~ %g +/- %g\n", time, ref, T.E, T.SE);
+ if(nfails)
+ printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N);
+ CHK(eq_eps(T.E, ref, T.SE * 3));
+ }
+
+ /* Test in 2D for various time values. */
+ nu = (4 * H) / (RHO*CP);
+ Tinf = (H * (T0 + T1 + T2 + T3)) / (4 * H);
+ printf("Temperature of the square at (%g %g)\n", SPLIT2(pos));
+ FOR_EACH(i, 0, 5) {
+ double time = i ? (double) i / nu : INF;
+ ref = Tf_0 * exp(-nu * time) + Tinf * (1 - exp(-nu * time));
+
+ /* Solve in 2D */
+ CHK(sdis_solve_probe(square_scn, N, pos, time, 1.0, 0, 0, &estimator) == RES_OK);
+ CHK(sdis_estimator_get_realisation_count(estimator, &nreals) == RES_OK);
+ CHK(sdis_estimator_get_failure_count(estimator, &nfails) == RES_OK);
+ CHK(nfails + nreals == N);
+ CHK(sdis_estimator_get_temperature(estimator, &T) == RES_OK);
+ CHK(sdis_estimator_ref_put(estimator) == RES_OK);
+ printf(" t=%g : %g ~ %g +/- %g\n", time, ref, T.E, T.SE);
+ if(nfails)
+ printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N);
+ CHK(eq_eps(T.E, ref, T.SE * 3));
+ }
+
+ CHK(sdis_scene_ref_put(box_scn) == RES_OK);
+ CHK(sdis_scene_ref_put(square_scn) == RES_OK);
+ CHK(sdis_device_ref_put(dev) == RES_OK);
+
+ check_memory_allocator(&allocator);
+ mem_shutdown_proxy_allocator(&allocator);
+ CHK(mem_allocated_size() == 0);
+ return 0;
+}
+
diff --git a/src/test_sdis_convection_non_uniform.c b/src/test_sdis_convection_non_uniform.c
@@ -0,0 +1,325 @@
+/* Copyright (C) 2016-2018 |Meso|Star> (contact@meso-star.com)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. */
+
+#include "sdis.h"
+#include "test_sdis_utils.h"
+
+#include <rsys/double3.h>
+#include <rsys/math.h>
+
+/*
+ * The scene is composed of an unit fluid cube/square whose temperature is
+ * unknown. The convection coefficient with the surrounding solid is H
+ * everywhere the temperature of the -/+X, -/+Y and -/+Z faces are fixed to T0
+ * and T1, T2, T3, T4 and T5, respectively. This test computes the temperature
+ * of the fluid Tf at an observation time t. This temperature is equal to:
+ *
+ * Tf(t) = Tf(0) * e^(-nu*t) + Tinf*(1-e^(-nu*t))
+ *
+ * nu = (Sum_{i=0..5}(H*Si)) / (RHO*CP*V)
+ * Tinf = (Sum_{i=0..5}(H*Si*Ti) / (Sum_{i=0..5}(H*Si));
+ *
+ * with Si surface of the faces (i.e. one), V the volume of the cube (i.e.
+ * one), RHO the volumic mass of the fluid and CP its calorific capacity.
+ *
+ * 3D 2D
+ *
+ * (1,1,1) (1,1)
+ * +---------+ +-----T3----+
+ * /' T3 /|T4 | |
+ * +---------+ | | H _\ |
+ * | ' H _\ |T1 T0 / / T1
+ * |T0 / / | | | \__/ |
+ * | +..\__/.|.+ | |
+ * T5|, T2 |/ +-----------+
+ * +---------+ (0,0)
+ * (0,0,0)
+ */
+
+#define UNKNOWN_TEMPERATURE -1
+#define N 100000 /* #realisations */
+
+#define Tf_0 280.0
+
+#define T0 300.0
+#define T1 310.0
+#define T2 320.0
+#define T3 330.0
+#define T4 340.0
+#define T5 350.0
+
+#define HC0 100.0
+#define HC1 30.0
+#define HC2 3020.0
+#define HC3 7300.0
+#define HC4 3400.0
+#define HC5 50.0
+
+#define RHO 25.0
+#define CP 2.0
+
+/*******************************************************************************
+ * Media
+ ******************************************************************************/
+static double
+fluid_get_temperature
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ (void)data;
+ CHK(vtx != NULL);
+ return vtx->time <= 0 ? Tf_0 : UNKNOWN_TEMPERATURE;
+}
+
+static double
+fluid_get_volumic_mass
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ (void)data;
+ CHK(vtx != NULL);
+ return RHO;
+}
+
+static double
+fluid_get_calorific_capacity
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ (void)data;
+ CHK(vtx != NULL);
+ return CP;
+}
+
+/*******************************************************************************
+ * Interface
+ ******************************************************************************/
+struct interf {
+ double temperature;
+ double hc;
+};
+
+static double
+interface_get_temperature
+ (const struct sdis_interface_fragment* frag, struct sdis_data* data)
+{
+ const struct interf* interf = sdis_data_cget(data);
+ CHK(frag && data);
+ return interf->temperature;
+}
+
+static double
+interface_get_convection_coef
+ (const struct sdis_interface_fragment* frag, struct sdis_data* data)
+{
+ const struct interf* interf = sdis_data_cget(data);
+ CHK(frag && data);
+ return interf->hc;
+}
+
+static double
+interface_get_emissivity
+ (const struct sdis_interface_fragment* frag, struct sdis_data* data)
+{
+ CHK(frag && data);
+ return 0;
+}
+
+static double
+interface_get_specular_fraction
+ (const struct sdis_interface_fragment* frag, struct sdis_data* data)
+{
+ CHK(frag && data);
+ return 0;
+}
+
+static struct sdis_interface*
+create_interface
+ (struct sdis_device* dev,
+ struct sdis_medium* front,
+ struct sdis_medium* back,
+ const struct sdis_interface_shader* interf_shader,
+ const double temperature,
+ const double hc)
+{
+ struct sdis_data* data;
+ struct sdis_interface* interf;
+ struct interf* interf_props;
+
+ CHK(sdis_data_create
+ (dev, sizeof(struct interf), ALIGNOF(struct interf), NULL, &data) == RES_OK);
+ interf_props = sdis_data_get(data);
+ interf_props->temperature = temperature;
+ interf_props->hc = hc;
+ CHK(sdis_interface_create
+ (dev, front, back, interf_shader, data, &interf) == RES_OK);
+ CHK(sdis_data_ref_put(data) == RES_OK);
+ return interf;
+}
+
+/*******************************************************************************
+ * 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* fluid = NULL;
+ struct sdis_medium* solid = NULL;
+ struct sdis_interface* interf_T0 = NULL;
+ struct sdis_interface* interf_T1 = NULL;
+ struct sdis_interface* interf_T2 = NULL;
+ struct sdis_interface* interf_T3 = NULL;
+ struct sdis_interface* interf_T4 = NULL;
+ struct sdis_interface* interf_T5 = NULL;
+ struct sdis_scene* box_scn = NULL;
+ struct sdis_scene* square_scn = NULL;
+ struct sdis_estimator* estimator = NULL;
+ struct sdis_fluid_shader fluid_shader = DUMMY_FLUID_SHADER;
+ struct sdis_solid_shader solid_shader = DUMMY_SOLID_SHADER;
+ struct sdis_interface_shader interf_shader = DUMMY_INTERFACE_SHADER;
+ struct sdis_interface* box_interfaces[12/*#triangles*/];
+ struct sdis_interface* square_interfaces[4/*#segments*/];
+ double pos[3];
+ double ref;
+ double Tinf;
+ double nu;
+ size_t nreals;
+ size_t nfails;
+ int 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 = fluid_get_temperature;
+ fluid_shader.calorific_capacity = fluid_get_calorific_capacity;
+ fluid_shader.volumic_mass = fluid_get_volumic_mass;
+ CHK(sdis_fluid_create(dev, &fluid_shader, NULL, &fluid) == RES_OK);
+
+ /* Create the solid_medium */
+ CHK(sdis_solid_create(dev, &solid_shader, NULL, &solid) == RES_OK);
+
+ /* Setup the interface shader */
+ interf_shader.convection_coef = interface_get_convection_coef;
+ interf_shader.front.temperature = interface_get_temperature;
+ interf_shader.front.emissivity = interface_get_emissivity;
+ interf_shader.front.specular_fraction = interface_get_specular_fraction;
+
+ /* Create the interfaces */
+ interf_shader.convection_coef_upper_bound = HC0;
+ interf_T0 = create_interface(dev, fluid, solid, &interf_shader, T0, HC0);
+ interf_shader.convection_coef_upper_bound = HC1;
+ interf_T1 = create_interface(dev, fluid, solid, &interf_shader, T1, HC1);
+ interf_shader.convection_coef_upper_bound = HC2;
+ interf_T2 = create_interface(dev, fluid, solid, &interf_shader, T2, HC2);
+ interf_shader.convection_coef_upper_bound = HC3;
+ interf_T3 = create_interface(dev, fluid, solid, &interf_shader, T3, HC3);
+ interf_shader.convection_coef_upper_bound = HC4;
+ interf_T4 = create_interface(dev, fluid, solid, &interf_shader, T4, HC4);
+ interf_shader.convection_coef_upper_bound = HC5;
+ interf_T5 = create_interface(dev, fluid, solid, &interf_shader, T5, HC5);
+
+ /* Release the media */
+ CHK(sdis_medium_ref_put(solid) == RES_OK);
+ CHK(sdis_medium_ref_put(fluid) == RES_OK);
+
+ /* Map the interfaces to their box triangles */
+ box_interfaces[0] = box_interfaces[1] = interf_T5; /* Front */
+ box_interfaces[2] = box_interfaces[3] = interf_T0; /* Left */
+ box_interfaces[4] = box_interfaces[5] = interf_T4; /* Back */
+ box_interfaces[6] = box_interfaces[7] = interf_T1; /* Right */
+ box_interfaces[8] = box_interfaces[9] = interf_T3; /* Top */
+ box_interfaces[10]= box_interfaces[11]= interf_T2; /* Bottom */
+
+ /* Map the interfaces to their square segments */
+ square_interfaces[0] = interf_T2; /* Bottom */
+ square_interfaces[1] = interf_T0; /* Left */
+ square_interfaces[2] = interf_T3; /* Top */
+ square_interfaces[3] = interf_T1; /* Right */
+
+ /* Create the box scene */
+ CHK(sdis_scene_create(dev, box_ntriangles, box_get_indices,
+ box_get_interface, box_nvertices, box_get_position, box_interfaces,
+ &box_scn) == RES_OK);
+
+ /* Create the square scene */
+ CHK(sdis_scene_2d_create(dev, square_nsegments, square_get_indices,
+ square_get_interface, square_nvertices, square_get_position,
+ square_interfaces, &square_scn) == RES_OK);
+
+ /* Release the interfaces */
+ CHK(sdis_interface_ref_put(interf_T0) == RES_OK);
+ CHK(sdis_interface_ref_put(interf_T1) == RES_OK);
+ CHK(sdis_interface_ref_put(interf_T2) == RES_OK);
+ CHK(sdis_interface_ref_put(interf_T3) == RES_OK);
+ CHK(sdis_interface_ref_put(interf_T4) == RES_OK);
+ CHK(sdis_interface_ref_put(interf_T5) == RES_OK);
+
+ d3_splat(pos, 0.25);
+
+ /* Test in 3D for various time values. */
+ nu = (HC0 + HC1 + HC2 + HC3 + HC4 + HC5) / (RHO * CP);
+ Tinf = (HC0 * T0 + HC1 * T1 + HC2 * T2 + HC3 * T3 + HC4 * T4 + HC5 * T5)
+ / (HC0 + HC1 + HC2 + HC3 + HC4 + HC5);
+ printf("Temperature of the box at (%g %g %g)\n", SPLIT3(pos));
+ FOR_EACH(i, 0, 5) {
+ double time = i ? (double) i / nu : INF;
+ ref = Tf_0 * exp(-nu * time) + Tinf * (1 - exp(-nu * time));
+
+ /* Solve in 3D */
+ CHK(sdis_solve_probe(box_scn, N, pos, time, 1.0, 0, 0, &estimator) == RES_OK);
+ CHK(sdis_estimator_get_realisation_count(estimator, &nreals) == RES_OK);
+ CHK(sdis_estimator_get_failure_count(estimator, &nfails) == RES_OK);
+ CHK(nfails + nreals == N);
+ CHK(sdis_estimator_get_temperature(estimator, &T) == RES_OK);
+ CHK(sdis_estimator_ref_put(estimator) == RES_OK);
+ printf(" t=%g : %g ~ %g +/- %g\n", time, ref, T.E, T.SE);
+ if(nfails)
+ printf("#failures = %lu/%lu\n", (unsigned long)nfails,(unsigned long)N);
+ CHK(eq_eps(T.E, ref, T.SE * 3));
+ }
+
+ /* Test in 2D for various time values. */
+ nu = (HC0 + HC1 + HC2 + HC3) / (RHO * CP);
+ Tinf = (HC0 * T0 + HC1 * T1 + HC2 * T2 + HC3 * T3) / (HC0 + HC1 + HC2 + HC3);
+ printf("Temperature of the square at (%g %g)\n", SPLIT2(pos));
+ FOR_EACH(i, 0, 5) {
+ double time = i ? (double) i / nu : INF;
+ ref = Tf_0 * exp(-nu * time) + Tinf * (1 - exp(-nu * time));
+
+ CHK(sdis_solve_probe(square_scn, N, pos, time, 1.0, 0, 0, &estimator) == RES_OK);
+ CHK(sdis_estimator_get_realisation_count(estimator, &nreals) == RES_OK);
+ CHK(sdis_estimator_get_failure_count(estimator, &nfails) == RES_OK);
+ CHK(nfails + nreals == N);
+ CHK(sdis_estimator_get_temperature(estimator, &T) == RES_OK);
+ CHK(sdis_estimator_ref_put(estimator) == RES_OK);
+ printf(" t=%g : %g ~ %g +/- %g\n", time, ref, T.E, T.SE);
+ if (nfails)
+ printf("#failures = %lu/%lu\n", (unsigned long)nfails,(unsigned long)N);
+ CHK(eq_eps(T.E, ref, T.SE * 3));
+ }
+
+ CHK(sdis_scene_ref_put(box_scn) == RES_OK);
+ CHK(sdis_scene_ref_put(square_scn) == RES_OK);
+ CHK(sdis_device_ref_put(dev) == RES_OK);
+
+ check_memory_allocator(&allocator);
+ mem_shutdown_proxy_allocator(&allocator);
+ CHK(mem_allocated_size() == 0);
+ return 0;
+}
+
diff --git a/src/test_sdis_flux.c b/src/test_sdis_flux.c
@@ -49,68 +49,6 @@
#define LAMBDA 0.1
/*******************************************************************************
- * Geometry 3D
- ******************************************************************************/
-static void
-box_get_indices(const size_t itri, size_t ids[3], void* context)
-{
- (void)context;
- CHK(ids);
- ids[0] = box_indices[itri*3+0];
- ids[1] = box_indices[itri*3+1];
- ids[2] = box_indices[itri*3+2];
-}
-
-static void
-box_get_position(const size_t ivert, double pos[3], void* context)
-{
- (void)context;
- CHK(pos);
- pos[0] = box_vertices[ivert*3+0];
- pos[1] = box_vertices[ivert*3+1];
- pos[2] = box_vertices[ivert*3+2];
-}
-
-static void
-box_get_interface(const size_t itri, struct sdis_interface** bound, void* context)
-{
- struct sdis_interface** interfaces = context;
- CHK(context && bound);
- *bound = interfaces[itri];
-}
-
-/*******************************************************************************
- * Geometry 2D
- ******************************************************************************/
-static void
-square_get_indices(const size_t iseg, size_t ids[2], void* context)
-{
- (void)context;
- CHK(ids);
- ids[0] = square_indices[iseg*2+0];
- ids[1] = square_indices[iseg*2+1];
-}
-
-static void
-square_get_position(const size_t ivert, double pos[2], void* context)
-{
- (void)context;
- CHK(pos);
- pos[0] = square_vertices[ivert*2+0];
- pos[1] = square_vertices[ivert*2+1];
-}
-
-static void
-square_get_interface
- (const size_t iseg, struct sdis_interface** bound, void* context)
-{
- struct sdis_interface** interfaces = context;
- CHK(context && bound);
- *bound = interfaces[iseg];
-}
-
-
-/*******************************************************************************
* Media
******************************************************************************/
static double
diff --git a/src/test_sdis_interface.c b/src/test_sdis_interface.c
@@ -103,6 +103,12 @@ main(int argc, char** argv)
shader.front.specular_fraction = dummy_interface_getter;
CHK(CREATE(dev, solid, fluid, &shader, NULL, &interf) == RES_OK); /* Warning */
CHK(sdis_interface_ref_put(interf) == RES_OK);
+ shader.front.specular_fraction = NULL;
+ shader.convection_coef_upper_bound = -1;
+ CHK(CREATE(dev, solid, solid, &shader, NULL, &interf) == RES_OK); /* Warning */
+ CHK(sdis_interface_ref_put(interf) == RES_OK);
+ CHK(CREATE(dev, solid, fluid, &shader, NULL, &interf) == RES_BAD_ARG);
+ shader.convection_coef_upper_bound = 0;
#undef CREATE
CHK(sdis_device_ref_put(dev) == RES_OK);
diff --git a/src/test_sdis_solve_probe_boundary.c b/src/test_sdis_solve_probe_boundary.c
@@ -51,68 +51,6 @@
#define LAMBDA 0.1
/*******************************************************************************
- * Geometry 3D
- ******************************************************************************/
-static void
-box_get_indices(const size_t itri, size_t ids[3], void* context)
-{
- (void)context;
- CHK(ids);
- ids[0] = box_indices[itri*3+0];
- ids[1] = box_indices[itri*3+1];
- ids[2] = box_indices[itri*3+2];
-}
-
-static void
-box_get_position(const size_t ivert, double pos[3], void* context)
-{
- (void)context;
- CHK(pos);
- pos[0] = box_vertices[ivert*3+0];
- pos[1] = box_vertices[ivert*3+1];
- pos[2] = box_vertices[ivert*3+2];
-}
-
-static void
-box_get_interface(const size_t itri, struct sdis_interface** bound, void* context)
-{
- struct sdis_interface** interfaces = context;
- CHK(context && bound);
- *bound = interfaces[itri];
-}
-
-/*******************************************************************************
- * Geometry 2D
- ******************************************************************************/
-static void
-square_get_indices(const size_t iseg, size_t ids[2], void* context)
-{
- (void)context;
- CHK(ids);
- ids[0] = square_indices[iseg*2+0];
- ids[1] = square_indices[iseg*2+1];
-}
-
-static void
-square_get_position(const size_t ivert, double pos[2], void* context)
-{
- (void)context;
- CHK(pos);
- pos[0] = square_vertices[ivert*2+0];
- pos[1] = square_vertices[ivert*2+1];
-}
-
-static void
-square_get_interface
- (const size_t iseg, struct sdis_interface** bound, void* context)
-{
- struct sdis_interface** interfaces = context;
- CHK(context && bound);
- *bound = interfaces[iseg];
-}
-
-
-/*******************************************************************************
* Media
******************************************************************************/
static double
diff --git a/src/test_sdis_utils.h b/src/test_sdis_utils.h
@@ -58,6 +58,37 @@ 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]);
+static INLINE void
+box_get_indices(const size_t itri, size_t ids[3], void* context)
+{
+ (void)context;
+ CHK(ids);
+ CHK(itri < box_ntriangles);
+ ids[0] = box_indices[itri*3+0];
+ ids[1] = box_indices[itri*3+1];
+ ids[2] = box_indices[itri*3+2];
+}
+
+static INLINE void
+box_get_position(const size_t ivert, double pos[3], void* context)
+{
+ (void)context;
+ CHK(pos);
+ CHK(ivert < box_nvertices);
+ pos[0] = box_vertices[ivert*3+0];
+ pos[1] = box_vertices[ivert*3+1];
+ pos[2] = box_vertices[ivert*3+2];
+}
+
+static INLINE void
+box_get_interface(const size_t itri, struct sdis_interface** bound, void* context)
+{
+ struct sdis_interface** interfaces = context;
+ CHK(context && bound);
+ CHK(itri < box_ntriangles);
+ *bound = interfaces[itri];
+}
+
/*******************************************************************************
* Square geometry
******************************************************************************/
@@ -77,6 +108,36 @@ static const size_t square_indices[4/*#segments*/*2/*#indices per segment*/]= {
};
static const size_t square_nsegments = sizeof(square_indices)/sizeof(size_t[2]);
+static INLINE void
+square_get_indices(const size_t iseg, size_t ids[2], void* context)
+{
+ (void)context;
+ CHK(ids);
+ CHK(iseg < square_nsegments);
+ ids[0] = square_indices[iseg*2+0];
+ ids[1] = square_indices[iseg*2+1];
+}
+
+static INLINE void
+square_get_position(const size_t ivert, double pos[2], void* context)
+{
+ (void)context;
+ CHK(pos);
+ CHK(ivert < square_nvertices);
+ pos[0] = square_vertices[ivert*2+0];
+ pos[1] = square_vertices[ivert*2+1];
+}
+
+static INLINE void
+square_get_interface
+ (const size_t iseg, struct sdis_interface** bound, void* context)
+{
+ struct sdis_interface** interfaces = context;
+ CHK(context && bound);
+ CHK(iseg < square_nsegments);
+ *bound = interfaces[iseg];
+}
+
/*******************************************************************************
* Medium & interface
******************************************************************************/
@@ -122,6 +183,7 @@ static const struct sdis_fluid_shader DUMMY_FLUID_SHADER = {
}
static const struct sdis_interface_shader DUMMY_INTERFACE_SHADER = {
dummy_interface_getter,
+ 0,
DUMMY_INTERFACE_SIDE_SHADER__,
DUMMY_INTERFACE_SIDE_SHADER__
};
diff --git a/src/test_sdis_volumic_power.c b/src/test_sdis_volumic_power.c
@@ -52,67 +52,6 @@
#define DELTA 1.0/20.0
/*******************************************************************************
- * Geometry 3D
- ******************************************************************************/
-static void
-box_get_indices(const size_t itri, size_t ids[3], void* context)
-{
- (void)context;
- CHK(ids);
- ids[0] = box_indices[itri*3+0];
- ids[1] = box_indices[itri*3+1];
- ids[2] = box_indices[itri*3+2];
-}
-
-static void
-box_get_position(const size_t ivert, double pos[3], void* context)
-{
- (void)context;
- CHK(pos);
- pos[0] = box_vertices[ivert*3+0];
- pos[1] = box_vertices[ivert*3+1];
- pos[2] = box_vertices[ivert*3+2];
-}
-
-static void
-box_get_interface(const size_t itri, struct sdis_interface** bound, void* context)
-{
- struct sdis_interface** interfaces = context;
- CHK(context && bound);
- *bound = interfaces[itri];
-}
-
-/*******************************************************************************
- * Geometry 2D
- ******************************************************************************/
-static void
-square_get_indices(const size_t iseg, size_t ids[2], void* context)
-{
- (void)context;
- CHK(ids);
- ids[0] = square_indices[iseg*2+0];
- ids[1] = square_indices[iseg*2+1];
-}
-
-static void
-square_get_position(const size_t ivert, double pos[2], void* context)
-{
- (void)context;
- CHK(pos);
- pos[0] = square_vertices[ivert*2+0];
- pos[1] = square_vertices[ivert*2+1];
-}
-
-static void
-square_get_interface
- (const size_t iseg, struct sdis_interface** bound, void* context)
-{
- struct sdis_interface** interfaces = context;
- CHK(context && bound);
- *bound = interfaces[iseg];
-}
-
-/*******************************************************************************
* Media
******************************************************************************/
struct solid {
@@ -364,3 +303,4 @@ main(int argc, char** argv)
CHK(mem_allocated_size() == 0);
return 0;
}
+
