stardis-solver

Solve coupled heat transfers
git clone git://git.meso-star.fr/stardis-solver.git
Log | Files | Refs | README | LICENSE

commit 22ac22d2334f5a38325f4ea8b71af885511b1d09
parent 1e166a8821fc4f8eb864c0bfd415b690025ab123
Author: Vincent Forest <vincent.forest@meso-star.com>
Date:   Wed, 18 Apr 2018 16:04:33 +0200

Merge branch 'release_0.3'

Diffstat:
MREADME.md | 16++++++++++++++++
Mcmake/CMakeLists.txt | 3++-
Msrc/sdis.h | 165+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++----------------
Msrc/sdis_interface.c | 68++++++++++++++++++++++++++++++++++++++++++++++++++------------------
Msrc/sdis_interface_c.h | 60++++++++++++++++++++++++++++++++++++++++++++++++------------
Msrc/sdis_medium.c | 4++--
Msrc/sdis_medium_c.h | 24++++++++++++------------
Msrc/sdis_scene.c | 122+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Msrc/sdis_solve.c | 9+++++----
Msrc/sdis_solve_Xd.h | 132+++++++++++++++++++++++++++++++++++++++++++++++++++++++++----------------------
Msrc/test_sdis_conducto_radiative.c | 19++++++++++++++-----
Msrc/test_sdis_conducto_radiative_2d.c | 105++++++++++++++++++++++++++++++++++++++++++++++++++++++-------------------------
Asrc/test_sdis_flux.c | 342+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Msrc/test_sdis_interface.c | 29++++++++++++++++++-----------
Msrc/test_sdis_scene.c | 207+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++----------
Msrc/test_sdis_solve_camera.c | 16+++++++++++-----
Msrc/test_sdis_solve_probe.c | 11++++++-----
Msrc/test_sdis_solve_probe2.c | 14+++-----------
Msrc/test_sdis_solve_probe2_2d.c | 18+++++++-----------
Msrc/test_sdis_solve_probe3.c | 19++++---------------
Msrc/test_sdis_solve_probe3_2d.c | 20++++++--------------
Msrc/test_sdis_solve_probe_2d.c | 13++-----------
Msrc/test_sdis_solve_probe_boundary.c | 47+++++++++++++++--------------------------------
Msrc/test_sdis_utils.h | 12+++++++++---
Msrc/test_sdis_volumic_power.c | 26++++----------------------
25 files changed, 1180 insertions(+), 321 deletions(-)

