commit d851233f46080a6f79827ea5a173a131ce52b3b2
parent 3a660ba7a83168ebd18642fb645f55225f7438df
Author: Christophe Coustet <christophe.coustet@meso-star.com>
Date: Wed, 10 Mar 2021 16:35:50 +0100
Merge branch 'feature_thermal_contact_resistance' into develop
Diffstat:
9 files changed, 1213 insertions(+), 14 deletions(-)
diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt
@@ -178,6 +178,8 @@ if(NOT NO_TEST)
new_test(test_sdis_compute_power)
new_test(test_sdis_conducto_radiative)
new_test(test_sdis_conducto_radiative_2d)
+ new_test(test_sdis_contact_resistance)
+ new_test(test_sdis_contact_resistance_2)
new_test(test_sdis_convection)
new_test(test_sdis_convection_non_uniform)
new_test(test_sdis_data)
diff --git a/src/sdis.h b/src/sdis.h
@@ -221,11 +221,16 @@ struct sdis_interface_shader {
* [0 convection_coef_upper_bound] */
double convection_coef_upper_bound;
+ /* May be NULL for solid/fluid or if the thermal contact resistance is 0 onto
+ * the whole interface. */
+ sdis_interface_getter_T thermal_contact_resistance; /* In K.m^2.W^-1 */
+
struct sdis_interface_side_shader front;
struct sdis_interface_side_shader back;
};
#define SDIS_INTERFACE_SHADER_NULL__ \
- {NULL, 0, SDIS_INTERFACE_SIDE_SHADER_NULL__, SDIS_INTERFACE_SIDE_SHADER_NULL__}
+ {NULL, 0, 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__;
diff --git a/src/sdis_heat_path_boundary_Xd.h b/src/sdis_heat_path_boundary_Xd.h
@@ -478,6 +478,7 @@ XD(solid_solid_boundary_path)
double tmp;
double r;
double power;
+ double tcr;
float dir0[DIM], dir1[DIM], dir2[DIM], dir3[DIM];
float dir_front[DIM], dir_back[DIM];
float* dir;
@@ -501,6 +502,9 @@ XD(solid_solid_boundary_path)
ASSERT(solid_front->type == SDIS_SOLID);
ASSERT(solid_back->type == SDIS_SOLID);
+ /* Retrieve the thermal contact resistance */
+ tcr = interface_get_thermal_contact_resistance(interf, frag);
+
/* Fetch the properties of the media */
lambda_front = solid_get_thermal_conductivity(solid_front, &rwalk->vtx);
lambda_back = solid_get_thermal_conductivity(solid_back, &rwalk->vtx);
@@ -560,19 +564,44 @@ XD(solid_solid_boundary_path)
goto error;
}
- /* Define the reinjection side. Note that the proba should be :
- * Lf/Df' / (Lf/Df' + Lb/Db')
- *
- * with L<f|b> the lambda of the <front|back> side and D<f|b>' the adjusted
- * delta of the <front|back> side, i.e. :
- * D<f|b>' = reinject_dst_<front|back> / sqrt(DIM)
- *
- * Anyway, one can avoid to compute the adjusted delta by directly using the
- * adjusted reinjection distance since the resulting proba is strictly the
- * same; sqrt(DIM) can be simplified. */
r = ssp_rng_canonical(rng);
- proba = (lambda_front/reinject_dst_front)
- / (lambda_front/reinject_dst_front + lambda_back/reinject_dst_back);
+ if(tcr == 0) { /* No thermal contact resistance */
+ /* Define the reinjection side. Note that the proba should be : Lf/Df' /
+ * (Lf/Df' + Lb/Db')
+ *
+ * with L<f|b> the lambda of the <front|back> side and D<f|b>' the adjusted
+ * delta of the <front|back> side, i.e. : D<f|b>' =
+ * reinject_dst_<front|back> / sqrt(DIM)
+ *
+ * Anyway, one can avoid to compute the adjusted delta by directly using the
+ * adjusted reinjection distance since the resulting proba is strictly the
+ * same; sqrt(DIM) can be simplified. */
+ proba = (lambda_front/reinject_dst_front)
+ / (lambda_front/reinject_dst_front + lambda_back/reinject_dst_back);
+ } else {
+ const double df = reinject_dst_front/sqrt(DIM);
+ const double db = reinject_dst_back/sqrt(DIM);
+ const double tmp_front = lambda_front/df;
+ const double tmp_back = lambda_back/db;
+ const double tmp_r = tcr*tmp_front*tmp_back;
+ switch(rwalk->hit_side) {
+ case SDIS_BACK:
+ /* When coming from the BACK side, the probability to be reinjected on
+ * the FRONT side depends on the thermal contact resistance: it decrases
+ * when the TCR increases (and tends to 0 when TCR -> +\infty) */
+ proba = (tmp_front) / (tmp_front + tmp_back + tmp_r);
+ break;
+ case SDIS_FRONT:
+ /* Same thing when coming from the FRONT side: the probability of
+ * reinjection on the FRONT side depends on the thermal contact
+ * resistance: it increases when the TCR increases (and tends to 1 when
+ * the TCR -> +\infty) */
+ proba = (tmp_front + tmp_r) / (tmp_front + tmp_back + tmp_r);
+ break;
+ default: FATAL("Unreachable code.\n"); break;
+ }
+ }
+
if(r < proba) { /* Reinject in front */
dir = dir_front;
hit = &hit0;
diff --git a/src/sdis_interface.c b/src/sdis_interface.c
@@ -56,7 +56,6 @@ check_interface_shader
" shader's pointer function for this attribute should be NULL.\n",
caller_name);
}
-
if(shader->convection_coef_upper_bound < 0) {
log_warn(dev,
"%s: Invalid upper bound for convection coefficient (%g).\n",
@@ -64,6 +63,14 @@ check_interface_shader
if(type[0] == SDIS_FLUID || type[1] == SDIS_FLUID) return 0;
}
+ if((type[0] != SDIS_SOLID || type[1] != SDIS_SOLID)
+ && shader->thermal_contact_resistance) {
+ log_warn(dev,
+ "%s: only solid/solid interface can have a thermal contact resistance. The "
+ " shader's pointer function for this attribute should be NULL.\n",
+ caller_name);
+ }
+
FOR_EACH(i, 0, 2) {
switch(type[i]) {
case SDIS_SOLID:
diff --git a/src/sdis_interface_c.h b/src/sdis_interface_c.h
@@ -89,6 +89,16 @@ interface_get_convection_coef
}
static INLINE double
+interface_get_thermal_contact_resistance
+ (const struct sdis_interface* interf,
+ const struct sdis_interface_fragment* frag)
+{
+ ASSERT(interf && frag);
+ return interf->shader.thermal_contact_resistance
+ ? interf->shader.thermal_contact_resistance(frag, interf->data) : 0;
+}
+
+static INLINE double
interface_get_convection_coef_upper_bound
(const struct sdis_interface* interf)
{
diff --git a/src/test_sdis_contact_resistance.c b/src/test_sdis_contact_resistance.c
@@ -0,0 +1,444 @@
+/* Copyright (C) 2016-2020 |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 "test_sdis_contact_resistance.h"
+
+#include <rsys/clock_time.h>
+#include <rsys/mem_allocator.h>
+#include <rsys/double3.h>
+#include <rsys/math.h>
+#include <star/ssp.h>
+
+/*
+ * The scene is composed of a solid cube/square of size L whose temperature is
+ * unknown. The cube is made of 2 solids that meet at x=e in ]0 L[ and the
+ * thermal contact resistance is R, thus T(X0-) differs from T(X0+).
+ * The faces are adiabatic exept at x=0 where T(0)=T0 and at x=L where T(L)=TL.
+ * At steady state:
+ *
+ * T(X0-) = (T0 * (R * LAMBDA1 / X0) * (1 + R * LAMBDA2 / (L - X0))
+ * + TL * (R * LAMBDA2 / (L - X0)))
+ * / ((1 + R * LAMBDA1 / X0) * (1 + R * LAMBDA2 / (L - X0)) - 1)
+ *
+ * T(X0+) = T(X0-) * (1 + r * LAMBDA1 / X0) - T0 * r * LAMBDA1 / X0
+ *
+ * T(x) is linear between T(0) and T(X0-) if x in [0 X0[
+ * T(x) is linear between T(X0+) and T(L) if x in ]X0 L]
+ *
+ * 3D 2D
+ *
+ * /////////(L,L,L) /////////(L,L)
+ * +-------+ +-------+
+ * /' / /| | ! |
+ * +-------+ TL T0 r TL
+ * | | ! | | | ! |
+ * T0 +.r...|.+ +-------+
+ * |, ! |/ (0,0)///x=X0///
+ * +-------+
+ * (0,0,0)///x=X0///
+ */
+
+#define UNKNOWN_TEMPERATURE -1
+#define N 10000 /* #realisations */
+
+#define T0 0.0
+#define LAMBDA1 0.1
+
+#define TL 100.0
+#define LAMBDA2 0.2
+
+#define DELTA1 X0/25.0
+#define DELTA2 (L-X0)/25.0
+
+/*******************************************************************************
+ * Media
+ ******************************************************************************/
+struct solid {
+ double lambda;
+ double rho;
+ double cp;
+ double delta;
+};
+
+static double
+fluid_get_temperature
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ (void)data;
+ CHK(vtx);
+ return UNKNOWN_TEMPERATURE;
+}
+
+static double
+solid_get_calorific_capacity
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ CHK(vtx && data);
+ return ((struct solid*)sdis_data_cget(data))->cp;
+}
+
+static double
+solid_get_thermal_conductivity
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ CHK(vtx && data);
+ return ((struct solid*)sdis_data_cget(data))->lambda;
+}
+
+static double
+solid_get_volumic_mass
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ CHK(vtx && data);
+ return ((struct solid*)sdis_data_cget(data))->rho;
+}
+
+static double
+solid_get_delta
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ CHK(vtx && data);
+ return ((struct solid*)sdis_data_cget(data))->delta;
+}
+
+static double
+solid_get_temperature
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ CHK(vtx && data);
+ return UNKNOWN_TEMPERATURE;
+}
+
+/*******************************************************************************
+ * Interfaces
+ ******************************************************************************/
+struct interf {
+ double temperature;
+ double resistance;
+};
+
+static double
+interface_get_temperature
+ (const struct sdis_interface_fragment* frag, struct sdis_data* data)
+{
+ const struct interf* interf = sdis_data_cget(data);
+ CHK(frag && data);
+ return interf->temperature;
+}
+
+static double
+interface_get_convection_coef
+ (const struct sdis_interface_fragment* frag, struct sdis_data* data)
+{
+ CHK(frag && data);
+ return 0;
+}
+
+static double
+interface_get_contact_resistance
+ (const struct sdis_interface_fragment* frag, struct sdis_data* data)
+{
+ const struct interf* interf = sdis_data_cget(data);
+ CHK(frag && data);
+ return interf->resistance;
+}
+
+/*******************************************************************************
+ * Helper functions
+ ******************************************************************************/
+static void
+solve
+ (struct sdis_scene* scn,
+ struct interf* interf_props,
+ struct ssp_rng* rng)
+{
+ char dump[128];
+ struct time t0, t1;
+ struct sdis_estimator* estimator;
+ struct sdis_solve_probe_args solve_args = SDIS_SOLVE_PROBE_ARGS_DEFAULT;
+ struct sdis_mc T = SDIS_MC_NULL;
+ struct sdis_mc time = SDIS_MC_NULL;
+ size_t nreals;
+ size_t nfails;
+ double ref_L, ref_R;
+ enum sdis_scene_dimension dim;
+ const int nsimuls = 8;
+ int isimul;
+ ASSERT(scn && interf_props && rng);
+
+ OK(sdis_scene_get_dimension(scn, &dim));
+
+ FOR_EACH(isimul, 0, nsimuls) {
+ double x, ref;
+ double r = pow(10, ssp_rng_uniform_double(rng, -2, 2));
+
+ interf_props->resistance = r;
+
+ ref_L = (
+ T0 * (r * LAMBDA1 / X0) * (1 + r * LAMBDA2 / (L - X0))
+ + TL * (r * LAMBDA2 / (L - X0))
+ )
+ / ((1 + r * LAMBDA1 / X0) * (1 + r * LAMBDA2 / (L - X0)) - 1);
+
+ ref_R = ref_L * (1 + r * LAMBDA1 / X0) - T0 * r * LAMBDA1 / X0;
+
+ if(isimul % 2) { /* In solid 1 */
+ x = ssp_rng_uniform_double(rng, 0.05 * X0, 0.95 * X0);
+ ref = T0 * (1 - x / X0) + ref_L * x / X0;
+ } else { /* In solid 2 */
+ x = X0 + ssp_rng_uniform_double(rng, 0.05 * (L - X0), 0.95 * (L - X0));
+ ref = ref_R * (1 - (x - X0) / (L - X0)) + TL * (x - X0) / (L - X0);
+ }
+
+ solve_args.position[0] = x;
+ solve_args.position[1] = ssp_rng_uniform_double(rng, 0.05 * L, 0.95 * L);
+ solve_args.position[2] = (dim == SDIS_SCENE_2D)
+ ? 0 : ssp_rng_uniform_double(rng, 0.05 * L, 0.95 * L);
+
+ solve_args.nrealisations = N;
+ solve_args.time_range[0] = solve_args.time_range[1] = INF;
+
+ time_current(&t0);
+ OK(sdis_solve_probe(scn, &solve_args, &estimator));
+ time_sub(&t0, time_current(&t1), &t0);
+ time_dump(&t0, TIME_ALL, NULL, dump, sizeof(dump));
+
+ OK(sdis_estimator_get_realisation_count(estimator, &nreals));
+ OK(sdis_estimator_get_failure_count(estimator, &nfails));
+ OK(sdis_estimator_get_temperature(estimator, &T));
+ OK(sdis_estimator_get_realisation_time(estimator, &time));
+
+ switch(dim) {
+ case SDIS_SCENE_2D:
+ printf("Steady temperature at (%g, %g) with R=%g = %g ~ %g +/- %g\n",
+ SPLIT2(solve_args.position), r, ref, T.E, T.SE);
+ break;
+ case SDIS_SCENE_3D:
+ printf("Steady temperature at (%g, %g, %g) with R=%g = %g ~ %g +/- %g\n",
+ SPLIT3(solve_args.position), r, ref, T.E, T.SE);
+ break;
+ default: FATAL("Unreachable code.\n"); break;
+ }
+ printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N);
+ printf("Elapsed time = %s\n", dump);
+ printf("Time per realisation (in usec) = %g +/- %g\n\n", time.E, time.SE);
+
+ CHK(nfails + nreals == N);
+ CHK(nfails <= N/1000);
+ CHK(eq_eps(T.E, ref, T.SE * 3));
+
+ OK(sdis_estimator_ref_put(estimator));
+ }
+}
+
+/*******************************************************************************
+ * Test
+ ******************************************************************************/
+int
+main(int argc, char** argv)
+{
+ struct mem_allocator allocator;
+ struct sdis_data* data = NULL;
+ struct sdis_device* dev = NULL;
+ struct sdis_medium* fluid = NULL;
+ struct sdis_medium* solid1 = NULL;
+ struct sdis_medium* solid2 = NULL;
+ struct sdis_interface* interf_adiabatic1 = NULL;
+ struct sdis_interface* interf_adiabatic2 = NULL;
+ struct sdis_interface* interf_T0 = NULL;
+ struct sdis_interface* interf_TL = NULL;
+ struct sdis_interface* interf_R = NULL;
+ struct sdis_scene* box_scn = NULL;
+ struct sdis_scene* square_scn = NULL;
+ struct sdis_scene_create_args scn_args = SDIS_SCENE_CREATE_ARGS_DEFAULT;
+ 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* model3d_interfaces[22 /*#triangles*/];
+ struct sdis_interface* model2d_interfaces[7/*#segments*/];
+ struct interf* interf_props = NULL;
+ struct solid* solid_props = NULL;
+ struct ssp_rng* rng = NULL;
+ (void)argc, (void)argv;
+
+ OK(mem_init_proxy_allocator(&allocator, &mem_default_allocator));
+ OK(sdis_device_create(NULL, &allocator, SDIS_NTHREADS_DEFAULT, 1, &dev));
+
+ fluid_shader.temperature = fluid_get_temperature;
+ OK(sdis_fluid_create(dev, &fluid_shader, NULL, &fluid));
+
+ /* Setup the solid shader */
+ solid_shader.calorific_capacity = solid_get_calorific_capacity;
+ solid_shader.thermal_conductivity = solid_get_thermal_conductivity;
+ solid_shader.volumic_mass = solid_get_volumic_mass;
+ solid_shader.delta_solid = solid_get_delta;
+ solid_shader.temperature = solid_get_temperature;
+
+ /* Create the solid medium #1 */
+ OK(sdis_data_create(dev, sizeof(struct solid), 16, NULL, &data));
+ solid_props = sdis_data_get(data);
+ solid_props->lambda = LAMBDA1;
+ solid_props->cp = 2;
+ solid_props->rho = 25;
+ solid_props->delta = DELTA1;
+ OK(sdis_solid_create(dev, &solid_shader, data, &solid1));
+ OK(sdis_data_ref_put(data));
+
+ /* Create the solid medium #2 */
+ OK(sdis_data_create(dev, sizeof(struct solid), 16, NULL, &data));
+ solid_props = sdis_data_get(data);
+ solid_props->lambda = LAMBDA2;
+ solid_props->cp = 2;
+ solid_props->rho = 25;
+ solid_props->delta = DELTA2;
+ OK(sdis_solid_create(dev, &solid_shader, data, &solid2));
+ OK(sdis_data_ref_put(data));
+
+ /* Setup the interface shader */
+ interf_shader.front.temperature = interface_get_temperature;
+ interf_shader.back.temperature = interface_get_temperature;
+ interf_shader.convection_coef = interface_get_convection_coef;
+
+ /* Create the adiabatic interfaces */
+ OK(sdis_data_create(dev, sizeof(struct interf), 16, NULL, &data));
+ interf_props = sdis_data_get(data);
+ interf_props->temperature = UNKNOWN_TEMPERATURE;
+ OK(sdis_interface_create
+ (dev, solid1, fluid, &interf_shader, data, &interf_adiabatic1));
+ OK(sdis_interface_create
+ (dev, solid2, fluid, &interf_shader, data, &interf_adiabatic2));
+ OK(sdis_data_ref_put(data));
+
+ /* Create the T0 interface */
+ OK(sdis_data_create(dev, sizeof(struct interf), 16, NULL, &data));
+ interf_props = sdis_data_get(data);
+ interf_props->temperature = T0;
+ OK(sdis_interface_create
+ (dev, solid1, fluid, &interf_shader, data, &interf_T0));
+ OK(sdis_data_ref_put(data));
+
+ /* Create the TL interface */
+ OK(sdis_data_create(dev, sizeof(struct interf), 16, NULL, &data));
+ interf_props = sdis_data_get(data);
+ interf_props->temperature = TL;
+ OK(sdis_interface_create
+ (dev, solid2, fluid, &interf_shader, data, &interf_TL));
+ OK(sdis_data_ref_put(data));
+
+ /* Create the solid1-solid2 interface */
+ interf_shader.convection_coef = NULL;
+ interf_shader.thermal_contact_resistance = interface_get_contact_resistance;
+ OK(sdis_data_create(dev, sizeof(struct interf), 16, NULL, &data));
+ interf_props = sdis_data_get(data);
+ interf_props->temperature = UNKNOWN_TEMPERATURE;
+ OK(sdis_interface_create
+ (dev, solid1, solid2, &interf_shader, data, &interf_R));
+ OK(sdis_data_ref_put(data));
+
+ /* Release the media */
+ OK(sdis_medium_ref_put(solid1));
+ OK(sdis_medium_ref_put(solid2));
+ OK(sdis_medium_ref_put(fluid));
+
+ /* Map the interfaces to their box triangles */
+ /* Front */
+ model3d_interfaces[0] = interf_adiabatic1;
+ model3d_interfaces[1] = interf_adiabatic1;
+ model3d_interfaces[2] = interf_adiabatic2;
+ model3d_interfaces[3] = interf_adiabatic2;
+ /* Left */
+ model3d_interfaces[4] = interf_T0;
+ model3d_interfaces[5] = interf_T0;
+ /* Back */
+ model3d_interfaces[6] = interf_adiabatic1;
+ model3d_interfaces[7] = interf_adiabatic1;
+ model3d_interfaces[8] = interf_adiabatic2;
+ model3d_interfaces[9] = interf_adiabatic2;
+ /* Right */
+ model3d_interfaces[10] = interf_TL;
+ model3d_interfaces[11] = interf_TL;
+ /* Top */
+ model3d_interfaces[12] = interf_adiabatic1;
+ model3d_interfaces[13] = interf_adiabatic1;
+ model3d_interfaces[14] = interf_adiabatic2;
+ model3d_interfaces[15] = interf_adiabatic2;
+ /* Bottom */
+ model3d_interfaces[16] = interf_adiabatic1;
+ model3d_interfaces[17] = interf_adiabatic1;
+ model3d_interfaces[18] = interf_adiabatic2;
+ model3d_interfaces[19] = interf_adiabatic2;
+ /* Inside */
+ model3d_interfaces[20] = interf_R;
+ model3d_interfaces[21] = interf_R;
+
+ /* Map the interfaces to their square segments */
+ /* Bottom */
+ model2d_interfaces[0] = interf_adiabatic2;
+ model2d_interfaces[1] = interf_adiabatic1;
+ /* Left */
+ model2d_interfaces[2] = interf_T0;
+ /* Top */
+ model2d_interfaces[3] = interf_adiabatic1;
+ model2d_interfaces[4] = interf_adiabatic2;
+ /* Right */
+ model2d_interfaces[5] = interf_TL;
+ /* Contact */
+ model2d_interfaces[6] = interf_R;
+
+ /* Create the box scene */
+ scn_args.get_indices = model3d_get_indices;
+ scn_args.get_interface = model3d_get_interface;
+ scn_args.get_position = model3d_get_position;
+ scn_args.nprimitives = model3d_ntriangles;
+ scn_args.nvertices = model3d_nvertices;
+ scn_args.context = model3d_interfaces;
+ OK(sdis_scene_create(dev, &scn_args, &box_scn));
+
+ /* Create the square scene */
+ scn_args.get_indices = model2d_get_indices;
+ scn_args.get_interface = model2d_get_interface;
+ scn_args.get_position = model2d_get_position;
+ scn_args.nprimitives = model2d_nsegments;
+ scn_args.nvertices = model2d_nvertices;
+ scn_args.context = model2d_interfaces;
+ OK(sdis_scene_2d_create(dev, &scn_args, &square_scn));
+
+ /* Release the interfaces */
+ OK(sdis_interface_ref_put(interf_adiabatic1));
+ OK(sdis_interface_ref_put(interf_adiabatic2));
+ OK(sdis_interface_ref_put(interf_T0));
+ OK(sdis_interface_ref_put(interf_TL));
+ OK(sdis_interface_ref_put(interf_R));
+
+ /* Solve */
+ OK(ssp_rng_create(&allocator, &ssp_rng_kiss, &rng));
+ printf(">> Box scene\n");
+ solve(box_scn, interf_props, rng);
+ printf("\n>> Square scene\n");
+ solve(square_scn, interf_props, rng);
+
+ OK(sdis_scene_ref_put(box_scn));
+ OK(sdis_scene_ref_put(square_scn));
+ OK(sdis_device_ref_put(dev));
+ OK(ssp_rng_ref_put(rng));
+
+ check_memory_allocator(&allocator);
+ mem_shutdown_proxy_allocator(&allocator);
+ CHK(mem_allocated_size() == 0);
+ return 0;
+}
+
diff --git a/src/test_sdis_contact_resistance.h b/src/test_sdis_contact_resistance.h
@@ -0,0 +1,161 @@
+/* Copyright (C) 2016-2020 |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/>. */
+
+
+/*
+ * The scene is composed of a solid cube/square of size L. The cube/square
+ * is made of 2 solids that meet at x=e in ]0 L[.
+ *
+ * 3D 2D
+ *
+ * /////////(L,L,L) /////////(L,L)
+ * +-------+ +-------+
+ * /' / /| | ! |
+ * +-------+ TL T0 r TL
+ * | | ! | | | ! |
+ * T0 +.r...|.+ +-------+
+ * |, ! |/ (0,0)///x=X0///
+ * +-------+
+ * (0,0,0)///x=X0///
+ */
+
+#define L 4.0
+#define X0 3.0
+
+ /*******************************************************************************
+ * Box geometry
+ ******************************************************************************/
+static const double model3d_vertices[12/*#vertices*/ * 3/*#coords per vertex*/]
+= {
+ 0, 0, 0,
+ X0, 0, 0,
+ L, 0, 0,
+ 0, L, 0,
+ X0, L, 0,
+ L, L, 0,
+ 0, 0, L,
+ X0, 0, L,
+ L, 0, L,
+ 0, L, L,
+ X0, L, L,
+ L, L, L
+};
+static const size_t model3d_nvertices = sizeof(model3d_vertices) / sizeof(double[3]);
+
+/* The following array lists the indices toward the 3D vertices of each
+ * triangle.
+ * ,3---,4---,5 ,3----4----5 ,4
+ * ,' | ,' | ,'/| ,'/| \ | \ | ,'/|
+ * 9----10---11 / | 9' / | \ | \ | 10 / | Y
+ * |', |', | / ,2 | / ,0---,1---,2 | / ,1 |
+ * | ',| ',|/,' |/,' | ,' | ,' |/,' o--X
+ * 6----7----8' 6----7'---8' 7 /
+ * Front, right Back, left and Internal Z
+ * and Top faces bottom faces face */
+static const size_t model3d_indices[22/*#triangles*/ * 3/*#indices per triangle*/]
+= {
+ 0, 3, 1, 1, 3, 4, 1, 4, 2, 2, 4, 5, /* -Z */
+ 0, 6, 3, 3, 6, 9, /* -X */
+ 6, 7, 9, 9, 7, 10, 7, 8, 10, 10, 8, 11, /* +Z */
+ 5, 11, 8, 8, 2, 5, /* +X */
+ 3, 9, 10, 10, 4, 3, 4, 10, 11, 11, 5, 4, /* +Y */
+ 0, 1, 7, 7, 6, 0, 1, 2, 8, 8, 7, 1, /* -Y */
+ 4, 10, 7, 7, 1, 4 /* Inside */
+};
+static const size_t model3d_ntriangles = sizeof(model3d_indices) / sizeof(size_t[3]);
+
+static INLINE void
+model3d_get_indices(const size_t itri, size_t ids[3], void* context)
+{
+ (void)context;
+ CHK(ids);
+ CHK(itri < model3d_ntriangles);
+ ids[0] = model3d_indices[itri * 3 + 0];
+ ids[1] = model3d_indices[itri * 3 + 1];
+ ids[2] = model3d_indices[itri * 3 + 2];
+}
+
+static INLINE void
+model3d_get_position(const size_t ivert, double pos[3], void* context)
+{
+ (void)context;
+ CHK(pos);
+ CHK(ivert < model3d_nvertices);
+ pos[0] = model3d_vertices[ivert * 3 + 0];
+ pos[1] = model3d_vertices[ivert * 3 + 1];
+ pos[2] = model3d_vertices[ivert * 3 + 2];
+}
+
+static INLINE void
+model3d_get_interface(const size_t itri, struct sdis_interface** bound, void* context)
+{
+ struct sdis_interface** interfaces = context;
+ CHK(context && bound);
+ CHK(itri < model3d_ntriangles);
+ *bound = interfaces[itri];
+}
+
+/*******************************************************************************
+ * Square geometry
+ ******************************************************************************/
+static const double model2d_vertices[6/*#vertices*/ * 2/*#coords per vertex*/] = {
+ L, 0,
+ X0, 0,
+ 0, 0,
+ 0, L,
+ X0, L,
+ L, L
+};
+static const size_t model2d_nvertices = sizeof(model2d_vertices) / sizeof(double[2]);
+
+static const size_t model2d_indices[7/*#segments*/ * 2/*#indices per segment*/] = {
+ 0, 1, 1, 2, /* Bottom */
+ 2, 3, /* Left */
+ 3, 4, 4, 5, /* Top */
+ 5, 0, /* Right */
+ 4, 1 /* Inside */
+};
+static const size_t model2d_nsegments = sizeof(model2d_indices) / sizeof(size_t[2]);
+
+
+static INLINE void
+model2d_get_indices(const size_t iseg, size_t ids[2], void* context)
+{
+ (void)context;
+ CHK(ids);
+ CHK(iseg < model2d_nsegments);
+ ids[0] = model2d_indices[iseg * 2 + 0];
+ ids[1] = model2d_indices[iseg * 2 + 1];
+}
+
+static INLINE void
+model2d_get_position(const size_t ivert, double pos[2], void* context)
+{
+ (void)context;
+ CHK(pos);
+ CHK(ivert < model2d_nvertices);
+ pos[0] = model2d_vertices[ivert * 2 + 0];
+ pos[1] = model2d_vertices[ivert * 2 + 1];
+}
+
+static INLINE void
+model2d_get_interface
+(const size_t iseg, struct sdis_interface** bound, void* context)
+{
+ struct sdis_interface** interfaces = context;
+ CHK(context && bound);
+ CHK(iseg < model2d_nsegments);
+ *bound = interfaces[iseg];
+}
diff --git a/src/test_sdis_contact_resistance_2.c b/src/test_sdis_contact_resistance_2.c
@@ -0,0 +1,540 @@
+/* Copyright (C) 2016-2020 |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 "test_sdis_contact_resistance.h"
+
+#include <rsys/clock_time.h>
+#include <rsys/mem_allocator.h>
+#include <rsys/double3.h>
+#include <rsys/math.h>
+#include <star/ssp.h>
+
+/*
+ * The scene is composed of a solid cube/square of size L whose temperature is
+ * unknown. The cube is made of 2 solids that meet at x=e in ]0 L[ and the
+ * thermal contact resistance is R, thus T(X0-) differs from T(X0+).
+ * The faces are adiabatic exept at x=0 where T(0)=T0 and at x=L where T(L)=TL.
+ * At steady state:
+ *
+ * Flux(x0) = (T(x0+) - T(x0-)) / R
+ *
+ * 3D 2D
+ *
+ * /////////(L,L,L) /////////(L,L)
+ * +-------+ +-------+
+ * /' / /| | ! |
+ * +-------+ TL T0 r TL
+ * | | ! | | | ! |
+ * T0 +.r...|.+ +-------+
+ * |, ! |/ (0,0)///x=X0///
+ * +-------+
+ * (0,0,0)///x=X0///
+ */
+
+#define UNKNOWN_TEMPERATURE -1
+#define N 10000 /* #realisations */
+
+#define T0 0.0
+#define LAMBDA1 0.1
+
+#define TL 100.0
+#define LAMBDA2 0.2
+
+#define DELTA1 X0/25.0
+#define DELTA2 (L-X0)/25.0
+
+/*******************************************************************************
+ * Media
+ ******************************************************************************/
+struct solid {
+ double lambda;
+ double rho;
+ double cp;
+ double delta;
+};
+
+static double
+fluid_get_temperature
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ (void)data;
+ CHK(vtx);
+ return UNKNOWN_TEMPERATURE;
+}
+
+static double
+solid_get_calorific_capacity
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ CHK(vtx && data);
+ return ((struct solid*)sdis_data_cget(data))->cp;
+}
+
+static double
+solid_get_thermal_conductivity
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ CHK(vtx && data);
+ return ((struct solid*)sdis_data_cget(data))->lambda;
+}
+
+static double
+solid_get_volumic_mass
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ CHK(vtx && data);
+ return ((struct solid*)sdis_data_cget(data))->rho;
+}
+
+static double
+solid_get_delta
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ CHK(vtx && data);
+ return ((struct solid*)sdis_data_cget(data))->delta;
+}
+
+static double
+solid_get_temperature
+ (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data)
+{
+ CHK(vtx && data);
+ return UNKNOWN_TEMPERATURE;
+}
+
+/*******************************************************************************
+ * Interfaces
+ ******************************************************************************/
+struct interf {
+ double temperature;
+ double resistance;
+};
+
+static double
+interface_get_temperature
+ (const struct sdis_interface_fragment* frag, struct sdis_data* data)
+{
+ const struct interf* interf = sdis_data_cget(data);
+ CHK(frag && data);
+ return interf->temperature;
+}
+
+static double
+interface_get_convection_coef
+ (const struct sdis_interface_fragment* frag, struct sdis_data* data)
+{
+ CHK(frag && data);
+ return 0;
+}
+
+static double
+interface_get_contact_resistance
+ (const struct sdis_interface_fragment* frag, struct sdis_data* data)
+{
+ const struct interf* interf = sdis_data_cget(data);
+ CHK(frag && data);
+ return interf->resistance;
+}
+
+/*******************************************************************************
+ * Helper functions
+ ******************************************************************************/
+static void
+solve_probe
+ (struct sdis_scene* scn,
+ struct interf* interf_props,
+ struct ssp_rng* rng)
+{
+ char dump[128];
+ struct time t0, t1;
+ struct sdis_estimator* estimator;
+ struct sdis_solve_probe_boundary_args solve_args
+ = SDIS_SOLVE_PROBE_BOUNDARY_ARGS_DEFAULT;
+ struct sdis_mc T = SDIS_MC_NULL;
+ struct sdis_mc time = SDIS_MC_NULL;
+ size_t nreals;
+ size_t nfails;
+ double ref_L, ref_R;
+ enum sdis_scene_dimension dim;
+ int nsimuls;
+ int isimul;
+ ASSERT(scn && interf_props && rng);
+
+ OK(sdis_scene_get_dimension(scn, &dim));
+
+ nsimuls = (dim == SDIS_SCENE_2D) ? 2 : 4;
+ FOR_EACH(isimul, 0, nsimuls) {
+ double ref;
+ double r = pow(10, ssp_rng_uniform_double(rng, -2, 2));
+
+ interf_props->resistance = r;
+
+ ref_L = (
+ T0 * (r * LAMBDA1 / X0) * (1 + r * LAMBDA2 / (L - X0))
+ + TL * (r * LAMBDA2 / (L - X0))
+ )
+ / ((1 + r * LAMBDA1 / X0) * (1 + r * LAMBDA2 / (L - X0)) - 1);
+
+ ref_R = ref_L * (1 + r * LAMBDA1 / X0) - T0 * r * LAMBDA1 / X0;
+
+ if(dim == SDIS_SCENE_2D)
+ /* last segment */
+ solve_args.iprim = model2d_nsegments - 1;
+ else
+ /* last 2 triangles */
+ solve_args.iprim = model3d_ntriangles - ((isimul % 2) ? 2 : 1);
+
+ solve_args.uv[0] = ssp_rng_canonical(rng);
+ solve_args.uv[1] = (dim == SDIS_SCENE_2D)
+ ? 0 : ssp_rng_uniform_double(rng, 0, 1 - solve_args.uv[0]);
+
+ if(isimul < nsimuls / 2) {
+ solve_args.side = SDIS_FRONT;
+ ref = ref_L;
+ } else {
+ solve_args.side = SDIS_BACK;
+ ref = ref_R;
+ }
+
+ solve_args.nrealisations = N;
+ solve_args.time_range[0] = solve_args.time_range[1] = INF;
+
+ time_current(&t0);
+ OK(sdis_solve_probe_boundary(scn, &solve_args, &estimator));
+ time_sub(&t0, time_current(&t1), &t0);
+ time_dump(&t0, TIME_ALL, NULL, dump, sizeof(dump));
+
+ OK(sdis_estimator_get_realisation_count(estimator, &nreals));
+ OK(sdis_estimator_get_failure_count(estimator, &nfails));
+ OK(sdis_estimator_get_temperature(estimator, &T));
+ OK(sdis_estimator_get_realisation_time(estimator, &time));
+
+ switch(dim) {
+ case SDIS_SCENE_2D:
+ printf("Steady temperature at (%lu/%s/%g) with R=%g = %g ~ %g +/- %g\n",
+ (unsigned long)solve_args.iprim,
+ (solve_args.side == SDIS_FRONT ? "front" : "back"),
+ solve_args.uv[0],
+ r, ref, T.E, T.SE);
+ break;
+ case SDIS_SCENE_3D:
+ printf("Steady temperature at (%lu/%s/%g,%g) with R=%g = %g ~ %g +/- %g\n",
+ (unsigned long)solve_args.iprim,
+ (solve_args.side == SDIS_FRONT ? "front" : "back"),
+ SPLIT2(solve_args.uv),
+ r, ref, T.E, T.SE);
+ break;
+ default: FATAL("Unreachable code.\n"); break;
+ }
+ printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N);
+ printf("Elapsed time = %s\n", dump);
+ printf("Time per realisation (in usec) = %g +/- %g\n\n", time.E, time.SE);
+
+ CHK(nfails + nreals == N);
+ CHK(nfails <= N/1000);
+ CHK(eq_eps(T.E, ref, T.SE * 3));
+
+ OK(sdis_estimator_ref_put(estimator));
+ }
+}
+static void
+solve
+ (struct sdis_scene* scn,
+ struct interf* interf_props,
+ struct ssp_rng* rng)
+{
+ char dump[128];
+ struct time t0, t1;
+ struct sdis_estimator* estimator;
+ struct sdis_solve_boundary_args solve_args = SDIS_SOLVE_BOUNDARY_ARGS_DEFAULT;
+ struct sdis_mc T = SDIS_MC_NULL;
+ struct sdis_mc time = SDIS_MC_NULL;
+ const enum sdis_side all_front[] = { SDIS_FRONT, SDIS_FRONT };
+ const enum sdis_side all_back[] = { SDIS_BACK, SDIS_BACK };
+ size_t plist[2];
+ size_t nreals;
+ size_t nfails;
+ double ref_L, ref_R;
+ enum sdis_scene_dimension dim;
+ int nsimuls;
+ int isimul;
+ ASSERT(scn && interf_props && rng);
+
+ OK(sdis_scene_get_dimension(scn, &dim));
+
+ nsimuls = (dim == SDIS_SCENE_2D) ? 2 : 4;
+ FOR_EACH(isimul, 0, nsimuls) {
+ double ref;
+ double r = pow(10, ssp_rng_uniform_double(rng, -2, 2));
+
+ interf_props->resistance = r;
+
+ ref_L = (
+ T0 * (r * LAMBDA1 / X0) * (1 + r * LAMBDA2 / (L - X0))
+ + TL * (r * LAMBDA2 / (L - X0))
+ )
+ / ((1 + r * LAMBDA1 / X0) * (1 + r * LAMBDA2 / (L - X0)) - 1);
+
+ ref_R = ref_L * (1 + r * LAMBDA1 / X0) - T0 * r * LAMBDA1 / X0;
+
+ if(dim == SDIS_SCENE_2D) {
+ /* last segment */
+ solve_args.nprimitives = 1;
+ plist[0] = model2d_nsegments - 1;
+ } else {
+ /* last 2 triangles */
+ solve_args.nprimitives = 2;
+ plist[0] = model3d_ntriangles - 2;
+ plist[1] = model3d_ntriangles - 1;
+ }
+ solve_args.primitives = plist;
+
+ if(isimul < nsimuls / 2) {
+ solve_args.sides = all_front;
+ ref = ref_L;
+ } else {
+ solve_args.sides = all_back;
+ ref = ref_R;
+ }
+
+ solve_args.nrealisations = N;
+ solve_args.time_range[0] = solve_args.time_range[1] = INF;
+
+ time_current(&t0);
+ OK(sdis_solve_boundary(scn, &solve_args, &estimator));
+ time_sub(&t0, time_current(&t1), &t0);
+ time_dump(&t0, TIME_ALL, NULL, dump, sizeof(dump));
+
+ OK(sdis_estimator_get_realisation_count(estimator, &nreals));
+ OK(sdis_estimator_get_failure_count(estimator, &nfails));
+ OK(sdis_estimator_get_temperature(estimator, &T));
+ OK(sdis_estimator_get_realisation_time(estimator, &time));
+
+ printf("Steady temperature at the %s side with R=%g = %g ~ %g +/- %g\n",
+ (solve_args.sides[0] == SDIS_FRONT ? "front" : "back"),
+ r, ref, T.E, T.SE);
+ printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N);
+ printf("Elapsed time = %s\n", dump);
+ printf("Time per realisation (in usec) = %g +/- %g\n\n", time.E, time.SE);
+
+ CHK(nfails + nreals == N);
+ CHK(nfails <= N / 1000);
+ CHK(eq_eps(T.E, ref, T.SE * 3));
+
+ OK(sdis_estimator_ref_put(estimator));
+ }
+}
+
+/*******************************************************************************
+ * Test
+ ******************************************************************************/
+int
+main(int argc, char** argv)
+{
+ struct mem_allocator allocator;
+ struct sdis_data* data = NULL;
+ struct sdis_device* dev = NULL;
+ struct sdis_medium* fluid = NULL;
+ struct sdis_medium* solid1 = NULL;
+ struct sdis_medium* solid2 = NULL;
+ struct sdis_interface* interf_adiabatic1 = NULL;
+ struct sdis_interface* interf_adiabatic2 = NULL;
+ struct sdis_interface* interf_T0 = NULL;
+ struct sdis_interface* interf_TL = NULL;
+ struct sdis_interface* interf_R = NULL;
+ struct sdis_scene* box_scn = NULL;
+ struct sdis_scene* square_scn = NULL;
+ struct sdis_scene_create_args scn_args = SDIS_SCENE_CREATE_ARGS_DEFAULT;
+ 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* model3d_interfaces[22 /*#triangles*/];
+ struct sdis_interface* model2d_interfaces[7/*#segments*/];
+ struct interf* interf_props = NULL;
+ struct solid* solid_props = NULL;
+ struct ssp_rng* rng = NULL;
+ (void)argc, (void)argv;
+
+ OK(mem_init_proxy_allocator(&allocator, &mem_default_allocator));
+ OK(sdis_device_create(NULL, &allocator, SDIS_NTHREADS_DEFAULT, 1, &dev));
+
+ fluid_shader.temperature = fluid_get_temperature;
+ OK(sdis_fluid_create(dev, &fluid_shader, NULL, &fluid));
+
+ /* Setup the solid shader */
+ solid_shader.calorific_capacity = solid_get_calorific_capacity;
+ solid_shader.thermal_conductivity = solid_get_thermal_conductivity;
+ solid_shader.volumic_mass = solid_get_volumic_mass;
+ solid_shader.delta_solid = solid_get_delta;
+ solid_shader.temperature = solid_get_temperature;
+
+ /* Create the solid medium #1 */
+ OK(sdis_data_create(dev, sizeof(struct solid), 16, NULL, &data));
+ solid_props = sdis_data_get(data);
+ solid_props->lambda = LAMBDA1;
+ solid_props->cp = 2;
+ solid_props->rho = 25;
+ solid_props->delta = DELTA1;
+ OK(sdis_solid_create(dev, &solid_shader, data, &solid1));
+ OK(sdis_data_ref_put(data));
+
+ /* Create the solid medium #2 */
+ OK(sdis_data_create(dev, sizeof(struct solid), 16, NULL, &data));
+ solid_props = sdis_data_get(data);
+ solid_props->lambda = LAMBDA2;
+ solid_props->cp = 2;
+ solid_props->rho = 25;
+ solid_props->delta = DELTA2;
+ OK(sdis_solid_create(dev, &solid_shader, data, &solid2));
+ OK(sdis_data_ref_put(data));
+
+ /* Setup the interface shader */
+ interf_shader.front.temperature = interface_get_temperature;
+ interf_shader.back.temperature = interface_get_temperature;
+ interf_shader.convection_coef = interface_get_convection_coef;
+
+ /* Create the adiabatic interfaces */
+ OK(sdis_data_create(dev, sizeof(struct interf), 16, NULL, &data));
+ interf_props = sdis_data_get(data);
+ interf_props->temperature = UNKNOWN_TEMPERATURE;
+ OK(sdis_interface_create
+ (dev, solid1, fluid, &interf_shader, data, &interf_adiabatic1));
+ OK(sdis_interface_create
+ (dev, solid2, fluid, &interf_shader, data, &interf_adiabatic2));
+ OK(sdis_data_ref_put(data));
+
+ /* Create the T0 interface */
+ OK(sdis_data_create(dev, sizeof(struct interf), 16, NULL, &data));
+ interf_props = sdis_data_get(data);
+ interf_props->temperature = T0;
+ OK(sdis_interface_create
+ (dev, solid1, fluid, &interf_shader, data, &interf_T0));
+ OK(sdis_data_ref_put(data));
+
+ /* Create the TL interface */
+ OK(sdis_data_create(dev, sizeof(struct interf), 16, NULL, &data));
+ interf_props = sdis_data_get(data);
+ interf_props->temperature = TL;
+ OK(sdis_interface_create
+ (dev, solid2, fluid, &interf_shader, data, &interf_TL));
+ OK(sdis_data_ref_put(data));
+
+ /* Create the solid1-solid2 interface */
+ interf_shader.convection_coef = NULL;
+ interf_shader.thermal_contact_resistance = interface_get_contact_resistance;
+ OK(sdis_data_create(dev, sizeof(struct interf), 16, NULL, &data));
+ interf_props = sdis_data_get(data);
+ interf_props->temperature = UNKNOWN_TEMPERATURE;
+ OK(sdis_interface_create
+ (dev, solid1, solid2, &interf_shader, data, &interf_R));
+ OK(sdis_data_ref_put(data));
+
+ /* Release the media */
+ OK(sdis_medium_ref_put(solid1));
+ OK(sdis_medium_ref_put(solid2));
+ OK(sdis_medium_ref_put(fluid));
+
+ /* Map the interfaces to their box triangles */
+ /* Front */
+ model3d_interfaces[0] = interf_adiabatic1;
+ model3d_interfaces[1] = interf_adiabatic1;
+ model3d_interfaces[2] = interf_adiabatic2;
+ model3d_interfaces[3] = interf_adiabatic2;
+ /* Left */
+ model3d_interfaces[4] = interf_T0;
+ model3d_interfaces[5] = interf_T0;
+ /* Back */
+ model3d_interfaces[6] = interf_adiabatic1;
+ model3d_interfaces[7] = interf_adiabatic1;
+ model3d_interfaces[8] = interf_adiabatic2;
+ model3d_interfaces[9] = interf_adiabatic2;
+ /* Right */
+ model3d_interfaces[10] = interf_TL;
+ model3d_interfaces[11] = interf_TL;
+ /* Top */
+ model3d_interfaces[12] = interf_adiabatic1;
+ model3d_interfaces[13] = interf_adiabatic1;
+ model3d_interfaces[14] = interf_adiabatic2;
+ model3d_interfaces[15] = interf_adiabatic2;
+ /* Bottom */
+ model3d_interfaces[16] = interf_adiabatic1;
+ model3d_interfaces[17] = interf_adiabatic1;
+ model3d_interfaces[18] = interf_adiabatic2;
+ model3d_interfaces[19] = interf_adiabatic2;
+ /* Inside */
+ model3d_interfaces[20] = interf_R;
+ model3d_interfaces[21] = interf_R;
+
+ /* Map the interfaces to their square segments */
+ /* Bottom */
+ model2d_interfaces[0] = interf_adiabatic2;
+ model2d_interfaces[1] = interf_adiabatic1;
+ /* Left */
+ model2d_interfaces[2] = interf_T0;
+ /* Top */
+ model2d_interfaces[3] = interf_adiabatic1;
+ model2d_interfaces[4] = interf_adiabatic2;
+ /* Right */
+ model2d_interfaces[5] = interf_TL;
+ /* Contact */
+ model2d_interfaces[6] = interf_R;
+
+ /* Create the box scene */
+ scn_args.get_indices = model3d_get_indices;
+ scn_args.get_interface = model3d_get_interface;
+ scn_args.get_position = model3d_get_position;
+ scn_args.nprimitives = model3d_ntriangles;
+ scn_args.nvertices = model3d_nvertices;
+ scn_args.context = model3d_interfaces;
+ OK(sdis_scene_create(dev, &scn_args, &box_scn));
+
+ /* Create the square scene */
+ scn_args.get_indices = model2d_get_indices;
+ scn_args.get_interface = model2d_get_interface;
+ scn_args.get_position = model2d_get_position;
+ scn_args.nprimitives = model2d_nsegments;
+ scn_args.nvertices = model2d_nvertices;
+ scn_args.context = model2d_interfaces;
+ OK(sdis_scene_2d_create(dev, &scn_args, &square_scn));
+
+ /* Release the interfaces */
+ OK(sdis_interface_ref_put(interf_adiabatic1));
+ OK(sdis_interface_ref_put(interf_adiabatic2));
+ OK(sdis_interface_ref_put(interf_T0));
+ OK(sdis_interface_ref_put(interf_TL));
+ OK(sdis_interface_ref_put(interf_R));
+
+ /* Solve */
+ OK(ssp_rng_create(&allocator, &ssp_rng_kiss, &rng));
+ printf(">> Box scene\n");
+ solve_probe(box_scn, interf_props, rng);
+ solve(box_scn, interf_props, rng);
+ printf("\n>> Square scene\n");
+ solve_probe(square_scn, interf_props, rng);
+ solve(square_scn, interf_props, rng);
+
+ OK(sdis_scene_ref_put(box_scn));
+ OK(sdis_scene_ref_put(square_scn));
+ OK(sdis_device_ref_put(dev));
+ OK(ssp_rng_ref_put(rng));
+
+ check_memory_allocator(&allocator);
+ mem_shutdown_proxy_allocator(&allocator);
+ CHK(mem_allocated_size() == 0);
+ return 0;
+}
+
diff --git a/src/test_sdis_utils.h b/src/test_sdis_utils.h
@@ -191,6 +191,7 @@ static const struct sdis_fluid_shader DUMMY_FLUID_SHADER = {
static const struct sdis_interface_shader DUMMY_INTERFACE_SHADER = {
dummy_interface_getter,
0,
+ dummy_interface_getter,
DUMMY_INTERFACE_SIDE_SHADER__,
DUMMY_INTERFACE_SIDE_SHADER__
};
