commit 1e166a8821fc4f8eb864c0bfd415b690025ab123 parent ef013176835b8e33141e87b7cb4c437454629824 Author: Vincent Forest <vincent.forest@meso-star.com> Date: Fri, 13 Apr 2018 12:11:51 +0200 Merge branch 'release_0.2' Diffstat:
38 files changed, 1122 insertions(+), 49 deletions(-)
diff --git a/README.md b/README.md @@ -23,13 +23,29 @@ variable the install directories of its dependencies. ## Release notes +### Version 0.2 + +- Add the support of volumic power to solid media: add the `volumic_power` + functor to the `sdis_solid_shader` data structure that, once defined, should + return the volumic power of the solid at a specific position and time. On + solve invocation, the conductive random walks take into account this + spatio-temporal volumic power in the computation of the solid temperature. +- Add the `sdis_solve_probe_boundary` function: it computes the temperature at + a given position and time onto a geometric primitive. The probe position is + defined by the index of the primitive and a parametric coordinates onto it. +- Add the `sdis_scene_get_boundary_position` function: it computes a world + space position from the index of a geometric primitive and a parametric + coordinate onto it. +- Fix how the `sdis_solve_probe` was parallelised. The submitted `threads_hint` + parameter was not correctly handled. + ### Version 0.1 - Add the support of radiative temperature. -- Add the `sdis_camera` API : it defines a pinhole camera into the scene. -- Add the `sdis_accum_buffer` API : it is a pool of MC accumulators, i.e. a sum +- Add the `sdis_camera` API: it defines a pinhole camera into the scene. +- Add the `sdis_accum_buffer` API: it is a pool of MC accumulators, i.e. a sum of MC weights and square weights. -- Add the `sdis_solve_camera` function : it relies on a `sdis_camera` and a +- Add the `sdis_solve_camera` function: it relies on a `sdis_camera` and a `sdis_accum_buffer` to compute the radiative temperature that reaches each pixel of an image whose definition is defined by the caller. Note that actually this function uses the same underlying MC algorithm behind the @@ -46,7 +62,7 @@ First version and implementation of the Stardis solver API. ## License -Stardis is Copyright (C) |Meso|Star> 2016-2018 (<contact@meso-star.com>). It is +Stardis is Copyright (C) 2016-2018 |Meso|Star> (<contact@meso-star.com>). It is free software released under the GPLv3+ license. You are welcome to redistribute it under certain conditions; refer to the COPYING files for details. diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt @@ -1,4 +1,4 @@ -# Copyright (C) |Meso|Star> 2016-2018 +# Copyright (C) 2016-2018 |Meso|Star> # # 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 @@ -40,7 +40,7 @@ rcmake_append_runtime_dirs(_runtime_dirs RSys Star3D StarSP) # Configure and define targets ################################################################################ set(VERSION_MAJOR 0) -set(VERSION_MINOR 1) +set(VERSION_MINOR 2) set(VERSION_PATCH 0) set(VERSION ${VERSION_MAJOR}.${VERSION_MINOR}.${VERSION_PATCH}) @@ -134,6 +134,8 @@ 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_volumic_power) target_link_libraries(test_sdis_solve_probe3 Star3DUT) target_link_libraries(test_sdis_solve_camera Star3DUT) diff --git a/src/sdis.h b/src/sdis.h @@ -1,4 +1,4 @@ -/* Copyright (C) |Meso|Star> 2016-2018 (contact@meso-star.com) +/* 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 @@ -129,11 +129,13 @@ struct sdis_solid_shader { sdis_medium_getter_T delta_solid; sdis_medium_getter_T delta_boundary; + sdis_medium_getter_T volumic_power; /* May be NULL <=> no volumic power */ + /* Initial/limit condition. A temperature < 0 means that the temperature is * unknown for the submitted random walk vertex. */ sdis_medium_getter_T temperature; }; -#define SDIS_SOLID_SHADER_NULL__ {NULL, NULL, NULL, NULL, NULL, NULL} +#define SDIS_SOLID_SHADER_NULL__ {NULL, NULL, NULL, NULL, NULL, NULL, NULL} static const struct sdis_solid_shader SDIS_SOLID_SHADER_NULL = SDIS_SOLID_SHADER_NULL__; @@ -398,6 +400,13 @@ sdis_scene_get_aabb double lower[3], double upper[3]); +SDIS_API res_T +sdis_scene_get_boundary_position + (const struct sdis_scene* scn, + const size_t iprim, /* Primitive index */ + const double uv[2], /* Parametric coordinate onto the pimitive */ + double pos[3]); /* World space position */ + /******************************************************************************* * An estimator stores the state of a simulation ******************************************************************************/ @@ -439,6 +448,18 @@ sdis_solve_probe struct sdis_estimator** estimator); SDIS_API res_T +sdis_solve_probe_boundary + (struct sdis_scene* scn, + const size_t nrealisations, /* #realisations */ + const size_t iprim, /* Identifier of the primitive on which the probe lies */ + const double uv[2], /* Parametric coordinates of the probe onto the primitve */ + const double time, /* Observation time */ + const double fp_to_meter, /* Scale from floating point units to meters */ + const double ambient_radiative_temperature, /* In Kelvin */ + const double reference_temperature, /* In Kelvin */ + struct sdis_estimator** estimator); + +SDIS_API res_T sdis_solve_camera (struct sdis_scene* scn, const struct sdis_camera* cam, /* Point of view */ diff --git a/src/sdis_accum_buffer.c b/src/sdis_accum_buffer.c @@ -1,4 +1,4 @@ -/* Copyright (C) |Meso|Star> 2016-2018 (contact@meso-star.com) +/* 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 diff --git a/src/sdis_camera.c b/src/sdis_camera.c @@ -1,4 +1,4 @@ -/* Copyright (C) |Meso|Star> 2016-2018 (contact@meso-star.com) +/* 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 diff --git a/src/sdis_camera.h b/src/sdis_camera.h @@ -1,4 +1,4 @@ -/* Copyright (C) |Meso|Star> 2016-2018 (contact@meso-star.com) +/* 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 diff --git a/src/sdis_data.c b/src/sdis_data.c @@ -1,4 +1,4 @@ -/* Copyright (C) |Meso|Star> 2016-2018 (contact@meso-star.com) +/* 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 diff --git a/src/sdis_device.c b/src/sdis_device.c @@ -1,4 +1,4 @@ -/* Copyright (C) |Meso|Star> 2016-2018 (contact@meso-star.com) +/* 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 diff --git a/src/sdis_device_c.h b/src/sdis_device_c.h @@ -1,4 +1,4 @@ -/* Copyright (C) |Meso|Star> 2016-2018 (contact@meso-star.com) +/* 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 diff --git a/src/sdis_estimator.c b/src/sdis_estimator.c @@ -1,4 +1,4 @@ -/* Copyright (C) |Meso|Star> 2016-2018 (contact@meso-star.com) +/* 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 diff --git a/src/sdis_estimator_c.h b/src/sdis_estimator_c.h @@ -1,4 +1,4 @@ -/* Copyright (C) |Meso|Star> 2016-2018 (contact@meso-star.com) +/* 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 diff --git a/src/sdis_interface.c b/src/sdis_interface.c @@ -1,4 +1,4 @@ -/* Copyright (C) |Meso|Star> 2016-2018 (contact@meso-star.com) +/* 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 diff --git a/src/sdis_interface_c.h b/src/sdis_interface_c.h @@ -1,4 +1,4 @@ -/* Copyright (C) |Meso|Star> 2016-2018 (contact@meso-star.com) +/* 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 diff --git a/src/sdis_medium.c b/src/sdis_medium.c @@ -1,4 +1,4 @@ -/* Copyright (C) |Meso|Star> 2016-2018 (contact@meso-star.com) +/* 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 diff --git a/src/sdis_medium_c.h b/src/sdis_medium_c.h @@ -1,4 +1,4 @@ -/* Copyright (C) |Meso|Star> 2016-2018 (contact@meso-star.com) +/* 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 @@ -102,6 +102,16 @@ solid_get_delta_boundary } static INLINE double +solid_get_volumic_power + (const struct sdis_medium* mdm, const struct sdis_rwalk_vertex* vtx) +{ + ASSERT(mdm && mdm->type == SDIS_MEDIUM_SOLID); + return mdm->shader.solid.volumic_power + ? mdm->shader.solid.volumic_power(vtx, mdm->data) + : 0; +} + +static INLINE double solid_get_temperature (const struct sdis_medium* mdm, const struct sdis_rwalk_vertex* vtx) { diff --git a/src/sdis_scene.c b/src/sdis_scene.c @@ -1,4 +1,4 @@ -/* Copyright (C) |Meso|Star> 2016-2018 (contact@meso-star.com) +/* 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 @@ -608,6 +608,37 @@ sdis_scene_get_aabb return RES_OK; } +res_T +sdis_scene_get_boundary_position + (const struct sdis_scene* scn, + const size_t iprim, + const double uv[], + double pos[]) +{ + if(!scn || !uv || !pos) return RES_BAD_ARG; + if(iprim >= scene_get_primitives_count(scn)) return RES_BAD_ARG; + + if(scene_is_2d(scn)) { + struct s2d_primitive prim; + struct s2d_attrib attr; + float s = (float)uv[0]; + + S2D(scene_view_get_primitive(scn->s2d_view, (unsigned int)iprim, &prim)); + S2D(primitive_get_attrib(&prim, S2D_POSITION, s, &attr)); + d2_set_f2(pos, attr.value); + } else { + struct s3d_primitive prim; + struct s3d_attrib attr; + float st[2]; + + f2_set_d2(st, uv); + S3D(scene_view_get_primitive(scn->s3d_view, (unsigned int)iprim, &prim)); + S3D(primitive_get_attrib(&prim, S3D_POSITION, st, &attr)); + d3_set_f3(pos, attr.value); + } + return RES_OK; +} + /******************************************************************************* * Local miscellaneous function ******************************************************************************/ diff --git a/src/sdis_scene_c.h b/src/sdis_scene_c.h @@ -1,4 +1,4 @@ -/* Copyright (C) |Meso|Star> 2016-2018 (contact@meso-star.com) +/* 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 @@ -46,6 +46,13 @@ struct sdis_scene { struct sdis_device* dev; }; +static FINLINE size_t +scene_get_primitives_count(const struct sdis_scene* scn) +{ + ASSERT(scn); + return darray_interf_size_get(&scn->prim_interfaces); +} + extern LOCAL_SYM const struct sdis_interface* scene_get_interface (const struct sdis_scene* scene, diff --git a/src/sdis_solve.c b/src/sdis_solve.c @@ -1,4 +1,4 @@ -/* Copyright (C) |Meso|Star> 2016-2018 (contact@meso-star.com) +/* 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 @@ -211,6 +211,7 @@ sdis_solve_probe if(res != RES_OK) goto error; /* Here we go! Launch the Monte Carlo estimation */ + omp_set_num_threads((int)scn->dev->nthreads); #pragma omp parallel for schedule(static) reduction(+:weight,sqr_weight,N) for(irealisation = 0; irealisation < nrealisations; ++irealisation) { res_T res_local; @@ -266,6 +267,146 @@ error: } res_T +sdis_solve_probe_boundary + (struct sdis_scene* scn, + const size_t nrealisations, /* #realisations */ + const size_t iprim, /* Identifier of the primitive on which the probe lies */ + const double uv[2], /* Parametric coordinates of the probe onto the primitve */ + 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 sdis_estimator* estimator = NULL; + struct ssp_rng_proxy* rng_proxy = NULL; + struct ssp_rng** rngs = NULL; + double weight = 0; + double sqr_weight = 0; + size_t irealisation = 0; + size_t N = 0; /* #realisations that do not fail */ + size_t i; + res_T res = RES_OK; + + if(!scn || !nrealisations || !uv || time < 0 || fp_to_meter <= 0 + || Tref < 0 || !out_estimator) { + res = RES_BAD_ARG; + goto error; + } + + /* Check the primitive identifier */ + if(iprim >= scene_get_primitives_count(scn)) { + log_err(scn->dev, +"%s: invalid primitive identifier `%lu'. It must be less than %lu.\n", + FUNC_NAME, + (unsigned long)iprim, + (unsigned long)scene_get_primitives_count(scn)); + res = RES_BAD_ARG; + goto error; + } + + /* Check parametric coordinates */ + if(scene_is_2d(scn)) { + const double v = CLAMP(1.0 - uv[0], 0, 1); + if(uv[0] < 0 || uv[0] > 1 || !eq_eps(uv[0] + v, 1, 1.e-6)) { + log_err(scn->dev, +"%s: invalid parametric coordinates %g.\n" +"u + (1-u) must be equal to 1 with u [0, 1].\n", + FUNC_NAME, uv[0]); + res = RES_BAD_ARG; + goto error; + } + } else { + const double w = CLAMP(1 - uv[0] - uv[1], 0, 1); + if(uv[0] < 0 || uv[1] < 0 || uv[0] > 1 || uv[1] > 1 + || !eq_eps(w + uv[0] + uv[1], 1, 1.e-6)) { + log_err(scn->dev, +"%s: invalid parametric coordinates [%g, %g].\n" +"u + v + (1-u-v) must be equal to 1 with u and v in [0, 1].\n", + FUNC_NAME, uv[0], uv[1]); + res = RES_BAD_ARG; + goto error; + } + } + + /* Create the proxy RNG */ + res = ssp_rng_proxy_create(scn->dev->allocator, &ssp_rng_mt19937_64, + scn->dev->nthreads, &rng_proxy); + if(res != RES_OK) goto error; + + /* Create the per thread RNG */ + rngs = MEM_CALLOC + (scn->dev->allocator, scn->dev->nthreads, sizeof(struct ssp_rng*)); + if(!rngs) { + res = RES_MEM_ERR; + goto error; + } + FOR_EACH(i, 0, scn->dev->nthreads) { + res = ssp_rng_proxy_create_rng(rng_proxy, i, rngs+i); + if(res != RES_OK) goto error; + } + + /* Create the estimator */ + res = estimator_create(scn->dev, &estimator); + if(res != RES_OK) goto error; + + /* Here we go! Launch the Monte Carlo estimation */ + omp_set_num_threads((int)scn->dev->nthreads); + #pragma omp parallel for schedule(static) reduction(+:weight,sqr_weight,N) + for(irealisation = 0; irealisation < nrealisations; ++irealisation) { + res_T res_local; + double w; + const int ithread = omp_get_thread_num(); + struct ssp_rng* rng = rngs[ithread]; + + if(ATOMIC_GET(&res) != RES_OK) continue; /* An error occured */ + + if(scene_is_2d(scn)) { + res_local = boundary_realisation_2d + (scn, rng, iprim, uv, time, fp_to_meter, Tarad, Tref, &w); + } else { + res_local = boundary_realisation_3d + (scn, rng, iprim, uv, time, fp_to_meter, Tarad, Tref, &w); + } + if(res_local != RES_OK) { + if(res_local != RES_BAD_OP) { + ATOMIC_SET(&res, res_local); + continue; + } + } else { + weight += w; + sqr_weight += w*w; + ++N; + } + } + + estimator->nrealisations = N; + estimator->nfailures = nrealisations - N; + estimator->temperature.E = weight / (double)N; + estimator->temperature.V = + sqr_weight / (double)N + - estimator->temperature.E * estimator->temperature.E; + estimator->temperature.SE = sqrt(estimator->temperature.V / (double)N); + +exit: + if(rngs) { + FOR_EACH(i, 0, scn->dev->nthreads) { + if(rngs[i]) SSP(rng_ref_put(rngs[i])); + } + MEM_RM(scn->dev->allocator, rngs); + } + if(rng_proxy) SSP(rng_proxy_ref_put(rng_proxy)); + if(out_estimator) *out_estimator = estimator; + return res; +error: + if(estimator) { + SDIS(estimator_ref_put(estimator)); + estimator = NULL; + } + goto exit; +} + +res_T sdis_solve_camera (struct sdis_scene* scn, const struct sdis_camera* cam, diff --git a/src/sdis_solve_Xd.h b/src/sdis_solve_Xd.h @@ -1,4 +1,4 @@ -/* Copyright (C) |Meso|Star> 2016-2018 (contact@meso-star.com) +/* 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 @@ -76,6 +76,8 @@ reflect(float res[3], const float V[3], const float N[3]) #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_POSITION CONCAT(CONCAT(S, DIM), D_POSITION) +#define SXD_GEOMETRY_NORMAL CONCAT(CONCAT(S, DIM), D_GEOMETRY_NORMAL) #define SXD CONCAT(CONCAT(S, DIM), D) /* Vector macros generic to SDIS_SOLVE_DIMENSION */ @@ -149,7 +151,6 @@ XD(radiative_temperature) /******************************************************************************* * Helper functions ******************************************************************************/ - static FINLINE void XD(move_pos)(double pos[DIM], const float dir[DIM], const float delta) { @@ -270,6 +271,7 @@ XD(trace_radiative_path) r = ssp_rng_canonical(rng); if(r < epsilon) { T->func = XD(boundary_temperature); + rwalk->mdm = NULL; /* The random walk is at an interface between 2 media */ break; } @@ -541,6 +543,7 @@ XD(boundary_temperature) const struct sdis_medium* mdm_back = NULL; double tmp; ASSERT(scn && fp_to_meter > 0 && ctx && rwalk && rng && T); + ASSERT(rwalk->mdm == NULL); ASSERT(!SXD_HIT_NONE(&rwalk->hit)); XD(setup_interface_fragment)(&frag, &rwalk->vtx, &rwalk->hit); @@ -602,6 +605,7 @@ XD(solid_temperature) double cp; /* Calorific capacity */ double tau, mu; double tmp; + double power; float delta, delta_solid; /* Random walk numerical parameter */ float range[2]; float dir0[DIM], dir1[DIM]; @@ -659,6 +663,14 @@ XD(solid_temperature) return RES_OK; } + /* Add the volumic power density to the measured temperature */ + power = solid_get_volumic_power(mdm, &rwalk->vtx); + if(power > 0) { + const double delta_in_meter = delta * fp_to_meter; + tmp = power * delta_in_meter * delta_in_meter / (2.0 * DIM * lambda); + T->value += tmp; + } + /* Define if the random walk hits something along dir0 */ rwalk->hit = hit0.distance > delta ? SXD_HIT_NULL : hit0; @@ -696,6 +708,7 @@ XD(solid_temperature) } while(SXD_HIT_NONE(&rwalk->hit)); T->func = XD(boundary_temperature); + rwalk->mdm = NULL; /* The random walk is at an interface between 2 media */ return RES_OK; } @@ -776,6 +789,70 @@ XD(probe_realisation) return RES_OK; } +static res_T +XD(boundary_realisation) + (struct sdis_scene* scn, + struct ssp_rng* rng, + const size_t iprim, + const double uv[DIM], + const double time, + const double fp_to_meter, + const double Tarad, + const double Tref, + double* weight) +{ + struct rwalk_context ctx; + struct XD(rwalk) rwalk = XD(RWALK_NULL); + struct XD(temperature) T = XD(TEMPERATURE_NULL); + struct sXd(attrib) attr; +#if SDIS_SOLVE_DIMENSION == 2 + float st; +#else + float st[2]; +#endif + res_T res = RES_OK; + ASSERT(uv && fp_to_meter > 0 && weight && time >= 0); + + T.func = XD(boundary_temperature); + + rwalk.hit.distance = 0; + rwalk.vtx.time = time; + rwalk.mdm = NULL; /* The random walk is at an interface between 2 media */ + +#if SDIS_SOLVE_DIMENSION == 2 + st = (float)uv[0]; +#else + f2_set_d2(st, uv); +#endif + + /* Fetch the primitive */ + SXD(scene_view_get_primitive + (scn->sXd(view), (unsigned int)iprim, &rwalk.hit.prim)); + + /* Retrieve the world space position of the probe onto the primitive */ + SXD(primitive_get_attrib(&rwalk.hit.prim, SXD_POSITION, st, &attr)); + dX_set_fX(rwalk.vtx.P, attr.value); + + /* Retrieve the primitive normal */ + SXD(primitive_get_attrib(&rwalk.hit.prim, SXD_GEOMETRY_NORMAL, st, &attr)); + fX(set)(rwalk.hit.normal, attr.value); + +#if SDIS_SOLVE_DIMENSION==2 + rwalk.hit.u = st; +#else + f2_set(rwalk.hit.uv, st); +#endif + + ctx.Tarad = Tarad; + ctx.Tref3 = Tref*Tref*Tref; + + res = XD(compute_temperature)(scn, fp_to_meter, &ctx, &rwalk, rng, &T); + if(res != RES_OK) return res; + + *weight = T.value; + return RES_OK; +} + #if SDIS_SOLVE_DIMENSION == 3 static res_T XD(ray_realisation) @@ -831,6 +908,8 @@ error: #undef SXD_HIT_NONE #undef SXD_HIT_NULL #undef SXD_HIT_NULL__ +#undef SXD_POSITION +#undef SXD_GEOMETRY_NORMAL #undef SXD #undef dX #undef fX diff --git a/src/test_sdis_accum_buffer.c b/src/test_sdis_accum_buffer.c @@ -1,4 +1,4 @@ -/* Copyright (C) |Meso|Star> 2016-2018 (contact@meso-star.com) +/* 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 diff --git a/src/test_sdis_camera.c b/src/test_sdis_camera.c @@ -1,4 +1,4 @@ -/* Copyright (C) |Meso|Star> 2016-2018 (contact@meso-star.com) +/* 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 diff --git a/src/test_sdis_conducto_radiative.c b/src/test_sdis_conducto_radiative.c @@ -1,4 +1,4 @@ -/* Copyright (C) |Meso|Star> 2016-2018 (contact@meso-star.com) +/* 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 diff --git a/src/test_sdis_conducto_radiative_2d.c b/src/test_sdis_conducto_radiative_2d.c @@ -1,4 +1,4 @@ -/* Copyright (C) |Meso|Star> 2016-2018 (contact@meso-star.com) +/* 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 diff --git a/src/test_sdis_data.c b/src/test_sdis_data.c @@ -1,4 +1,4 @@ -/* Copyright (C) |Meso|Star> 2016-2018 (contact@meso-star.com) +/* 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 diff --git a/src/test_sdis_device.c b/src/test_sdis_device.c @@ -1,4 +1,4 @@ -/* Copyright (C) |Meso|Star> 2016-2018 (contact@meso-star.com) +/* 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 diff --git a/src/test_sdis_interface.c b/src/test_sdis_interface.c @@ -1,4 +1,4 @@ -/* Copyright (C) |Meso|Star> 2016-2018 (contact@meso-star.com) +/* 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 diff --git a/src/test_sdis_medium.c b/src/test_sdis_medium.c @@ -1,4 +1,4 @@ -/* Copyright (C) |Meso|Star> 2016-2018 (contact@meso-star.com) +/* 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 diff --git a/src/test_sdis_scene.c b/src/test_sdis_scene.c @@ -1,4 +1,4 @@ -/* Copyright (C) |Meso|Star> 2016-2018 (contact@meso-star.com) +/* 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 diff --git a/src/test_sdis_solve_camera.c b/src/test_sdis_solve_camera.c @@ -1,4 +1,4 @@ -/* Copyright (C) |Meso|Star> 2016-2018 (contact@meso-star.com) +/* 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 diff --git a/src/test_sdis_solve_probe.c b/src/test_sdis_solve_probe.c @@ -1,4 +1,4 @@ -/* Copyright (C) |Meso|Star> 2016-2018 (contact@meso-star.com) +/* 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 diff --git a/src/test_sdis_solve_probe2.c b/src/test_sdis_solve_probe2.c @@ -1,4 +1,4 @@ -/* Copyright (C) |Meso|Star> 2016-2018 (contact@meso-star.com) +/* 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 diff --git a/src/test_sdis_solve_probe2_2d.c b/src/test_sdis_solve_probe2_2d.c @@ -1,4 +1,4 @@ -/* Copyright (C) |Meso|Star> 2016-2018 (contact@meso-star.com) +/* 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 diff --git a/src/test_sdis_solve_probe3.c b/src/test_sdis_solve_probe3.c @@ -1,4 +1,4 @@ -/* Copyright (C) |Meso|Star> 2016-2018 (contact@meso-star.com) +/* 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 diff --git a/src/test_sdis_solve_probe3_2d.c b/src/test_sdis_solve_probe3_2d.c @@ -1,4 +1,4 @@ -/* Copyright (C) |Meso|Star> 2016-2018 (contact@meso-star.com) +/* 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 diff --git a/src/test_sdis_solve_probe_2d.c b/src/test_sdis_solve_probe_2d.c @@ -1,4 +1,4 @@ -/* Copyright (C) |Meso|Star> 2016-2018 (contact@meso-star.com) +/* 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 diff --git a/src/test_sdis_solve_probe_boundary.c b/src/test_sdis_solve_probe_boundary.c @@ -0,0 +1,398 @@ +/* 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 + +/******************************************************************************* + * 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 +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_delta_boundary + (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) +{ + (void)data; + CHK(vtx != NULL); + return 2.1/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; +} + +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; +} + +/******************************************************************************* + * 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 = DUMMY_INTERFACE_SHADER; + 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, 0, &dev) == RES_OK); + + /* Create the fluid medium */ + fluid_shader.temperature = fluid_get_temperature; + CHK(sdis_fluid_create(dev, &fluid_shader, NULL, &fluid) == RES_OK); + + /* Create the solid_medium */ + solid_shader.calorific_capacity = solid_get_calorific_capacity; + solid_shader.thermal_conductivity = solid_get_thermal_conductivity; + solid_shader.volumic_mass = solid_get_volumic_mass; + solid_shader.delta_solid = solid_get_delta; + solid_shader.delta_boundary = solid_get_delta_boundary; + solid_shader.temperature = solid_get_temperature; + CHK(sdis_solid_create(dev, &solid_shader, NULL, &solid) == RES_OK); + + /* Setup the interface shader */ + interf_shader.temperature = interface_get_temperature; + interf_shader.convection_coef = interface_get_convection_coef; + interf_shader.emissivity = interface_get_emissivity; + interf_shader.specular_fraction = interface_get_specular_fraction; + + /* 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 + CHK(SOLVE(NULL, N, iprim, uv, INF, 1.0, 0, 0, &estimator) == RES_BAD_ARG); + CHK(SOLVE(box_scn, 0, iprim, uv, INF, 1.0, 0, 0, &estimator) == RES_BAD_ARG); + CHK(SOLVE(box_scn, N, 12, uv, INF, 1.0, 0, 0, &estimator) == RES_BAD_ARG); + CHK(SOLVE(box_scn, N, iprim, NULL, INF, 1.0, 0, 0, &estimator) == RES_BAD_ARG); + CHK(SOLVE(box_scn, N, iprim, uv, -1, 1.0, 0, 0, &estimator) == RES_BAD_ARG); + CHK(SOLVE(box_scn, N, iprim, uv, INF, 1.0, 0, 0, NULL) == RES_BAD_ARG); + CHK(SOLVE(box_scn, N, iprim, uv, INF, 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); + + CHK(sdis_scene_get_boundary_position(NULL, iprim, uv, pos) == RES_BAD_ARG); + CHK(sdis_scene_get_boundary_position(box_scn, 12, uv, pos) == RES_BAD_ARG); + CHK(sdis_scene_get_boundary_position(box_scn, iprim, NULL, pos) == RES_BAD_ARG); + CHK(sdis_scene_get_boundary_position(box_scn, iprim, uv, NULL) == RES_BAD_ARG); + CHK(sdis_scene_get_boundary_position(box_scn, iprim, uv, pos) == RES_OK); + + ref = (H*Tf + LAMBDA * Tb) / (H + LAMBDA); + + 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(eq_eps(T.E, ref, T.SE)); + + uv[0] = 0.5; + iprim = 3; + CHK(SOLVE(square_scn, N, iprim, uv, INF, 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); + + 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(eq_eps(T.E, ref, 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; +} + diff --git a/src/test_sdis_utils.h b/src/test_sdis_utils.h @@ -1,4 +1,4 @@ -/* Copyright (C) |Meso|Star> 2016-2018 (contact@meso-star.com) +/* 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 @@ -40,13 +40,13 @@ static const size_t box_nvertices = sizeof(box_vertices) / sizeof(double[3]); /* The following array lists the indices toward the 3D vertices of each * triangle. - * ,6---,7 ,6----7 - * ,' | ,'/| ,' | \ | - * 2----3' / | 2', | \ | - * |', | / ,5 | ',4---,5 - * | ',|/,' | ,' | ,' - * 0----1' 0----1' - * Front, right Back, left and + * ,2---,3 ,2----3 + * ,' | ,'/| ,'/| \ | + * 6----7' / | 6' / | \ | Y + * |', | / ,1 | / ,0---,1 | + * | ',|/,' |/,' | ,' o--X + * 4----5' 4----5' / + * Front, right Back, left and Z * and Top faces bottom faces */ static const size_t box_indices[12/*#triangles*/*3/*#indices per triangle*/] = { 0, 2, 1, 1, 2, 3, /* Front face */ @@ -104,6 +104,7 @@ static const struct sdis_solid_shader DUMMY_SOLID_SHADER = { dummy_medium_getter, dummy_medium_getter, dummy_medium_getter, + dummy_medium_getter, dummy_medium_getter }; diff --git a/src/test_sdis_volumic_power.c b/src/test_sdis_volumic_power.c @@ -0,0 +1,367 @@ +/* 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 a solid cube/square whose temperature is unknown. + * The convection coefficient with the surrounding fluid is null everywhere The + * Temperature of the +/- X faces are fixed to T0, and the solid has a volumic + * power of P0. This test computes the temperature of a probe position pos into + * the solid and check that it is is equal to: + * + * T(pos) = P0 / (2*LAMBDA) * (A^2/4 - (pos-0.5)^2) + T0 + * + * with LAMBDA the conductivity of the cube and A the size of the cube/square, + * i.e. 1. + * + * 3D 2D + * + * ///// (1,1,1) ///// (1,1) + * +-------+ +-------+ + * /' /| | | + * +-------+ T0 T0 T0 + * T0 +.....|.+ | | + * |, |/ +-------+ + * +-------+ (0,0) ///// + * (0,0,0) ///// + */ + +#define UNKNOWN_TEMPERATURE -1 +#define N 10000 /* #realisations */ + +#define T0 320 +#define LAMBDA 0.1 +#define P0 10 + +/******************************************************************************* + * 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 +fluid_get_temperature + (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) +{ + (void)data; + CHK(vtx != NULL); + return UNKNOWN_TEMPERATURE; +} + +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_delta_boundary + (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) +{ + (void)data; + CHK(vtx != NULL); + return 2.1/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; +} + +static double +solid_get_volumic_power + (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) +{ + (void)data; + CHK(vtx != NULL); + return P0; +} + +/******************************************************************************* + * Interfaces + ******************************************************************************/ +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 0; +} + +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; +} + +/******************************************************************************* + * 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_T0 = 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*/]; + struct interf* interf_props = NULL; + double pos[3]; + double x; + double ref; + size_t nreals; + size_t nfails; + (void)argc, (void)argv; + + CHK(mem_init_proxy_allocator(&allocator, &mem_default_allocator) == RES_OK); + CHK(sdis_device_create + (NULL, &allocator, SDIS_NTHREADS_DEFAULT, 0, &dev) == RES_OK); + + /* Create the fluid medium */ + fluid_shader.temperature = fluid_get_temperature; + CHK(sdis_fluid_create(dev, &fluid_shader, NULL, &fluid) == RES_OK); + + /* Create the solid_medium */ + solid_shader.calorific_capacity = solid_get_calorific_capacity; + solid_shader.thermal_conductivity = solid_get_thermal_conductivity; + solid_shader.volumic_mass = solid_get_volumic_mass; + solid_shader.delta_solid = solid_get_delta; + solid_shader.delta_boundary = solid_get_delta_boundary; + solid_shader.temperature = solid_get_temperature; + solid_shader.volumic_power = solid_get_volumic_power; + CHK(sdis_solid_create(dev, &solid_shader, NULL, &solid) == RES_OK); + + /* Setup the interface shader */ + interf_shader.temperature = interface_get_temperature; + interf_shader.convection_coef = interface_get_convection_coef; + interf_shader.emissivity = interface_get_emissivity; + interf_shader.specular_fraction = interface_get_specular_fraction; + + /* 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->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 T0 interface */ + CHK(sdis_data_create(dev, sizeof(struct interf), 16, NULL, &data) == RES_OK); + interf_props = sdis_data_get(data); + interf_props->temperature = T0; + CHK(sdis_interface_create + (dev, solid, fluid, &interf_shader, data, &interf_T0) == 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_T0; /* Left */ + box_interfaces[4] = box_interfaces[5] = interf_adiabatic; /* Back */ + box_interfaces[6] = box_interfaces[7] = interf_T0; /* 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_T0; /* Lef */ + square_interfaces[2] = interf_adiabatic; /* Top */ + square_interfaces[3] = interf_T0; /* 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_T0) == RES_OK); + + d3_splat(pos, 0.25); + x = pos[0] - 0.5; + ref = P0 / (2*LAMBDA) * (1.0/4.0 - x*x) + T0; + + /* Solve in 3D */ + CHK(sdis_solve_probe(box_scn, N, pos, INF, 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("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(eq_eps(T.E, ref, T.SE*2)); + + /* Solve in 2D */ + CHK(sdis_solve_probe(square_scn, N, pos, INF, 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("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(eq_eps(T.E, ref, T.SE*2.0)); + + 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; +}