commit 29f02068eb7805f0e4793c2046ab06cd58f5a99c
parent 6fdf0552279f3b1d8e296dda218ce00f4b1e4226
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Wed, 17 Jan 2018 10:15:28 +0100
Add support of 2D scenes in the sdis_solve_probe function
Diffstat:
7 files changed, 185 insertions(+), 108 deletions(-)
diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt
@@ -26,7 +26,7 @@ option(NO_TEST "Do not build tests" OFF)
find_package(RCMake 0.3 REQUIRED)
find_package(Star2D 0.1 REQUIRED)
find_package(Star3D 0.4 REQUIRED)
-find_package(StarSP 0.5 REQUIRED)
+find_package(StarSP 0.7 REQUIRED)
find_package(RSys 0.6 REQUIRED)
find_package(OpenMP 1.2 REQUIRED)
diff --git a/src/sdis_interface.c b/src/sdis_interface.c
@@ -21,6 +21,7 @@
#include <rsys/double3.h>
#include <rsys/mem_allocator.h>
+#include <star/s2d.h>
#include <star/s3d.h>
/*******************************************************************************
@@ -176,7 +177,22 @@ interface_get_id(const struct sdis_interface* interf)
}
void
-setup_interface_fragment
+setup_interface_fragment_2d
+ (struct sdis_interface_fragment* frag,
+ const struct sdis_rwalk_vertex* vertex,
+ const struct s2d_hit* hit)
+{
+ ASSERT(frag && vertex && hit && !S2D_HIT_NONE(hit));
+ d2_set(frag->P, vertex->P);
+ frag->P[2] = 0;
+ d2_normalize(frag->Ng, d2_set_f2(frag->Ng, hit->normal));
+ frag->Ng[2] = 0;
+ frag->uv[0] = hit->u;
+ frag->time = vertex->time;
+}
+
+void
+setup_interface_fragment_3d
(struct sdis_interface_fragment* frag,
const struct sdis_rwalk_vertex* vertex,
const struct s3d_hit* hit)
diff --git a/src/sdis_interface_c.h b/src/sdis_interface_c.h
@@ -22,6 +22,7 @@
#include <float.h>
/* Forward declaration of external type */
+struct s2d_hit;
struct s3d_hit;
struct sdis_interface {
@@ -45,7 +46,13 @@ interface_get_id
(const struct sdis_interface* interf);
extern LOCAL_SYM void
-setup_interface_fragment
+setup_interface_fragment_2d
+ (struct sdis_interface_fragment* frag,
+ const struct sdis_rwalk_vertex* vertex,
+ const struct s2d_hit* hit);
+
+extern LOCAL_SYM void
+setup_interface_fragment_3d
(struct sdis_interface_fragment* frag,
const struct sdis_rwalk_vertex* vertex,
const struct s3d_hit* hit);
diff --git a/src/sdis_scene.c b/src/sdis_scene.c
@@ -39,6 +39,22 @@ struct geometry_context {
/*******************************************************************************
* Helper function
******************************************************************************/
+/* Check that `hit' roughly lies on a vertex. For segments, a simple but
+ * approximative way is to test that its position have at least one barycentric
+ * coordinate roughly equal to 0 or 1. */
+static FINLINE int
+hit_on_vertex(const struct s2d_hit* hit)
+{
+ const float on_vertex_eps = 1.e-4f;
+ float v;
+ ASSERT(hit && !S2D_HIT_NONE(hit));
+ v = 1.f - hit->u;
+ return eq_epsf(hit->u, 0.f, on_vertex_eps)
+ || eq_epsf(hit->u, 1.f, on_vertex_eps)
+ || eq_epsf(v, 0.f, on_vertex_eps)
+ || eq_epsf(v, 1.f, on_vertex_eps);
+}
+
/* Check that `hit' roughly lies on an edge. For triangular primitives, a
* simple but approximative way is to test that its position have at least one
* barycentric coordinate roughly equal to 0 or 1. */
@@ -58,7 +74,7 @@ hit_on_edge(const struct s3d_hit* hit)
}
static int
-hit_filter_function
+hit_filter_function_3d
(const struct s3d_hit* hit,
const float org[3],
const float dir[3],
@@ -81,22 +97,6 @@ hit_filter_function
return 0;
}
-/* 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);
-}
-
static int
hit_filter_function_2d
(const struct s2d_hit* hit,
@@ -122,49 +122,49 @@ hit_filter_function_2d
}
static void
-get_indices(const unsigned itri, unsigned out_ids[3], void* data)
+get_indices_2d(const unsigned iseg, unsigned out_ids[2], void* data)
{
struct geometry_context* ctx = data;
- size_t ids[3];
+ size_t ids[2];
ASSERT(ctx);
- ctx->indices(itri, ids, ctx->data);
+ ctx->indices(iseg, ids, ctx->data);
out_ids[0] = (unsigned)ids[0];
out_ids[1] = (unsigned)ids[1];
- out_ids[2] = (unsigned)ids[2];
}
static void
-get_indices_2d(const unsigned iseg, unsigned out_ids[2], void* data)
+get_indices_3d(const unsigned itri, unsigned out_ids[3], void* data)
{
struct geometry_context* ctx = data;
- size_t ids[2];
+ size_t ids[3];
ASSERT(ctx);
- ctx->indices(iseg, ids, ctx->data);
+ 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(const unsigned ivert, float out_pos[3], void* data)
+get_position_2d(const unsigned ivert, float out_pos[2], void* data)
{
struct geometry_context* ctx = data;
- double pos[3];
+ double pos[2];
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
-get_position_2d(const unsigned ivert, float out_pos[2], void* data)
+get_position_3d(const unsigned ivert, float out_pos[3], void* data)
{
struct geometry_context* ctx = data;
- double pos[2];
+ 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
@@ -229,20 +229,20 @@ error:
}
static res_T
-setup_geometry
+setup_geometry_2d
(struct sdis_scene* scn,
- const size_t ntris, /* #triangles */
- void (*indices)(const size_t itri, size_t ids[3], void*),
+ 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[3], void* ctx),
+ void (*position)(const size_t ivert, double pos[2], 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;
+ 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 && ntris && indices && nverts && position);
+ ASSERT(scn && nsegs && indices && nverts && position);
/* Setup the intermediary geometry context */
context.indices = indices;
@@ -250,50 +250,51 @@ setup_geometry
context.data = ctx;
/* Setup the vertex data */
- vdata.usage = S3D_POSITION;
- vdata.type = S3D_FLOAT3;
- vdata.get = get_position;
+ vdata.usage = S2D_POSITION;
+ vdata.type = S2D_FLOAT2;
+ vdata.get = get_position_2d;
- /* Create the Star-3D geometry */
- res = s3d_scene_create(scn->dev->s3d, &s3d_scn);
+ /* Create the Star-2D geometry */
+ res = s2d_scene_create(scn->dev->s2d, &s2d_scn);
if(res != RES_OK) goto error;
- res = s3d_shape_create_mesh(scn->dev->s3d, &s3d_msh);
+ res = s2d_shape_create_line_segments(scn->dev->s2d, &s2d_msh);
if(res != RES_OK) goto error;
- res = s3d_mesh_set_hit_filter_function(s3d_msh, hit_filter_function, NULL);
+ res = s2d_line_segments_set_hit_filter_function
+ (s2d_msh, hit_filter_function_2d, NULL);
if(res != RES_OK) goto error;
- res = s3d_scene_attach_shape(s3d_scn, s3d_msh);
+ res = s2d_scene_attach_shape(s2d_scn, s2d_msh);
if(res != RES_OK) goto error;
- res = s3d_mesh_setup_indexed_vertices(s3d_msh, (unsigned)ntris, get_indices,
- (unsigned)nverts, &vdata, 1, &context);
+ res = 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 = s3d_scene_view_create(s3d_scn, S3D_SAMPLE|S3D_TRACE|S3D_GET_PRIMITIVE,
- &scn->s3d_view);
+ 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(s3d_msh) S3D(shape_ref_put(s3d_msh));
- if(s3d_scn) S3D(scene_ref_put(s3d_scn));
+ if(s2d_msh) S2D(shape_ref_put(s2d_msh));
+ if(s2d_scn) S2D(scene_ref_put(s2d_scn));
return res;
error:
- if(scn->s3d_view) S3D(scene_view_ref_put(scn->s3d_view));
+ if(scn->s2d_view) S2D(scene_view_ref_put(scn->s2d_view));
goto exit;
}
static res_T
-setup_geometry_2d
+setup_geometry_3d
(struct sdis_scene* scn,
- const size_t nsegs, /* #segments */
- void (*indices)(const size_t itri, size_t ids[2], void*),
+ 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[2], void* ctx),
+ void (*position)(const size_t ivert, double pos[3], 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;
+ 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 && nsegs && indices && nverts && position);
+ ASSERT(scn && ntris && indices && nverts && position);
/* Setup the intermediary geometry context */
context.indices = indices;
@@ -301,33 +302,32 @@ setup_geometry_2d
context.data = ctx;
/* Setup the vertex data */
- vdata.usage = S2D_POSITION;
- vdata.type = S2D_FLOAT2;
- vdata.get = get_position_2d;
+ vdata.usage = S3D_POSITION;
+ vdata.type = S3D_FLOAT3;
+ vdata.get = get_position_3d;
- /* Create the Star-2D geometry */
- res = s2d_scene_create(scn->dev->s2d, &s2d_scn);
+ /* Create the Star-3D geometry */
+ res = s3d_scene_create(scn->dev->s3d, &s3d_scn);
if(res != RES_OK) goto error;
- res = s2d_shape_create_line_segments(scn->dev->s2d, &s2d_msh);
+ res = s3d_shape_create_mesh(scn->dev->s3d, &s3d_msh);
if(res != RES_OK) goto error;
- res = s2d_line_segments_set_hit_filter_function
- (s2d_msh, hit_filter_function_2d, NULL);
+ res = s3d_mesh_set_hit_filter_function(s3d_msh, hit_filter_function_3d, NULL);
if(res != RES_OK) goto error;
- res = s2d_scene_attach_shape(s2d_scn, s2d_msh);
+ res = s3d_scene_attach_shape(s3d_scn, s3d_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);
+ 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 = s2d_scene_view_create(s2d_scn, S2D_SAMPLE|S2D_TRACE|S2D_GET_PRIMITIVE,
- &scn->s2d_view);
+ 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(s2d_msh) S2D(shape_ref_put(s2d_msh));
- if(s2d_scn) S2D(scene_ref_put(s2d_scn));
+ if(s3d_msh) S3D(shape_ref_put(s3d_msh));
+ if(s3d_scn) S3D(scene_ref_put(s3d_scn));
return res;
error:
- if(scn->s2d_view) S2D(scene_view_ref_put(scn->s2d_view));
+ if(scn->s3d_view) S3D(scene_view_ref_put(scn->s3d_view));
goto exit;
}
@@ -373,7 +373,7 @@ scene_create
if(is_2d) {
res = setup_geometry_2d(scn, nprims, indices, nverts, position, ctx);
} else {
- res = setup_geometry(scn, nprims, indices, nverts, position, ctx);
+ 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);
@@ -623,7 +623,6 @@ scene_get_medium
const double pos[],
const struct sdis_medium** out_medium)
{
-
return scene_is_2d(scn)
? scene_get_medium_2d(scn, pos, out_medium)
: scene_get_medium_3d(scn, pos, out_medium);
diff --git a/src/sdis_solve_probe.c b/src/sdis_solve_probe.c
@@ -18,6 +18,10 @@
#include "sdis_estimator_c.h"
#include "sdis_solve_probe_Xd.h"
+/* Generate the 2D solver */
+#define SDIS_SOLVE_PROBE_DIMENSION 2
+#include "sdis_solve_probe_Xd.h"
+
/* Generate the 3D solver */
#define SDIS_SOLVE_PROBE_DIMENSION 3
#include "sdis_solve_probe_Xd.h"
@@ -47,12 +51,6 @@ sdis_solve_probe
goto error;
}
- if(scene_is_2d(scn)) {
- log_err(scn->dev, "%s: 2D scene are not supported yet.\n", FUNC_NAME);
- res = RES_BAD_ARG;
- goto error;
- }
-
res = scene_get_medium(scn, position, &medium);
if(res != RES_OK) goto error;
@@ -62,23 +60,15 @@ sdis_solve_probe
if(res != RES_OK) goto error;
FOR_EACH(irealisation, 0, nrealisations) {
- struct rwalk_3d rwalk = RWALK_NULL_3d;
- struct temperature_3d T = TEMPERATURE_NULL_3d;
+ double w;
- switch(medium->type) {
- case SDIS_MEDIUM_FLUID: T.func = fluid_temperature_3d; break;
- case SDIS_MEDIUM_SOLID: T.func = solid_temperature_3d; break;
- default: FATAL("Unreachable code\n"); break;
+ if(scene_is_2d(scn)) {
+ res = probe_realisation_2d
+ (scn, rng, medium, position, time, fp_to_meter, &w);
+ } else {
+ res = probe_realisation_3d
+ (scn, rng, medium, position, time, fp_to_meter, &w);
}
-
- rwalk.vtx.P[0] = position[0];
- rwalk.vtx.P[1] = position[1];
- rwalk.vtx.P[2] = position[2];
- rwalk.vtx.time = time;
- rwalk.hit = S3D_HIT_NULL;
- rwalk.mdm = medium;
-
- res = compute_temperature_3d(scn, fp_to_meter, &rwalk, rng, &T);
if(res != RES_OK) {
if(res == RES_BAD_OP) {
++estimator->nfailures;
@@ -86,8 +76,8 @@ sdis_solve_probe
goto error;
}
} else {
- weight += T.value;
- sqr_weight += T.value*T.value;
+ weight += w;
+ sqr_weight += w*w;
++estimator->nrealisations;
}
}
diff --git a/src/sdis_solve_probe_Xd.h b/src/sdis_solve_probe_Xd.h
@@ -46,16 +46,23 @@
#error "Invalid dimension "STR(SDIS_SOLVE_PROBE_DIMENSION)
#endif
+/* Syntactic sugar */
#define DIM SDIS_SOLVE_PROBE_DIMENSION
+
+/* Star-XD macros generic to SDIS_SOLVE_PROBE_DIMENSION */
#define sXd(Name) CONCAT(CONCAT(CONCAT(s, DIM), d_), Name)
#define SXD_HIT_NONE CONCAT(CONCAT(S,DIM), D_HIT_NONE)
#define SXD_HIT_NULL CONCAT(CONCAT(S,DIM), D_HIT_NULL)
#define SXD_HIT_NULL__ CONCAT(CONCAT(S, DIM), D_HIT_NULL__)
#define SXD CONCAT(CONCAT(S, DIM), D)
+
+/* Vector macros generic to SDIS_SOLVE_PROBE_DIMENSION */
#define dX(Func) CONCAT(CONCAT(CONCAT(d, DIM), _), Func)
#define fX(Func) CONCAT(CONCAT(CONCAT(f, DIM), _), Func)
#define fX_set_dX CONCAT(CONCAT(CONCAT(f, DIM), _set_d), DIM)
#define dX_set_fX CONCAT(CONCAT(CONCAT(d, DIM), _set_f), DIM)
+
+/* Macro making generic its subimitted name to SDIS_SOLVE_PROBE_DIMENSION */
#define XD(Name) CONCAT(CONCAT(CONCAT(Name, _), DIM), d)
/* Current state of the random walk */
@@ -322,7 +329,7 @@ XD(boundary_temperature)
ASSERT(scn && fp_to_meter > 0 && rwalk && rng && T);
ASSERT(!SXD_HIT_NONE(&rwalk->hit));
- setup_interface_fragment(&frag, &rwalk->vtx, &rwalk->hit);
+ XD(setup_interface_fragment)(&frag, &rwalk->vtx, &rwalk->hit);
/* Retrieve the current interface */
interf = scene_get_interface(scn, rwalk->hit.prim.prim_id);
@@ -396,14 +403,19 @@ XD(solid_temperature)
rho = solid_get_volumic_mass(mdm, &rwalk->vtx);
cp = solid_get_calorific_capacity(mdm, &rwalk->vtx);
+#if (SDIS_SOLVE_PROBE_DIMENSION == 2)
+ /* Sample a direction around 2PI */
+ ssp_ran_circle_uniform_float(rng, dir0, NULL);
+#else
/* Sample a direction around 4PI */
ssp_ran_sphere_uniform_float(rng, dir0, NULL);
+#endif
/* Trace a ray along the sampled direction and its opposite to check if a
* surface is hit in [0, delta_solid]. */
fX_set_dX(org, rwalk->vtx.P);
fX(minus)(dir1, dir0);
- hit0 = hit1 = S3D_HIT_NULL;
+ hit0 = hit1 = SXD_HIT_NULL;
range[0] = 0.f, range[1] = delta_solid*RAY_RANGE_MAX_SCALE;
SXD(scene_view_trace_ray(scn->sXd(view), org, dir0, range, NULL, &hit0));
SXD(scene_view_trace_ray(scn->sXd(view), org, dir1, range, NULL, &hit1));
@@ -417,7 +429,7 @@ XD(solid_temperature)
}
/* Sample the time */
- mu = (2 * DIM * lambda) / (rho * cp * delta * fp_to_meter * delta * fp_to_meter);
+ mu = (2*DIM*lambda) / (rho*cp*delta*fp_to_meter*delta*fp_to_meter);
tau = ssp_ran_exp(rng, mu);
rwalk->vtx.time -= tau;
@@ -503,6 +515,40 @@ error:
goto exit;
}
+static res_T
+XD(probe_realisation)
+ (struct sdis_scene* scn,
+ struct ssp_rng* rng,
+ const struct sdis_medium* medium,
+ const double position[],
+ const double time,
+ const double fp_to_meter,/* Scale factor from floating point unit to meter */
+ double* weight)
+{
+ struct XD(rwalk) rwalk = XD(RWALK_NULL);
+ struct XD(temperature) T = XD(TEMPERATURE_NULL);
+ res_T res = RES_OK;
+ ASSERT(medium && position && fp_to_meter > 0 && weight && time >= 0);
+
+ switch(medium->type) {
+ case SDIS_MEDIUM_FLUID: T.func = XD(fluid_temperature); break;
+ case SDIS_MEDIUM_SOLID: T.func = XD(solid_temperature); break;
+ default: FATAL("Unreachable code\n"); break;
+ }
+
+ dX(set)(rwalk.vtx.P, position);
+ rwalk.vtx.time = time;
+ rwalk.hit = SXD_HIT_NULL;
+ rwalk.mdm = medium;
+
+ res = XD(compute_temperature)(scn, fp_to_meter, &rwalk, rng, &T);
+ if(res != RES_OK) return res;
+
+ *weight = T.value;
+ return RES_OK;
+}
+
+
#undef SDIS_SOLVE_PROBE_DIMENSION
#undef DIM
#undef sXd
diff --git a/src/test_sdis_utils.h b/src/test_sdis_utils.h
@@ -20,7 +20,7 @@
#include <stdio.h>
/*******************************************************************************
- * Geometry
+ * Box geometry
******************************************************************************/
static const double box_vertices[8/*#vertices*/*3/*#coords per vertex*/] = {
0.0, 0.0, 0.0,
@@ -55,6 +55,25 @@ static const size_t box_indices[12/*#triangles*/*3/*#indices per triangle*/] = {
static const size_t box_ntriangles = sizeof(box_indices) / sizeof(size_t[3]);
/*******************************************************************************
+ * Square geometry
+ ******************************************************************************/
+static const double square_vertices[4/*#vertices*/*2/*#coords per vertex*/] = {
+ 1.0, 0.0,
+ 0.0, 0.0,
+ 0.0, 1.0,
+ 1.0, 1.0
+};
+static const size_t square_nvertices = sizeof(square_vertices)/sizeof(double[2]);
+
+static const size_t square_indices[4/*#triangles*/*2/*#indices per segment*/]= {
+ 0, 1, /* Bottom */
+ 1, 2, /* Left */
+ 2, 3, /* Top */
+ 3, 0 /* Right */
+};
+static const size_t square_nsegments = sizeof(square_indices)/sizeof(size_t[2]);
+
+/*******************************************************************************
* Medium & interface
******************************************************************************/
static INLINE double
