commit 257b7e805d12062dbcb3a34a1cc613515b27b80d parent 2bfab77721b5c4060317343582ec2d0caa6842a9 Author: Vincent Forest <vincent.forest@meso-star.com> Date: Mon, 9 Apr 2018 15:29:26 +0200 Merge branch 'feature_boundary_temperature' into develop Diffstat:
38 files changed, 1102 insertions(+), 45 deletions(-)
diff --git a/README.md b/README.md @@ -46,7 +46,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 @@ -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; +}