commit 56a09c32a1e1d4a45476f7d553c3a882ad6e7172
parent 500a3d18b2eca9176573552f48a22fe16c055c20
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Tue, 16 Oct 2018 10:17:11 +0200
Merge branch 'feature_solve_boundary' into develop
Diffstat:
7 files changed, 659 insertions(+), 328 deletions(-)
diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt
@@ -29,9 +29,9 @@ CMAKE_DEPENDENT_OPTION(ALL_TESTS
# Check dependencies
################################################################################
find_package(RCMake 0.4 REQUIRED)
-find_package(Star2D 0.1.1 REQUIRED)
+find_package(Star2D 0.2 REQUIRED)
find_package(Star3D 0.5 REQUIRED)
-find_package(StarSP 0.7 REQUIRED)
+find_package(StarSP 0.8 REQUIRED)
find_package(StarEnc 0.2.0 REQUIRED)
find_package(StarEnc2D 0.2.0 REQUIRED)
find_package(RSys 0.6.1 REQUIRED)
@@ -155,7 +155,7 @@ if(NOT NO_TEST)
new_test(test_sdis_solve_probe_2d)
new_test(test_sdis_solve_probe2_2d)
new_test(test_sdis_solve_probe3_2d)
- new_test(test_sdis_solve_probe_boundary)
+ new_test(test_sdis_solve_boundary)
new_test(test_sdis_volumic_power)
# Additionnal tests
diff --git a/src/sdis.h b/src/sdis.h
@@ -403,7 +403,7 @@ sdis_interface_ref_put
/*******************************************************************************
* A scene is a collection of primitives. Each primitive is the geometric
- * support of the interface between 2 mediums.
+ * support of the interface between 2 media.
******************************************************************************/
/* Create a 3D scene. The geometry of the scene is defined by an indexed
* triangular mesh: each triangle is composed of 3 indices where each index
@@ -582,6 +582,19 @@ sdis_solve_camera
sdis_write_accums_T writer,
void* writer_data);
+SDIS_API res_T
+sdis_solve_boundary
+ (struct sdis_scene* scn,
+ const size_t nrealisations, /* #realisations */
+ const size_t primitives[], /* List of boundary primitives to handle */
+ const enum sdis_side sides[], /* Per primitive side to consider */
+ const size_t nprimitives, /* #primitives */
+ const double time, /* Observation time */
+ const double fp_to_meter, /* Scale from floating point units to meters */
+ const double ambient_radiative_temperature, /* In Kelvin */
+ const double reference_temperature, /* In Kelvin */
+ struct sdis_estimator** estimator);
+
END_DECLS
#endif /* SDIS_H */
diff --git a/src/sdis_scene.c b/src/sdis_scene.c
@@ -231,8 +231,8 @@ sdis_scene_boundary_project_position
/* Retrieve the segment vertices */
S2D(scene_view_get_primitive(scn->s2d_view, (unsigned int)iprim, &prim));
- S2D(primitive_get_attrib(&prim, S2D_POSITION, 0, &a)); d2_set_f2(V[0], a.value);
- S2D(primitive_get_attrib(&prim, S2D_POSITION, 1, &a)); d2_set_f2(V[1], a.value);
+ S2D(segment_get_vertex_attrib(&prim, 0, S2D_POSITION, &a)); d2_set_f2(V[0], a.value);
+ S2D(segment_get_vertex_attrib(&prim, 1, S2D_POSITION, &a)); d2_set_f2(V[1], a.value);
/* Compute the parametric coordinate of the project of `pos' onto the
* segment.*/
diff --git a/src/sdis_solve.c b/src/sdis_solve.c
@@ -245,13 +245,7 @@ sdis_solve_probe
}
if(res != RES_OK) goto error;
- estimator->nrealisations = N;
- estimator->nfailures = nrealisations - N;
- estimator->temperature.E = weight / (double)N;
- estimator->temperature.V =
- sqr_weight / (double)N
- - estimator->temperature.E * estimator->temperature.E;
- estimator->temperature.SE = sqrt(estimator->temperature.V / (double)N);
+ setup_estimator(estimator, nrealisations, N, weight, sqr_weight);
exit:
if(rngs) {
@@ -296,8 +290,8 @@ sdis_solve_probe_boundary
ATOMIC res = RES_OK;
if(!scn || !nrealisations || nrealisations > INT64_MAX || !uv || time < 0
- || fp_to_meter <= 0 || Tref < 0 || (side != SDIS_FRONT && side != SDIS_BACK)
- || !out_estimator) {
+ || fp_to_meter <= 0 || Tref < 0 || (side != SDIS_FRONT && side != SDIS_BACK)
+ || !out_estimator) {
res = RES_BAD_ARG;
goto error;
}
@@ -367,7 +361,7 @@ sdis_solve_probe_boundary
const int ithread = omp_get_thread_num();
struct ssp_rng* rng = rngs[ithread];
- if(ATOMIC_GET(&res) != RES_OK) continue; /* An error occured */
+ if(ATOMIC_GET(&res) != RES_OK) continue; /* An error occurred */
if(scene_is_2d(scn)) {
res_local = boundary_realisation_2d
@@ -389,13 +383,7 @@ sdis_solve_probe_boundary
}
if(res != RES_OK) goto error;
- estimator->nrealisations = N;
- estimator->nfailures = nrealisations - N;
- estimator->temperature.E = weight / (double)N;
- estimator->temperature.V =
- sqr_weight / (double)N
- - estimator->temperature.E * estimator->temperature.E;
- estimator->temperature.SE = sqrt(estimator->temperature.V / (double)N);
+ setup_estimator(estimator, nrealisations, N, weight, sqr_weight);
exit:
if(rngs) {
@@ -553,3 +541,28 @@ error:
goto exit;
}
+res_T
+sdis_solve_boundary
+ (struct sdis_scene* scn,
+ const size_t nrealisations, /* #realisations */
+ const size_t primitives[], /* List of boundary primitives to handle */
+ const enum sdis_side sides[], /* Per primitive side to consider */
+ const size_t nprimitives, /* #primitives */
+ const double time, /* Observation time */
+ const double fp_to_meter, /* Scale from floating point units to meters */
+ const double Tarad, /* In Kelvin */
+ const double Tref, /* In Kelvin */
+ struct sdis_estimator** out_estimator)
+{
+ res_T res = RES_OK;
+ if(!scn) return RES_BAD_ARG;
+ if(scene_is_2d(scn)) {
+ res = solve_boundary_2d(scn, nrealisations, primitives, sides, nprimitives,
+ time, fp_to_meter, Tarad, Tref, out_estimator);
+ } else {
+ res = solve_boundary_3d(scn, nrealisations, primitives, sides, nprimitives,
+ time, fp_to_meter, Tarad, Tref, out_estimator);
+ }
+ return res;
+}
+
diff --git a/src/sdis_solve_Xd.h b/src/sdis_solve_Xd.h
@@ -27,6 +27,7 @@
#include <rsys/stretchy_array.h>
#include <star/ssp.h>
+#include <omp.h>
/* Define a new result code from RES_BAD_OP saying that the bad operation is
* definitive, i.e. in the current state, the realisation will inevitably fail.
@@ -80,6 +81,25 @@ reflect_3d(float res[3], const float V[3], const float N[3])
return res;
}
+static INLINE void
+setup_estimator
+ (struct sdis_estimator* estimator,
+ const size_t nrealisations,
+ const size_t nsuccesses,
+ const double accum_weights,
+ const double accum_sqr_weights)
+{
+ ASSERT(estimator && nrealisations && nsuccesses);
+ estimator->nrealisations = nsuccesses;
+ estimator->nfailures = nrealisations - nsuccesses;
+ estimator->temperature.E = accum_weights / (double)nsuccesses;
+ estimator->temperature.V =
+ accum_sqr_weights / (double)nsuccesses
+ - estimator->temperature.E * estimator->temperature.E;
+ estimator->temperature.V = MMAX(estimator->temperature.V, 0);
+ estimator->temperature.SE = sqrt(estimator->temperature.V / (double)nsuccesses);
+}
+
#endif /* SDIS_SOLVE_XD_H */
#else
@@ -106,7 +126,12 @@ reflect_3d(float res[3], const float V[3], const float N[3])
#define SXD_HIT_NULL__ CONCAT(CONCAT(S, DIM), D_HIT_NULL__)
#define SXD_POSITION CONCAT(CONCAT(S, DIM), D_POSITION)
#define SXD_GEOMETRY_NORMAL CONCAT(CONCAT(S, DIM), D_GEOMETRY_NORMAL)
+#define SXD_VERTEX_DATA_NULL CONCAT(CONCAT(S, DIM), D_VERTEX_DATA_NULL)
#define SXD CONCAT(CONCAT(S, DIM), D)
+#define SXD_FLOAT2 CONCAT(CONCAT(S, DIM), D_FLOAT2)
+#define SXD_FLOAT3 CONCAT(CONCAT(S, DIM), D_FLOAT3)
+#define SXD_SAMPLE CONCAT(CONCAT(S, DIM), D_SAMPLE)
+#define sXd_dev CONCAT(CONCAT(s, DIM), d)
/* Vector macros generic to SDIS_SOLVE_DIMENSION */
#define dX(Func) CONCAT(CONCAT(CONCAT(d, DIM), _), Func)
@@ -128,6 +153,14 @@ static const struct XD(rwalk) XD(RWALK_NULL) = {
SDIS_RWALK_VERTEX_NULL__, NULL, SXD_HIT_NULL__, SDIS_SIDE_NULL__
};
+struct XD(boundary_context) {
+ struct sXd(scene_view)* view;
+ const size_t* primitives;
+};
+static const struct XD(boundary_context) XD(BOUNDARY_CONTEXT_NULL) = {
+ NULL, NULL
+};
+
struct XD(temperature) {
res_T (*func)/* Next function to invoke in order to compute the temperature */
(struct sdis_scene* scn,
@@ -180,6 +213,38 @@ XD(radiative_temperature)
/*******************************************************************************
* Helper functions
******************************************************************************/
+static INLINE void
+XD(boundary_get_indices)(const unsigned iprim, unsigned ids[DIM], void* context)
+{
+ unsigned i;
+ (void)context;
+ ASSERT(ids);
+ FOR_EACH(i, 0, DIM) ids[i] = iprim*DIM + i;
+}
+
+static INLINE void
+XD(boundary_get_position)(const unsigned ivert, float pos[DIM], void* context)
+{
+ struct XD(boundary_context)* ctx = context;
+ struct sXd(primitive) prim;
+ struct sXd(attrib) attr;
+ const unsigned iprim_id = ivert / DIM;
+ const unsigned iprim_vert = ivert % DIM;
+ unsigned iprim;
+ ASSERT(pos && context);
+
+ iprim = (unsigned)ctx->primitives[iprim_id];
+ SXD(scene_view_get_primitive(ctx->view, iprim, &prim));
+#if DIM == 2
+ SXD(segment_get_vertex_attrib(&prim, iprim_vert, SXD_POSITION, &attr));
+ ASSERT(attr.type == SXD_FLOAT2);
+#else
+ SXD(triangle_get_vertex_attrib(&prim, iprim_vert, SXD_POSITION, &attr));
+ ASSERT(attr.type == SXD_FLOAT3);
+#endif
+ fX(set)(pos, attr.value);
+}
+
static FINLINE void
XD(move_pos)(double pos[DIM], const float dir[DIM], const float delta)
{
@@ -1504,14 +1569,186 @@ error:
}
#endif /* SDIS_SOLVE_DIMENSION == 3 */
+static res_T
+XD(solve_boundary)
+ (struct sdis_scene* scn,
+ const size_t nrealisations, /* #realisations */
+ const size_t primitives[], /* List of boundary primitives to handle */
+ const enum sdis_side sides[], /* Per primitive side to consider */
+ const size_t nprimitives, /* #primitives */
+ const double time, /* Observation time */
+ const double fp_to_meter, /* Scale from floating point units to meters */
+ const double Tarad, /* In Kelvin */
+ const double Tref, /* In Kelvin */
+ struct sdis_estimator** out_estimator)
+{
+ struct XD(boundary_context) ctx = XD(BOUNDARY_CONTEXT_NULL);
+ struct sXd(vertex_data) vdata = SXD_VERTEX_DATA_NULL;
+ struct sXd(scene)* scene = NULL;
+ struct sXd(shape)* shape = NULL;
+ struct sXd(scene_view)* view = NULL;
+ struct sdis_estimator* estimator = NULL;
+ struct ssp_rng_proxy* rng_proxy = NULL;
+ struct ssp_rng** rngs = NULL;
+ size_t i;
+ size_t N = 0; /* #realisations that do not fail */
+ size_t view_nprims;
+ double weight=0, sqr_weight=0;
+ int64_t irealisation;
+ ATOMIC res = RES_OK;
+
+ if(!scn || !nrealisations || nrealisations > INT64_MAX || !primitives
+ || !sides || !nprimitives || time < 0 || fp_to_meter < 0 || Tref < 0
+ || !out_estimator) {
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ SXD(scene_view_primitives_count(scn->sXd(view), &view_nprims));
+ FOR_EACH(i, 0, nprimitives) {
+ if(primitives[i] >= view_nprims) {
+ res = RES_BAD_ARG;
+ goto error;
+ }
+ }
+
+ /* Create the Star-XD shape of the boundary */
+#if DIM == 2
+ res = sXd(shape_create_line_segments)(scn->dev->sXd_dev, &shape);
+#else
+ res = sXd(shape_create_mesh)(scn->dev->sXd_dev, &shape);
+#endif
+ if(res != RES_OK) goto error;
+
+ /* Initialise the boundary shape with the triangles/segments of the
+ * submitted primitives */
+ ctx.primitives = primitives;
+ ctx.view = scn->sXd(view);
+ vdata.usage = SXD_POSITION;
+ vdata.type = DIM == 2 ? SXD_FLOAT2 : SXD_FLOAT3;
+ vdata.get = XD(boundary_get_position);
+#if DIM == 2
+ res = sXd(line_segments_setup_indexed_vertices)(shape, (unsigned)nprimitives,
+ XD(boundary_get_indices), (unsigned)(nprimitives*DIM), &vdata, 1, &ctx);
+#else /* DIM == 3 */
+ res = sXd(mesh_setup_indexed_vertices)(shape, (unsigned)nprimitives,
+ XD(boundary_get_indices), (unsigned)(nprimitives*DIM), &vdata, 1, &ctx);
+#endif
+ if(res != RES_OK) goto error;
+
+ /* Create and setup the boundary Star-XD scene */
+ res = sXd(scene_create)(scn->dev->sXd_dev, &scene);
+ if(res != RES_OK) goto error;
+ res = sXd(scene_attach_shape)(scene, shape);
+ if(res != RES_OK) goto error;
+ res = sXd(scene_view_create)(scene, SXD_SAMPLE, &view);
+ if(res != RES_OK) goto error;
+
+ /* Create the proxy RNG */
+ res = ssp_rng_proxy_create(scn->dev->allocator, &ssp_rng_mt19937_64,
+ scn->dev->nthreads, &rng_proxy);
+ if(res != RES_OK) goto error;
+
+ /* Create the per thread RNG */
+ rngs = MEM_CALLOC(scn->dev->allocator, scn->dev->nthreads, sizeof(*rngs));
+ if(!rngs) { res = RES_MEM_ERR; goto error; }
+ FOR_EACH(i, 0, scn->dev->nthreads) {
+ res = ssp_rng_proxy_create_rng(rng_proxy, i, rngs+i);
+ if(res != RES_OK) goto error;
+ }
+
+ /* Create the estimator */
+ res = estimator_create(scn->dev, &estimator);
+ if(res != RES_OK) goto error;
+
+ omp_set_num_threads((int)scn->dev->nthreads);
+ #pragma omp parallel for schedule(static) reduction(+:weight,sqr_weight,N)
+ for(irealisation=0; irealisation<(int64_t)nrealisations; ++irealisation) {
+ const int ithread = omp_get_thread_num();
+ struct sXd(primitive) prim;
+ struct ssp_rng* rng = rngs[ithread];
+ enum sdis_side side;
+ size_t iprim;
+ double w = NaN;
+ double uv[DIM-1];
+ float st[DIM-1];
+ res_T res_local = RES_OK;
+
+ if(ATOMIC_GET(&res) != RES_OK) continue; /* An error occurred */
+
+ /* Sample a position onto the boundary */
+#if DIM == 2
+ res_local = s2d_scene_view_sample
+ (view,
+ ssp_rng_canonical_float(rng),
+ ssp_rng_canonical_float(rng),
+ &prim, st);
+ uv[0] = (double)st[0];
+#else
+ res_local = s3d_scene_view_sample
+ (view,
+ ssp_rng_canonical_float(rng),
+ ssp_rng_canonical_float(rng),
+ ssp_rng_canonical_float(rng),
+ &prim, st);
+ d2_set_f2(uv, st);
+#endif
+ if(res_local != RES_OK) { ATOMIC_SET(&res, res_local); continue; }
+
+ /* Map from boundary scene to sdis scene */
+ ASSERT(prim.prim_id < nprimitives);
+ iprim = primitives[prim.prim_id];
+ side = sides[prim.prim_id];
+
+ /* Invoke the boundary realisation */
+ res_local = XD(boundary_realisation)
+ (scn, rng, iprim, uv, time, side, fp_to_meter, Tarad, Tref, &w);
+
+ /* Update the MC accumulators */
+ if(res_local == RES_OK) {
+ weight += w;
+ sqr_weight += w*w;
+ ++N;
+ } else if(res_local != RES_BAD_OP) {
+ ATOMIC_SET(&res, res_local);
+ continue;
+ }
+ }
+
+ setup_estimator(estimator, nrealisations, N, weight, sqr_weight);
+
+exit:
+ if(scene) SXD(scene_ref_put(scene));
+ if(shape) SXD(shape_ref_put(shape));
+ if(view) SXD(scene_view_ref_put(view));
+ if(rng_proxy) SSP(rng_proxy_ref_put(rng_proxy));
+ if(out_estimator) *out_estimator = estimator;
+ if(rngs) {
+ FOR_EACH(i, 0, scn->dev->nthreads) {if(rngs[i]) SSP(rng_ref_put(rngs[i]));}
+ MEM_RM(scn->dev->allocator, rngs);
+ }
+ return (res_T)res;
+error:
+ if(estimator) {
+ SDIS(estimator_ref_put(estimator));
+ estimator = NULL;
+ }
+ goto exit;
+}
+
#undef SDIS_SOLVE_DIMENSION
#undef DIM
#undef sXd
+#undef sXd_dev
#undef SXD_HIT_NONE
#undef SXD_HIT_NULL
#undef SXD_HIT_NULL__
#undef SXD_POSITION
#undef SXD_GEOMETRY_NORMAL
+#undef SXD_VERTEX_DATA_NULL
+#undef SXD_FLOAT2
+#undef SXD_FLOAT3
+#undef SXD_SAMPLE
#undef SXD
#undef dX
#undef fX
diff --git a/src/test_sdis_solve_boundary.c b/src/test_sdis_solve_boundary.c
@@ -0,0 +1,373 @@
+/* Copyright (C) 2016-2018 |Meso|Star> (contact@meso-star.com)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. */
+
+#include "sdis.h"
+#include "test_sdis_utils.h"
+
+#include <rsys/math.h>
+
+/*
+ * The scene is composed of a solid cube/square whose temperature is unknown.
+ * The convection coefficient with the surrounding fluid is null exepted for
+ * the +X face whose value is 'H'. The Temperature of the -X face is fixed to
+ * Tb. This test computes the temperature on the +X face and check that it is
+ * equal to:
+ *
+ * T = (H*Tf + LAMBDA/A * Tb) / (H+LAMBDA/A)
+ *
+ * with Tf the temperature of the surrounding fluid, lambda the conductivity of
+ * the cube and A the size of the cube/square, i.e. 1.
+ *
+ * 3D 2D
+ *
+ * ///// (1,1,1) ///// (1,1)
+ * +-------+ +-------+
+ * /' /| _\ | | _\
+ * +-------+ | / / Tf Tb | / / Tf
+ * Tb +.....|.+ \__/ | | \__/
+ * |, |/ +-------+
+ * +-------+ (0,0) /////
+ * (0,0,0) /////
+ */
+
+#define UNKNOWN_TEMPERATURE -1
+#define N 10000 /* #realisations */
+
+#define Tf 310.0
+#define Tb 300.0
+#define H 0.5
+#define LAMBDA 0.1
+
+/*******************************************************************************
+ * Media
+ ******************************************************************************/
+static double
+fluid_get_temperature
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ (void)data;
+ CHK(vtx != NULL);
+ return Tf;
+}
+
+static double
+solid_get_calorific_capacity
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ (void)data;
+ CHK(vtx != NULL);
+ return 2.0;
+}
+
+static double
+solid_get_thermal_conductivity
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ (void)data;
+ CHK(vtx != NULL);
+ return LAMBDA;
+}
+
+static double
+solid_get_volumic_mass
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ (void)data;
+ CHK(vtx != NULL);
+ return 25.0;
+}
+
+static double
+solid_get_delta
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ (void)data;
+ CHK(vtx != NULL);
+ return 1.0/20.0;
+}
+
+static double
+solid_get_temperature
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ (void)data;
+ CHK(vtx != NULL);
+ return UNKNOWN_TEMPERATURE;
+}
+
+/*******************************************************************************
+ * Interfaces
+ ******************************************************************************/
+struct interf {
+ double temperature;
+ double 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;
+}
+
+/*******************************************************************************
+ * Helper function
+ ******************************************************************************/
+static void
+check_estimator
+ (const struct sdis_estimator* estimator,
+ const size_t nrealisations, /* #realisations */
+ const double ref)
+{
+ struct sdis_mc T = SDIS_MC_NULL;
+ size_t nreals;
+ size_t nfails;
+ CHK(estimator && nrealisations);
+
+ CHK(sdis_estimator_get_temperature(estimator, &T) == RES_OK);
+ CHK(sdis_estimator_get_realisation_count(estimator, &nreals) == RES_OK);
+ CHK(sdis_estimator_get_failure_count(estimator, &nfails) == RES_OK);
+ printf("%g ~ %g +/- %g\n", ref, T.E, T.SE);
+ printf("#failures = %lu/%lu\n",
+ (unsigned long)nfails, (unsigned long)nrealisations);
+ CHK(nfails + nreals == nrealisations);
+ CHK(nfails < N/1000);
+ CHK(eq_eps(T.E, ref, 3*T.SE));
+}
+
+/*******************************************************************************
+ * Test
+ ******************************************************************************/
+int
+main(int argc, char** argv)
+{
+ struct mem_allocator allocator;
+ struct sdis_data* data = NULL;
+ struct sdis_device* dev = NULL;
+ struct sdis_medium* fluid = NULL;
+ struct sdis_medium* solid = NULL;
+ struct sdis_interface* interf_adiabatic = NULL;
+ struct sdis_interface* interf_Tb = NULL;
+ struct sdis_interface* interf_H = NULL;
+ struct sdis_scene* box_scn = NULL;
+ struct sdis_scene* square_scn = NULL;
+ struct sdis_estimator* estimator = NULL;
+ struct sdis_fluid_shader fluid_shader = DUMMY_FLUID_SHADER;
+ struct sdis_solid_shader solid_shader = DUMMY_SOLID_SHADER;
+ struct sdis_interface_shader interf_shader = SDIS_INTERFACE_SHADER_NULL;
+ struct sdis_interface* box_interfaces[12 /*#triangles*/];
+ struct sdis_interface* square_interfaces[4/*#segments*/];
+ struct interf* interf_props = NULL;
+ double uv[2];
+ double pos[3];
+ double ref;
+ size_t prims[4];
+ enum sdis_side sides[4];
+ size_t iprim;
+ (void)argc, (void)argv;
+
+ CHK(mem_init_proxy_allocator(&allocator, &mem_default_allocator) == RES_OK);
+ CHK(sdis_device_create
+ (NULL, &allocator, SDIS_NTHREADS_DEFAULT, 1, &dev) == RES_OK);
+
+ /* Create the fluid medium */
+ fluid_shader.temperature = fluid_get_temperature;
+ CHK(sdis_fluid_create(dev, &fluid_shader, NULL, &fluid) == RES_OK);
+
+ /* Create the solid_medium */
+ solid_shader.calorific_capacity = solid_get_calorific_capacity;
+ solid_shader.thermal_conductivity = solid_get_thermal_conductivity;
+ solid_shader.volumic_mass = solid_get_volumic_mass;
+ solid_shader.delta_solid = solid_get_delta;
+ solid_shader.temperature = solid_get_temperature;
+ CHK(sdis_solid_create(dev, &solid_shader, NULL, &solid) == RES_OK);
+
+ /* Setup the interface shader */
+ interf_shader.convection_coef = interface_get_convection_coef;
+ interf_shader.front.temperature = interface_get_temperature;
+ interf_shader.front.emissivity = NULL;
+ interf_shader.front.specular_fraction = NULL;
+ interf_shader.back = SDIS_INTERFACE_SIDE_SHADER_NULL;
+
+ /* Create the adiabatic interface */
+ CHK(sdis_data_create(dev, sizeof(struct interf), 16, NULL, &data) == RES_OK);
+ interf_props = sdis_data_get(data);
+ interf_props->hc = 0;
+ interf_props->temperature = UNKNOWN_TEMPERATURE;
+ CHK(sdis_interface_create
+ (dev, solid, fluid, &interf_shader, data, &interf_adiabatic) == RES_OK);
+ CHK(sdis_data_ref_put(data) == RES_OK);
+
+ /* Create the Tb interface */
+ CHK(sdis_data_create(dev, sizeof(struct interf), 16, NULL, &data) == RES_OK);
+ interf_props = sdis_data_get(data);
+ interf_props->hc = 0;
+ interf_props->temperature = Tb;
+ CHK(sdis_interface_create
+ (dev, solid, fluid, &interf_shader, data, &interf_Tb) == RES_OK);
+ CHK(sdis_data_ref_put(data) == RES_OK);
+
+ /* Create the H interface */
+ CHK(sdis_data_create(dev, sizeof(struct interf), 16, NULL, &data) == RES_OK);
+ interf_props = sdis_data_get(data);
+ interf_props->hc = H;
+ interf_props->temperature = UNKNOWN_TEMPERATURE;
+ CHK(sdis_interface_create
+ (dev, solid, fluid, &interf_shader, data, &interf_H) == RES_OK);
+ CHK(sdis_data_ref_put(data) == RES_OK);
+
+ /* Release the media */
+ CHK(sdis_medium_ref_put(solid) == RES_OK);
+ CHK(sdis_medium_ref_put(fluid) == RES_OK);
+
+ /* Map the interfaces to their box triangles */
+ box_interfaces[0] = box_interfaces[1] = interf_adiabatic; /* Front */
+ box_interfaces[2] = box_interfaces[3] = interf_Tb; /* Left */
+ box_interfaces[4] = box_interfaces[5] = interf_adiabatic; /* Back */
+ box_interfaces[6] = box_interfaces[7] = interf_H; /* Right */
+ box_interfaces[8] = box_interfaces[9] = interf_adiabatic; /* Top */
+ box_interfaces[10]= box_interfaces[11]= interf_adiabatic; /* Bottom */
+
+ /* Map the interfaces to their square segments */
+ square_interfaces[0] = interf_adiabatic; /* Bottom */
+ square_interfaces[1] = interf_Tb; /* Lef */
+ square_interfaces[2] = interf_adiabatic; /* Top */
+ square_interfaces[3] = interf_H; /* Right */
+
+ /* Create the box scene */
+ CHK(sdis_scene_create(dev, box_ntriangles, box_get_indices,
+ box_get_interface, box_nvertices, box_get_position, box_interfaces,
+ &box_scn) == RES_OK);
+
+ /* Create the square scene */
+ CHK(sdis_scene_2d_create(dev, square_nsegments, square_get_indices,
+ square_get_interface, square_nvertices, square_get_position,
+ square_interfaces, &square_scn) == RES_OK);
+
+ /* Release the interfaces */
+ CHK(sdis_interface_ref_put(interf_adiabatic) == RES_OK);
+ CHK(sdis_interface_ref_put(interf_Tb) == RES_OK);
+ CHK(sdis_interface_ref_put(interf_H) == RES_OK);
+
+ ref = (H*Tf + LAMBDA * Tb) / (H + LAMBDA);
+
+ #define SOLVE sdis_solve_probe_boundary
+ #define F SDIS_FRONT
+ uv[0] = 0.3;
+ uv[1] = 0.3;
+ iprim = 6;
+
+ CHK(SOLVE(NULL, N, iprim, uv, INF, F, 1.0, 0, 0, &estimator) == RES_BAD_ARG);
+ CHK(SOLVE(box_scn, 0, iprim, uv, INF, F, 1.0, 0, 0, &estimator) == RES_BAD_ARG);
+ CHK(SOLVE(box_scn, N, 12, uv, INF, F, 1.0, 0, 0, &estimator) == RES_BAD_ARG);
+ CHK(SOLVE(box_scn, N, iprim, NULL, INF, F, 1.0, 0, 0, &estimator) == RES_BAD_ARG);
+ CHK(SOLVE(box_scn, N, iprim, uv, -1, F, 1.0, 0, 0, &estimator) == RES_BAD_ARG);
+ CHK(SOLVE(box_scn, N, iprim, uv, INF, -1, 1.0, 0, 0, &estimator) == RES_BAD_ARG);
+ CHK(SOLVE(box_scn, N, iprim, uv, INF, F, 1.0, 0, 0, NULL) == RES_BAD_ARG);
+
+ CHK(SOLVE(box_scn, N, iprim, uv, INF, F, 1.0, 0, 0, &estimator) == RES_OK);
+ CHK(sdis_scene_get_boundary_position(box_scn, iprim, uv, pos) == RES_OK);
+ printf("Boundary temperature of the box at (%g %g %g) = ", SPLIT3(pos));
+ check_estimator(estimator, N, ref);
+ CHK(sdis_estimator_ref_put(estimator) == RES_OK);
+
+ uv[0] = 0.5;
+ iprim = 3;
+ CHK(SOLVE(square_scn, N, iprim, uv, INF, F, 1.0, 0, 0, &estimator) == RES_OK);
+ CHK(sdis_scene_get_boundary_position(square_scn, iprim, uv, pos) == RES_OK);
+ printf("Boundary temperature of the square at (%g %g) = ", SPLIT2(pos));
+ check_estimator(estimator, N, ref);
+ CHK(sdis_estimator_ref_put(estimator) == RES_OK);
+ #undef F
+ #undef SOLVE
+
+ sides[0] = SDIS_FRONT;
+ sides[1] = SDIS_FRONT;
+ sides[2] = SDIS_FRONT;
+ sides[3] = SDIS_FRONT;
+
+ #define SOLVE sdis_solve_boundary
+ prims[0] = 6;
+ prims[1] = 7;
+ CHK(SOLVE(NULL, N, prims, sides, 2, INF, 1.0, 0, 0, &estimator) == RES_BAD_ARG);
+ CHK(SOLVE(box_scn, 0, prims, sides, 2, INF, 1.0, 0, 0, &estimator) == RES_BAD_ARG);
+ CHK(SOLVE(box_scn, N, NULL, sides, 2, INF, 1.0, 0, 0, &estimator) == RES_BAD_ARG);
+ CHK(SOLVE(box_scn, N, prims, NULL, 2, INF, 1.0, 0, 0, &estimator) == RES_BAD_ARG);
+ CHK(SOLVE(box_scn, N, prims, sides, 0, INF, 1.0, 0, 0, &estimator) == RES_BAD_ARG);
+ CHK(SOLVE(box_scn, N, prims, sides, 2, -1, 1.0, 0, 0, &estimator) == RES_BAD_ARG);
+ CHK(SOLVE(box_scn, N, prims, sides, 2, INF, 1.0, 0, 0, NULL) == RES_BAD_ARG);
+
+ /* Average temperature on the right side of the box */
+ CHK(SOLVE(box_scn, N, prims, sides, 2, INF, 1.0, 0, 0, &estimator) == RES_OK);
+ printf("Average temperature of the right side of the box = ");
+ check_estimator(estimator, N, ref);
+ CHK(sdis_estimator_ref_put(estimator) == RES_OK);
+
+ /* Average temperature on the right side of the square */
+ prims[0] = 3;
+ sides[0] = SDIS_FRONT;
+ CHK(SOLVE(square_scn, N, prims, sides, 1, INF, 1.0, 0, 0, &estimator) == RES_OK);
+ printf("Average temperature of the right side of the square = ");
+ check_estimator(estimator, N, ref);
+ CHK(sdis_estimator_ref_put(estimator) == RES_OK);
+
+ /* Check out of bound prims */
+ prims[0] = 12;
+ CHK(SOLVE(box_scn, N, prims, sides, 2, INF, 1.0, 0, 0, &estimator) == RES_BAD_ARG);
+ prims[0] = 4;
+ CHK(SOLVE(square_scn, N, prims, sides, 1, INF, 1.0, 0, 0, &estimator) == RES_BAD_ARG);
+
+ ref = (ref + Tb) / 2;
+
+ /* Average temperature on the left/right side of the box */
+ prims[0] = 2;
+ prims[1] = 3;
+ prims[2] = 6;
+ prims[3] = 7;
+ CHK(SOLVE(box_scn, N, prims, sides, 4, INF, 1.0, 0, 0, &estimator) == RES_OK);
+ printf("Average temperature of the right/left side of the box = ");
+ check_estimator(estimator, N, ref);
+ CHK(sdis_estimator_ref_put(estimator) == RES_OK);
+
+ /* Average temperature on the left/right side of the square */
+ prims[0] = 1;
+ prims[1] = 3;
+ CHK(SOLVE(square_scn, N, prims, sides, 2, INF, 1.0, 0, 0, &estimator) == RES_OK);
+ printf("Average temperature of the right/left side of the square = ");
+ check_estimator(estimator, N, ref);
+ CHK(sdis_estimator_ref_put(estimator) == RES_OK);
+ #undef sdis_solve_boundary
+
+ 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_solve_probe_boundary.c b/src/test_sdis_solve_probe_boundary.c
@@ -1,305 +0,0 @@
-/* Copyright (C) 2016-2018 |Meso|Star> (contact@meso-star.com)
- *
- * This program is free software: you can redistribute it and/or modify
- * it under the terms of the GNU General Public License as published by
- * the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * This program is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU General Public License
- * along with this program. If not, see <http://www.gnu.org/licenses/>. */
-
-#include "sdis.h"
-#include "test_sdis_utils.h"
-
-#include <rsys/math.h>
-
-/*
- * The scene is composed of a solid cube/square whose temperature is unknown.
- * The convection coefficient with the surrounding fluid is null exepted for
- * the +X face whose value is 'H'. The Temperature of the -X face is fixed to
- * Tb. This test computes the temperature on the +X face and check that it is
- * equal to:
- *
- * T = (H*Tf + LAMBDA/A * Tb) / (H+LAMBDA/A)
- *
- * with Tf the temperature of the surrounding fluid, lambda the conductivity of
- * the cube and A the size of the cube/square, i.e. 1.
- *
- * 3D 2D
- *
- * ///// (1,1,1) ///// (1,1)
- * +-------+ +-------+
- * /' /| _\ | | _\
- * +-------+ | / / Tf Tb | / / Tf
- * Tb +.....|.+ \__/ | | \__/
- * |, |/ +-------+
- * +-------+ (0,0) /////
- * (0,0,0) /////
- */
-
-#define UNKNOWN_TEMPERATURE -1
-#define N 10000 /* #realisations */
-
-#define Tf 310
-#define Tb 300
-#define H 0.5
-#define LAMBDA 0.1
-
-/*******************************************************************************
- * Media
- ******************************************************************************/
-static double
-fluid_get_temperature
- (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
-{
- (void)data;
- CHK(vtx != NULL);
- return Tf;
-}
-
-static double
-solid_get_calorific_capacity
- (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
-{
- (void)data;
- CHK(vtx != NULL);
- return 2.0;
-}
-
-static double
-solid_get_thermal_conductivity
- (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
-{
- (void)data;
- CHK(vtx != NULL);
- return LAMBDA;
-}
-
-static double
-solid_get_volumic_mass
- (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
-{
- (void)data;
- CHK(vtx != NULL);
- return 25.0;
-}
-
-static double
-solid_get_delta
- (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
-{
- (void)data;
- CHK(vtx != NULL);
- return 1.0/20.0;
-}
-
-static double
-solid_get_temperature
- (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
-{
- (void)data;
- CHK(vtx != NULL);
- return UNKNOWN_TEMPERATURE;
-}
-
-/*******************************************************************************
- * Interfaces
- ******************************************************************************/
-struct interf {
- double temperature;
- double 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;
-}
-
-/*******************************************************************************
- * Test
- ******************************************************************************/
-int
-main(int argc, char** argv)
-{
- struct mem_allocator allocator;
- struct sdis_mc T = SDIS_MC_NULL;
- struct sdis_data* data = NULL;
- struct sdis_device* dev = NULL;
- struct sdis_medium* fluid = NULL;
- struct sdis_medium* solid = NULL;
- struct sdis_interface* interf_adiabatic = NULL;
- struct sdis_interface* interf_Tb = NULL;
- struct sdis_interface* interf_H = NULL;
- struct sdis_scene* box_scn = NULL;
- struct sdis_scene* square_scn = NULL;
- struct sdis_estimator* estimator = NULL;
- struct sdis_fluid_shader fluid_shader = DUMMY_FLUID_SHADER;
- struct sdis_solid_shader solid_shader = DUMMY_SOLID_SHADER;
- struct sdis_interface_shader interf_shader = SDIS_INTERFACE_SHADER_NULL;
- struct sdis_interface* box_interfaces[12 /*#triangles*/];
- struct sdis_interface* square_interfaces[4/*#segments*/];
- struct interf* interf_props = NULL;
- double uv[2];
- double pos[3];
- double ref;
- size_t iprim;
- size_t nreals;
- size_t nfails;
- (void)argc, (void)argv;
-
- CHK(mem_init_proxy_allocator(&allocator, &mem_default_allocator) == RES_OK);
- CHK(sdis_device_create
- (NULL, &allocator, SDIS_NTHREADS_DEFAULT, 1, &dev) == RES_OK);
-
- /* Create the fluid medium */
- fluid_shader.temperature = fluid_get_temperature;
- CHK(sdis_fluid_create(dev, &fluid_shader, NULL, &fluid) == RES_OK);
-
- /* Create the solid_medium */
- solid_shader.calorific_capacity = solid_get_calorific_capacity;
- solid_shader.thermal_conductivity = solid_get_thermal_conductivity;
- solid_shader.volumic_mass = solid_get_volumic_mass;
- solid_shader.delta_solid = solid_get_delta;
- solid_shader.temperature = solid_get_temperature;
- CHK(sdis_solid_create(dev, &solid_shader, NULL, &solid) == RES_OK);
-
- /* Setup the interface shader */
- interf_shader.convection_coef = interface_get_convection_coef;
- interf_shader.front.temperature = interface_get_temperature;
- interf_shader.front.emissivity = NULL;
- interf_shader.front.specular_fraction = NULL;
- interf_shader.back = SDIS_INTERFACE_SIDE_SHADER_NULL;
-
- /* Create the adiabatic interface */
- CHK(sdis_data_create(dev, sizeof(struct interf), 16, NULL, &data) == RES_OK);
- interf_props = sdis_data_get(data);
- interf_props->hc = 0;
- interf_props->temperature = UNKNOWN_TEMPERATURE;
- CHK(sdis_interface_create
- (dev, solid, fluid, &interf_shader, data, &interf_adiabatic) == RES_OK);
- CHK(sdis_data_ref_put(data) == RES_OK);
-
- /* Create the Tb interface */
- CHK(sdis_data_create(dev, sizeof(struct interf), 16, NULL, &data) == RES_OK);
- interf_props = sdis_data_get(data);
- interf_props->hc = 0;
- interf_props->temperature = Tb;
- CHK(sdis_interface_create
- (dev, solid, fluid, &interf_shader, data, &interf_Tb) == RES_OK);
- CHK(sdis_data_ref_put(data) == RES_OK);
-
- /* Create the H interface */
- CHK(sdis_data_create(dev, sizeof(struct interf), 16, NULL, &data) == RES_OK);
- interf_props = sdis_data_get(data);
- interf_props->hc = H;
- interf_props->temperature = UNKNOWN_TEMPERATURE;
- CHK(sdis_interface_create
- (dev, solid, fluid, &interf_shader, data, &interf_H) == RES_OK);
- CHK(sdis_data_ref_put(data) == RES_OK);
-
- /* Release the media */
- CHK(sdis_medium_ref_put(solid) == RES_OK);
- CHK(sdis_medium_ref_put(fluid) == RES_OK);
-
- /* Map the interfaces to their box triangles */
- box_interfaces[0] = box_interfaces[1] = interf_adiabatic; /* Front */
- box_interfaces[2] = box_interfaces[3] = interf_Tb; /* Left */
- box_interfaces[4] = box_interfaces[5] = interf_adiabatic; /* Back */
- box_interfaces[6] = box_interfaces[7] = interf_H; /* Right */
- box_interfaces[8] = box_interfaces[9] = interf_adiabatic; /* Top */
- box_interfaces[10]= box_interfaces[11]= interf_adiabatic; /* Bottom */
-
- /* Map the interfaces to their square segments */
- square_interfaces[0] = interf_adiabatic; /* Bottom */
- square_interfaces[1] = interf_Tb; /* Lef */
- square_interfaces[2] = interf_adiabatic; /* Top */
- square_interfaces[3] = interf_H; /* Right */
-
- /* Create the box scene */
- CHK(sdis_scene_create(dev, box_ntriangles, box_get_indices,
- box_get_interface, box_nvertices, box_get_position, box_interfaces,
- &box_scn) == RES_OK);
-
- /* Create the square scene */
- CHK(sdis_scene_2d_create(dev, square_nsegments, square_get_indices,
- square_get_interface, square_nvertices, square_get_position,
- square_interfaces, &square_scn) == RES_OK);
-
- /* Release the interfaces */
- CHK(sdis_interface_ref_put(interf_adiabatic) == RES_OK);
- CHK(sdis_interface_ref_put(interf_Tb) == RES_OK);
- CHK(sdis_interface_ref_put(interf_H) == RES_OK);
-
- uv[0] = 0.3;
- uv[1] = 0.3;
- iprim = 6;
-
- #define SOLVE sdis_solve_probe_boundary
- #define F SDIS_FRONT
- CHK(SOLVE(NULL, N, iprim, uv, INF, F, 1.0, 0, 0, &estimator) == RES_BAD_ARG);
- CHK(SOLVE(box_scn, 0, iprim, uv, INF, F, 1.0, 0, 0, &estimator) == RES_BAD_ARG);
- CHK(SOLVE(box_scn, N, 12, uv, INF, F, 1.0, 0, 0, &estimator) == RES_BAD_ARG);
- CHK(SOLVE(box_scn, N, iprim, NULL, INF, F, 1.0, 0, 0, &estimator) == RES_BAD_ARG);
- CHK(SOLVE(box_scn, N, iprim, uv, -1, F, 1.0, 0, 0, &estimator) == RES_BAD_ARG);
- CHK(SOLVE(box_scn, N, iprim, uv, INF, -1, 1.0, 0, 0, &estimator) == RES_BAD_ARG);
- CHK(SOLVE(box_scn, N, iprim, uv, INF, F, 1.0, 0, 0, NULL) == RES_BAD_ARG);
- CHK(SOLVE(box_scn, N, iprim, uv, INF, F, 1.0, 0, 0, &estimator) == RES_OK);
-
- ref = (H*Tf + LAMBDA * Tb) / (H + LAMBDA);
-
- CHK(sdis_estimator_get_realisation_count(estimator, &nreals) == RES_OK);
- CHK(sdis_estimator_get_failure_count(estimator, &nfails) == RES_OK);
- CHK(sdis_estimator_get_temperature(estimator, &T) == RES_OK);
- CHK(sdis_estimator_ref_put(estimator) == RES_OK);
- CHK(sdis_scene_get_boundary_position(box_scn, iprim, uv, pos) == RES_OK);
- printf("Boundary temperature of the box at (%g %g %g) = %g ~ %g +/- %g\n",
- SPLIT3(pos), ref, T.E, T.SE);
- printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N);
- CHK(nfails + nreals == N);
- CHK(nfails < N/1000);
- CHK(eq_eps(T.E, ref, 3*T.SE));
-
- uv[0] = 0.5;
- iprim = 3;
- CHK(SOLVE(square_scn, N, iprim, uv, INF, F, 1.0, 0, 0, &estimator) == RES_OK);
- CHK(sdis_estimator_get_realisation_count(estimator, &nreals) == RES_OK);
- CHK(sdis_estimator_get_failure_count(estimator, &nfails) == RES_OK);
- CHK(sdis_estimator_get_temperature(estimator, &T) == RES_OK);
- CHK(sdis_estimator_ref_put(estimator) == RES_OK);
- CHK(sdis_scene_get_boundary_position(square_scn, iprim, uv, pos) == RES_OK);
- printf("Boundary temperature of the square at (%g %g) = %g ~ %g +/- %g\n",
- SPLIT2(pos), ref, T.E, T.SE);
- printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N);
- CHK(nfails + nreals == N);
- CHK(nfails < N/1000);
- CHK(eq_eps(T.E, ref, 3*T.SE));
- #undef SOLVE
-
- CHK(sdis_scene_ref_put(box_scn) == RES_OK);
- 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;
-}
-