diff --git a/README.md b/README.md @@ -23,6 +23,22 @@ variable the install directories of its dependencies. ## Release notes +### Version 0.3 + +- Some interface properties become double sided: the temperature, emissivity + and specular fraction is defined for each side of the interface. Actually, + only the convection coefficient is shared by the 2 sides of the interface. + The per side interface properties are grouped into the new `struct + sdis_interface_side_shader` data structure. +- Add the support of fixed flux: the flux is per side interface property. + Currently, the flux can be fixed only for the interface sides facing a solid + medium. +- Update the default comportment of the interface shader when a function is not + set. +- Rename the `SDIS_MEDIUM_<FLUID|SOLID>` constants in `SDIS_<FLUID|SOLID>`. +- Rename the `enum sdis_side_flag` enumerate in `enum sdis_side` and update its + values. + ### Version 0.2 - Add the support of volumic power to solid media: add the `volumic_power` diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt @@ -40,7 +40,7 @@ rcmake_append_runtime_dirs(_runtime_dirs RSys Star3D StarSP) # Configure and define targets ################################################################################ set(VERSION_MAJOR 0) -set(VERSION_MINOR 2) +set(VERSION_MINOR 3) set(VERSION_PATCH 0) set(VERSION ${VERSION_MAJOR}.${VERSION_MINOR}.${VERSION_PATCH}) @@ -124,6 +124,7 @@ if(NOT NO_TEST) new_test(test_sdis_conducto_radiative_2d) new_test(test_sdis_data) new_test(test_sdis_device) + new_test(test_sdis_flux) new_test(test_sdis_interface) new_test(test_sdis_medium) new_test(test_sdis_scene) diff --git a/src/sdis.h b/src/sdis.h @@ -40,6 +40,9 @@ * as CPU cores */ #define SDIS_NTHREADS_DEFAULT (~0u) +#define SDIS_VOLUMIC_POWER_NONE DBL_MAX /* <=> No volumic power */ +#define SDIS_FLUX_NONE DBL_MAX /* <=> No flux */ + /* Forward declaration of external opaque data types */ struct logger; struct mem_allocator; @@ -60,15 +63,15 @@ struct sdis_interface; struct sdis_medium; struct sdis_scene; -enum sdis_side_flag { - SDIS_FRONT = BIT(0), - SDIS_BACK = BIT(1), - SDIS_SIDE_NULL__ = BIT(2) +enum sdis_side { + SDIS_FRONT, + SDIS_BACK, + SDIS_SIDE_NULL__ }; enum sdis_medium_type { - SDIS_MEDIUM_FLUID, - SDIS_MEDIUM_SOLID, + SDIS_FLUID, + SDIS_SOLID, SDIS_MEDIUM_TYPES_COUNT__ }; @@ -82,20 +85,24 @@ struct sdis_rwalk_vertex { static const struct sdis_rwalk_vertex SDIS_RWALK_VERTEX_NULL = SDIS_RWALK_VERTEX_NULL__; -/* Spatiotemporal position onto an interface */ +/* Spatiotemporal position onto an interface. As a random walk vertex, it + * stores the position and time of the random walk, but since it lies onto an + * interface, it has additionnal parameters as the normal of the interface and + * the parametric coordinate of the position onto the interface */ struct sdis_interface_fragment { double P[3]; /* World space position */ double Ng[3]; /* Normalized world space geometry normal at the interface */ double uv[2]; /* Parametric coordinates of the interface */ double time; /* Current time */ + enum sdis_side side; }; -#define SDIS_INTERFACE_FRAGMENT_NULL__ {{0}, {0}, {0}, -1} +#define SDIS_INTERFACE_FRAGMENT_NULL__ {{0}, {0}, {0}, -1, SDIS_SIDE_NULL__} static const struct sdis_interface_fragment SDIS_INTERFACE_FRAGMENT_NULL = SDIS_INTERFACE_FRAGMENT_NULL__; /* Monte-Carlo accumulator */ struct sdis_accum { - double sum_weights; /* Sum of Monte-Carlo weight */ + double sum_weights; /* Sum of Monte-Carlo weights */ double sum_weights_sqr; /* Sum of Monte-Carlo square weights */ size_t nweights; /* #accumulated weights */ }; @@ -109,27 +116,33 @@ struct sdis_mc { #define SDIS_MC_NULL__ {0, 0, 0} static const struct sdis_mc SDIS_MC_NULL = SDIS_MC_NULL__; -/* Functor type to retrieve the medium properties. */ +/* Functor type used to retrieve the spatio temporal physical properties of a + * medium. */ typedef double (*sdis_medium_getter_T) (const struct sdis_rwalk_vertex* vert, struct sdis_data* data); -/* Functor type to retrieve the interface properties. */ +/* Functor type used to retrieve the spatio temporal physical properties of an + * interface. */ typedef double (*sdis_interface_getter_T) (const struct sdis_interface_fragment* frag, struct sdis_data* data); +/* Define the physical properties of a solid */ struct sdis_solid_shader { /* Properties */ - sdis_medium_getter_T calorific_capacity; - sdis_medium_getter_T thermal_conductivity; - sdis_medium_getter_T volumic_mass; + sdis_medium_getter_T calorific_capacity; /* In J.K^-1.kg^-1 */ + sdis_medium_getter_T thermal_conductivity; /* In W.m^-1.K^-1 */ + sdis_medium_getter_T volumic_mass; /* In kg.m^-3 */ sdis_medium_getter_T delta_solid; sdis_medium_getter_T delta_boundary; - sdis_medium_getter_T volumic_power; /* May be NULL <=> no volumic power */ + /* May be NULL if there is no volumic power. One can also return + * SDIS_VOLUMIC_POWER_NONE to define that there is no volumic power at the + * submitted position and time */ + sdis_medium_getter_T volumic_power; /* In W.m^-3 */ /* Initial/limit condition. A temperature < 0 means that the temperature is * unknown for the submitted random walk vertex. */ @@ -139,10 +152,11 @@ struct sdis_solid_shader { static const struct sdis_solid_shader SDIS_SOLID_SHADER_NULL = SDIS_SOLID_SHADER_NULL__; +/* Define the physical properties of a fluid */ struct sdis_fluid_shader { /* Properties */ - sdis_medium_getter_T calorific_capacity; - sdis_medium_getter_T volumic_mass; + sdis_medium_getter_T calorific_capacity; /* In J.K^-1.kg^-1 */ + sdis_medium_getter_T volumic_mass; /* In kg.m^-3 */ /* Initial/limit condition. A temperature < 0 means that the temperature is * unknown for the submitted position and time. */ @@ -152,15 +166,33 @@ struct sdis_fluid_shader { static const struct sdis_fluid_shader SDIS_FLUID_SHADER_NULL = SDIS_FLUID_SHADER_NULL__; +/* Define the physical properties of one side of an interface. */ +struct sdis_interface_side_shader { + /* Fixed temperature/flux. May be NULL if the temperature/flux is unknown + * onto the whole interface */ + sdis_interface_getter_T temperature; /* In Kelvin. < 0 <=> Unknown temp */ + sdis_interface_getter_T flux; /* In W.m^-2. SDIS_FLUX_NONE <=> no flux */ + + /* Control the emissivity of the interface. May be NULL for solid/solid + * interface or if the emissivity is 0 onto the whole interface. */ + sdis_interface_getter_T emissivity; /* Overall emissivity. */ + sdis_interface_getter_T specular_fraction; /* Specular part in [0,1] */ +}; +#define SDIS_INTERFACE_SIDE_SHADER_NULL__ { NULL, NULL, NULL, NULL } +static const struct sdis_interface_side_shader SDIS_INTERFACE_SIDE_SHADER_NULL = + SDIS_INTERFACE_SIDE_SHADER_NULL__; + +/* Define the physical properties of an interface between 2 media .*/ struct sdis_interface_shader { - sdis_interface_getter_T temperature; /* Limit condition. NULL <=> Unknown */ - sdis_interface_getter_T convection_coef; /* May be NULL for solid/solid */ + /* May be NULL for solid/solid or if the convection coefficient is 0 onto + * the whole interface. */ + sdis_interface_getter_T convection_coef; /* In W.K^-1.m^-2 */ - /* Interface emssivity. May be NULL for solid/solid interface */ - sdis_interface_getter_T emissivity; /* Overall emissivity */ - sdis_interface_getter_T specular_fraction; /* Specular fraction in [0, 1] */ + struct sdis_interface_side_shader front; + struct sdis_interface_side_shader back; }; -#define SDIS_INTERFACE_SHADER_NULL__ {NULL, NULL, NULL, NULL} +#define SDIS_INTERFACE_SHADER_NULL__ \ + {NULL, SDIS_INTERFACE_SIDE_SHADER_NULL__, SDIS_INTERFACE_SIDE_SHADER_NULL__} static const struct sdis_interface_shader SDIS_INTERFACE_SHADER_NULL = SDIS_INTERFACE_SHADER_NULL__; @@ -204,7 +236,8 @@ sdis_device_ref_put /******************************************************************************* * A data stores in the Stardis memory space a set of user defined data. It can - * be seen as a ref counted memory space allocated by Stardis. + * be seen as a ref counted memory space allocated by Stardis. It is used to + * attach user data to the media and to the interfaces. ******************************************************************************/ SDIS_API res_T sdis_data_create @@ -265,13 +298,15 @@ sdis_camera_look_at const double up[3]); /******************************************************************************* - * A buffer of accumulations + * A buffer of accumulations is a 2D array whose each cell stores an + * Monte-Carlo accumulation, i.e. a sum of MC weights, the sum of their square + * and the overall number of summed weights (see struct sdis_accum) ******************************************************************************/ SDIS_API res_T sdis_accum_buffer_create (struct sdis_device* dev, - const size_t width, - const size_t height, + const size_t width, /* #cells in X */ + const size_t height, /* #cells in Y */ struct sdis_accum_buffer** buf); SDIS_API res_T @@ -287,6 +322,7 @@ sdis_accum_buffer_get_layout (const struct sdis_accum_buffer* buf, struct sdis_accum_buffer_layout* layout); +/* Get a read only pointer toward the memory space of the accum buffer */ SDIS_API res_T sdis_accum_buffer_map (const struct sdis_accum_buffer* buf, @@ -296,7 +332,11 @@ SDIS_API res_T sdis_accum_buffer_unmap (const struct sdis_accum_buffer* buf); -/* Helper function that matches the `sdis_write_accums_T' functor type */ +/* Helper function that matches the `sdis_write_accums_T' functor type. On can + * send this function directly to the sdis_solve_camera function, to fill the + * accum buffer with the estimation of the radiative temperature that reaches + * each pixel of an image whose definition matches the definition of the accum + * buffer. */ SDIS_API res_T sdis_accum_buffer_write (void* buf, /* User data */ @@ -334,13 +374,13 @@ sdis_medium_get_type (const struct sdis_medium* medium); /******************************************************************************* - * An interface is the boundary between 2 mediums. + * An interface is the boundary between 2 media. ******************************************************************************/ SDIS_API res_T sdis_interface_create (struct sdis_device* dev, - struct sdis_medium* front, - struct sdis_medium* back, + struct sdis_medium* front, /* Medium on the front side of the geometry */ + struct sdis_medium* back, /* Medium on the back side of the geometry */ const struct sdis_interface_shader* shader, struct sdis_data* data, /* Data sent to the shader. May be NULL */ struct sdis_interface** interf); @@ -357,6 +397,18 @@ sdis_interface_ref_put * A scene is a collection of primitives. Each primitive is the geometric * support of the interface between 2 mediums. ******************************************************************************/ +/* Create a 3D scene. The geometry of the scene is defined by an indexed + * triangular mesh: each triangle is composed of 3 indices where each index + * references an absolute 3D position. The physical properties of an interface + * is defined by the interface of the triangle. + * + * Note that each triangle has 2 sides: a front and a back side. By convention, + * the front side of a triangle is the side where its vertices are clock wise + * ordered. The back side of a triangle is the exact opposite: it is the side + * where the triangle vertices are counter-clock wise ordered. The front and + * back media of a triangle interface directly refer to this convention and + * thus one has to take care of how the triangle vertices are defined to ensure + * that the front and the back media are correctly defined wrt the geometry. */ SDIS_API res_T sdis_scene_create (struct sdis_device* dev, @@ -371,6 +423,18 @@ sdis_scene_create void* ctx, /* Client side data sent as input of the previous callbacks */ struct sdis_scene** scn); +/* Create a 2D scene. The geometry of the 2D scene is defined by an indexed + * line segments: each segment is composed of 2 indices where each index + * references an absolute 2D position. The physical properties of an interface + * is defined by the interface of the segment. + * + * Note that each segment has 2 sides: a front and a back side. By convention, + * the front side of a segment is the side where its vertices are clock wise + * ordered. The back side of a segment is the exact opposite: it is the side + * where the segment vertices are counter-clock wise ordered. The front and + * back media of a segment interface directly refer to this convention and + * thus one has to take care of how the segment vertices are defined to ensure + * that the front and the back media are correctly defined wrt the geometry. */ SDIS_API res_T sdis_scene_2d_create (struct sdis_device* dev, @@ -400,13 +464,49 @@ sdis_scene_get_aabb double lower[3], double upper[3]); +/* Define the world space position of a point onto the primitive `iprim' whose + * parametric coordinate is uv. */ 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 */ + const double uv[2], /* Parametric coordinate onto the primitive */ double pos[3]); /* World space position */ +/* Project a world space position onto a primitive wrt its normal and compute + * the parametric coordinates of the projected point onto the primitive. This + * function may help to define the probe position onto a boundary as expected + * by the sdis_solve_probe_boundary function. + * + * Note that the projected point can lie outside the submitted primitive. In + * this case, the parametric coordinates are clamped against the primitive + * boundaries in order to ensure that the returned parametric coordinates are + * valid according to the primitive. To ensure this, in 2D, the parametric + * coordinate is simply clamped to [0, 1]. In 3D, the `uv' coordinates are + * clamped against the triangle edges. For instance, let the + * following triangle whose vertices are `a', `b' and `c': + * , , + * , B , + * , , + * b E1 + * E0 / \ ,P + * / \,*' + * / \ + * ....a-------c...... + * ' ' + * A ' E2 ' C + * ' ' + * The projected point `P' is orthogonally wrapped to the edge `ab', `bc' or + * `ca' if it lies in the `E0', `E1' or `E2' region, respectively. If `P' is in + * the `A', `B' or `C' region, then it is taken back to the `a', `b' or `c' + * vertex, respectively. */ +SDIS_API res_T +sdis_scene_boundary_project_position + (const struct sdis_scene* scn, + const size_t iprim, + const double pos[3], + double uv[]); + /******************************************************************************* * An estimator stores the state of a simulation ******************************************************************************/ @@ -454,6 +554,7 @@ sdis_solve_probe_boundary 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 enum sdis_side side, /* Side of iprim on which the probe lies */ 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 */ diff --git a/src/sdis_interface.c b/src/sdis_interface.c @@ -29,26 +29,52 @@ ******************************************************************************/ static int check_interface_shader - (const struct sdis_interface_shader* shader, + (struct sdis_device* dev, + const char* caller_name, + const struct sdis_interface_shader* shader, const struct sdis_medium* front, const struct sdis_medium* back) { - enum sdis_medium_type type0; - enum sdis_medium_type type1; - ASSERT(shader && front && back); + enum sdis_medium_type type[2]; + const struct sdis_interface_side_shader* shaders[2]; + int i; + ASSERT(dev && caller_name && shader && front && back); - type0 = sdis_medium_get_type(front); - type1 = sdis_medium_get_type(back); + type[0] = sdis_medium_get_type(front); + type[1] = sdis_medium_get_type(back); + shaders[0] = &shader->front; + shaders[1] = &shader->back; /* Fluid<->solid interface */ - if(type0 != type1) { - if(shader->convection_coef == NULL - || shader->emissivity == NULL - || shader->specular_fraction == NULL) { - return 0; - } + if(type[0] == SDIS_SOLID + && type[1] == SDIS_SOLID + && shader->convection_coef) { + log_warn(dev, + "%s: a solid/solid interface can't have a convection coefficient. This " + "function of the interface shader should be NULL.\n", caller_name); } + FOR_EACH(i, 0, 2) { + switch(type[i]) { + case SDIS_SOLID: + if(shaders[i]->emissivity || shaders[i]->specular_fraction) { + log_warn(dev, + "%s: the interface side toward a solid can't have the emissivity " + "and specular_fraction properties. The shader functions that return " + "these attributes should be NULL.\n", caller_name); + } + break; + case SDIS_FLUID: + if(shaders[i]->flux) { + log_warn(dev, + "%s: the interface side toward a fluid can't have a flux property. " + "The shader function that returns this attribute should be NULL.\n", + caller_name); + } + break; + default: FATAL("Unreachable code.\n"); break; + } + } return 1; } @@ -88,14 +114,14 @@ sdis_interface_create goto error; } - if(sdis_medium_get_type(front) == SDIS_MEDIUM_FLUID - && sdis_medium_get_type(back) == SDIS_MEDIUM_FLUID) { + if(sdis_medium_get_type(front) == SDIS_FLUID + && sdis_medium_get_type(back) == SDIS_FLUID) { log_err(dev, "%s: invalid fluid<->fluid interface.\n", FUNC_NAME); res = RES_BAD_ARG; goto error; } - if(!check_interface_shader(shader, front, back)) { + if(!check_interface_shader(dev, FUNC_NAME, shader, front, back)) { log_err(dev, "%s: invalid interface shader.\n", FUNC_NAME); res = RES_BAD_ARG; goto error; @@ -154,7 +180,7 @@ sdis_interface_ref_put(struct sdis_interface* interf) ******************************************************************************/ const struct sdis_medium* interface_get_medium - (const struct sdis_interface* interf, const enum sdis_side_flag side) + (const struct sdis_interface* interf, const enum sdis_side side) { struct sdis_medium* mdm = NULL; ASSERT(interf); @@ -177,27 +203,33 @@ void setup_interface_fragment_2d (struct sdis_interface_fragment* frag, const struct sdis_rwalk_vertex* vertex, - const struct s2d_hit* hit) + const struct s2d_hit* hit, + const enum sdis_side side) { ASSERT(frag && vertex && hit && !S2D_HIT_NONE(hit)); + ASSERT(side == SDIS_FRONT || side == SDIS_BACK); d2_set(frag->P, vertex->P); frag->P[2] = 0; d2_normalize(frag->Ng, d2_set_f2(frag->Ng, hit->normal)); frag->Ng[2] = 0; frag->uv[0] = hit->u; frag->time = vertex->time; + frag->side = side; } void setup_interface_fragment_3d (struct sdis_interface_fragment* frag, const struct sdis_rwalk_vertex* vertex, - const struct s3d_hit* hit) + const struct s3d_hit* hit, + const enum sdis_side side) { ASSERT(frag && vertex && hit && !S3D_HIT_NONE(hit)); + ASSERT(side == SDIS_FRONT || side == SDIS_BACK); d3_set(frag->P, vertex->P); d3_normalize(frag->Ng, d3_set_f3(frag->Ng, hit->normal)); d2_set_f2(frag->uv, hit->uv); frag->time = vertex->time; + frag->side = side; } diff --git a/src/sdis_interface_c.h b/src/sdis_interface_c.h @@ -39,7 +39,7 @@ struct sdis_interface { extern LOCAL_SYM const struct sdis_medium* interface_get_medium (const struct sdis_interface* interf, - const enum sdis_side_flag side); + const enum sdis_side side); extern LOCAL_SYM unsigned interface_get_id @@ -49,49 +49,85 @@ extern LOCAL_SYM void setup_interface_fragment_2d (struct sdis_interface_fragment* frag, const struct sdis_rwalk_vertex* vertex, - const struct s2d_hit* hit); + const struct s2d_hit* hit, + const enum sdis_side side); extern LOCAL_SYM void setup_interface_fragment_3d (struct sdis_interface_fragment* frag, const struct sdis_rwalk_vertex* vertex, - const struct s3d_hit* hit); + const struct s3d_hit* hit, + const enum sdis_side side); static INLINE double -interface_get_temperature +interface_get_convection_coef (const struct sdis_interface* interf, const struct sdis_interface_fragment* frag) { ASSERT(interf && frag); - if(!interf->shader.temperature) return -DBL_MAX; - return interf->shader.temperature(frag, interf->data); + return interf->shader.convection_coef + ? interf->shader.convection_coef(frag, interf->data) : 0; } static INLINE double -interface_get_convection_coef +interface_side_get_temperature + (const struct sdis_interface* interf, + const struct sdis_interface_fragment* frag) +{ + const struct sdis_interface_side_shader* shader; + ASSERT(interf && frag); + switch(frag->side) { + case SDIS_FRONT: shader = &interf->shader.front; break; + case SDIS_BACK: shader = &interf->shader.back; break; + default: FATAL("Unreachable code.\n"); + } + return shader->temperature ? shader->temperature(frag, interf->data) : -1; +} + +static INLINE double +interface_side_get_flux (const struct sdis_interface* interf, const struct sdis_interface_fragment* frag) { + const struct sdis_interface_side_shader* shader; ASSERT(interf && frag); - return interf->shader.convection_coef(frag, interf->data); + switch(frag->side) { + case SDIS_FRONT: shader = &interf->shader.front; break; + case SDIS_BACK: shader = &interf->shader.back; break; + default: FATAL("Unreachable code.\n"); + } + return shader->flux ? shader->flux(frag, interf->data) : SDIS_FLUX_NONE; } static INLINE double -interface_get_emissivity +interface_side_get_emissivity (const struct sdis_interface* interf, const struct sdis_interface_fragment* frag) { + const struct sdis_interface_side_shader* shader; ASSERT(interf && frag); - return interf->shader.emissivity(frag, interf->data); + switch(frag->side) { + case SDIS_FRONT: shader = &interf->shader.front; break; + case SDIS_BACK: shader = &interf->shader.back; break; + default: FATAL("Unreachable code\n"); break; + } + return shader->emissivity ? shader->emissivity(frag, interf->data) : 0; } static INLINE double -interface_get_specular_fraction +interface_side_get_specular_fraction (const struct sdis_interface* interf, const struct sdis_interface_fragment* frag) { + const struct sdis_interface_side_shader* shader; ASSERT(interf && frag); - return interf->shader.specular_fraction(frag, interf->data); + switch(frag->side) { + case SDIS_FRONT: shader = &interf->shader.front; break; + case SDIS_BACK: shader = &interf->shader.back; break; + default: FATAL("Unreachable code\n"); break; + } + return shader->specular_fraction + ? shader->specular_fraction(frag, interf->data) : 0; } #endif /* SDIS_INTERFACE_C_H */ diff --git a/src/sdis_medium.c b/src/sdis_medium.c @@ -114,7 +114,7 @@ sdis_fluid_create goto error; } - res = medium_create(dev, &medium, SDIS_MEDIUM_FLUID); + res = medium_create(dev, &medium, SDIS_FLUID); if(res != RES_OK) { log_err(dev, "%s: could not create the fluid medium.\n", FUNC_NAME); goto error; @@ -159,7 +159,7 @@ sdis_solid_create goto error; } - res = medium_create(dev, &medium, SDIS_MEDIUM_SOLID); + res = medium_create(dev, &medium, SDIS_SOLID); if(res != RES_OK) { log_err(dev, "%s: could not create the solid medium.\n", FUNC_NAME); goto error; diff --git a/src/sdis_medium_c.h b/src/sdis_medium_c.h @@ -38,7 +38,7 @@ static INLINE double fluid_get_calorific_capacity (const struct sdis_medium* mdm, const struct sdis_rwalk_vertex* vtx) { - ASSERT(mdm && mdm->type == SDIS_MEDIUM_FLUID); + ASSERT(mdm && mdm->type == SDIS_FLUID); return mdm->shader.fluid.calorific_capacity(vtx, mdm->data); } @@ -46,7 +46,7 @@ static INLINE double fluid_get_volumic_mass (const struct sdis_medium* mdm, const struct sdis_rwalk_vertex* vtx) { - ASSERT(mdm && mdm->type == SDIS_MEDIUM_FLUID); + ASSERT(mdm && mdm->type == SDIS_FLUID); return mdm->shader.fluid.volumic_mass(vtx, mdm->data); } @@ -54,7 +54,7 @@ static INLINE double fluid_get_temperature (const struct sdis_medium* mdm, const struct sdis_rwalk_vertex* vtx) { - ASSERT(mdm && mdm->type == SDIS_MEDIUM_FLUID); + ASSERT(mdm && mdm->type == SDIS_FLUID); return mdm->shader.fluid.temperature(vtx, mdm->data); } @@ -65,7 +65,7 @@ static INLINE double solid_get_calorific_capacity (const struct sdis_medium* mdm, const struct sdis_rwalk_vertex* vtx) { - ASSERT(mdm && mdm->type == SDIS_MEDIUM_SOLID); + ASSERT(mdm && mdm->type == SDIS_SOLID); return mdm->shader.solid.calorific_capacity(vtx, mdm->data); } @@ -73,7 +73,7 @@ static INLINE double solid_get_thermal_conductivity (const struct sdis_medium* mdm, const struct sdis_rwalk_vertex* vtx) { - ASSERT(mdm && mdm->type == SDIS_MEDIUM_SOLID); + ASSERT(mdm && mdm->type == SDIS_SOLID); return mdm->shader.solid.thermal_conductivity(vtx, mdm->data); } @@ -81,7 +81,7 @@ static INLINE double solid_get_volumic_mass (const struct sdis_medium* mdm, const struct sdis_rwalk_vertex* vtx) { - ASSERT(mdm && mdm->type == SDIS_MEDIUM_SOLID); + ASSERT(mdm && mdm->type == SDIS_SOLID); return mdm->shader.solid.volumic_mass(vtx, mdm->data); } @@ -89,7 +89,7 @@ static INLINE double solid_get_delta (const struct sdis_medium* mdm, const struct sdis_rwalk_vertex* vtx) { - ASSERT(mdm && mdm->type == SDIS_MEDIUM_SOLID); + ASSERT(mdm && mdm->type == SDIS_SOLID); return mdm->shader.solid.delta_solid(vtx, mdm->data); } @@ -97,7 +97,7 @@ static INLINE double solid_get_delta_boundary (const struct sdis_medium* mdm, const struct sdis_rwalk_vertex* vtx) { - ASSERT(mdm && mdm->type == SDIS_MEDIUM_SOLID); + ASSERT(mdm && mdm->type == SDIS_SOLID); return mdm->shader.solid.delta_boundary(vtx, mdm->data); } @@ -105,17 +105,17 @@ 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 + ASSERT(mdm && mdm->type == SDIS_SOLID); + return mdm->shader.solid.volumic_power ? mdm->shader.solid.volumic_power(vtx, mdm->data) - : 0; + : SDIS_VOLUMIC_POWER_NONE; } static INLINE double solid_get_temperature (const struct sdis_medium* mdm, const struct sdis_rwalk_vertex* vtx) { - ASSERT(mdm && mdm->type == SDIS_MEDIUM_SOLID); + ASSERT(mdm && mdm->type == SDIS_SOLID); return mdm->shader.solid.temperature(vtx, mdm->data); } diff --git a/src/sdis_scene.c b/src/sdis_scene.c @@ -39,6 +39,65 @@ struct geometry_context { /******************************************************************************* * Helper function ******************************************************************************/ +static void +project_position + (const double V0[3], + const double E0[3], + const double N[3], + const double NxE1[3], + const double rcp_det, + const double pos[3], + double uvw[3]) +{ + double T[3], Q[3], k; + ASSERT(V0 && E0 && N && NxE1 && pos && uvw); + + /* Use Moller/Trumbore intersection test the compute the parametric + * coordinates of the intersection between the triangle and the ray + * `r = pos + N*d' */ + d3_sub(T, pos, V0); + uvw[0] = d3_dot(T, NxE1) * rcp_det; + d3_cross(Q, T, E0); + uvw[1] = d3_dot(Q, N) * rcp_det; + uvw[2] = 1.0 - uvw[0] - uvw[1]; + + if(uvw[0] >= 0 && uvw[1] >= 0 && uvw[2] >= 0) {/* The ray hits the triangle */ + ASSERT(eq_eps(uvw[0] + uvw[1] + uvw[2], 1.0, 1.e-6)); + return; + } + + /* Clamp barycentric coordinates to triangle edges */ + if(uvw[0] >= 0) { + if(uvw[1] >= 0) { + k = 1.0 / (uvw[0] + uvw[1]); + uvw[0] *= k; + uvw[1] *= k; + uvw[2] = 0; + } else if( uvw[2] >= 0) { + k = 1.0 / (uvw[0] + uvw[2]); + uvw[0] *= k; + uvw[1] = 0; + uvw[2] *= k; + } else { + ASSERT(uvw[0] >= 1.f); + d3(uvw, 1, 0, 0); + } + } else if(uvw[1] >= 0) { + if(uvw[2] >= 0) { + k = 1.0 / (uvw[1] + uvw[2]); + uvw[0] = 0; + uvw[1] *= k; + uvw[2] *= k; + } else { + ASSERT(uvw[1] >= 1); + d3(uvw, 0, 1, 0); + } + } else { + ASSERT(uvw[2] >= 1); + d3(uvw, 0, 0, 1); + } +} + /* Check that `hit' roughly lies on a vertex. For segments, a simple but * approximative way is to test that its position have at least one barycentric * coordinate roughly equal to 0 or 1. */ @@ -639,6 +698,69 @@ sdis_scene_get_boundary_position return RES_OK; } +res_T +sdis_scene_boundary_project_position + (const struct sdis_scene* scn, + const size_t iprim, + const double pos[], + double uv[]) +{ + if(!scn || !pos || !uv) 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 a; + double V[2][2]; /* Vertices */ + double E[2][3]; /* V0->V1 and V0->pos */ + double proj; + + /* Retrieve the segment vertices */ + S2D(scene_view_get_primitive(scn->s2d_view, (unsigned int)iprim, &prim)); + S2D(primitive_get_attrib(&prim, S2D_POSITION, 0, &a)); d2_set_f2(V[0], a.value); + S2D(primitive_get_attrib(&prim, S2D_POSITION, 1, &a)); d2_set_f2(V[1], a.value); + + /* Compute the parametric coordinate of the project of `pos' onto the + * segment.*/ + d2_sub(E[0], V[1], V[0]); + d2_normalize(E[0], E[0]); + d2_sub(E[1], pos, V[0]); + proj = d2_dot(E[0], E[1]); + + uv[0] = CLAMP(proj, 0, 1); /* Clamp the parametric coordinate in [0, 1] */ + + } else { + struct s3d_primitive prim; + struct s3d_attrib a; + double V[3][3]; /* Vertices */ + double E[2][3]; /* V0->V1 and V0->V2 edges */ + double N[3]; /* Normal */ + double NxE1[3], rcp_det; /* Muller/Trumboer triangle parameters */ + double uvw[3]; + + S3D(scene_view_get_primitive(scn->s3d_view, (unsigned int)iprim, &prim)); + S3D(triangle_get_vertex_attrib(&prim, 0, S3D_POSITION, &a)); d3_set_f3(V[0], a.value); + S3D(triangle_get_vertex_attrib(&prim, 1, S3D_POSITION, &a)); d3_set_f3(V[1], a.value); + S3D(triangle_get_vertex_attrib(&prim, 2, S3D_POSITION, &a)); d3_set_f3(V[2], a.value); + d3_sub(E[0], V[1], V[0]); + d3_sub(E[1], V[2], V[0]); + d3_cross(N, E[0], E[1]); + + /* Muller/Trumbore triangle parameters */ + d3_cross(NxE1, N, E[1]); + rcp_det = 1.0 / d3_dot(NxE1, E[0]); + + /* Use the Muller/Trumbore intersection test to project `pos' onto the + * triangle and to retrieve the parametric coordinates of the projection + * point */ + project_position(V[0], E[0], N, NxE1, rcp_det, pos, uvw); + + uv[0] = uvw[2]; + uv[1] = uvw[0]; + } + return RES_OK; +} + /******************************************************************************* * Local miscellaneous function ******************************************************************************/ diff --git a/src/sdis_solve.c b/src/sdis_solve.c @@ -273,6 +273,7 @@ sdis_solve_probe_boundary 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 enum sdis_side side, /* Side of iprim on which the probe lies */ const double fp_to_meter, /* Scale from floating point units to meters */ const double Tarad, /* In Kelvin */ const double Tref, /* In Kelvin */ @@ -289,7 +290,7 @@ sdis_solve_probe_boundary res_T res = RES_OK; if(!scn || !nrealisations || !uv || time < 0 || fp_to_meter <= 0 - || Tref < 0 || !out_estimator) { + || Tref < 0 || (side != SDIS_FRONT && side != SDIS_BACK) || !out_estimator) { res = RES_BAD_ARG; goto error; } @@ -363,10 +364,10 @@ sdis_solve_probe_boundary if(scene_is_2d(scn)) { res_local = boundary_realisation_2d - (scn, rng, iprim, uv, time, fp_to_meter, Tarad, Tref, &w); + (scn, rng, iprim, uv, time, side, fp_to_meter, Tarad, Tref, &w); } else { res_local = boundary_realisation_3d - (scn, rng, iprim, uv, time, fp_to_meter, Tarad, Tref, &w); + (scn, rng, iprim, uv, time, side, fp_to_meter, Tarad, Tref, &w); } if(res_local != RES_OK) { if(res_local != RES_BAD_OP) { @@ -448,7 +449,7 @@ sdis_solve_camera res = scene_get_medium(scn, cam->position, &medium); if(res != RES_OK) goto error; - if(medium->type != SDIS_MEDIUM_FLUID) { + if(medium->type != SDIS_FLUID) { log_err(scn->dev, "%s: the camera position `%g %g %g' is not in a fluid.\n", FUNC_NAME, SPLIT3(cam->position)); res = RES_BAD_ARG; diff --git a/src/sdis_solve_Xd.h b/src/sdis_solve_Xd.h @@ -94,9 +94,10 @@ struct XD(rwalk) { struct sdis_rwalk_vertex vtx; /* Position and time of the Random walk */ const struct sdis_medium* mdm; /* Medium in which the random walk lies */ struct sXd(hit) hit; /* Hit of the random walk */ + enum sdis_side hit_side; }; static const struct XD(rwalk) XD(RWALK_NULL) = { - SDIS_RWALK_VERTEX_NULL__, NULL, SXD_HIT_NULL__ + SDIS_RWALK_VERTEX_NULL__, NULL, SXD_HIT_NULL__, SDIS_SIDE_NULL__ }; struct XD(temperature) { @@ -231,17 +232,19 @@ XD(trace_radiative_path) (scn->sXd(view), pos, dir, range, &rwalk->hit, &rwalk->hit)); #endif if(SXD_HIT_NONE(&rwalk->hit)) { /* Fetch the ambient radiative temperature */ + rwalk->hit_side = SDIS_SIDE_NULL__; if(ctx->Tarad >= 0) { T->value += ctx->Tarad; T->done = 1; break; } else { log_err(scn->dev, -"%s: the random walk reaches an invalid ambient radiative temperature of `%gK'\n" -"at position `%g %g %g'. This may be due to numerical inaccuracies or to\n" -"inconsistency in the simulated system (eg: unclosed geometry). For systems\n" -"where the random walks can reach such temperature, one has to setup a valid\n" -"ambient radiative temperature, i.e. it must be greater or equal to 0.\n", + "%s: the random walk reaches an invalid ambient radiative temperature " + "of `%gK' at position `%g %g %g'. This may be due to numerical " + "inaccuracies or to inconsistency in the simulated system (eg: " + "unclosed geometry). For systems where the random walks can reach " + "such temperature, one has to setup a valid ambient radiative " + "temperature, i.e. it must be greater or equal to 0.\n", FUNC_NAME, ctx->Tarad, SPLIT3(rwalk->vtx.P)); @@ -250,16 +253,20 @@ XD(trace_radiative_path) } } + /* Define the hit side */ + rwalk->hit_side = fX(dot)(dir, rwalk->hit.normal) < 0 + ? SDIS_FRONT : SDIS_BACK; + /* Move the random walk to the hit position */ XD(move_pos)(rwalk->vtx.P, dir, rwalk->hit.distance); /* Fetch the new interface and setup the hit fragment */ interf = scene_get_interface(scn, rwalk->hit.prim.prim_id); - XD(setup_interface_fragment)(&frag, &rwalk->vtx, &rwalk->hit); + XD(setup_interface_fragment)(&frag, &rwalk->vtx, &rwalk->hit, rwalk->hit_side); /* Fetch the interface emissivity */ - epsilon = interface_get_emissivity(interf, &frag); - if(epsilon > 1 && epsilon >= 0) { + epsilon = interface_side_get_emissivity(interf, &frag); + if(epsilon > 1 || epsilon < 0) { log_err(scn->dev, "%s: invalid overall emissivity `%g' at position `%g %g %g'.\n", FUNC_NAME, epsilon, SPLIT3(rwalk->vtx.P)); @@ -278,7 +285,7 @@ XD(trace_radiative_path) /* Normalize the normal of the interface and ensure that it points toward the * current medium */ fX(normalize)(N, rwalk->hit.normal); - if(f3_dot(N, dir) > 0) { + if(rwalk->hit_side == SDIS_BACK){ chk_mdm = interf->medium_back; fX(minus)(N, N); } else { @@ -288,10 +295,10 @@ XD(trace_radiative_path) if(chk_mdm != rwalk->mdm) { log_err(scn->dev, "%s: inconsistent medium definition at `%g %g %g'.\n", FUNC_NAME, SPLIT3(rwalk->vtx.P)); - res = RES_BAD_ARG; + res = RES_BAD_OP; goto error; } - alpha = interface_get_specular_fraction(interf, &frag); + alpha = interface_side_get_specular_fraction(interf, &frag); r = ssp_rng_canonical(rng); if(r < alpha) { /* Sample specular part */ reflect(dir, f3_minus(dir, dir), N); @@ -315,8 +322,6 @@ XD(radiative_temperature) struct ssp_rng* rng, struct XD(temperature)* T) { - const struct sdis_interface* interf; - /* The radiative random walk is always perform in 3D. In 2D, the geometry are * assumed to be extruded to the infinty along the Z dimension. */ float N[3] = {0, 0, 0}; @@ -327,13 +332,10 @@ XD(radiative_temperature) ASSERT(!SXD_HIT_NONE(&rwalk->hit)); (void)fp_to_meter; - /* Fetch the current interface */ - interf = scene_get_interface(scn, rwalk->hit.prim.prim_id); - /* Normalize the normal of the interface and ensure that it points toward the * current medium */ fX(normalize(N, rwalk->hit.normal)); - if(interf->medium_back == rwalk->mdm) { + if(rwalk->hit_side == SDIS_BACK) { fX(minus(N, N)); } @@ -362,7 +364,7 @@ XD(fluid_temperature) double tmp; (void)rng, (void)fp_to_meter, (void)ctx; ASSERT(scn && fp_to_meter > 0 && ctx && rwalk && rng && T); - ASSERT(rwalk->mdm->type == SDIS_MEDIUM_FLUID); + ASSERT(rwalk->mdm->type == SDIS_FLUID); tmp = fluid_get_temperature(rwalk->mdm, &rwalk->vtx); if(tmp < 0) { @@ -404,8 +406,8 @@ XD(solid_solid_boundary_temperature) interf = scene_get_interface(scn, rwalk->hit.prim.prim_id); solid_front = interface_get_medium(interf, SDIS_FRONT); solid_back = interface_get_medium(interf, SDIS_BACK); - ASSERT(solid_front->type == SDIS_MEDIUM_SOLID); - ASSERT(solid_back->type == SDIS_MEDIUM_SOLID); + ASSERT(solid_front->type == SDIS_SOLID); + ASSERT(solid_back->type == SDIS_SOLID); /* Fetch the properties of the media */ lambda_front = solid_get_thermal_conductivity(solid_front, &rwalk->vtx); @@ -442,6 +444,8 @@ XD(solid_solid_boundary_temperature) /* Switch in solid random walk */ T->func = XD(solid_temperature); + rwalk->hit = SXD_HIT_NULL; + rwalk->hit_side = SDIS_SIDE_NULL__; } static void @@ -459,6 +463,7 @@ XD(solid_fluid_boundary_temperature) const struct sdis_medium* mdm_back = NULL; const struct sdis_medium* solid = NULL; const struct sdis_medium* fluid = NULL; + struct sdis_interface_fragment frag_fluid; double hc; double hr; double epsilon; /* Interface emissivity */ @@ -478,12 +483,16 @@ XD(solid_fluid_boundary_temperature) mdm_front = interface_get_medium(interf, SDIS_FRONT); mdm_back = interface_get_medium(interf, SDIS_BACK); ASSERT(mdm_front->type != mdm_back->type); - if(mdm_front->type == SDIS_MEDIUM_SOLID) { + + frag_fluid = *frag; + if(mdm_front->type == SDIS_SOLID) { solid = mdm_front; fluid = mdm_back; + frag_fluid.side = SDIS_BACK; } else { solid = mdm_back; fluid = mdm_front; + frag_fluid.side = SDIS_FRONT; } /* Fetch the solid properties */ @@ -491,7 +500,7 @@ XD(solid_fluid_boundary_temperature) delta_boundary = solid_get_delta_boundary(solid, &rwalk->vtx); /* Fetch the boundary properties */ - epsilon = interface_get_emissivity(interf, frag); + epsilon = interface_side_get_emissivity(interf, &frag_fluid); hc = interface_get_convection_coef(interf, frag); /* Compute the radiative coefficient */ @@ -505,11 +514,13 @@ XD(solid_fluid_boundary_temperature) r = ssp_rng_canonical(rng); if(r < radia_proba) { /* Switch in radiative random walk */ - rwalk->mdm = fluid; T->func = XD(radiative_temperature); - } else if(r < fluid_proba + radia_proba) { /* Switch to fluid random walk */ rwalk->mdm = fluid; + rwalk->hit_side = rwalk->mdm == mdm_front ? SDIS_FRONT : SDIS_BACK; + } else if(r < fluid_proba + radia_proba) { /* Switch to fluid random walk */ T->func = XD(fluid_temperature); + rwalk->mdm = fluid; + rwalk->hit_side = rwalk->mdm == mdm_front ? SDIS_FRONT : SDIS_BACK; } else { /* Solid random walk */ rwalk->mdm = solid; fX(normalize)(dir, rwalk->hit.normal); @@ -525,6 +536,8 @@ XD(solid_fluid_boundary_temperature) /* Switch in solid random walk */ T->func = XD(solid_temperature); + rwalk->hit = SXD_HIT_NULL; + rwalk->hit_side = SDIS_SIDE_NULL__; } } @@ -541,24 +554,62 @@ XD(boundary_temperature) const struct sdis_interface* interf = NULL; const struct sdis_medium* mdm_front = NULL; const struct sdis_medium* mdm_back = NULL; + const struct sdis_medium* mdm = 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); + XD(setup_interface_fragment)(&frag, &rwalk->vtx, &rwalk->hit, rwalk->hit_side); /* Retrieve the current interface */ interf = scene_get_interface(scn, rwalk->hit.prim.prim_id); - /* Check if the boundary condition is known */ - tmp = interface_get_temperature(interf, &frag); + /* Check if the boundary temperature is known */ + tmp = interface_side_get_temperature(interf, &frag); if(tmp >= 0) { T->value += tmp; T->done = 1; return RES_OK; } + /* Check if the boundary flux is known. Note that actually, only solid media + * can have a flux as limit condition */ + mdm = interface_get_medium(interf, frag.side); + if(sdis_medium_get_type(mdm) == SDIS_SOLID) { + const double phi = interface_side_get_flux(interf, &frag); + + if(phi != SDIS_FLUX_NONE) { + double lambda = solid_get_thermal_conductivity(mdm, &rwalk->vtx); + double delta_b = solid_get_delta_boundary(mdm, &rwalk->vtx); + double delta_b_in_meter = delta_b * fp_to_meter; + float pos[3]; + float range[2]; + float dir[3]; + + /* Update the temperature */ + T->value += phi * delta_b_in_meter / lambda; + + /* Ensuure that the normal points toward the solid */ + fX(normalize)(dir, rwalk->hit.normal); + if(frag.side == SDIS_BACK) fX(minus)(dir, dir); + + /* "Reinject" the random walk into the solid */ + fX_set_dX(pos, rwalk->vtx.P); + range[0] = 0, range[1] = (float)delta_b*RAY_RANGE_MAX_SCALE; + SXD(scene_view_trace_ray + (scn->sXd(view), pos, dir, range, &rwalk->hit, &rwalk->hit)); + if(!SXD_HIT_NONE(&rwalk->hit)) delta_b = rwalk->hit.distance * 0.5; + XD(move_pos)(rwalk->vtx.P, dir, (float)delta_b); + + /* Switch in solid random walk */ + T->func = XD(solid_temperature); + rwalk->hit = SXD_HIT_NULL; + rwalk->hit_side = SDIS_SIDE_NULL__; + rwalk->mdm = mdm; + return RES_OK; + } + } mdm_front = interface_get_medium(interf, SDIS_FRONT); mdm_back = interface_get_medium(interf, SDIS_BACK); @@ -584,7 +635,7 @@ XD(solid_temperature) double position_start[DIM]; const struct sdis_medium* mdm; ASSERT(scn && fp_to_meter > 0 && rwalk && rng && T); - ASSERT(rwalk->mdm->type == SDIS_MEDIUM_SOLID); + ASSERT(rwalk->mdm->type == SDIS_SOLID); (void)ctx; /* Check the random walk consistency */ @@ -665,14 +716,20 @@ XD(solid_temperature) /* Add the volumic power density to the measured temperature */ power = solid_get_volumic_power(mdm, &rwalk->vtx); - if(power > 0) { + if(power != SDIS_VOLUMIC_POWER_NONE) { 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; + if(hit0.distance > delta) { + rwalk->hit = SXD_HIT_NULL; + rwalk->hit_side = SDIS_SIDE_NULL__; + } else { + rwalk->hit = hit0; + rwalk->hit_side = fX(dot)(hit0.normal, dir0) < 0 ? SDIS_FRONT : SDIS_BACK; + } /* Update the random walk position */ XD(move_pos)(rwalk->vtx.P, dir0, delta); @@ -683,9 +740,7 @@ XD(solid_temperature) } else { const struct sdis_interface* interf; interf = scene_get_interface(scn, rwalk->hit.prim.prim_id); - mdm = interface_get_medium - (interf, - fX(dot(rwalk->hit.normal, dir0)) < 0 ? SDIS_FRONT : SDIS_BACK); + mdm = interface_get_medium(interf, rwalk->hit_side); } /* Check random walk consistency */ @@ -766,8 +821,8 @@ XD(probe_realisation) ASSERT(medium && position && fp_to_meter > 0 && weight && time >= 0); switch(medium->type) { - case SDIS_MEDIUM_FLUID: T.func = XD(fluid_temperature); break; - case SDIS_MEDIUM_SOLID: T.func = XD(solid_temperature); break; + case SDIS_FLUID: T.func = XD(fluid_temperature); break; + case SDIS_SOLID: T.func = XD(solid_temperature); break; default: FATAL("Unreachable code\n"); break; } @@ -796,6 +851,7 @@ XD(boundary_realisation) const size_t iprim, const double uv[DIM], const double time, + const enum sdis_side side, const double fp_to_meter, const double Tarad, const double Tref, @@ -815,6 +871,7 @@ XD(boundary_realisation) T.func = XD(boundary_temperature); + rwalk.hit_side = side; rwalk.hit.distance = 0; rwalk.vtx.time = time; rwalk.mdm = NULL; /* The random walk is at an interface between 2 media */ @@ -873,11 +930,12 @@ XD(ray_realisation) float dir[3]; res_T res = RES_OK; ASSERT(scn && position && direction && time>=0 && fp_to_meter>0 && weight); - ASSERT(medium && medium->type == SDIS_MEDIUM_FLUID); + ASSERT(medium && medium->type == SDIS_FLUID); dX(set)(rwalk.vtx.P, position); rwalk.vtx.time = time; rwalk.hit = SXD_HIT_NULL; + rwalk.hit_side = SDIS_SIDE_NULL__; rwalk.mdm = medium; ctx.Tarad = Tarad; diff --git a/src/test_sdis_conducto_radiative.c b/src/test_sdis_conducto_radiative.c @@ -229,15 +229,24 @@ create_interface const struct interface* interf, struct sdis_interface** out_interf) { - struct sdis_interface_shader shader = DUMMY_INTERFACE_SHADER; + struct sdis_interface_shader shader = SDIS_INTERFACE_SHADER_NULL; struct sdis_data* data = NULL; CHK(interf != NULL); - shader.temperature = interface_get_temperature; - shader.convection_coef = interface_get_convection_coef; - shader.emissivity = interface_get_emissivity; - shader.specular_fraction = interface_get_specular_fraction; + shader.front.temperature = interface_get_temperature; + shader.back.temperature = interface_get_temperature; + if(sdis_medium_get_type(front) != sdis_medium_get_type(back)) { + shader.convection_coef = interface_get_convection_coef; + } + if(sdis_medium_get_type(front) == SDIS_FLUID) { + shader.front.emissivity = interface_get_emissivity; + shader.front.specular_fraction = interface_get_specular_fraction; + } + if(sdis_medium_get_type(back) == SDIS_FLUID) { + shader.back.emissivity = interface_get_emissivity; + shader.back.specular_fraction = interface_get_specular_fraction; + } CHK(sdis_data_create(dev, sizeof(struct interface), ALIGNOF(struct interface), NULL, &data) == RES_OK); diff --git a/src/test_sdis_conducto_radiative_2d.c b/src/test_sdis_conducto_radiative_2d.c @@ -161,42 +161,74 @@ solid_get_delta_boundary * Interface ******************************************************************************/ struct interface { - double temperature; double convection_coef; - double emissivity; - double specular_fraction; + struct { + double temperature; + double emissivity; + double specular_fraction; + } front, back; +}; + +static const struct interface INTERFACE_NULL = { + 0, {-1, -1, -1}, {-1, -1, -1} }; static double interface_get_temperature (const struct sdis_interface_fragment* frag, struct sdis_data* data) { + const struct interface* interf; + double T = -1; CHK(data != NULL && frag != NULL); - return ((const struct interface*)sdis_data_cget(data))->temperature; + interf = sdis_data_cget(data); + switch(frag->side) { + case SDIS_FRONT: T = interf->front.temperature; break; + case SDIS_BACK: T = interf->back.temperature; break; + default: FATAL("Unreachable code.\n"); break; + } + return T; } static double interface_get_convection_coef (const struct sdis_interface_fragment* frag, struct sdis_data* data) { + const struct interface* interf; CHK(data != NULL && frag != NULL); - return ((const struct interface*)sdis_data_cget(data))->convection_coef; + interf = sdis_data_cget(data); + return interf->convection_coef; } static double interface_get_emissivity (const struct sdis_interface_fragment* frag, struct sdis_data* data) { + const struct interface* interf; + double e = -1; CHK(data != NULL && frag != NULL); - return ((const struct interface*)sdis_data_cget(data))->emissivity; + interf = sdis_data_cget(data); + switch(frag->side) { + case SDIS_FRONT: e = interf->front.emissivity; break; + case SDIS_BACK: e = interf->back.emissivity; break; + default: FATAL("Unreachable code.\n"); break; + } + return e; } static double interface_get_specular_fraction (const struct sdis_interface_fragment* frag, struct sdis_data* data) { + const struct interface* interf; + double f = -1; CHK(data != NULL && frag != NULL); - return ((const struct interface*)sdis_data_cget(data))->specular_fraction; + interf = sdis_data_cget(data); + switch(frag->side) { + case SDIS_FRONT: f = interf->front.specular_fraction; break; + case SDIS_BACK: f = interf->back.specular_fraction; break; + default: FATAL("Unreachable code.\n"); break; + } + return f; } /******************************************************************************* @@ -210,16 +242,27 @@ create_interface const struct interface* interf, struct sdis_interface** out_interf) { - struct sdis_interface_shader shader = DUMMY_INTERFACE_SHADER; + struct sdis_interface_shader shader = SDIS_INTERFACE_SHADER_NULL; struct sdis_data* data = NULL; + const enum sdis_medium_type type_f = sdis_medium_get_type(front); + const enum sdis_medium_type type_b = sdis_medium_get_type(back); CHK(interf != NULL); - shader.temperature = interface_get_temperature; - shader.convection_coef = interface_get_convection_coef; - shader.emissivity = interface_get_emissivity; - shader.specular_fraction = interface_get_specular_fraction; + shader.back.temperature = interface_get_temperature; + shader.front.temperature = interface_get_temperature; + if(type_f != type_b) { + shader.convection_coef = interface_get_convection_coef; + } + if(type_f == SDIS_FLUID) { + shader.front.emissivity = interface_get_emissivity; + shader.front.specular_fraction = interface_get_specular_fraction; + } + if(type_b == SDIS_FLUID) { + shader.back.emissivity = interface_get_emissivity; + shader.back.specular_fraction = interface_get_specular_fraction; + } CHK(sdis_data_create(dev, sizeof(struct interface), ALIGNOF(struct interface), NULL, &data) == RES_OK); *((struct interface*)sdis_data_get(data)) = *interf; @@ -228,7 +271,6 @@ create_interface CHK(sdis_data_ref_put(data) == RES_OK); } - /******************************************************************************* * Test ******************************************************************************/ @@ -293,38 +335,35 @@ main(int argc, char** argv) CHK(sdis_data_ref_put(data) == RES_OK); /* Create the interface that forces to keep in conduction */ - interf.temperature = UNKNOWN_TEMPERATURE; - interf.convection_coef = -1; - interf.emissivity = -1; - interf.specular_fraction = -1; + interf = INTERFACE_NULL; create_interface(dev, solid, solid2, &interf, interfaces+0); /* Create the interface that emits radiative heat from the solid */ - interf.temperature = UNKNOWN_TEMPERATURE; - interf.convection_coef = 0; - interf.emissivity = emissivity; - interf.specular_fraction = -1; + interf = INTERFACE_NULL; + interf.back.temperature = UNKNOWN_TEMPERATURE; + interf.back.emissivity = emissivity; + interf.back.specular_fraction = -1; /* Should not be fetched */ create_interface(dev, solid, fluid, &interf, interfaces+1); /* Create the interface that forces the radiative heat to bounce */ - interf.temperature = UNKNOWN_TEMPERATURE; - interf.convection_coef = 0; - interf.emissivity = 0; - interf.specular_fraction = 1; + interf = INTERFACE_NULL; + interf.front.temperature = UNKNOWN_TEMPERATURE; + interf.front.emissivity = 0; + interf.front.specular_fraction = 1; create_interface(dev, fluid, solid2, &interf, interfaces+2); /* Create the interface with a limit condition of T0 Kelvin */ - interf.temperature = T0; - interf.convection_coef = 0; - interf.emissivity = 1; - interf.specular_fraction = 1; + interf = INTERFACE_NULL; + interf.front.temperature = T0; + interf.front.emissivity = 1; + interf.front.specular_fraction = 1; create_interface(dev, fluid, solid2, &interf, interfaces+3); /* Create the interface with a limit condition of T1 Kelvin */ - interf.temperature = T1; - interf.convection_coef = 0; - interf.emissivity = 1; - interf.specular_fraction = 1; + interf = INTERFACE_NULL; + interf.front.temperature = T1; + interf.front.emissivity = 1; + interf.front.specular_fraction = 1; create_interface(dev, fluid, solid2, &interf, interfaces+4); /* Setup the per primitive interface of the solid medium */ diff --git a/src/test_sdis_flux.c b/src/test_sdis_flux.c @@ -0,0 +1,342 @@ +/* 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> + +/* + * The scene is composed of a solid cube/square whose temperature is unknown. + * The temperature is fixed at T0 on the +X face. The Flux of the -X face is + * fixed to PHI. The flux on the other faces is null (i.e. adiabatic). This + * test computes the temperature of a probe position pos into the solid and + * check that it is equal to: + * + * T(pos) = T0 + (A-pos) * PHI/LAMBDA + * + * with LAMBDA the conductivity of the solid and A the size of cube/square. + * + * 3D 2D + * + * ///// (1,1,1) ///// (1,1) + * +-------+ +-------+ + * /' /| | | + * +-------+ T0 PHI T0 + * PHI +.....|.+ | | + * |, |/ +-------+ + * +-------+ (0,0) ///// + * (0,0,0) ///// + */ + +#define UNKNOWN_TEMPERATURE -1 +#define N 10000 + +#define PHI 10.0 +#define T0 320.0 +#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 +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 phi; +}; + +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_flux + (const struct sdis_interface_fragment* frag, struct sdis_data* data) +{ + const struct interf* interf = sdis_data_cget(data); + CHK(frag && data); + return interf->phi; +} + +/******************************************************************************* + * 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_interface* interf_phi = NULL; + struct sdis_scene* box_scn = NULL; + struct sdis_scene* square_scn = NULL; + struct sdis_estimator* estimator = NULL; + struct sdis_fluid_shader fluid_shader = DUMMY_FLUID_SHADER; + struct sdis_solid_shader solid_shader = DUMMY_SOLID_SHADER; + struct sdis_interface_shader interf_shader = SDIS_INTERFACE_SHADER_NULL; + struct sdis_interface* box_interfaces[12 /*#triangles*/]; + struct sdis_interface* square_interfaces[4/*#segments*/]; + struct interf* interf_props = NULL; + double pos[3]; + 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 dummy fluid medium */ + 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.front.temperature = interface_get_temperature; + interf_shader.front.flux = interface_get_flux; + + /* 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; + interf_props->phi = 0; + 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; + interf_props->phi = 0; /* Unused */ + CHK(sdis_interface_create + (dev, solid, fluid, &interf_shader, data, &interf_T0) == RES_OK); + CHK(sdis_data_ref_put(data) == RES_OK); + + /* Create the PHI 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; + interf_props->phi = PHI; + CHK(sdis_interface_create + (dev, solid, fluid, &interf_shader, data, &interf_phi) == 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_phi; /* 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_phi; /* Left */ + 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); + CHK(sdis_interface_ref_put(interf_phi) == RES_OK); + + d3_splat(pos, 0.25); + ref = T0 + (1 - pos[0]) * PHI/LAMBDA; + + /* 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; + +} diff --git a/src/test_sdis_interface.c b/src/test_sdis_interface.c @@ -31,11 +31,14 @@ main(int argc, char** 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); + (NULL, &allocator, SDIS_NTHREADS_DEFAULT, 1, &dev) == RES_OK); CHK(sdis_fluid_create(dev, &fluid_shader, NULL, &fluid) == RES_OK); CHK(sdis_solid_create(dev, &solid_shader, NULL, &solid) == RES_OK); + + shader = SDIS_INTERFACE_SHADER_NULL; + #define CREATE sdis_interface_create CHK(CREATE(NULL, NULL, NULL, NULL, NULL, NULL) == RES_BAD_ARG); CHK(CREATE(dev, NULL, NULL, NULL, NULL, NULL) == RES_BAD_ARG); @@ -78,24 +81,28 @@ main(int argc, char** argv) CHK(CREATE(dev, solid, solid, &shader, NULL, &interf) == RES_OK); CHK(sdis_interface_ref_put(interf) == RES_OK); - shader.convection_coef = NULL; - shader.specular_fraction = NULL; - shader.emissivity = NULL; + shader = SDIS_INTERFACE_SHADER_NULL; CHK(CREATE(dev, solid, solid, &shader, NULL, &interf) == RES_OK); CHK(sdis_interface_ref_put(interf) == RES_OK); - shader.temperature = NULL; + shader.front.temperature = dummy_interface_getter; CHK(CREATE(dev, solid, solid, &shader, NULL, &interf) == RES_OK); CHK(sdis_interface_ref_put(interf) == RES_OK); - CHK(CREATE(dev, solid, fluid, &shader, NULL, &interf) == RES_BAD_ARG); - shader.convection_coef = DUMMY_INTERFACE_SHADER.convection_coef; - CHK(CREATE(dev, solid, fluid, &shader, NULL, &interf) == RES_BAD_ARG); - shader.emissivity = DUMMY_INTERFACE_SHADER.emissivity; - CHK(CREATE(dev, solid, fluid, &shader, NULL, &interf) == RES_BAD_ARG); - shader.specular_fraction = DUMMY_INTERFACE_SHADER.specular_fraction; + shader.back.emissivity = dummy_interface_getter; + CHK(CREATE(dev, solid, fluid, &shader, NULL, &interf) == RES_OK); + CHK(sdis_interface_ref_put(interf) == RES_OK); + shader.back.specular_fraction = dummy_interface_getter; CHK(CREATE(dev, solid, fluid, &shader, NULL, &interf) == RES_OK); CHK(sdis_interface_ref_put(interf) == RES_OK); + shader.back = SDIS_INTERFACE_SIDE_SHADER_NULL; + shader.front.emissivity = dummy_interface_getter; + CHK(CREATE(dev, solid, fluid, &shader, NULL, &interf) == RES_OK); /* Warning */ + CHK(sdis_interface_ref_put(interf) == RES_OK); + shader.front.emissivity = NULL; + shader.front.specular_fraction = dummy_interface_getter; + CHK(CREATE(dev, solid, fluid, &shader, NULL, &interf) == RES_OK); /* Warning */ + CHK(sdis_interface_ref_put(interf) == RES_OK); #undef CREATE CHK(sdis_device_ref_put(dev) == RES_OK); diff --git a/src/test_sdis_scene.c b/src/test_sdis_scene.c @@ -15,6 +15,9 @@ #include "sdis.h" #include "test_sdis_utils.h" + +#include <rsys/double2.h> +#include <rsys/double3.h> #include <rsys/math.h> struct context { @@ -23,8 +26,14 @@ struct context { struct sdis_interface* interf; }; +static INLINE double +rand_canonic() +{ + return (double)rand()/(double)((double)RAND_MAX+1); +} + static INLINE void -get_indices(const size_t itri, size_t ids[3], void* context) +get_indices_3d(const size_t itri, size_t ids[3], void* context) { struct context* ctx = context; CHK(ctx != NULL); @@ -35,7 +44,17 @@ get_indices(const size_t itri, size_t ids[3], void* context) } static INLINE void -get_position(const size_t ivert, double pos[3], void* context) +get_indices_2d(const size_t iseg, size_t ids[2], void* context) +{ + struct context* ctx = context; + CHK(ctx != NULL); + CHK(iseg < square_nsegments); + ids[0] = ctx->indices[iseg*2+0]; + ids[1] = ctx->indices[iseg*2+1]; +} + +static INLINE void +get_position_3d(const size_t ivert, double pos[3], void* context) { struct context* ctx = context; CHK(ctx != NULL); @@ -46,6 +65,16 @@ get_position(const size_t ivert, double pos[3], void* context) } static INLINE void +get_position_2d(const size_t ivert, double pos[2], void* context) +{ + struct context* ctx = context; + CHK(ctx != NULL); + CHK(ivert < square_nvertices); + pos[0] = ctx->positions[ivert*2+0]; + pos[1] = ctx->positions[ivert*2+1]; +} + +static INLINE void get_interface(const size_t itri, struct sdis_interface** bound, void* context) { struct context* ctx = context; @@ -55,30 +84,15 @@ get_interface(const size_t itri, struct sdis_interface** bound, void* context) *bound = ctx->interf; } -int -main(int argc, char** argv) +static void +test_scene_3d(struct sdis_device* dev, struct sdis_interface* interf) { - struct mem_allocator allocator; - struct sdis_device* dev = NULL; - struct sdis_medium* solid = NULL; - struct sdis_medium* fluid = NULL; - struct sdis_interface* interf = NULL; struct sdis_scene* scn = NULL; - struct sdis_fluid_shader fluid_shader = DUMMY_FLUID_SHADER; - struct sdis_solid_shader solid_shader = DUMMY_SOLID_SHADER; - struct sdis_interface_shader interface_shader = DUMMY_INTERFACE_SHADER; double lower[3], upper[3]; + double uv0[2], uv1[2], pos[3], pos1[3]; struct context ctx; size_t ntris, npos; - (void)argc, (void)argv; - - CHK(mem_init_proxy_allocator(&allocator, &mem_default_allocator) == RES_OK); - CHK(sdis_device_create(NULL, &allocator, 1, 0, &dev) == RES_OK); - - CHK(sdis_fluid_create(dev, &fluid_shader, NULL, &fluid) == RES_OK); - CHK(sdis_solid_create(dev, &solid_shader, NULL, &solid) == RES_OK); - CHK(sdis_interface_create - (dev, solid, fluid, &interface_shader, NULL, &interf) == RES_OK); + size_t i; ctx.positions = box_vertices; ctx.indices = box_indices; @@ -87,8 +101,8 @@ main(int argc, char** argv) npos = box_nvertices; #define CREATE sdis_scene_create - #define IDS get_indices - #define POS get_position + #define IDS get_indices_3d + #define POS get_position_3d #define IFA get_interface CHK(CREATE(NULL, 0, NULL, NULL, 0, NULL, &ctx, NULL) == RES_BAD_ARG); @@ -115,17 +129,160 @@ main(int argc, char** argv) CHK(eq_eps(upper[1], 1, 1.e-6)); CHK(eq_eps(upper[2], 1, 1.e-6)); + uv0[0] = 0.3; + uv0[1] = 0.3; + + CHK(sdis_scene_get_boundary_position(NULL, 6, uv0, pos) == RES_BAD_ARG); + CHK(sdis_scene_get_boundary_position(scn, 12, uv0, pos) == RES_BAD_ARG); + CHK(sdis_scene_get_boundary_position(scn, 6, NULL, pos) == RES_BAD_ARG); + CHK(sdis_scene_get_boundary_position(scn, 6, uv0, NULL) == RES_BAD_ARG); + CHK(sdis_scene_get_boundary_position(scn, 6, uv0, pos) == RES_OK); + + CHK(sdis_scene_boundary_project_position(NULL, 6, pos, uv1) == RES_BAD_ARG); + CHK(sdis_scene_boundary_project_position(scn, 12, pos, uv1) == RES_BAD_ARG); + CHK(sdis_scene_boundary_project_position(scn, 6, NULL, uv1) == RES_BAD_ARG); + CHK(sdis_scene_boundary_project_position(scn, 6, pos, NULL) == RES_BAD_ARG); + CHK(sdis_scene_boundary_project_position(scn, 6, pos, uv1) == RES_OK); + + CHK(d2_eq_eps(uv0, uv1, 1.e-6)); + + FOR_EACH(i, 0, 64) { + uv0[0] = rand_canonic(); + uv0[1] = rand_canonic() * (1 - uv0[0]); + + CHK(sdis_scene_get_boundary_position(scn, 4, uv0, pos) == RES_OK); + CHK(sdis_scene_boundary_project_position(scn, 4, pos, uv1) == RES_OK); + CHK(d2_eq_eps(uv0, uv1, 1.e-6)); + } + + pos[0] = 10; + pos[1] = 0.1; + pos[2] = 0.5; + CHK(sdis_scene_boundary_project_position(scn, 6, pos, uv1) == RES_OK); + CHK(sdis_scene_get_boundary_position(scn, 6, uv1, pos1) == RES_OK); + CHK(!d3_eq_eps(pos1, pos, 1.e-6)); + CHK(sdis_scene_ref_get(NULL) == RES_BAD_ARG); CHK(sdis_scene_ref_get(scn) == RES_OK); CHK(sdis_scene_ref_put(NULL) == RES_BAD_ARG); CHK(sdis_scene_ref_put(scn) == RES_OK); CHK(sdis_scene_ref_put(scn) == RES_OK); +} + +static void +test_scene_2d(struct sdis_device* dev, struct sdis_interface* interf) +{ + struct sdis_scene* scn = NULL; + double lower[2], upper[2]; + double u0, u1, pos[2]; + struct context ctx; + size_t nsegs, npos; + size_t i; + + ctx.positions = square_vertices; + ctx.indices = square_indices; + ctx.interf = interf; + nsegs = square_nsegments; + npos = square_nvertices; + + #define CREATE sdis_scene_2d_create + #define IDS get_indices_2d + #define POS get_position_2d + #define IFA get_interface + + CHK(CREATE(NULL, 0, NULL, NULL, 0, NULL, &ctx, NULL) == RES_BAD_ARG); + CHK(CREATE(dev, 0, IDS, IFA, npos, POS, &ctx, &scn) == RES_BAD_ARG); + CHK(CREATE(dev, nsegs, NULL, IFA, npos, POS, &ctx, &scn) == RES_BAD_ARG); + CHK(CREATE(dev, nsegs, IDS, NULL, npos, POS, &ctx, &scn) == RES_BAD_ARG); + CHK(CREATE(dev, nsegs, IDS, IFA, 0, POS, &ctx, &scn) == RES_BAD_ARG); + CHK(CREATE(dev, nsegs, IDS, IFA, npos, NULL, &ctx, &scn) == RES_BAD_ARG); + CHK(CREATE(dev, nsegs, IDS, IFA, npos, POS, &ctx, &scn) == RES_OK); + + #undef CREATE + #undef IDS + #undef POS + #undef IFA + + CHK(sdis_scene_get_aabb(NULL, lower, upper) == RES_BAD_ARG); + CHK(sdis_scene_get_aabb(scn, NULL, upper) == RES_BAD_ARG); + CHK(sdis_scene_get_aabb(scn, lower, NULL) == RES_BAD_ARG); + CHK(sdis_scene_get_aabb(scn, lower, upper) == RES_OK); + CHK(eq_eps(lower[0], 0, 1.e-6)); + CHK(eq_eps(lower[1], 0, 1.e-6)); + CHK(eq_eps(upper[0], 1, 1.e-6)); + CHK(eq_eps(upper[1], 1, 1.e-6)); + + u0 = 0.5; + + CHK(sdis_scene_get_boundary_position(NULL, 1, &u0, pos) == RES_BAD_ARG); + CHK(sdis_scene_get_boundary_position(scn, 4, &u0, pos) == RES_BAD_ARG); + CHK(sdis_scene_get_boundary_position(scn, 1, NULL, pos) == RES_BAD_ARG); + CHK(sdis_scene_get_boundary_position(scn, 1, &u0, NULL) == RES_BAD_ARG); + CHK(sdis_scene_get_boundary_position(scn, 1, &u0, pos) == RES_OK); + + CHK(sdis_scene_boundary_project_position(NULL, 1, pos, &u1) == RES_BAD_ARG); + CHK(sdis_scene_boundary_project_position(scn, 4, pos, &u1) == RES_BAD_ARG); + CHK(sdis_scene_boundary_project_position(scn, 1, NULL, &u1) == RES_BAD_ARG); + CHK(sdis_scene_boundary_project_position(scn, 1, pos, NULL) == RES_BAD_ARG); + CHK(sdis_scene_boundary_project_position(scn, 1, pos, &u1) == RES_OK); + + CHK(eq_eps(u0, u1, 1.e-6)); + + FOR_EACH(i, 0, 64) { + u0 = rand_canonic(); + + CHK(sdis_scene_get_boundary_position(scn, 2, &u0, pos) == RES_OK); + CHK(sdis_scene_boundary_project_position(scn, 2, pos, &u1) == RES_OK); + CHK(eq_eps(u0, u1, 1.e-6)); + } + + d2(pos, 5, 0.5); + CHK(sdis_scene_boundary_project_position(scn, 3, pos, &u0) == RES_OK); + CHK(eq_eps(u0, 0.5, 1.e-6)); + + d2(pos, 1, 2); + CHK(sdis_scene_boundary_project_position(scn, 3, pos, &u0) == RES_OK); + CHK(eq_eps(u0, 0, 1.e-6)); + + d2(pos, 1, -1); + CHK(sdis_scene_boundary_project_position(scn, 3, pos, &u0) == RES_OK); + CHK(eq_eps(u0, 1, 1.e-6)); + + CHK(sdis_scene_ref_put(scn) == RES_OK); +} + +int +main(int argc, char** argv) +{ + struct mem_allocator allocator; + struct sdis_device* dev = NULL; + struct sdis_medium* solid = NULL; + struct sdis_medium* fluid = NULL; + struct sdis_interface* interf = NULL; + struct sdis_fluid_shader fluid_shader = DUMMY_FLUID_SHADER; + struct sdis_solid_shader solid_shader = DUMMY_SOLID_SHADER; + struct sdis_interface_shader interface_shader = SDIS_INTERFACE_SHADER_NULL; + (void)argc, (void)argv; + + interface_shader.convection_coef = DUMMY_INTERFACE_SHADER.convection_coef; + + CHK(mem_init_proxy_allocator(&allocator, &mem_default_allocator) == RES_OK); + CHK(sdis_device_create(NULL, &allocator, 1, 0, &dev) == RES_OK); + + CHK(sdis_fluid_create(dev, &fluid_shader, NULL, &fluid) == RES_OK); + CHK(sdis_solid_create(dev, &solid_shader, NULL, &solid) == RES_OK); + CHK(sdis_interface_create + (dev, solid, fluid, &interface_shader, NULL, &interf) == RES_OK); - CHK(sdis_device_ref_put(dev) == RES_OK); - CHK(sdis_interface_ref_put(interf) == RES_OK); CHK(sdis_medium_ref_put(solid) == RES_OK); CHK(sdis_medium_ref_put(fluid) == RES_OK); + test_scene_3d(dev, interf); + test_scene_2d(dev, interf); + + CHK(sdis_device_ref_put(dev) == RES_OK); + CHK(sdis_interface_ref_put(interf) == RES_OK); + check_memory_allocator(&allocator); mem_shutdown_proxy_allocator(&allocator); CHK(mem_allocated_size() == 0); diff --git a/src/test_sdis_solve_camera.c b/src/test_sdis_solve_camera.c @@ -361,7 +361,7 @@ create_interface { struct sdis_data* data = NULL; struct interf* interface_param = NULL; - struct sdis_interface_shader interface_shader = DUMMY_INTERFACE_SHADER; + struct sdis_interface_shader interface_shader = SDIS_INTERFACE_SHADER_NULL; CHK(mdm_front != NULL); CHK(mdm_back != NULL); @@ -376,10 +376,16 @@ create_interface /* Setup the interface shader */ interface_shader.convection_coef = interface_get_convection_coef; - interface_shader.temperature = interface_get_temperature; - interface_shader.emissivity = interface_get_emissivity; - interface_shader.specular_fraction = interface_get_specular_fraction; - + interface_shader.front.temperature = interface_get_temperature; + interface_shader.back.temperature = interface_get_temperature; + if(sdis_medium_get_type(mdm_front) == SDIS_FLUID) { + interface_shader.front.emissivity = interface_get_emissivity; + interface_shader.front.specular_fraction = interface_get_specular_fraction; + } + if(sdis_medium_get_type(mdm_back) == SDIS_FLUID) { + interface_shader.back.emissivity = interface_get_emissivity; + interface_shader.back.specular_fraction = interface_get_specular_fraction; + } /* Create the interface */ CHK(sdis_interface_create (dev, mdm_front, mdm_back, &interface_shader, data, interf) == RES_OK); diff --git a/src/test_sdis_solve_probe.c b/src/test_sdis_solve_probe.c @@ -191,7 +191,7 @@ main(int argc, char** argv) 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 interface_shader = DUMMY_INTERFACE_SHADER; + struct sdis_interface_shader interface_shader = SDIS_INTERFACE_SHADER_NULL; struct context ctx; struct fluid* fluid_param; struct solid* solid_param; @@ -206,7 +206,7 @@ main(int argc, char** 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); + (NULL, &allocator, SDIS_NTHREADS_DEFAULT, 1, &dev) == RES_OK); /* Create the fluid medium */ CHK(sdis_data_create @@ -243,9 +243,10 @@ main(int argc, char** argv) interface_param->epsilon = 0; interface_param->specular_fraction = 0; interface_shader.convection_coef = interface_get_convection_coef; - interface_shader.temperature = NULL; - interface_shader.emissivity = interface_get_emissivity; - interface_shader.specular_fraction = interface_get_specular_fraction; + interface_shader.front = SDIS_INTERFACE_SIDE_SHADER_NULL; + interface_shader.back.temperature = NULL; + interface_shader.back.emissivity = interface_get_emissivity; + interface_shader.back.specular_fraction = interface_get_specular_fraction; CHK(sdis_interface_create (dev, solid, fluid, &interface_shader, data, &interf) == RES_OK); CHK(sdis_data_ref_put(data) == RES_OK); diff --git a/src/test_sdis_solve_probe2.c b/src/test_sdis_solve_probe2.c @@ -197,9 +197,8 @@ main(int argc, char** argv) /* Create the fluid/solid interface with no limit conidition */ interface_shader.convection_coef = null_interface_value; - interface_shader.temperature = NULL; - interface_shader.emissivity = null_interface_value; - interface_shader.specular_fraction = null_interface_value; + interface_shader.front = SDIS_INTERFACE_SIDE_SHADER_NULL; + interface_shader.back = SDIS_INTERFACE_SIDE_SHADER_NULL; CHK(sdis_interface_create (dev, solid, fluid, &interface_shader, NULL, &Tnone) == RES_OK); @@ -208,10 +207,7 @@ main(int argc, char** argv) ALIGNOF(struct interf), NULL, &data) == RES_OK); interface_param = sdis_data_get(data); interface_param->temperature = 300; - interface_shader.convection_coef = null_interface_value; - interface_shader.temperature = interface_get_temperature; - interface_shader.emissivity = null_interface_value; - interface_shader.specular_fraction = null_interface_value; + interface_shader.front.temperature = interface_get_temperature; CHK(sdis_interface_create (dev, solid, fluid, &interface_shader, data, &T300) == RES_OK); CHK(sdis_data_ref_put(data) == RES_OK); @@ -221,10 +217,6 @@ main(int argc, char** argv) ALIGNOF(struct interf), NULL, &data) == RES_OK); interface_param = sdis_data_get(data); interface_param->temperature = 350; - interface_shader.convection_coef = null_interface_value; - interface_shader.temperature = interface_get_temperature; - interface_shader.emissivity = null_interface_value; - interface_shader.specular_fraction = null_interface_value; CHK(sdis_interface_create (dev, solid, fluid, &interface_shader, data, &T350) == RES_OK); CHK(sdis_data_ref_put(data) == RES_OK); diff --git a/src/test_sdis_solve_probe2_2d.c b/src/test_sdis_solve_probe2_2d.c @@ -194,21 +194,21 @@ main(int argc, char** argv) /* Create the fluid/solid interface with no limit conidition */ interface_shader.convection_coef = null_interface_value; - interface_shader.temperature = NULL; - interface_shader.emissivity = null_interface_value; - interface_shader.specular_fraction = null_interface_value; + interface_shader.front = SDIS_INTERFACE_SIDE_SHADER_NULL; + interface_shader.back = SDIS_INTERFACE_SIDE_SHADER_NULL; CHK(sdis_interface_create (dev, solid, fluid, &interface_shader, NULL, &Tnone) == RES_OK); + interface_shader.convection_coef = null_interface_value; + interface_shader.front = SDIS_INTERFACE_SIDE_SHADER_NULL; + interface_shader.back = SDIS_INTERFACE_SIDE_SHADER_NULL; + interface_shader.front.temperature = interface_get_temperature; + /* Create the fluid/solid interface with a fixed temperature of 300K */ CHK(sdis_data_create(dev, sizeof(struct interf), ALIGNOF(struct interf), NULL, &data) == RES_OK); interface_param = sdis_data_get(data); interface_param->temperature = 300; - interface_shader.convection_coef = null_interface_value; - interface_shader.temperature = interface_get_temperature; - interface_shader.emissivity = null_interface_value; - interface_shader.specular_fraction = null_interface_value; CHK(sdis_interface_create (dev, solid, fluid, &interface_shader, data, &T300) == RES_OK); CHK(sdis_data_ref_put(data) == RES_OK); @@ -218,10 +218,6 @@ main(int argc, char** argv) ALIGNOF(struct interf), NULL, &data) == RES_OK); interface_param = sdis_data_get(data); interface_param->temperature = 350; - interface_shader.convection_coef = null_interface_value; - interface_shader.temperature = interface_get_temperature; - interface_shader.emissivity = null_interface_value; - interface_shader.specular_fraction = null_interface_value; CHK(sdis_interface_create (dev, solid, fluid, &interface_shader, data, &T350) == RES_OK); CHK(sdis_data_ref_put(data) == RES_OK); diff --git a/src/test_sdis_solve_probe3.c b/src/test_sdis_solve_probe3.c @@ -224,9 +224,8 @@ main(int argc, char** argv) /* Create the fluid/solid interface with no limit conidition */ interface_shader.convection_coef = null_interface_value; - interface_shader.temperature = NULL; - interface_shader.emissivity = null_interface_value; - interface_shader.specular_fraction = null_interface_value; + interface_shader.front = SDIS_INTERFACE_SIDE_SHADER_NULL; + interface_shader.back = SDIS_INTERFACE_SIDE_SHADER_NULL; CHK(sdis_interface_create (dev, solid, fluid, &interface_shader, NULL, &Tnone) == RES_OK); @@ -235,10 +234,7 @@ main(int argc, char** argv) ALIGNOF(struct interf), NULL, &data) == RES_OK); interface_param = sdis_data_get(data); interface_param->temperature = 300; - interface_shader.convection_coef = null_interface_value; - interface_shader.temperature = interface_get_temperature; - interface_shader.emissivity = null_interface_value; - interface_shader.specular_fraction = null_interface_value; + interface_shader.front.temperature = interface_get_temperature; CHK(sdis_interface_create (dev, solid, fluid, &interface_shader, data, &T300) == RES_OK); CHK(sdis_data_ref_put(data) == RES_OK); @@ -248,19 +244,12 @@ main(int argc, char** argv) ALIGNOF(struct interf), NULL, &data) == RES_OK); interface_param = sdis_data_get(data); interface_param->temperature = 350; - interface_shader.convection_coef = null_interface_value; - interface_shader.temperature = interface_get_temperature; - interface_shader.emissivity = null_interface_value; - interface_shader.specular_fraction = null_interface_value; CHK(sdis_interface_create (dev, solid, fluid, &interface_shader, data, &T350) == RES_OK); CHK(sdis_data_ref_put(data) == RES_OK); /* Create the solid/solid interface */ - interface_shader.convection_coef = NULL; - interface_shader.temperature = NULL; - interface_shader.specular_fraction = NULL; - interface_shader.emissivity = NULL; + interface_shader = SDIS_INTERFACE_SHADER_NULL; CHK(sdis_interface_create (dev, solid, solid, &interface_shader, NULL, &solid_solid) == RES_OK); diff --git a/src/test_sdis_solve_probe3_2d.c b/src/test_sdis_solve_probe3_2d.c @@ -218,10 +218,7 @@ main(int argc, char** argv) CHK(sdis_solid_create(dev, &solid_shader, NULL, &solid) == RES_OK); /* Create the fluid/solid interface with no limit conidition */ - interface_shader.convection_coef = null_interface_value; - interface_shader.temperature = NULL; - interface_shader.emissivity = null_interface_value; - interface_shader.specular_fraction = null_interface_value; + interface_shader = SDIS_INTERFACE_SHADER_NULL; CHK(sdis_interface_create (dev, solid, fluid, &interface_shader, NULL, &Tnone) == RES_OK); @@ -231,9 +228,8 @@ main(int argc, char** argv) interface_param = sdis_data_get(data); interface_param->temperature = 300; interface_shader.convection_coef = null_interface_value; - interface_shader.temperature = interface_get_temperature; - interface_shader.emissivity = null_interface_value; - interface_shader.specular_fraction = null_interface_value; + interface_shader.front.temperature = interface_get_temperature; + interface_shader.back = SDIS_INTERFACE_SIDE_SHADER_NULL; CHK(sdis_interface_create (dev, solid, fluid, &interface_shader, data, &T300) == RES_OK); CHK(sdis_data_ref_put(data) == RES_OK); @@ -244,18 +240,14 @@ main(int argc, char** argv) interface_param = sdis_data_get(data); interface_param->temperature = 350; interface_shader.convection_coef = null_interface_value; - interface_shader.temperature = interface_get_temperature; - interface_shader.emissivity = null_interface_value; - interface_shader.specular_fraction = null_interface_value; + interface_shader.front.temperature = interface_get_temperature; + interface_shader.back = SDIS_INTERFACE_SIDE_SHADER_NULL; CHK(sdis_interface_create (dev, solid, fluid, &interface_shader, data, &T350) == RES_OK); CHK(sdis_data_ref_put(data) == RES_OK); /* Create the solid/solid interface */ - interface_shader.convection_coef = NULL; - interface_shader.temperature = NULL; - interface_shader.specular_fraction = NULL; - interface_shader.emissivity = NULL; + interface_shader = SDIS_INTERFACE_SHADER_NULL; CHK(sdis_interface_create (dev, solid, solid, &interface_shader, NULL, &solid_solid) == RES_OK); diff --git a/src/test_sdis_solve_probe_2d.c b/src/test_sdis_solve_probe_2d.c @@ -131,14 +131,6 @@ interface_get_convection_coef return 0.5; } -static double -interface_null_reflectivity - (const struct sdis_interface_fragment* frag, struct sdis_data* data) -{ - (void)frag, (void)data; - return 0; -} - /******************************************************************************* * Main test ******************************************************************************/ @@ -184,9 +176,8 @@ main(int argc, char** argv) /* Create the solid/fluid interface */ interface_shader.convection_coef = interface_get_convection_coef; - interface_shader.temperature = NULL; - interface_shader.emissivity = interface_null_reflectivity; - interface_shader.specular_fraction = interface_null_reflectivity; + interface_shader.front = SDIS_INTERFACE_SIDE_SHADER_NULL; + interface_shader.back = SDIS_INTERFACE_SIDE_SHADER_NULL; CHK(sdis_interface_create (dev, solid, fluid, &interface_shader, NULL, &interf) == RES_OK); diff --git a/src/test_sdis_solve_probe_boundary.c b/src/test_sdis_solve_probe_boundary.c @@ -204,22 +204,6 @@ interface_get_convection_coef 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 ******************************************************************************/ @@ -240,7 +224,7 @@ main(int argc, char** argv) 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_shader interf_shader = SDIS_INTERFACE_SHADER_NULL; struct sdis_interface* box_interfaces[12 /*#triangles*/]; struct sdis_interface* square_interfaces[4/*#segments*/]; struct interf* interf_props = NULL; @@ -270,10 +254,11 @@ main(int argc, char** argv) 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; + interf_shader.front.temperature = interface_get_temperature; + interf_shader.front.emissivity = NULL; + interf_shader.front.specular_fraction = NULL; + interf_shader.back = SDIS_INTERFACE_SIDE_SHADER_NULL; /* Create the adiabatic interface */ CHK(sdis_data_create(dev, sizeof(struct interf), 16, NULL, &data) == RES_OK); @@ -340,13 +325,15 @@ main(int argc, char** argv) 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); + #define F SDIS_FRONT + CHK(SOLVE(NULL, N, iprim, uv, INF, F, 1.0, 0, 0, &estimator) == RES_BAD_ARG); + CHK(SOLVE(box_scn, 0, iprim, uv, INF, F, 1.0, 0, 0, &estimator) == RES_BAD_ARG); + CHK(SOLVE(box_scn, N, 12, uv, INF, F, 1.0, 0, 0, &estimator) == RES_BAD_ARG); + CHK(SOLVE(box_scn, N, iprim, NULL, INF, F, 1.0, 0, 0, &estimator) == RES_BAD_ARG); + CHK(SOLVE(box_scn, N, iprim, uv, -1, F, 1.0, 0, 0, &estimator) == RES_BAD_ARG); + CHK(SOLVE(box_scn, N, iprim, uv, INF, -1, 1.0, 0, 0, &estimator) == RES_BAD_ARG); + CHK(SOLVE(box_scn, N, iprim, uv, INF, F, 1.0, 0, 0, NULL) == RES_BAD_ARG); + CHK(SOLVE(box_scn, N, iprim, uv, INF, F, 1.0, 0, 0, &estimator) == RES_OK); CHK(sdis_estimator_get_realisation_count(estimator, &nreals) == RES_OK); CHK(sdis_estimator_get_failure_count(estimator, &nfails) == RES_OK); @@ -355,10 +342,6 @@ main(int argc, char** argv) 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); @@ -370,7 +353,7 @@ main(int argc, char** argv) uv[0] = 0.5; iprim = 3; - CHK(SOLVE(square_scn, N, iprim, uv, INF, 1.0, 0, 0, &estimator) == RES_OK); + CHK(SOLVE(square_scn, N, iprim, uv, INF, F, 1.0, 0, 0, &estimator) == RES_OK); CHK(sdis_estimator_get_realisation_count(estimator, &nreals) == RES_OK); CHK(sdis_estimator_get_failure_count(estimator, &nfails) == RES_OK); CHK(nfails + nreals == N); diff --git a/src/test_sdis_utils.h b/src/test_sdis_utils.h @@ -114,11 +114,17 @@ static const struct sdis_fluid_shader DUMMY_FLUID_SHADER = { dummy_medium_getter }; + +#define DUMMY_INTERFACE_SIDE_SHADER__ { \ + dummy_interface_getter, \ + dummy_interface_getter, \ + dummy_interface_getter, \ + dummy_interface_getter \ +} static const struct sdis_interface_shader DUMMY_INTERFACE_SHADER = { dummy_interface_getter, - dummy_interface_getter, - dummy_interface_getter, - dummy_interface_getter + DUMMY_INTERFACE_SIDE_SHADER__, + DUMMY_INTERFACE_SIDE_SHADER__ }; /******************************************************************************* diff --git a/src/test_sdis_volumic_power.c b/src/test_sdis_volumic_power.c @@ -28,7 +28,7 @@ * * 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, + * with LAMBDA the conductivity of the solid and A the size of the cube/square, * i.e. 1. * * 3D 2D @@ -210,22 +210,6 @@ interface_get_convection_coef 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 ******************************************************************************/ @@ -245,7 +229,7 @@ main(int argc, char** argv) 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_shader interf_shader = SDIS_INTERFACE_SHADER_NULL; struct sdis_interface* box_interfaces[12 /*#triangles*/]; struct sdis_interface* square_interfaces[4/*#segments*/]; struct interf* interf_props = NULL; @@ -275,10 +259,8 @@ main(int argc, char** argv) 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; + interf_shader.front.temperature = interface_get_temperature; /* Create the adiabatic interface */ CHK(sdis_data_create(dev, sizeof(struct interf), 16, NULL, &data) == RES_OK); @@ -310,7 +292,7 @@ main(int argc, char** argv) /* Map the interfaces to their square segments */ square_interfaces[0] = interf_adiabatic; /* Bottom */ - square_interfaces[1] = interf_T0; /* Lef */ + square_interfaces[1] = interf_T0; /* Left */ square_interfaces[2] = interf_adiabatic; /* Top */ square_interfaces[3] = interf_T0; /* Right */