commit c28a073d6d70d8e90354d73172a699f0cb64088d
parent aa1eed1b81eeb2d121639cd61454d16a1bf4bfff
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Mon, 24 Jun 2024 10:33:50 +0200
Merge branch 'feature_custom_solid_path_sampling' into develop
Diffstat:
21 files changed, 2610 insertions(+), 69 deletions(-)
diff --git a/Makefile b/Makefile
@@ -43,6 +43,7 @@ SRC =\
src/sdis_log.c\
src/sdis_medium.c\
src/sdis_misc.c\
+ src/sdis_primkey.c\
src/sdis_radiative_env.c\
src/sdis_realisation.c\
src/sdis_scene.c\
@@ -192,6 +193,8 @@ TEST_SRC =\
src/test_sdis_interface.c\
src/test_sdis_medium.c\
src/test_sdis_picard.c\
+ src/test_sdis_primkey.c\
+ src/test_sdis_primkey_2d.c\
src/test_sdis_radiative_env.c\
src/test_sdis_scene.c\
src/test_sdis_solid_random_walk_robustness.c\
@@ -216,6 +219,8 @@ TEST_SRC_LONG =\
TEST_SRC_MPI =\
src/test_sdis.c\
src/test_sdis_compute_power.c\
+ src/test_sdis_custom_solid_path_sampling.c\
+ src/test_sdis_custom_solid_path_sampling_2d.c\
src/test_sdis_device.c\
src/test_sdis_external_flux.c\
src/test_sdis_solve_camera.c\
@@ -396,6 +401,7 @@ test_sdis_volumic_power4 \
################################################################################
src/test_sdis_draw_external_flux.d \
src/test_sdis_external_flux_with_diffuse_radiance.d \
+src/test_sdis_primkey.d \
src/test_sdis_solid_random_walk_robustness.d \
src/test_sdis_solve_probe3.d \
src/test_sdis_unsteady_analytic_profile.d \
@@ -404,6 +410,7 @@ src/test_sdis_unsteady_analytic_profile.d \
src/test_sdis_draw_external_flux.o \
src/test_sdis_external_flux_with_diffuse_radiance.o \
+src/test_sdis_primkey.o \
src/test_sdis_solid_random_walk_robustness.o \
src/test_sdis_solve_probe3.o \
src/test_sdis_unsteady_analytic_profile.o \
@@ -412,6 +419,7 @@ src/test_sdis_unsteady_analytic_profile.o \
test_sdis_draw_external_flux \
test_sdis_external_flux_with_diffuse_radiance \
+test_sdis_primkey \
test_sdis_solid_random_walk_robustness \
test_sdis_solve_probe3 \
test_sdis_unsteady_analytic_profile \
@@ -419,6 +427,36 @@ test_sdis_unsteady_analytic_profile \
$(CC) $(TEST_CFLAGS) $(S3DUT_CFLAGS) -o $@ src/$@.o $(TEST_LIBS) $(S3DUT_LIBS)
################################################################################
+# Test based on Star-2D
+################################################################################
+src/test_sdis_primkey_2d.d \
+: config.mk sdis-local.pc
+ @$(CC) $(TEST_CFLAGS) $(S2D_CFLAGS) -MM -MT "$(@:.d=.o) $@" $(@:.d=.c) -MF $@
+
+src/test_sdis_primkey_2d.o \
+: config.mk sdis-local.pc
+ $(CC) $(TEST_CFLAGS) $(S2D_CFLAGS) -c $(@:.o=.c) -o $@
+
+test_sdis_primkey_2d \
+: config.mk sdis-local.pc $(LIBNAME) src/test_sdis_utils.o
+ $(CC) $(TEST_CFLAGS) $(S2D_CFLAGS) -o $@ src/$@.o $(TEST_LIBS) $(S2D_LIBS)
+
+################################################################################
+# Test based on Star-2D with (optional) MPI support
+################################################################################
+src/test_sdis_custom_solid_path_sampling_2d.d \
+: config.mk sdis-local.pc
+ @$(CC) $(TEST_CFLAGS_MPI) $(S2D_CFLAGS) -MM -MT "$(@:.d=.o) $@" $(@:.d=.c) -MF $@
+
+src/test_sdis_custom_solid_path_sampling_2d.o \
+: config.mk sdis-local.pc
+ $(CC) $(TEST_CFLAGS_MPI) $(S2D_CFLAGS) -c $(@:.o=.c) -o $@
+
+test_sdis_custom_solid_path_sampling_2d \
+: config.mk sdis-local.pc $(LIBNAME) src/test_sdis_utils.o
+ $(CC) $(TEST_CFLAGS_MPI) $(S2D_CFLAGS) -o $@ src/$@.o $(TEST_LIBS_MPI) $(S2D_LIBS)
+
+################################################################################
# Tests based on Star-Enclosures-<2|3>D
################################################################################
src/test_sdis_scene.d \
@@ -482,19 +520,24 @@ test_sdis_solve_probe_boundary_list \
$(CC) $(TEST_CFLAGS_MPI) $(S3DUT_CFLAGS) -o $@ src/$@.o $(TEST_LIBS_MPI) $(S3DUT_LIBS)
################################################################################
-# Tests based on Star-3D and Star-3DUT with (optional) MPI support
+# Tests based on Star-3D, Star-3DUT and Star-SP with (optional) MPI support
################################################################################
+src/test_sdis_custom_solid_path_sampling.d \
src/test_sdis_solve_probe_list.d \
: config.mk sdis-local.pc
- @$(CC) $(TEST_CFLAGS_MPI) $(S3D_CFLAGS) $(S3DUT_CFLAGS) -MM -MT "$(@:.d=.o) $@" $(@:.d=.c) -MF $@
+ @$(CC) $(TEST_CFLAGS_MPI) $(S3D_CFLAGS) $(S3DUT_CFLAGS) $(SSP_CFLAGS) \
+ -MM -MT "$(@:.d=.o) $@" $(@:.d=.c) -MF $@
+src/test_sdis_custom_solid_path_sampling.o \
src/test_sdis_solve_probe_list.o \
: config.mk sdis-local.pc
- $(CC) $(TEST_CFLAGS_MPI) $(S3D_CFLAGS) $(S3DUT_CFLAGS) -c $(@:.o=.c) -o $@
+ $(CC) $(TEST_CFLAGS_MPI) $(S3D_CFLAGS) $(S3DUT_CFLAGS) $(SSP_CFLAGS) -c $(@:.o=.c) -o $@
+test_sdis_custom_solid_path_sampling \
test_sdis_solve_probe_list \
: config.mk sdis-local.pc $(LIBNAME) src/test_sdis_utils.o
- $(CC) $(TEST_CFLAGS_MPI) $(S3D_CFLAGS) $(S3DUT_CFLAGS) -o $@ src/$@.o $(TEST_LIBS_MPI) $(S3D_LIBS) $(S3DUT_LIBS)
+ $(CC) $(TEST_CFLAGS_MPI) $(S3D_CFLAGS) $(S3DUT_CFLAGS) $(SSP_CFLAGS) \
+ -o $@ src/$@.o $(TEST_LIBS_MPI) $(S3D_LIBS) $(S3DUT_LIBS) $(SSP_CFLAGS)
################################################################################
# Tests based on Star-SP with (optional) MPI support
diff --git a/src/sdis.h b/src/sdis.h
@@ -220,6 +220,53 @@ struct sdis_scene_find_closest_point_args {
static const struct sdis_scene_find_closest_point_args
SDIS_SCENE_FIND_CLOSEST_POINT_ARGS_NULL = SDIS_SCENE_FIND_CLOSEST_POINT_ARGS_NULL__;
+/* A sampled path */
+struct sdis_path {
+ struct sdis_rwalk_vertex vtx; /* Current position and time */
+
+ /* Surface intersection. When defined, the path is on a border */
+ struct s2d_primitive prim_2d;
+ struct s3d_primitive prim_3d;
+
+ double elapsed_time; /* Time elapsed along the path */
+ double weight; /* Monte Carlo weight update along the path */
+
+ /* Define whether the path has reached a boundary condition in time/space */
+ int at_limit;
+};
+#define SDIS_PATH_NULL__ { \
+ SDIS_RWALK_VERTEX_NULL__, \
+ S2D_PRIMITIVE_NULL__, \
+ S3D_PRIMITIVE_NULL__, \
+ 0, /* Elapsed time */ \
+ 0, /* MC weight */ \
+ 0 /* At limit */ \
+}
+static const struct sdis_path SDIS_PATH_NULL = SDIS_PATH_NULL__;
+
+/* Type of functor used by the user to write the way in which the path is
+ * sampled. So it's no longer Stardis that samples the path, but the user
+ * through his own function. */
+typedef res_T
+(*sdis_sample_path_T)
+ (struct sdis_scene* scn,
+ struct ssp_rng* rng,
+ struct sdis_path* path,
+ struct sdis_data* data);
+
+/* Key to a geometric primitive, i.e its unique identifier. Its member variables
+ * must be treated as private variables, i.e. the caller must not access them
+ * directly but use the primkey API functions instead (see below) */
+struct sdis_primkey {
+ /* List of primitive nodes sorted in ascending order */
+ double nodes[9];
+
+ /* Overall number of coordinates (4 in 2D, 9 in 3D) */
+ unsigned ncoords;
+};
+#define SDIS_PRIMKEY_NULL__ {{0,0,0,0,0,0,0,0,0},0}
+static const struct sdis_primkey SDIS_PRIMKEY_NULL = SDIS_PRIMKEY_NULL__;
+
/*******************************************************************************
* Estimation data types
******************************************************************************/
@@ -294,12 +341,16 @@ struct sdis_solid_shader {
* This getter is always called at time >= t0 (see below). */
sdis_medium_getter_T temperature;
+ /* Function to be used to sample the path through the solid. If not defined,
+ * let stardis sample a conductive path */
+ sdis_sample_path_T sample_path;
+
/* The time until the initial condition is maintained for this solid.
* Can be negative or set to +/- infinity to simulate a system that is always
* in the initial state or never reaches it, respectively. */
double t0;
};
-#define SDIS_SOLID_SHADER_NULL__ {NULL, NULL, NULL, NULL, NULL, NULL, 0}
+#define SDIS_SOLID_SHADER_NULL__ {NULL, NULL, NULL, NULL, NULL, NULL, NULL, 0}
static const struct sdis_solid_shader SDIS_SOLID_SHADER_NULL =
SDIS_SOLID_SHADER_NULL__;
@@ -1311,18 +1362,30 @@ sdis_scene_boundary_project_position
const double pos[],
double uv[]);
-/* Get the 2D scene's enclosures. Only defined for a 2D scene. */
+/* Get Star-Enclosure-2D scene. Defined on 2D scene only */
SDIS_API res_T
sdis_scene_get_senc2d_scene
(struct sdis_scene* scn,
struct senc2d_scene** senc2d_scn);
-/* Get the 3D scene's enclosures. Only defined for a 3D scene. */
+/* Get Star-Enclosure-3D scene. Defined on 3D scene only */
SDIS_API res_T
sdis_scene_get_senc3d_scene
(struct sdis_scene* scn,
struct senc3d_scene** senc3d_scn);
+/* Get Star-2D scene view. Defined on 2D scene only */
+SDIS_API res_T
+sdis_scene_get_s2d_scene_view
+ (struct sdis_scene* scn,
+ struct s2d_scene_view** s2d_view);
+
+/* Get Star-3D scene view. Defined on 3D scene only */
+SDIS_API res_T
+sdis_scene_get_s3d_scene_view
+ (struct sdis_scene* scn,
+ struct s3d_scene_view** s3d_view);
+
SDIS_API res_T
sdis_scene_get_dimension
(const struct sdis_scene* scn,
@@ -1353,6 +1416,20 @@ sdis_scene_get_radiative_env
/* The returned pointer can be NULL, i.e. there is no radiative environement*/
struct sdis_radiative_env** radenv);
+/* Get the internal Star-2D primitive corresponding to the primitive key */
+SDIS_API res_T
+sdis_scene_get_s2d_primitive
+ (struct sdis_scene* scn,
+ const struct sdis_primkey* key,
+ struct s2d_primitive* primitive);
+
+/* Get the internal Star-3D primitive corresponding to the primitive key */
+SDIS_API res_T
+sdis_scene_get_s3d_primitive
+ (struct sdis_scene* scn,
+ const struct sdis_primkey* key,
+ struct s3d_primitive* primitive);
+
/*******************************************************************************
* An estimator stores the state of a simulation
******************************************************************************/
@@ -1717,6 +1794,31 @@ SDIS_API res_T
sdis_get_info
(struct sdis_info* info);
+/*******************************************************************************
+ * Primitive identifier
+ ******************************************************************************/
+SDIS_API void
+sdis_primkey_setup
+ (struct sdis_primkey* key,
+ const double node0[3],
+ const double node1[3],
+ const double node2[3]);
+
+SDIS_API void
+sdis_primkey_2d_setup
+ (struct sdis_primkey* key,
+ const double node0[2],
+ const double node1[2]);
+
+SDIS_API size_t
+sdis_primkey_hash
+ (const struct sdis_primkey* key);
+
+SDIS_API char
+sdis_primkey_eq
+ (const struct sdis_primkey* key0,
+ const struct sdis_primkey* key1);
+
END_DECLS
#endif /* SDIS_H */
diff --git a/src/sdis_Xd_begin.h b/src/sdis_Xd_begin.h
@@ -65,6 +65,7 @@
#define SXD_GET_PRIMITIVE CONCAT(CONCAT(S, DIM), D_GET_PRIMITIVE)
#define SXD_SAMPLE CONCAT(CONCAT(S, DIM), D_SAMPLE)
#define SXD_TRACE CONCAT(CONCAT(S, DIM), D_TRACE)
+#define SXD_PRIMITIVE_NULL CONCAT(CONCAT(S, DIM), D_PRIMITIVE_NULL)
#define SXD_PRIMITIVE_EQ CONCAT(CONCAT(S, DIM), D_PRIMITIVE_EQ)
/* Vector macros generic to SDIS_XD_DIMENSION */
diff --git a/src/sdis_Xd_end.h b/src/sdis_Xd_end.h
@@ -34,6 +34,8 @@
#undef SXD_GET_PRIMITIVE
#undef SXD_SAMPLE
#undef SXD_TRACE
+#undef SXD_PRIMITIVE_NULL
+#undef SXD_PRIMITIVE_EQ
#undef dX
#undef fX
diff --git a/src/sdis_heat_path_conductive.c b/src/sdis_heat_path_conductive.c
@@ -13,10 +13,10 @@
* 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_heat_path.h"
-#include "sdis_heat_path_conductive_c.h"
#include "sdis_log.h"
+#include "sdis_heat_path_conductive_c.h"
#include "sdis_medium_c.h"
+#include "sdis_scene_c.h"
/*******************************************************************************
* Local function
@@ -71,54 +71,16 @@ error:
goto exit;
}
-res_T
-conductive_path_2d
- (struct sdis_scene* scn,
- struct rwalk_context* ctx,
- struct rwalk* rwalk,
- struct ssp_rng* rng,
- struct temperature* T)
-{
- res_T res = RES_OK;
- ASSERT(ctx);
-
- switch(ctx->diff_algo) {
- case SDIS_DIFFUSION_DELTA_SPHERE:
- res = conductive_path_delta_sphere_2d(scn, ctx, rwalk, rng, T);
- break;
- case SDIS_DIFFUSION_WOS:
- res = conductive_path_wos_2d(scn, ctx, rwalk, rng, T);
- break;
- default: FATAL("Unreachable code.\n"); break;
- }
- return res;
-}
-
-res_T
-conductive_path_3d
- (struct sdis_scene* scn,
- struct rwalk_context* ctx,
- struct rwalk* rwalk,
- struct ssp_rng* rng,
- struct temperature* T)
-{
- res_T res = RES_OK;
- ASSERT(ctx);
-
- switch(ctx->diff_algo) {
- case SDIS_DIFFUSION_DELTA_SPHERE:
- res = conductive_path_delta_sphere_3d(scn, ctx, rwalk, rng, T);
- break;
- case SDIS_DIFFUSION_WOS:
- res = conductive_path_wos_3d(scn, ctx, rwalk, rng, T);
- break;
- default: FATAL("Unreachable code.\n"); break;
- }
- return res;
-}
-
/* Generate the conductive path functions */
#define SDIS_XD_DIMENSION 2
+#include "sdis_heat_path_conductive_Xd.h"
+#define SDIS_XD_DIMENSION 3
+#include "sdis_heat_path_conductive_Xd.h"
+#define SDIS_XD_DIMENSION 2
+#include "sdis_heat_path_conductive_custom_Xd.h"
+#define SDIS_XD_DIMENSION 3
+#include "sdis_heat_path_conductive_custom_Xd.h"
+#define SDIS_XD_DIMENSION 2
#include "sdis_heat_path_conductive_delta_sphere_Xd.h"
#define SDIS_XD_DIMENSION 3
#include "sdis_heat_path_conductive_delta_sphere_Xd.h"
diff --git a/src/sdis_heat_path_conductive_Xd.h b/src/sdis_heat_path_conductive_Xd.h
@@ -0,0 +1,73 @@
+/* Copyright (C) 2016-2024 |Méso|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 "sdis_heat_path.h"
+#include "sdis_heat_path_conductive_c.h"
+#include "sdis_log.h"
+#include "sdis_medium_c.h"
+#include "sdis_scene_c.h"
+
+#include <rsys/cstr.h>
+
+#include "sdis_Xd_begin.h"
+
+/*******************************************************************************
+ * Local functions
+ ******************************************************************************/
+res_T
+XD(conductive_path)
+ (struct sdis_scene* scn,
+ struct rwalk_context* ctx,
+ struct rwalk* rwalk,
+ struct ssp_rng* rng,
+ struct temperature* T)
+{
+ struct sdis_medium* mdm = NULL;
+ unsigned enc_id = ENCLOSURE_ID_NULL;
+ res_T res = RES_OK;
+ ASSERT(ctx && rwalk);
+
+ res = scene_get_enclosure_id_in_closed_boundaries(scn, rwalk->vtx.P, &enc_id);
+ if(res != RES_OK) goto error;
+ res = scene_get_enclosure_medium(scn, scene_get_enclosure(scn, enc_id), &mdm);
+ if(res != RES_OK) goto error;
+ ASSERT(sdis_medium_get_type(mdm) == SDIS_SOLID);
+
+ /* The caller defined a custom function to sample the solid path */
+ if(mdm->shader.solid.sample_path != NULL) {
+ res = XD(conductive_path_custom)(scn, enc_id, mdm, rwalk, rng, T);
+ if(res != RES_OK) goto error;
+
+ } else {
+ switch(ctx->diff_algo) {
+ case SDIS_DIFFUSION_DELTA_SPHERE:
+ res = XD(conductive_path_delta_sphere)(scn, ctx, rwalk, rng, T);
+ break;
+ case SDIS_DIFFUSION_WOS:
+ res = XD(conductive_path_wos)(scn, ctx, rwalk, rng, T);
+ break;
+ default: FATAL("Unreachable code.\n"); break;
+ }
+ if(res != RES_OK) goto error;
+ }
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
+#include "sdis_Xd_end.h"
diff --git a/src/sdis_heat_path_conductive_c.h b/src/sdis_heat_path_conductive_c.h
@@ -22,6 +22,7 @@
struct rwalk;
struct rwalk_context;
struct sdis_device;
+struct sdis_medium;
struct sdis_scene;
struct solid_props;
struct ssp_rng;
@@ -36,6 +37,27 @@ check_solid_constant_properties
const struct solid_props* props);
/*******************************************************************************
+ * Conductive paths using custom user algorithm
+ ******************************************************************************/
+extern LOCAL_SYM res_T
+conductive_path_custom_2d
+ (struct sdis_scene* scn,
+ const unsigned enc_id, /* Enclosure in which path is sampled */
+ const struct sdis_medium* mdm, /* Medium in which path is sampled */
+ struct rwalk* rwalk,
+ struct ssp_rng* rng,
+ struct temperature* T);
+
+extern LOCAL_SYM res_T
+conductive_path_custom_3d
+ (struct sdis_scene* scn,
+ const unsigned enc_id, /* Enclosure in which path is sampled */
+ const struct sdis_medium* mdm, /* Medium in which path is sampled */
+ struct rwalk* rwalk,
+ struct ssp_rng* rng,
+ struct temperature* T);
+
+/*******************************************************************************
* Conductive paths using the delta sphere algorithm
******************************************************************************/
extern LOCAL_SYM res_T
diff --git a/src/sdis_heat_path_conductive_custom_Xd.h b/src/sdis_heat_path_conductive_custom_Xd.h
@@ -0,0 +1,194 @@
+/* Copyright (C) 2016-2024 |Méso|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 "sdis_log.h"
+#include "sdis_medium_c.h"
+#include "sdis_scene_c.h"
+
+#include <rsys/cstr.h>
+
+#include "sdis_Xd_begin.h"
+
+/*******************************************************************************
+ * Helper function
+ ******************************************************************************/
+static res_T
+XD(check_sampled_path)
+ (struct sdis_scene* scn,
+ const struct sdis_path* path)
+{
+ int null_prim = 0;
+ res_T res = RES_OK;
+ ASSERT(scn && path);
+
+ null_prim = SXD_PRIMITIVE_EQ(&path->XD(prim), &SXD_PRIMITIVE_NULL);
+
+ /* Check end of path */
+ if(!path->at_limit && null_prim) {
+ log_err(scn->dev,
+ "%s: the sampled path should have reached a limit condition or a boundary"
+ " -- pos=("FORMAT_VECX")\n",
+ FUNC_NAME, SPLITX(path->vtx.P));
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ if(!null_prim) {
+ struct sXd(primitive) prim;
+ const unsigned iprim = path->XD(prim).prim_id;
+
+ /* Check intersected primitive */
+ res = sXd(scene_view_get_primitive)(scn->sXd(view), iprim, &prim);
+ if(res != RES_OK || !SXD_PRIMITIVE_EQ(&path->XD(prim), &prim)) {
+ log_err(scn->dev,
+ "%s: invalid intersected primitive on sampled path -- %s\n",
+ FUNC_NAME, res_to_cstr(res));
+ goto error;
+ }
+ }
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
+static int
+XD(keep_only_one_primitive)
+ (const struct sXd(hit)* hit,
+ const float org[DIM],
+ const float dir[DIM],
+ const float range[2],
+ void* query_data,
+ void* filter_data)
+{
+ const struct sXd(primitive)* prim = query_data;
+ (void)org, (void)dir, (void)range, (void)filter_data;
+ return !SXD_PRIMITIVE_EQ(prim, &hit->prim);
+}
+
+static res_T
+XD(get_path_hit)
+ (struct sdis_scene* scn,
+ struct sdis_path* path,
+ const struct sdis_medium* mdm,
+ struct sXd(hit)* hit)
+{
+ struct hit_filter_data filter_data = HIT_FILTER_DATA_NULL;
+ float query_radius = 0;
+ float query_pos[DIM] = {0};
+ double delta = 0;
+ res_T res = RES_OK;
+ ASSERT(scn && path && hit);
+
+ filter_data.XD(custom_filter) = XD(keep_only_one_primitive);
+ filter_data.custom_filter_data = &path->XD(prim);
+
+ /* Search for the hit corresponding to the path position on a primitive. Search
+ * at a maximum distance from the delta of the medium, as this hit should be
+ * very close to the submitted position, since it should represent the same
+ * point. */
+ delta = solid_get_delta(mdm, &path->vtx);
+ fX_set_dX(query_pos, path->vtx.P);
+ query_radius = (float)delta;
+ SXD(scene_view_closest_point(scn->sXd(view), query_pos, query_radius,
+ &filter_data, hit));
+ ASSERT(SXD_PRIMITIVE_EQ(&hit->prim, &path->XD(prim)));
+
+ if(SXD_HIT_NONE(hit)) {
+ log_warn(scn->dev,
+ "%s: the position returned by custom sampling of the conductive path "
+ "is too far from the primitive it should be on -- search distance=%g\n",
+ FUNC_NAME, delta);
+ res = RES_BAD_OP;
+ goto error;
+ }
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
+/*******************************************************************************
+ * Local function
+ ******************************************************************************/
+res_T
+XD(conductive_path_custom)
+ (struct sdis_scene* scn,
+ const unsigned enc_id,
+ const struct sdis_medium* mdm,
+ struct rwalk* rwalk,
+ struct ssp_rng* rng,
+ struct temperature* T)
+{
+ struct sdis_path path = SDIS_PATH_NULL;
+ res_T res = RES_OK;
+
+ /* Check pre-conditions */
+ ASSERT(scn && rwalk && rng && T);
+ ASSERT(sdis_medium_get_type(mdm) == SDIS_SOLID);
+ ASSERT(mdm->shader.solid.sample_path);
+
+ /* Sample a conductive path */
+ path.vtx = rwalk->vtx;
+ res = mdm->shader.solid.sample_path(scn, rng, &path, mdm->data);
+ if(res != RES_OK) {
+ log_err(scn->dev,
+ "%s: error in customized sampling of a conductive path "
+ "from pos=("FORMAT_VECX") -- %s\n",
+ FUNC_NAME, SPLITX(rwalk->vtx.P), res_to_cstr(res));
+ goto error;
+ }
+
+ res = XD(check_sampled_path)(scn, &path);
+ if(res!= RES_OK) goto error;
+
+ res = XD(get_path_hit)(scn, &path, mdm, &rwalk->XD(hit));
+ if(res != RES_OK) goto error;
+
+ /* Update random walk position and time from sampled path */
+ rwalk->vtx = path.vtx;
+ rwalk->elapsed_time += path.elapsed_time;
+
+ /* The path reached a boundary */
+ if(!SXD_HIT_NONE(&rwalk->XD(hit))) {
+ unsigned enc_ids[2] = {ENCLOSURE_ID_NULL, ENCLOSURE_ID_NULL};
+
+ rwalk->enc_id = ENCLOSURE_ID_NULL; /* Not in an enclosure */
+
+ /* Define which side of the interface the path is on */
+ scene_get_enclosure_ids(scn, rwalk->XD(hit).prim.prim_id, enc_ids);
+ if(enc_id == enc_ids[SDIS_FRONT]) rwalk->hit_side = SDIS_FRONT;
+ else if(enc_id == enc_ids[SDIS_BACK]) rwalk->hit_side = SDIS_BACK;
+ else FATAL("Unreachable code.\n");
+ }
+
+ /* Update Monte Carlo weight */
+ T->value += path.weight;
+ T->done = path.at_limit;
+
+ /* The path either reaches a boundary condition and will be stopped, or is
+ * on a boundary and a boundary path must be sampled */
+ T->func = XD(boundary_path);
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
+#include "sdis_Xd_end.h"
diff --git a/src/sdis_log.h b/src/sdis_log.h
@@ -28,6 +28,8 @@
#define MSG_ERROR_PREFIX_PLAIN_TEXT "stardis-solver (error): "
#define MSG_WARNING_PREFIX_PLAIN_TEXT "stardis-solver (warning): "
+struct sdis_device;
+
extern LOCAL_SYM res_T
setup_log_default
(struct sdis_device* dev);
diff --git a/src/sdis_primkey.c b/src/sdis_primkey.c
@@ -0,0 +1,113 @@
+/* Copyright (C) 2016-2024 |Méso|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 <rsys/double2.h>
+#include <rsys/double3.h>
+
+/*******************************************************************************
+ * Helper functions
+ ******************************************************************************/
+static INLINE int
+cmp_dbl2(const double v0[2], const double v1[2])
+{
+ if(v0[0] < v1[0]) return -1;
+ if(v0[0] > v1[0]) return +1;
+ if(v0[1] < v1[1]) return -1;
+ if(v0[1] > v1[1]) return +1;
+ return 0;
+}
+
+static INLINE int
+cmp_dbl3(const double v0[3], const double v1[3])
+{
+ if(v0[0] < v1[0]) return -1;
+ if(v0[0] > v1[0]) return +1;
+ if(v0[1] < v1[1]) return -1;
+ if(v0[1] > v1[1]) return +1;
+ if(v0[2] < v1[2]) return -1;
+ if(v0[2] > v1[2]) return +1;
+ return 0;
+}
+
+/*******************************************************************************
+ * Exported functions
+ ******************************************************************************/
+void
+sdis_primkey_2d_setup
+ (struct sdis_primkey* key,
+ const double node0[2],
+ const double node1[2])
+{
+ double *v0, *v1;
+ ASSERT(key && node0 && node1);
+ v0 = d2_set(key->nodes+0, node0);
+ v1 = d2_set(key->nodes+2, node1);
+ if(cmp_dbl2(v0, v1) > 0) {
+ SWAP(double, v0[0], v1[0]);
+ SWAP(double, v0[1], v1[1]);
+ }
+ key->ncoords = 4;
+}
+
+void
+sdis_primkey_setup
+ (struct sdis_primkey* key,
+ const double node0[3],
+ const double node1[3],
+ const double node2[3])
+{
+ double *v0, *v1, *v2;
+ ASSERT(key && node0 && node1 && node2);
+ v0 = d3_set(key->nodes+0, node0);
+ v1 = d3_set(key->nodes+3, node1);
+ v2 = d3_set(key->nodes+6, node2);
+
+ /* Bubble sort */
+ #define SWAP_DBL3(V0, V1) { \
+ SWAP(double, (V0)[0], (V1)[0]); \
+ SWAP(double, (V0)[1], (V1)[1]); \
+ SWAP(double, (V0)[2], (V1)[2]); \
+ } (void)0
+ if(cmp_dbl3(v0, v1) > 0) SWAP_DBL3(v0, v1);
+ if(cmp_dbl3(v1, v2) > 0) SWAP_DBL3(v1, v2);
+ if(cmp_dbl3(v0, v1) > 0) SWAP_DBL3(v0, v1);
+ #undef SWAP_DBL3
+ key->ncoords = 9;
+}
+
+size_t
+sdis_primkey_hash(const struct sdis_primkey* key)
+{
+ return hash_fnv64(key->nodes, sizeof(double)*key->ncoords);
+}
+
+char
+sdis_primkey_eq
+ (const struct sdis_primkey* key0,
+ const struct sdis_primkey* key1)
+{
+ unsigned i = 0;
+ ASSERT(key0 && key1);
+
+ if(key0->ncoords != key1->ncoords) return 0;
+ if(key0->ncoords != 4 && key0->ncoords != 9) return 0;
+ if(key1->ncoords != 4 && key1->ncoords != 9) return 0;
+ FOR_EACH(i, 0, key0->ncoords) {
+ if(key0->nodes[i] != key1->nodes[i]) return 0;
+ }
+ return 1;
+}
diff --git a/src/sdis_scene.c b/src/sdis_scene.c
@@ -103,6 +103,8 @@ scene_release(ref_T * ref)
darray_prim_prop_release(&scn->prim_props);
htable_enclosure_release(&scn->enclosures);
htable_d_release(&scn->tmp_hc_ub);
+ htable_key2prim2d_release(&scn->key2prim2d);
+ htable_key2prim3d_release(&scn->key2prim3d);
if(scn->s2d_view) S2D(scene_view_ref_put(scn->s2d_view));
if(scn->s3d_view) S3D(scene_view_ref_put(scn->s3d_view));
if(scn->senc2d_scn) SENC2D(scene_ref_put(scn->senc2d_scn));
@@ -332,7 +334,6 @@ sdis_scene_get_senc2d_scene
{
if(!scn || !senc2d_scn) return RES_BAD_ARG;
if(!scn->senc2d_scn) return RES_BAD_ARG; /* Scene is 3D */
- SENC2D(scene_ref_get(scn->senc2d_scn));
*senc2d_scn = scn->senc2d_scn;
return RES_OK;
}
@@ -344,12 +345,33 @@ sdis_scene_get_senc3d_scene
{
if(!scn || !senc3d_scn) return RES_BAD_ARG;
if(!scn->senc3d_scn) return RES_BAD_ARG; /* Scene is 2D */
- SENC3D(scene_ref_get(scn->senc3d_scn));
*senc3d_scn = scn->senc3d_scn;
return RES_OK;
}
res_T
+sdis_scene_get_s2d_scene_view
+ (struct sdis_scene* scn,
+ struct s2d_scene_view** s2d_view)
+{
+ if(!scn || !s2d_view) return RES_BAD_ARG;
+ if(!scn->s2d_view) return RES_BAD_ARG; /* Scene is 3D */
+ *s2d_view = scn->s2d_view;
+ return RES_OK;
+}
+
+res_T
+sdis_scene_get_s3d_scene_view
+ (struct sdis_scene* scn,
+ struct s3d_scene_view** s3d_view)
+{
+ if(!scn || !s3d_view) return RES_BAD_ARG;
+ if(!scn->s3d_view) return RES_BAD_ARG; /* Scene is 2D */
+ *s3d_view = scn->s3d_view;
+ return RES_OK;
+}
+
+res_T
sdis_scene_get_dimension
(const struct sdis_scene* scn, enum sdis_scene_dimension* dim)
{
@@ -416,6 +438,38 @@ sdis_scene_get_radiative_env
return RES_OK;
}
+res_T
+sdis_scene_get_s2d_primitive
+ (struct sdis_scene* scn,
+ const struct sdis_primkey* key,
+ struct s2d_primitive* out_prim)
+{
+ struct s2d_primitive* prim = NULL;
+
+ if(!scn || !key || !out_prim || !scene_is_2d(scn)) return RES_BAD_ARG;
+
+ if((prim = htable_key2prim2d_find(&scn->key2prim2d, key)) == NULL)
+ return RES_BAD_ARG;
+ *out_prim = *prim;
+ return RES_OK;
+}
+
+res_T
+sdis_scene_get_s3d_primitive
+ (struct sdis_scene* scn,
+ const struct sdis_primkey* key,
+ struct s3d_primitive* out_prim)
+{
+ struct s3d_primitive* prim = NULL;
+
+ if(!scn || !key || !out_prim || scene_is_2d(scn)) return RES_BAD_ARG;
+
+ if((prim = htable_key2prim3d_find(&scn->key2prim3d, key)) == NULL)
+ return RES_BAD_ARG;
+ *out_prim = *prim;
+ return RES_OK;
+}
+
/*******************************************************************************
* Local miscellaneous function
******************************************************************************/
diff --git a/src/sdis_scene_Xd.h b/src/sdis_scene_Xd.h
@@ -938,6 +938,89 @@ error:
goto exit;
}
+#if DIM == 2
+static res_T
+setup_primitive_keys_2d(struct sdis_scene* scn, struct senc2d_scene* senc_scn)
+{
+ unsigned iprim = 0;
+ unsigned nprims = 0;
+ res_T res = RES_OK;
+ ASSERT(scn && senc_scn);
+
+ SENC2D(scene_get_primitives_count(senc_scn, &nprims));
+
+ FOR_EACH(iprim, 0, nprims) {
+ struct s2d_primitive prim = S2D_PRIMITIVE_NULL;
+ struct sdis_primkey key = SDIS_PRIMKEY_NULL;
+ unsigned ids[2] = {0,0};
+ double v0[2] = {0,0};
+ double v1[2] = {0,0};
+
+ /* Retrieve positions from Star-Enclosre, not Star-2D. Star-Enclosure keeps
+ * the positions submitted by the user as they are, without any
+ * transformation or conversion (Star-2D converts them to float). This
+ * ensures that the caller can construct the same key from his data */
+ SENC2D(scene_get_primitive(senc_scn, iprim, ids));
+ SENC2D(scene_get_vertex(senc_scn, ids[0], v0));
+ SENC2D(scene_get_vertex(senc_scn, ids[1], v1));
+ S2D(scene_view_get_primitive(scn->s2d_view, iprim, &prim));
+
+ sdis_primkey_2d_setup(&key, v0, v1);
+
+ res = htable_key2prim2d_set(&scn->key2prim2d, &key, &prim);
+ if(res != RES_OK) goto error;
+ }
+
+exit:
+ return res;
+error:
+ htable_key2prim2d_purge(&scn->key2prim2d);
+ goto exit;
+}
+
+#elif DIM == 3
+static res_T
+setup_primitive_keys_3d(struct sdis_scene* scn, struct senc3d_scene* senc_scn)
+{
+ unsigned iprim = 0;
+ unsigned nprims = 0;
+ res_T res = RES_OK;
+ ASSERT(scn && senc_scn);
+
+ SENC3D(scene_get_primitives_count(senc_scn, &nprims));
+
+ FOR_EACH(iprim, 0, nprims) {
+ struct s3d_primitive prim = S3D_PRIMITIVE_NULL;
+ struct sdis_primkey key = SDIS_PRIMKEY_NULL;
+ unsigned ids[3] = {0};
+ double v0[3] = {0};
+ double v1[3] = {0};
+ double v2[3] = {0};
+
+ /* Retrieve positions from Star-Enclosre, not Star-3D. Star-Enclosure keeps
+ * the positions submitted by the user as they are, without any
+ * transformation or conversion (Star-3D converts them to float). This
+ * ensures that the caller can construct the same key from his data */
+ SENC3D(scene_get_primitive(senc_scn, iprim, ids));
+ SENC3D(scene_get_vertex(senc_scn, ids[0], v0));
+ SENC3D(scene_get_vertex(senc_scn, ids[1], v1));
+ SENC3D(scene_get_vertex(senc_scn, ids[2], v2));
+ S3D(scene_view_get_primitive(scn->s3d_view, iprim, &prim));
+
+ sdis_primkey_setup(&key, v0, v1, v2);
+
+ res = htable_key2prim3d_set(&scn->key2prim3d, &key, &prim);
+ if(res != RES_OK) goto error;
+ }
+
+exit:
+ return res;
+error:
+ htable_key2prim3d_purge(&scn->key2prim3d);
+ goto exit;
+}
+#endif
+
/* Create a Stardis scene */
static res_T
XD(scene_create)
@@ -956,8 +1039,9 @@ XD(scene_create)
scn = MEM_CALLOC(dev->allocator, 1, sizeof(struct sdis_scene));
if(!scn) {
- log_err(dev, "%s: could not allocate the Stardis scene.\n", FUNC_NAME);
res = RES_MEM_ERR;
+ log_err(dev, "%s: unabale to allocate the scene -- %s\n",
+ FUNC_NAME, res_to_cstr(res));
goto error;
}
@@ -973,6 +1057,8 @@ XD(scene_create)
darray_prim_prop_init(dev->allocator, &scn->prim_props);
htable_enclosure_init(dev->allocator, &scn->enclosures);
htable_d_init(dev->allocator, &scn->tmp_hc_ub);
+ htable_key2prim2d_init(dev->allocator, &scn->key2prim2d);
+ htable_key2prim3d_init(dev->allocator, &scn->key2prim3d);
if(args->source) {
SDIS(source_ref_get(args->source));
@@ -994,23 +1080,32 @@ XD(scene_create)
args->context,
&senc_scn);
if(res != RES_OK) {
- log_err(dev, "%s: error during the scene analysis.\n", FUNC_NAME);
+ log_err(dev, "%s: unable to analyze the scene -- %s\n",
+ FUNC_NAME, res_to_cstr(res));
goto error;
}
res = XD(setup_properties)(scn, senc_scn, args->get_interface, args->context);
if(res != RES_OK) {
- log_err(dev, "%s: could not setup the scene interfaces and their media.\n",
- FUNC_NAME);
+ log_err(dev, "%s: unable to configure interfaces and media -- %s\n",
+ FUNC_NAME, res_to_cstr(res));
goto error;
}
res = XD(setup_scene_geometry)(scn, senc_scn);
if(res != RES_OK) {
- log_err(dev, "%s: could not setup the scene geometry.\n", FUNC_NAME);
+ log_err(dev, "%s: unable to configure scene geometry -- %s\n",
+ FUNC_NAME, res_to_cstr(res));
goto error;
}
res = XD(setup_enclosures)(scn, senc_scn);
if(res != RES_OK) {
- log_err(dev, "%s: could not setup the enclosures.\n", FUNC_NAME);
+ log_err(dev, "%s: unable to configure enclosures -- %s\n",
+ FUNC_NAME, res_to_cstr(res));
+ goto error;
+ }
+ res = XD(setup_primitive_keys)(scn, senc_scn);
+ if(res != RES_OK) {
+ log_err(dev, "%s: unable to configure primitive keys -- %s\n",
+ FUNC_NAME, res_to_cstr(res));
goto error;
}
scn->sencXd(scn) = senc_scn;
diff --git a/src/sdis_scene_c.h b/src/sdis_scene_c.h
@@ -159,6 +159,16 @@ enclosure_local2global_prim_id
return darray_uint_cdata_get(&enc->local2global)[local_prim_id];
}
+static INLINE void
+primkey_init
+ (const struct mem_allocator* allocator,
+ struct sdis_primkey* key)
+{
+ ASSERT(allocator && key);
+ (void)allocator;
+ *key = SDIS_PRIMKEY_NULL;
+}
+
/* Declare the array of interfaces */
#define DARRAY_NAME interf
#define DARRAY_DATA struct sdis_interface*
@@ -193,6 +203,24 @@ enclosure_local2global_prim_id
#define HTABLE_DATA double
#include <rsys/hash_table.h>
+/* Declare the hash table that maps the primitive key to its 2D primitve */
+#define HTABLE_NAME key2prim2d
+#define HTABLE_KEY struct sdis_primkey
+#define HTABLE_KEY_FUNCTOR_INIT primkey_init
+#define HTABLE_KEY_FUNCTOR_HASH sdis_primkey_hash
+#define HTABLE_KEY_FUNCTOR_EQ sdis_primkey_eq
+#define HTABLE_DATA struct s2d_primitive
+#include <rsys/hash_table.h>
+
+/* Declare the hash table that maps the primitive key to its 3D primitive */
+#define HTABLE_NAME key2prim3d
+#define HTABLE_KEY struct sdis_primkey
+#define HTABLE_KEY_FUNCTOR_INIT primkey_init
+#define HTABLE_KEY_FUNCTOR_HASH sdis_primkey_hash
+#define HTABLE_KEY_FUNCTOR_EQ sdis_primkey_eq
+#define HTABLE_DATA struct s3d_primitive
+#include <rsys/hash_table.h>
+
struct sdis_scene {
struct darray_interf interfaces; /* List of interfaces own by the scene */
struct darray_medium media; /* List of media own by the scene */
@@ -206,6 +234,10 @@ struct sdis_scene {
struct htable_enclosure enclosures; /* Map an enclosure id to its data */
unsigned outer_enclosure_id;
+ /* Map a primivei key to its Star-2D/Star-3D primitive */
+ struct htable_key2prim2d key2prim2d;
+ struct htable_key2prim3d key2prim3d;
+
double fp_to_meter;
double tmin; /* Minimum temperature of the system (In Kelvin) */
double tmax; /* Maximum temperature of the system (In Kelvin) */
diff --git a/src/test_sdis_custom_solid_path_sampling.c b/src/test_sdis_custom_solid_path_sampling.c
@@ -0,0 +1,583 @@
+/* Copyright (C) 2016-2024 |Méso|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_mesh.h"
+#include "test_sdis_utils.h"
+
+#include <star/s3d.h>
+#include <star/s3dut.h>
+#include <star/ssp.h>
+
+#include <rsys/double3.h>
+
+/*
+ * The system is a trilinear profile of the temperature at steady state, i.e. at
+ * each point of the system we can calculate the temperature analytically. Two
+ * forms are immersed in this temperature field: a super shape and a sphere
+ * included in the super shape. On the Monte Carlo side, the temperature is
+ * unknown everywhere except on the surface of the super shape whose
+ * temperature is defined from the aformentionned trilinear profile.
+ *
+ * We will estimate the temperature at the position of a probe in solids by
+ * providing a user-side function to sample the conductive path in the sphere.
+ * We should find the temperature of the trilinear profile at the probe position
+ * by Monte Carlo, independently of this coupling with an external path sampling
+ * routine.
+ *
+ *
+ * /\ <-- T(x,y,z)
+ * ___/ \___
+ * T(z) \ __ /
+ * | T(y) T=?/. \ /
+ * |/ / \__/ \
+ * o--- T(x) /_ __ _\
+ * \/ \/
+ */
+
+#define NREALISATIONS 10000
+#define SPHERE_RADIUS 1
+
+/*******************************************************************************
+ * Helper functions
+ ******************************************************************************/
+static double
+trilinear_profile(const double pos[3])
+{
+ /* Range in X, Y and Z in which the trilinear profile is defined */
+ const double lower = -4;
+ const double upper = +4;
+
+ /* Upper temperature limit in X, Y and Z [K]. Lower temperature limit is
+ * implicitly 0 */
+ const double a = 333; /* Upper temperature limit in X [K] */
+ const double b = 432; /* Upper temperature limit in Y [K] */
+ const double c = 579; /* Upper temperature limit in Z [K] */
+
+ double x, y, z;
+
+ /* Check pre-conditions */
+ CHK(pos);
+ CHK(lower <= pos[0] && pos[0] <= upper);
+ CHK(lower <= pos[1] && pos[1] <= upper);
+ CHK(lower <= pos[2] && pos[2] <= upper);
+
+ x = (pos[0] - lower) / (upper - lower);
+ y = (pos[1] - lower) / (upper - lower);
+ z = (pos[2] - lower) / (upper - lower);
+ return a*x + b*y + c*z;
+}
+
+/*******************************************************************************
+ * Scene view
+ ******************************************************************************/
+/* Parameter for get_inidices/get_vertices functions. Although its layout is the
+ * same as that of the mesh data structure, we have deliberately defined a new
+ * data type since mesh functions cannot be invoked. This is because the form's
+ * member variables are not necessarily extensible arrays, as is the case with
+ * mesh */
+struct shape {
+ double* pos;
+ size_t* ids;
+ size_t npos;
+ size_t ntri;
+};
+#define SHAPE_NULL__ {NULL, NULL, 0, 0}
+static const struct shape SHAPE_NULL = SHAPE_NULL__;
+
+static void
+get_position(const unsigned ivert, float pos[3], void* ctx)
+{
+ const struct shape* shape = ctx;
+ CHK(shape && pos && ivert < shape->npos);
+ f3_set_d3(pos, shape->pos + ivert*3);
+}
+
+static void
+get_indices(const unsigned itri, unsigned ids[3], void* ctx)
+{
+ const struct shape* shape = ctx;
+ CHK(shape && ids && itri < shape->ntri);
+ ids[0] = (unsigned)shape->ids[itri*3+0];
+ ids[1] = (unsigned)shape->ids[itri*3+1];
+ ids[2] = (unsigned)shape->ids[itri*3+2];
+}
+
+static struct s3d_scene_view*
+create_view(struct shape* shape_data)
+{
+ /* Star3D */
+ struct s3d_vertex_data vdata = S3D_VERTEX_DATA_NULL;
+ struct s3d_device* dev = NULL;
+ struct s3d_scene* scn = NULL;
+ struct s3d_shape* shape = NULL;
+ struct s3d_scene_view* view = NULL;
+
+ OK(s3d_device_create(NULL, NULL, 0, &dev));
+
+ /* Shape */
+ vdata.usage = S3D_POSITION;
+ vdata.type = S3D_FLOAT3;
+ vdata.get = get_position;
+ OK(s3d_shape_create_mesh(dev, &shape));
+ OK(s3d_mesh_setup_indexed_vertices(shape, (unsigned)shape_data->ntri,
+ get_indices, (unsigned)shape_data->npos, &vdata, 1, shape_data));
+
+ /* Scene view */
+ OK(s3d_scene_create(dev, &scn));
+ OK(s3d_scene_attach_shape(scn, shape));
+ OK(s3d_scene_view_create(scn, S3D_TRACE|S3D_GET_PRIMITIVE, &view));
+
+ /* Clean up */
+ OK(s3d_device_ref_put(dev));
+ OK(s3d_scene_ref_put(scn));
+ OK(s3d_shape_ref_put(shape));
+
+ return view;
+}
+
+/*******************************************************************************
+ * Mesh, i.e. supershape and sphere
+ ******************************************************************************/
+static void
+mesh_add_super_shape(struct mesh* mesh)
+{
+ struct s3dut_mesh* sshape = NULL;
+ struct s3dut_mesh_data sshape_data;
+ struct s3dut_super_formula f0 = S3DUT_SUPER_FORMULA_NULL;
+ struct s3dut_super_formula f1 = S3DUT_SUPER_FORMULA_NULL;
+ const double radius = 2;
+ const unsigned nslices = 256;
+
+ f0.A = 1.5; f0.B = 1; f0.M = 11.0; f0.N0 = 1; f0.N1 = 1; f0.N2 = 2.0;
+ f1.A = 1.0; f1.B = 2; f1.M = 3.6; f1.N0 = 1; f1.N1 = 2; f1.N2 = 0.7;
+ OK(s3dut_create_super_shape(NULL, &f0, &f1, radius, nslices, nslices/2, &sshape));
+ OK(s3dut_mesh_get_data(sshape, &sshape_data));
+ mesh_append(mesh, sshape_data.positions, sshape_data.nvertices,
+ sshape_data.indices, sshape_data.nprimitives, NULL);
+ OK(s3dut_mesh_ref_put(sshape));
+}
+
+static void
+mesh_add_sphere(struct mesh* mesh)
+{
+ struct s3dut_mesh* sphere = NULL;
+ struct s3dut_mesh_data sphere_data;
+ const double radius = SPHERE_RADIUS;
+ const unsigned nslices = 128;
+
+ OK(s3dut_create_sphere(NULL, radius, nslices, nslices/2, &sphere));
+ OK(s3dut_mesh_get_data(sphere, &sphere_data));
+ mesh_append(mesh, sphere_data.positions, sphere_data.nvertices,
+ sphere_data.indices, sphere_data.nprimitives, NULL);
+ OK(s3dut_mesh_ref_put(sphere));
+}
+
+/*******************************************************************************
+ * Custom conductive path
+ ******************************************************************************/
+struct custom_solid {
+ struct s3d_scene_view* view; /* Star-3D view of the shape */
+ const struct shape* shape; /* Raw data */
+};
+
+static void
+setup_solver_primitive
+ (struct sdis_scene* scn,
+ const struct shape* shape,
+ const struct s3d_hit* user_hit,
+ struct s3d_primitive* prim)
+{
+ struct sdis_primkey key = SDIS_PRIMKEY_NULL;
+ const double *v0, *v1, *v2;
+ float v0f[3], v1f[3], v2f[3];
+ struct s3d_attrib attr0, attr1, attr2;
+
+ v0 = shape->pos + shape->ids[user_hit->prim.prim_id*3+0]*3;
+ v1 = shape->pos + shape->ids[user_hit->prim.prim_id*3+1]*3;
+ v2 = shape->pos + shape->ids[user_hit->prim.prim_id*3+2]*3;
+ sdis_primkey_setup(&key, v0, v1, v2);
+ OK(sdis_scene_get_s3d_primitive(scn, &key, prim));
+
+ /* Check that the primitive on the solver side is the same as that on the
+ * user side. On the solver side, vertices are stored in simple precision in
+ * Star-3D view. We therefore need to take care of this conversion to check
+ * that the vertices are the same */
+ OK(s3d_triangle_get_vertex_attrib(prim, 0, S3D_POSITION, &attr0));
+ OK(s3d_triangle_get_vertex_attrib(prim, 1, S3D_POSITION, &attr1));
+ OK(s3d_triangle_get_vertex_attrib(prim, 2, S3D_POSITION, &attr2));
+ f3_set_d3(v0f, v0);
+ f3_set_d3(v1f, v1);
+ f3_set_d3(v2f, v2);
+
+ /* The vertices have been inverted on the user's side to reverse the normal
+ * orientation. Below it is taken into account */
+ CHK(f3_eq(v0f, attr0.value));
+ CHK(f3_eq(v1f, attr2.value));
+ CHK(f3_eq(v2f, attr1.value));
+}
+
+/* Implementation of a Walk on Sphere algorithm for an origin-centered spherical
+ * geometry of radius SPHERE_RADIUS */
+static res_T
+sample_steady_diffusive_path
+ (struct sdis_scene* scn,
+ struct ssp_rng* rng,
+ struct sdis_path* path,
+ struct sdis_data* data)
+{
+ struct custom_solid* solid = NULL;
+ struct s3d_hit hit = S3D_HIT_NULL;
+
+ double pos[3];
+ const double epsilon = SPHERE_RADIUS * 1.e-6; /* Epsilon shell */
+
+ CHK(scn && rng && path && data);
+
+ solid = sdis_data_get(data);
+
+ d3_set(pos, path->vtx.P);
+
+ do {
+ /* Distance from the geometry center to the current position */
+ const double dst = d3_len(pos);
+
+ double dir[3] = {0,0,0};
+ double r = 0; /* Radius */
+
+ r = SPHERE_RADIUS - dst;
+ CHK(dst > 0);
+
+ if(r > epsilon) {
+ /* Uniformly sample a new position on the surrounding sphere */
+ ssp_ran_sphere_uniform(rng, dir, NULL);
+
+ /* Move to the new position */
+ d3_muld(dir, dir, r);
+ d3_add(pos, pos, dir);
+
+ /* The current position is in the epsilon shell:
+ * move it to the nearest interface position */
+ } else {
+ float posf[3];
+
+ d3_set(dir, pos);
+ d3_normalize(dir, dir);
+ d3_muld(pos, dir, SPHERE_RADIUS);
+
+ /* Map the position to the sphere geometry */
+ f3_set_d3(posf, pos);
+ OK(s3d_scene_view_closest_point(solid->view, posf, (float)INF, NULL, &hit));
+ }
+
+ /* The calculation is performed in steady state, so the path necessarily stops
+ * at a boundary */
+ } while(S3D_HIT_NONE(&hit));
+
+ /* Setup the path state */
+ d3_set(path->vtx.P, pos);
+ path->weight = 0;
+ path->at_limit = 0;
+ path->prim_2d = S2D_PRIMITIVE_NULL;
+ path->elapsed_time = 0;
+ path->weight = 0;
+ path->at_limit = 0;
+ setup_solver_primitive(scn, solid->shape, &hit, &path->prim_3d);
+
+ return RES_OK;
+}
+
+/*******************************************************************************
+ * The solids, i.e. media of the super shape and the sphere
+ ******************************************************************************/
+#define SOLID_PROP(Prop, Val) \
+ static double \
+ solid_get_##Prop \
+ (const struct sdis_rwalk_vertex* vtx, \
+ struct sdis_data* data) \
+ { \
+ (void)vtx, (void)data; /* Avoid the "unused variable" warning */ \
+ return Val; \
+ }
+SOLID_PROP(calorific_capacity, 500.0) /* [J/K/Kg] */
+SOLID_PROP(thermal_conductivity, 25.0) /* [W/m/K] */
+SOLID_PROP(volumic_mass, 7500.0) /* [kg/m^3] */
+SOLID_PROP(temperature, SDIS_TEMPERATURE_NONE/*<=> unknown*/) /* [K] */
+SOLID_PROP(delta, 1.0/40.0) /* [m] */
+
+static struct sdis_medium*
+create_solid(struct sdis_device* sdis)
+{
+ struct sdis_solid_shader shader = SDIS_SOLID_SHADER_NULL;
+ struct sdis_medium* solid = NULL;
+
+ shader.calorific_capacity = solid_get_calorific_capacity;
+ shader.thermal_conductivity = solid_get_thermal_conductivity;
+ shader.volumic_mass = solid_get_volumic_mass;
+ shader.delta = solid_get_delta;
+ shader.temperature = solid_get_temperature;
+
+ OK(sdis_solid_create(sdis, &shader, NULL, &solid));
+ return solid;
+}
+
+static struct sdis_medium*
+create_custom
+ (struct sdis_device* sdis,
+ struct s3d_scene_view* view,
+ const struct shape* shape)
+{
+ /* Stardis variables */
+ struct sdis_solid_shader shader = SDIS_SOLID_SHADER_NULL;
+ struct sdis_medium* solid = NULL;
+ struct sdis_data* data = NULL;
+
+ /* Mesh variables */
+ struct custom_solid* custom_solid = NULL;
+ const size_t sz = sizeof(struct custom_solid);
+ const size_t al = ALIGNOF(struct custom_solid);
+
+ OK(sdis_data_create(sdis, sz, al, NULL, &data));
+ custom_solid = sdis_data_get(data);
+ custom_solid->view = view;
+ custom_solid->shape = shape;
+
+ shader.calorific_capacity = solid_get_calorific_capacity;
+ shader.thermal_conductivity = solid_get_thermal_conductivity;
+ shader.volumic_mass = solid_get_volumic_mass;
+ shader.delta = solid_get_delta;
+ shader.temperature = solid_get_temperature;
+ shader.sample_path = sample_steady_diffusive_path;
+
+ OK(sdis_solid_create(sdis, &shader, data, &solid));
+ OK(sdis_data_ref_put(data));
+
+ return solid;
+}
+
+/*******************************************************************************
+ * Dummy environment, i.e. environment surrounding the super shape. It is
+ * defined only for Stardis compliance: in Stardis, an interface must divide 2
+ * media.
+ ******************************************************************************/
+static struct sdis_medium*
+create_dummy(struct sdis_device* sdis)
+{
+ struct sdis_fluid_shader shader = SDIS_FLUID_SHADER_NULL;
+ struct sdis_medium* dummy = NULL;
+
+ shader.calorific_capacity = dummy_medium_getter;
+ shader.volumic_mass = dummy_medium_getter;
+ shader.temperature = dummy_medium_getter;
+ OK(sdis_fluid_create(sdis, &shader, NULL, &dummy));
+ return dummy;
+}
+
+/*******************************************************************************
+ * Interface: its temperature is fixed to the trilinear profile
+ ******************************************************************************/
+static double
+interface_get_temperature
+ (const struct sdis_interface_fragment* frag,
+ struct sdis_data* data)
+{
+ (void)data; /* Avoid the "unused variable" warning */
+ return trilinear_profile(frag->P);
+}
+
+static struct sdis_interface*
+create_interface
+ (struct sdis_device* sdis,
+ struct sdis_medium* front,
+ struct sdis_medium* back)
+{
+ struct sdis_interface* interf = NULL;
+ struct sdis_interface_shader shader = SDIS_INTERFACE_SHADER_NULL;
+
+ /* Fluid/solid interface: fix the temperature to the temperature profile */
+ if(sdis_medium_get_type(front) != sdis_medium_get_type(back)) {
+ shader.front.temperature = interface_get_temperature;
+ shader.back.temperature = interface_get_temperature;
+ }
+
+ OK(sdis_interface_create(sdis, front, back, &shader, NULL, &interf));
+ return interf;
+}
+
+/*******************************************************************************
+ * Scene, i.e. the system to simulate
+ ******************************************************************************/
+struct scene_context {
+ const struct mesh* mesh;
+ size_t sshape_end_id; /* Last triangle index of the super shape */
+ struct sdis_interface* sshape;
+ struct sdis_interface* sphere;
+};
+#define SCENE_CONTEXT_NULL__ {NULL, 0, 0, 0}
+static const struct scene_context SCENE_CONTEXT_NULL = SCENE_CONTEXT_NULL__;
+
+static void
+scene_get_indices(const size_t itri, size_t ids[3], void* ctx)
+{
+ struct scene_context* context = ctx;
+ CHK(ids && context && itri < mesh_ntriangles(context->mesh));
+ /* Flip the indices to ensure that the normal points into the super shape */
+ ids[0] = (unsigned)context->mesh->indices[itri*3+0];
+ ids[1] = (unsigned)context->mesh->indices[itri*3+2];
+ ids[2] = (unsigned)context->mesh->indices[itri*3+1];
+}
+
+static void
+scene_get_interface(const size_t itri, struct sdis_interface** interf, void* ctx)
+{
+ struct scene_context* context = ctx;
+ CHK(interf && context && itri < mesh_ntriangles(context->mesh));
+ if(itri < context->sshape_end_id) {
+ *interf = context->sshape;
+ } else {
+ *interf = context->sphere;
+ }
+}
+
+static void
+scene_get_position(const size_t ivert, double pos[3], void* ctx)
+{
+ struct scene_context* context = ctx;
+ CHK(pos && context && ivert < mesh_nvertices(context->mesh));
+ pos[0] = context->mesh->positions[ivert*3+0];
+ pos[1] = context->mesh->positions[ivert*3+1];
+ pos[2] = context->mesh->positions[ivert*3+2];
+}
+
+static struct sdis_scene*
+create_scene(struct sdis_device* sdis, struct scene_context* ctx)
+{
+ struct sdis_scene* scn = NULL;
+ struct sdis_scene_create_args scn_args = SDIS_SCENE_CREATE_ARGS_DEFAULT;
+
+ scn_args.get_indices = scene_get_indices;
+ scn_args.get_interface = scene_get_interface;
+ scn_args.get_position = scene_get_position;
+ scn_args.nprimitives = mesh_ntriangles(ctx->mesh);
+ scn_args.nvertices = mesh_nvertices(ctx->mesh);
+ scn_args.context = ctx;
+ OK(sdis_scene_create(sdis, &scn_args, &scn));
+ return scn;
+}
+
+/*******************************************************************************
+ * Validations
+ ******************************************************************************/
+static void
+check_probe(struct sdis_scene* scn, const int is_master_process)
+{
+ struct sdis_solve_probe_args args = SDIS_SOLVE_PROBE_ARGS_DEFAULT;
+ struct sdis_mc T = SDIS_MC_NULL;
+ struct sdis_estimator* estimator = NULL;
+ double ref = 0;
+
+ args.position[0] = SPHERE_RADIUS*0.125;
+ args.position[1] = SPHERE_RADIUS*0.250;
+ args.position[2] = SPHERE_RADIUS*0.375;
+ args.nrealisations = NREALISATIONS;
+
+ OK(sdis_solve_probe(scn, &args, &estimator));
+
+ if(!is_master_process) return;
+
+ OK(sdis_estimator_get_temperature(estimator, &T));
+
+ ref = trilinear_profile(args.position);
+
+ printf("T(%g, %g, %g) = %g ~ %g +/- %g\n",
+ SPLIT3(args.position), ref, T.E, T.SE);
+
+ CHK(eq_eps(ref, T.E, T.SE*3));
+ OK(sdis_estimator_ref_put(estimator));
+}
+
+/*******************************************************************************
+ * The test
+ ******************************************************************************/
+int
+main(int argc, char** argv)
+{
+ /* Stardis */
+ struct sdis_device* dev = NULL;
+ struct sdis_interface* solid_dummy = NULL;
+ struct sdis_interface* custom_solid = NULL;
+ struct sdis_medium* solid = NULL;
+ struct sdis_medium* custom = NULL;
+ struct sdis_medium* dummy = NULL; /* Medium surrounding the solid */
+ struct sdis_scene* scn = NULL;
+
+ /* Star3D */
+ struct shape shape = SHAPE_NULL;
+ struct s3d_scene_view* sphere_view = NULL;
+
+ /* Miscellaneous */
+ struct scene_context ctx = SCENE_CONTEXT_NULL;
+ struct mesh mesh = MESH_NULL;
+ size_t sshape_end_id = 0; /* Last index of the super shape */
+ int is_master_process = 1;
+ (void)argc, (void)argv;
+
+ create_default_device(&argc, &argv, &is_master_process, &dev);
+
+ /* Mesh */
+ mesh_init(&mesh);
+ mesh_add_super_shape(&mesh);
+ sshape_end_id = mesh_ntriangles(&mesh);
+ mesh_add_sphere(&mesh);
+
+ /* Create a view of the sphere's geometry. This will be used to couple custom
+ * solid path sampling to the solver */
+ shape.pos = mesh.positions;
+ shape.ids = mesh.indices + sshape_end_id*3;
+ shape.npos = mesh_nvertices(&mesh);
+ shape.ntri = mesh_ntriangles(&mesh) - sshape_end_id/* #sshape triangles*/;
+ sphere_view = create_view(&shape);
+
+ /* Physical properties */
+ dummy = create_dummy(dev);
+ solid = create_solid(dev);
+ custom = create_custom(dev, sphere_view, &shape);
+ solid_dummy = create_interface(dev, solid, dummy);
+ custom_solid = create_interface(dev, custom, solid);
+
+ /* Scene */
+ ctx.mesh = &mesh;
+ ctx.sshape_end_id = sshape_end_id;
+ ctx.sshape = solid_dummy;
+ ctx.sphere = custom_solid;
+ scn = create_scene(dev, &ctx);
+
+ check_probe(scn, is_master_process);
+
+ mesh_release(&mesh);
+
+ OK(s3d_scene_view_ref_put(sphere_view));
+ OK(sdis_interface_ref_put(solid_dummy));
+ OK(sdis_interface_ref_put(custom_solid));
+ OK(sdis_medium_ref_put(custom));
+ OK(sdis_medium_ref_put(dummy));
+ OK(sdis_medium_ref_put(solid));
+ OK(sdis_scene_ref_put(scn));
+
+ free_default_device(dev);
+
+ CHK(mem_allocated_size() == 0);
+ return 0;
+}
diff --git a/src/test_sdis_custom_solid_path_sampling_2d.c b/src/test_sdis_custom_solid_path_sampling_2d.c
@@ -0,0 +1,603 @@
+/* Copyright (C) 2016-2024 |Méso|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_mesh.h"
+#include "test_sdis_utils.h"
+
+#include <rsys/double2.h>
+#include <rsys/float2.h>
+
+#include <star/s2d.h>
+
+/*
+ * The system is a bilinear profile of the temperature at steady state, i.e. at
+ * each point of the system we can calculate the temperature analytically. Two
+ * forms are immersed in this temperature field: a super shape and a circle
+ * included in the super shape. On the Monte Carlo side, the temperature is
+ * unknown everywhere except on the surface of the super shape whose
+ * temperature is defined from the aformentionned bilinear profile.
+ *
+ * We will estimate the temperature at the position of a probe in solids by
+ * providing a user-side function to sample the conductive path in the circle.
+ * We should find the temperature of the bilinear profile at the probe position
+ * by Monte Carlo, independently of this coupling with an external path sampling
+ * routine.
+ *
+ *
+ * /\ <-- T(x,y,z)
+ * ___/ \___
+ * T(y) \ __ /
+ * | T(y) \ / \ /
+ * |/ T=? *__/ \
+ * o--- T(x) /_ __ _\
+ * \/ \/
+ */
+
+#define NREALISATIONS 10000
+#define CIRCLE_RADIUS 1
+
+/*******************************************************************************
+ * Helper functions
+ ******************************************************************************/
+static double
+bilinear_profile(const double pos[2])
+{
+ /* Range in X, Y in which the trilinear profile is defined */
+ const double lower = -4;
+ const double upper = +4;
+
+ /* Upper temperature limit in X, Y and Z [K]. Lower temperature limit is
+ * implicitly 0 */
+ const double a = 333; /* Upper temperature limit in X [K] */
+ const double b = 432; /* Upper temperature limit in Y [K] */
+
+ double x, y;
+
+ /* Check pre-conditions */
+ CHK(pos);
+ CHK(lower <= pos[0] && pos[0] <= upper);
+ CHK(lower <= pos[1] && pos[1] <= upper);
+
+ x = (pos[0] - lower) / (upper - lower);
+ y = (pos[1] - lower) / (upper - lower);
+ return a*x + b*y;
+}
+
+/*******************************************************************************
+ * Scene view
+ ******************************************************************************/
+/* Parameter for get_inidices/get_vertices functions. Although its layout is the
+ * same as that of the mesh data structure, we have deliberately defined a new
+ * data type since mesh functions cannot be invoked. This is because the form's
+ * member variables are not necessarily extensible arrays, as is the case with
+ * mesh */
+struct shape {
+ double* pos;
+ size_t* ids;
+ size_t npos;
+ size_t nseg;
+};
+#define SHAPE_NULL__ {NULL, NULL, 0, 0}
+static const struct shape SHAPE_NULL = SHAPE_NULL__;
+
+static void
+get_position(const unsigned ivert, float pos[2], void* ctx)
+{
+ const struct shape* shape = ctx;
+ CHK(shape && pos && ivert < shape->npos);
+ f2_set_d2(pos, shape->pos + ivert*2);
+}
+
+static void
+get_indices(const unsigned iseg, unsigned ids[2], void* ctx)
+{
+ const struct shape* shape = ctx;
+ CHK(shape && ids && iseg < shape->nseg);
+ ids[0] = (unsigned)shape->ids[iseg*2+0];
+ ids[1] = (unsigned)shape->ids[iseg*2+1];
+}
+
+static struct s2d_scene_view*
+create_view(struct shape* shape_data)
+{
+ /* Star2D */
+ struct s2d_vertex_data vdata = S2D_VERTEX_DATA_NULL;
+ struct s2d_device* dev = NULL;
+ struct s2d_scene* scn = NULL;
+ struct s2d_shape* shape = NULL;
+ struct s2d_scene_view* view = NULL;
+
+ OK(s2d_device_create(NULL, NULL, 0, &dev));
+
+ /* Shape */
+ vdata.usage = S2D_POSITION;
+ vdata.type = S2D_FLOAT2;
+ vdata.get = get_position;
+ OK(s2d_shape_create_line_segments(dev, &shape));
+ OK(s2d_line_segments_setup_indexed_vertices(shape, (unsigned)shape_data->nseg,
+ get_indices, (unsigned)shape_data->npos, &vdata, 1, shape_data));
+
+ /* Scene view */
+ OK(s2d_scene_create(dev, &scn));
+ OK(s2d_scene_attach_shape(scn, shape));
+ OK(s2d_scene_view_create(scn, S2D_TRACE|S2D_GET_PRIMITIVE, &view));
+
+ /* Clean up */
+ OK(s2d_device_ref_put(dev));
+ OK(s2d_scene_ref_put(scn));
+ OK(s2d_shape_ref_put(shape));
+
+ return view;
+}
+
+/*******************************************************************************
+ * Mesh, i.e. supershape and sphere
+ ******************************************************************************/
+static void
+mesh_add_super_shape(struct mesh* mesh)
+{
+ double* pos = NULL;
+ size_t* ids = NULL;
+
+ const unsigned nslices = 128;
+ const double a = 1.0;
+ const double b = 1.0;
+ const double n1 = 1.0;
+ const double n2 = 1.0;
+ const double n3 = 1.0;
+ const double m = 6.0;
+ size_t i = 0;
+
+ CHK(mesh);
+
+ FOR_EACH(i, 0, nslices) {
+ const double theta = (double)i * (2.0*PI / (double)nslices);
+ const double tmp0 = pow(fabs(1.0/a * cos(m/4.0*theta)), n2);
+ const double tmp1 = pow(fabs(1.0/b * sin(m/4.0*theta)), n3);
+ const double tmp2 = pow(tmp0 + tmp1, 1.0/n1);
+ const double r = 1.0 / tmp2;
+ const double x = cos(theta) * r * CIRCLE_RADIUS*2;
+ const double y = sin(theta) * r * CIRCLE_RADIUS*2;
+ const size_t j = (i + 1) % nslices;
+
+ sa_push(pos, x);
+ sa_push(pos, y);
+ sa_push(ids, i);
+ sa_push(ids, j);
+ }
+
+ mesh_2d_append(mesh, pos, sa_size(pos)/2, ids, sa_size(ids)/2, NULL);
+
+ sa_release(pos);
+ sa_release(ids);
+}
+
+static void
+mesh_add_circle(struct mesh* mesh)
+{
+ double* pos = NULL;
+ size_t* ids = NULL;
+
+ const size_t nverts = 64;
+ size_t i = 0;
+
+ CHK(mesh);
+
+ FOR_EACH(i, 0, nverts) {
+ const double theta = (double)i * (2*PI)/(double)nverts;
+ const double x = cos(theta)*CIRCLE_RADIUS;
+ const double y = sin(theta)*CIRCLE_RADIUS;
+ const size_t j = (i+1)%nverts;
+
+ sa_push(pos, x);
+ sa_push(pos, y);
+ sa_push(ids, i);
+ sa_push(ids, j);
+ }
+
+ mesh_2d_append(mesh, pos, sa_size(pos)/2, ids, sa_size(ids)/2, NULL);
+
+ sa_release(pos);
+ sa_release(ids);
+}
+
+/*******************************************************************************
+ * Custom conductive path
+ ******************************************************************************/
+struct custom_solid {
+ struct s2d_scene_view* view; /* Star-2D view of the shape */
+ const struct shape* shape; /* Raw data */
+};
+
+static void
+setup_solver_primitive
+ (struct sdis_scene* scn,
+ const struct shape* shape,
+ const struct s2d_hit* user_hit,
+ struct s2d_primitive* prim)
+{
+ struct sdis_primkey key = SDIS_PRIMKEY_NULL;
+ const double *v0, *v1;
+ float v0f[2], v1f[2];
+ struct s2d_attrib attr0, attr1;
+
+ v0 = shape->pos + shape->ids[user_hit->prim.prim_id*2+0]*2;
+ v1 = shape->pos + shape->ids[user_hit->prim.prim_id*2+1]*2;
+ sdis_primkey_2d_setup(&key, v0, v1);
+ OK(sdis_scene_get_s2d_primitive(scn, &key, prim));
+
+ /* Check that the primitive on the solver side is the same as that on the
+ * user side. On the solver side, vertices are stored in simple precision in
+ * Star-3D view. We therefore need to take care of this conversion to check
+ * that the vertices are the same */
+ OK(s2d_segment_get_vertex_attrib(prim, 0, S2D_POSITION, &attr0));
+ OK(s2d_segment_get_vertex_attrib(prim, 1, S2D_POSITION, &attr1));
+ f2_set_d2(v0f, v0);
+ f2_set_d2(v1f, v1);
+
+ /* The vertices have been inverted on the user's side to reverse the normal
+ * orientation. Below it is taken into account */
+ CHK(f2_eq(v0f, attr0.value));
+ CHK(f2_eq(v1f, attr1.value));
+}
+
+/* Implementation of a Walk on Sphere algorithm for an origin-centered spherical
+ * geometry of radius SPHERE_RADIUS */
+static res_T
+sample_steady_diffusive_path
+ (struct sdis_scene* scn,
+ struct ssp_rng* rng,
+ struct sdis_path* path,
+ struct sdis_data* data)
+{
+ struct custom_solid* solid = NULL;
+ struct s2d_hit hit = S2D_HIT_NULL;
+
+ double pos[2];
+ const double epsilon = CIRCLE_RADIUS * 1.e-6; /* Epsilon shell */
+
+ CHK(scn && rng && path && data);
+
+ solid = sdis_data_get(data);
+
+ d2_set(pos, path->vtx.P);
+
+ do {
+ /* Distance from the geometry center to the current position */
+ const double dst = d2_len(pos);
+
+ double dir[3] = {0,0,0};
+ double r = 0; /* Radius */
+
+ r = CIRCLE_RADIUS - dst;
+ CHK(dst > 0);
+
+ if(r > epsilon) {
+ /* Uniformly sample a new position on the surrounding sphere */
+ ssp_ran_circle_uniform(rng, dir, NULL);
+
+ /* Move to the new position */
+ d2_muld(dir, dir, r);
+ d2_add(pos, pos, dir);
+
+ /* The current position is in the epsilon shell:
+ * move it to the nearest interface position */
+ } else {
+ float posf[2];
+
+ d2_set(dir, pos);
+ d2_normalize(dir, dir);
+ d2_muld(pos, dir, CIRCLE_RADIUS);
+
+ /* Map the position to the sphere geometry */
+ f2_set_d2(posf, pos);
+ OK(s2d_scene_view_closest_point(solid->view, posf, (float)INF, NULL, &hit));
+ }
+
+ /* The calculation is performed in steady state, so the path necessarily stops
+ * at a boundary */
+ } while(S2D_HIT_NONE(&hit));
+
+ /* Setup the path state */
+ d2_set(path->vtx.P, pos);
+ path->weight = 0;
+ path->at_limit = 0;
+ path->prim_3d = S3D_PRIMITIVE_NULL;
+ path->elapsed_time = 0;
+ path->weight = 0;
+ path->at_limit = 0;
+ setup_solver_primitive(scn, solid->shape, &hit, &path->prim_2d);
+
+ return RES_OK;
+}
+/*******************************************************************************
+ * The solids, i.e. media of the super shape and the sphere
+ ******************************************************************************/
+#define SOLID_PROP(Prop, Val) \
+ static double \
+ solid_get_##Prop \
+ (const struct sdis_rwalk_vertex* vtx, \
+ struct sdis_data* data) \
+ { \
+ (void)vtx, (void)data; /* Avoid the "unused variable" warning */ \
+ return Val; \
+ }
+SOLID_PROP(calorific_capacity, 500.0) /* [J/K/Kg] */
+SOLID_PROP(thermal_conductivity, 25.0) /* [W/m/K] */
+SOLID_PROP(volumic_mass, 7500.0) /* [kg/m^3] */
+SOLID_PROP(temperature, SDIS_TEMPERATURE_NONE/*<=> unknown*/) /* [K] */
+SOLID_PROP(delta, 1.0/40.0) /* [m] */
+
+static struct sdis_medium*
+create_solid(struct sdis_device* sdis)
+{
+ struct sdis_solid_shader shader = SDIS_SOLID_SHADER_NULL;
+ struct sdis_medium* solid = NULL;
+
+ shader.calorific_capacity = solid_get_calorific_capacity;
+ shader.thermal_conductivity = solid_get_thermal_conductivity;
+ shader.volumic_mass = solid_get_volumic_mass;
+ shader.delta = solid_get_delta;
+ shader.temperature = solid_get_temperature;
+
+ OK(sdis_solid_create(sdis, &shader, NULL, &solid));
+ return solid;
+}
+
+static struct sdis_medium*
+create_custom
+ (struct sdis_device* sdis,
+ struct s2d_scene_view* view,
+ const struct shape* shape)
+{
+ /* Stardis variables */
+ struct sdis_solid_shader shader = SDIS_SOLID_SHADER_NULL;
+ struct sdis_medium* solid = NULL;
+ struct sdis_data* data = NULL;
+
+ /* Mesh variables */
+ struct custom_solid* custom_solid = NULL;
+ const size_t sz = sizeof(struct custom_solid);
+ const size_t al = ALIGNOF(struct custom_solid);
+
+ OK(sdis_data_create(sdis, sz, al, NULL, &data));
+ custom_solid = sdis_data_get(data);
+ custom_solid->view = view;
+ custom_solid->shape = shape;
+
+ shader.calorific_capacity = solid_get_calorific_capacity;
+ shader.thermal_conductivity = solid_get_thermal_conductivity;
+ shader.volumic_mass = solid_get_volumic_mass;
+ shader.delta = solid_get_delta;
+ shader.temperature = solid_get_temperature;
+ shader.sample_path = sample_steady_diffusive_path;
+
+ OK(sdis_solid_create(sdis, &shader, data, &solid));
+ OK(sdis_data_ref_put(data));
+
+ return solid;
+}
+
+/*******************************************************************************
+ * Dummy environment, i.e. environment surrounding the super shape. It is
+ * defined only for Stardis compliance: in Stardis, an interface must divide 2
+ * media.
+ ******************************************************************************/
+static struct sdis_medium*
+create_dummy(struct sdis_device* sdis)
+{
+ struct sdis_fluid_shader shader = SDIS_FLUID_SHADER_NULL;
+ struct sdis_medium* dummy = NULL;
+
+ shader.calorific_capacity = dummy_medium_getter;
+ shader.volumic_mass = dummy_medium_getter;
+ shader.temperature = dummy_medium_getter;
+ OK(sdis_fluid_create(sdis, &shader, NULL, &dummy));
+ return dummy;
+}
+
+/*******************************************************************************
+ * Interface: its temperature is fixed to the trilinear profile
+ ******************************************************************************/
+static double
+interface_get_temperature
+ (const struct sdis_interface_fragment* frag,
+ struct sdis_data* data)
+{
+ (void)data; /* Avoid the "unused variable" warning */
+ return bilinear_profile(frag->P);
+}
+
+static struct sdis_interface*
+create_interface
+ (struct sdis_device* sdis,
+ struct sdis_medium* front,
+ struct sdis_medium* back)
+{
+ struct sdis_interface* interf = NULL;
+ struct sdis_interface_shader shader = SDIS_INTERFACE_SHADER_NULL;
+
+ /* Fluid/solid interface: fix the temperature to the temperature profile */
+ if(sdis_medium_get_type(front) != sdis_medium_get_type(back)) {
+ shader.front.temperature = interface_get_temperature;
+ shader.back.temperature = interface_get_temperature;
+ }
+
+ OK(sdis_interface_create(sdis, front, back, &shader, NULL, &interf));
+ return interf;
+}
+
+/*******************************************************************************
+ * Scene, i.e. the system to simulate
+ ******************************************************************************/
+struct scene_context {
+ const struct mesh* mesh;
+ size_t sshape_end_id; /* Last segment index of the super shape */
+ struct sdis_interface* sshape;
+ struct sdis_interface* sphere;
+};
+#define SCENE_CONTEXT_NULL__ {NULL, 0, 0, 0}
+static const struct scene_context SCENE_CONTEXT_NULL = SCENE_CONTEXT_NULL__;
+
+static void
+scene_get_indices(const size_t iseg, size_t ids[2], void* ctx)
+{
+ struct scene_context* context = ctx;
+ CHK(ids && context && iseg < mesh_2d_nsegments(context->mesh));
+ ids[0] = (unsigned)context->mesh->indices[iseg*2+0];
+ ids[1] = (unsigned)context->mesh->indices[iseg*2+1];
+}
+
+static void
+scene_get_interface(const size_t iseg, struct sdis_interface** interf, void* ctx)
+{
+ struct scene_context* context = ctx;
+ CHK(interf && context && iseg < mesh_2d_nsegments(context->mesh));
+ if(iseg < context->sshape_end_id) {
+ *interf = context->sshape;
+ } else {
+ *interf = context->sphere;
+ }
+}
+
+static void
+scene_get_position(const size_t ivert, double pos[2], void* ctx)
+{
+ struct scene_context* context = ctx;
+ CHK(pos && context && ivert < mesh_2d_nvertices(context->mesh));
+ pos[0] = context->mesh->positions[ivert*2+0];
+ pos[1] = context->mesh->positions[ivert*2+1];
+}
+
+static struct sdis_scene*
+create_scene(struct sdis_device* sdis, struct scene_context* ctx)
+{
+ struct sdis_scene* scn = NULL;
+ struct sdis_scene_create_args scn_args = SDIS_SCENE_CREATE_ARGS_DEFAULT;
+
+ scn_args.get_indices = scene_get_indices;
+ scn_args.get_interface = scene_get_interface;
+ scn_args.get_position = scene_get_position;
+ scn_args.nprimitives = mesh_2d_nsegments(ctx->mesh);
+ scn_args.nvertices = mesh_2d_nvertices(ctx->mesh);
+ scn_args.context = ctx;
+ OK(sdis_scene_2d_create(sdis, &scn_args, &scn));
+ return scn;
+}
+
+/*******************************************************************************
+ * Validations
+ ******************************************************************************/
+static void
+check_probe(struct sdis_scene* scn, const int is_master_process)
+{
+ struct sdis_solve_probe_args args = SDIS_SOLVE_PROBE_ARGS_DEFAULT;
+ struct sdis_mc T = SDIS_MC_NULL;
+ struct sdis_estimator* estimator = NULL;
+ double ref = 0;
+
+ args.position[0] = CIRCLE_RADIUS*0.125;
+ args.position[1] = CIRCLE_RADIUS*0.250;
+ args.nrealisations = NREALISATIONS;
+
+ OK(sdis_solve_probe(scn, &args, &estimator));
+
+ if(!is_master_process) return;
+
+ OK(sdis_estimator_get_temperature(estimator, &T));
+
+ ref = bilinear_profile(args.position);
+
+ printf("T(%g, %g) = %g ~ %g +/- %g\n",
+ SPLIT2(args.position), ref, T.E, T.SE);
+
+ CHK(eq_eps(ref, T.E, T.SE*3));
+ OK(sdis_estimator_ref_put(estimator));
+}
+
+/*******************************************************************************
+ * The test
+ ******************************************************************************/
+int
+main(int argc, char** argv)
+{
+ /* Stardis */
+ struct sdis_device* dev = NULL;
+ struct sdis_interface* solid_dummy = NULL;
+ struct sdis_interface* custom_solid = NULL;
+ struct sdis_medium* solid = NULL;
+ struct sdis_medium* custom = NULL;
+ struct sdis_medium* dummy = NULL; /* Medium surrounding the solid */
+ struct sdis_scene* scn = NULL;
+
+ /* Star 2D */
+ struct shape shape = SHAPE_NULL;
+ struct s2d_scene_view* circle_view = NULL;
+
+ /* Miscellaneous */
+ struct scene_context ctx = SCENE_CONTEXT_NULL;
+ struct mesh mesh = MESH_NULL;
+ size_t sshape_end_id = 0; /* Last index of the super shape */
+ int is_master_process = 1;
+ (void)argc, (void)argv;
+
+ create_default_device(&argc, &argv, &is_master_process, &dev);
+
+ /* Mesh */
+ mesh_init(&mesh);
+ mesh_add_super_shape(&mesh);
+ sshape_end_id = mesh_2d_nsegments(&mesh);
+ mesh_add_circle(&mesh);
+
+ /* Create a view of the circle's geometry. This will be used to couple custom
+ * solid path sampling to the solver */
+ shape.pos = mesh.positions;
+ shape.ids = mesh.indices + sshape_end_id*2;
+ shape.npos = mesh_2d_nvertices(&mesh);
+ shape.nseg = mesh_2d_nsegments(&mesh) - sshape_end_id;
+ circle_view = create_view(&shape);
+
+ /* Physical properties */
+ dummy = create_dummy(dev);
+ solid = create_solid(dev);
+ custom = create_custom(dev, circle_view, &shape);
+ solid_dummy = create_interface(dev, dummy, solid);
+ custom_solid = create_interface(dev, solid, custom);
+
+ /* Scene */
+ ctx.mesh = &mesh;
+ ctx.sshape_end_id = sshape_end_id;
+ ctx.sshape = solid_dummy;
+ ctx.sphere = custom_solid;
+ scn = create_scene(dev, &ctx);
+
+ check_probe(scn, is_master_process);
+
+ mesh_release(&mesh);
+
+ OK(s2d_scene_view_ref_put(circle_view));
+ OK(sdis_interface_ref_put(solid_dummy));
+ OK(sdis_interface_ref_put(custom_solid));
+ OK(sdis_medium_ref_put(custom));
+ OK(sdis_medium_ref_put(dummy));
+ OK(sdis_medium_ref_put(solid));
+ OK(sdis_scene_ref_put(scn));
+
+ free_default_device(dev);
+
+ CHK(mem_allocated_size() == 0);
+ return 0;
+}
diff --git a/src/test_sdis_mesh.h b/src/test_sdis_mesh.h
@@ -49,6 +49,13 @@ mesh_nvertices(const struct mesh* mesh)
return sa_size(mesh->positions) / 3/* #coords per vertex */;
}
+static INLINE size_t
+mesh_2d_nvertices(const struct mesh* mesh)
+{
+ CHK(mesh);
+ return sa_size(mesh->positions) / 2/* #coords per vertex */;
+}
+
/* Number of triangles */
static INLINE size_t
mesh_ntriangles(const struct mesh* mesh)
@@ -57,6 +64,13 @@ mesh_ntriangles(const struct mesh* mesh)
return sa_size(mesh->indices) / 3/* #indices per triangle */;
}
+static INLINE size_t
+mesh_2d_nsegments(const struct mesh* mesh)
+{
+ CHK(mesh);
+ return sa_size(mesh->indices) / 2/* #indices per segment */;
+}
+
static INLINE void
mesh_append
(struct mesh* mesh,
@@ -97,6 +111,42 @@ mesh_append
}
static INLINE void
+mesh_2d_append
+ (struct mesh* mesh,
+ const double* in_positions,
+ const size_t in_nvertices,
+ const size_t* in_indices,
+ const size_t in_nsegments,
+ const double in_translate[2]) /* May be NULL */
+{
+ double translate[2] = {0, 0};
+ double* positions = NULL;
+ size_t* indices = NULL;
+ size_t ivert = 0;
+ size_t i = 0;
+ CHK(mesh != NULL);
+
+ ivert = mesh_2d_nvertices(mesh);
+ positions = sa_add(mesh->positions, in_nvertices*2);
+ indices = sa_add(mesh->indices, in_nsegments*2);
+
+ if(in_translate) {
+ translate[0] = in_translate[0];
+ translate[1] = in_translate[1];
+ }
+
+ FOR_EACH(i, 0, in_nvertices) {
+ positions[i*2 + 0] = in_positions[i*2 + 0] + translate[0];
+ positions[i*2 + 1] = in_positions[i*2 + 1] + translate[1];
+ }
+
+ FOR_EACH(i, 0, in_nsegments) {
+ indices[i*2 + 0] = in_indices[i*2 + 0] + ivert;
+ indices[i*2 + 1] = in_indices[i*2 + 1] + ivert;
+ }
+}
+
+static INLINE void
mesh_dump(const struct mesh* mesh, FILE* stream)
{
size_t i, n;
@@ -117,4 +167,24 @@ mesh_dump(const struct mesh* mesh, FILE* stream)
fflush(stream);
}
+static INLINE void
+mesh_2d_dump(const struct mesh* mesh, FILE* stream)
+{
+ size_t i, n;
+ CHK(mesh != NULL);
+
+ n = mesh_2d_nvertices(mesh);
+ FOR_EACH(i, 0, n) {
+ fprintf(stream, "v %g %g\n", SPLIT2(mesh->positions+i*2));
+ }
+
+ n = mesh_2d_nsegments(mesh);
+ FOR_EACH(i, 0, n) {
+ fprintf(stream, "l %lu %lu\n",
+ (unsigned long)(mesh->indices[i*2+0] + 1),
+ (unsigned long)(mesh->indices[i*2+1] + 1));
+ }
+ fflush(stream);
+}
+
#endif /* TEST_SDIS_MESH_H */
diff --git a/src/test_sdis_primkey.c b/src/test_sdis_primkey.c
@@ -0,0 +1,281 @@
+/* Copyright (C) 2016-2024 |Méso|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 <star/s3dut.h>
+#include <rsys/float3.h>
+
+/*
+ * The system is a solid supershape whose boundary temperature is set to a
+ * constant. The temperature of the solid is therefore this same temperature.
+ * This simplistic test case is not used to verify a Monte Carlo estimate, but
+ * to ensure that the caller can recover the internal representation of the
+ * geometric primitives from his own data.
+ *
+ * /\
+ * ___/ \___
+ * \ /
+ * /_ __ _\
+ * \/ \/
+ *
+ */
+
+/*******************************************************************************
+ * Super shape
+ ******************************************************************************/
+static struct s3dut_mesh*
+create_super_shape(void)
+{
+ struct s3dut_mesh* mesh = NULL;
+ struct s3dut_super_formula f0 = S3DUT_SUPER_FORMULA_NULL;
+ struct s3dut_super_formula f1 = S3DUT_SUPER_FORMULA_NULL;
+ const double radius = 1;
+ const unsigned nslices = 256;
+
+ f0.A = 1.5; f0.B = 1; f0.M = 11.0; f0.N0 = 1; f0.N1 = 1; f0.N2 = 2.0;
+ f1.A = 1.0; f1.B = 2; f1.M = 3.6; f1.N0 = 1; f1.N1 = 2; f1.N2 = 0.7;
+ OK(s3dut_create_super_shape(NULL, &f0, &f1, radius, nslices, nslices/2, &mesh));
+
+ return mesh;
+}
+
+/*******************************************************************************
+ * Solid, i.e. medium of the super shape
+ ******************************************************************************/
+#define SOLID_PROP(Prop, Val) \
+ static double \
+ solid_get_##Prop \
+ (const struct sdis_rwalk_vertex* vtx, \
+ struct sdis_data* data) \
+ { \
+ (void)vtx, (void)data; /* Avoid the "unused variable" warning */ \
+ return Val; \
+ }
+SOLID_PROP(calorific_capacity, 1) /* [J/K/Kg] */
+SOLID_PROP(thermal_conductivity, 1) /* [W/m/K] */
+SOLID_PROP(volumic_mass, 1) /* [kg/m^3] */
+SOLID_PROP(temperature, SDIS_TEMPERATURE_NONE) /* [K] */
+SOLID_PROP(delta, 1.0/20.0) /* [m/fp_to_meter] */
+#undef SOLID_PROP
+
+static struct sdis_medium*
+create_solid(struct sdis_device* sdis)
+{
+ struct sdis_solid_shader shader = SDIS_SOLID_SHADER_NULL;
+ struct sdis_medium* solid = NULL;
+
+ shader.calorific_capacity = solid_get_calorific_capacity;
+ shader.thermal_conductivity = solid_get_thermal_conductivity;
+ shader.volumic_mass = solid_get_volumic_mass;
+ shader.delta = solid_get_delta;
+ shader.temperature = solid_get_temperature;
+ OK(sdis_solid_create(sdis, &shader, NULL, &solid));
+ return solid;
+}
+
+/*******************************************************************************
+ * Dummy environment, i.e. environment surrounding the super shape. It is
+ * defined only for Stardis compliance: in Stardis, an interface must divide 2
+ * media.
+ ******************************************************************************/
+static struct sdis_medium*
+create_dummy(struct sdis_device* sdis)
+{
+ struct sdis_fluid_shader shader = SDIS_FLUID_SHADER_NULL;
+ struct sdis_medium* dummy = NULL;
+
+ shader.calorific_capacity = dummy_medium_getter;
+ shader.volumic_mass = dummy_medium_getter;
+ shader.temperature = dummy_medium_getter;
+ OK(sdis_fluid_create(sdis, &shader, NULL, &dummy));
+ return dummy;
+}
+
+/*******************************************************************************
+ * Interface
+ ******************************************************************************/
+static double
+interface_get_temperature
+ (const struct sdis_interface_fragment* frag,
+ struct sdis_data* data)
+{
+ (void)frag, (void)data; /* Avoid the "unused variable" warning */
+ return 300; /* [K] */
+}
+
+static struct sdis_interface*
+create_interface
+ (struct sdis_device* sdis,
+ struct sdis_medium* front,
+ struct sdis_medium* back)
+{
+ struct sdis_interface* interf = NULL;
+ struct sdis_interface_shader shader = SDIS_INTERFACE_SHADER_NULL;
+
+ shader.front.temperature = interface_get_temperature;
+ shader.back.temperature = interface_get_temperature;
+ OK(sdis_interface_create(sdis, front, back, &shader, NULL, &interf));
+ return interf;
+}
+
+/*******************************************************************************
+ * Scene, i.e. the system to simulate
+ ******************************************************************************/
+struct scene_context {
+ struct s3dut_mesh_data mesh_data;
+ struct sdis_interface* interf;
+};
+static const struct scene_context SCENE_CONTEXT_NULL = {{0}, NULL};
+
+static void
+scene_get_indices(const size_t itri, size_t ids[3], void* ctx)
+{
+ struct scene_context* context = ctx;
+ CHK(ids && context && itri < context->mesh_data.nprimitives);
+ /* Flip the indices to ensure that the normal points into the super shape */
+ ids[0] = (unsigned)context->mesh_data.indices[itri*3+0];
+ ids[1] = (unsigned)context->mesh_data.indices[itri*3+2];
+ ids[2] = (unsigned)context->mesh_data.indices[itri*3+1];
+}
+
+static void
+scene_get_interface(const size_t itri, struct sdis_interface** interf, void* ctx)
+{
+ struct scene_context* context = ctx;
+ CHK(interf && context && itri < context->mesh_data.nprimitives);
+ *interf = context->interf;
+}
+
+static void
+scene_get_position(const size_t ivert, double pos[3], void* ctx)
+{
+ struct scene_context* context = ctx;
+ CHK(pos && context && ivert < context->mesh_data.nvertices);
+ pos[0] = context->mesh_data.positions[ivert*3+0];
+ pos[1] = context->mesh_data.positions[ivert*3+1];
+ pos[2] = context->mesh_data.positions[ivert*3+2];
+}
+
+static struct sdis_scene*
+create_scene
+ (struct sdis_device* sdis,
+ const struct s3dut_mesh* mesh,
+ struct sdis_interface* interf)
+{
+ struct sdis_scene* scn = NULL;
+ struct sdis_scene_create_args scn_args = SDIS_SCENE_CREATE_ARGS_DEFAULT;
+ struct scene_context context = SCENE_CONTEXT_NULL;
+
+ OK(s3dut_mesh_get_data(mesh, &context.mesh_data));
+ context.interf = interf;
+
+ scn_args.get_indices = scene_get_indices;
+ scn_args.get_interface = scene_get_interface;
+ scn_args.get_position = scene_get_position;
+ scn_args.nprimitives = context.mesh_data.nprimitives;
+ scn_args.nvertices = context.mesh_data.nvertices;
+ scn_args.context = &context;
+ OK(sdis_scene_create(sdis, &scn_args, &scn));
+ return scn;
+}
+
+/*******************************************************************************
+ * Validation
+ ******************************************************************************/
+static void
+check(struct sdis_scene* scn, const struct s3dut_mesh* mesh)
+{
+ struct s3d_primitive prim = S3D_PRIMITIVE_NULL;
+ struct s3dut_mesh_data mesh_data;
+ struct sdis_primkey key = SDIS_PRIMKEY_NULL;
+ size_t iprim = 0;
+
+ OK(s3dut_mesh_get_data(mesh, &mesh_data));
+
+ BA(sdis_scene_get_s3d_primitive(NULL, &key, &prim));
+ BA(sdis_scene_get_s3d_primitive(scn, NULL, &prim));
+ BA(sdis_scene_get_s3d_primitive(scn, &key, NULL));
+ BA(sdis_scene_get_s3d_primitive(scn, &key, &prim));
+
+ FOR_EACH(iprim, 0, mesh_data.nprimitives) {
+ const double *v0, *v1, *v2;
+ float v0f[3], v1f[3], v2f[3];
+ struct s3d_attrib attr0, attr1, attr2;
+
+ /* Check that a primitive can be obtained from the key constructed on the
+ * user side */
+ v0 = mesh_data.positions + mesh_data.indices[iprim*3 + 0]*3;
+ v1 = mesh_data.positions + mesh_data.indices[iprim*3 + 1]*3;
+ v2 = mesh_data.positions + mesh_data.indices[iprim*3 + 2]*3;
+ sdis_primkey_setup(&key, v0, v1, v2);
+ OK(sdis_scene_get_s3d_primitive(scn, &key, &prim));
+
+ /* Check that the primitive on the solver side is the same as that on the
+ * user side. On the solver side, vertices are stored in simple precision in
+ * Star-3D view. We therefore need to take care of this conversion to check
+ * that the vertices are the same */
+ OK(s3d_triangle_get_vertex_attrib(&prim, 0, S3D_POSITION, &attr0));
+ OK(s3d_triangle_get_vertex_attrib(&prim, 1, S3D_POSITION, &attr1));
+ OK(s3d_triangle_get_vertex_attrib(&prim, 2, S3D_POSITION, &attr2));
+ f3_set_d3(v0f, v0);
+ f3_set_d3(v1f, v1);
+ f3_set_d3(v2f, v2);
+
+ /* The vertices have been inverted on the user's side to reverse the normal
+ * orientation. Below it is taken into account */
+ CHK(f3_eq(v0f, attr0.value));
+ CHK(f3_eq(v1f, attr2.value));
+ CHK(f3_eq(v2f, attr1.value));
+ }
+}
+
+/*******************************************************************************
+ * The test
+ ******************************************************************************/
+int
+main(int argc, char** argv)
+{
+ /* Stardis variables */
+ struct sdis_device* sdis = NULL;
+ struct sdis_interface* interf = NULL;
+ struct sdis_medium* solid = NULL;
+ struct sdis_medium* dummy = NULL;
+ struct sdis_scene* scn = NULL;
+
+ struct s3dut_mesh* super_shape = NULL;
+ (void)argc, (void)argv;
+
+ OK(sdis_device_create(&SDIS_DEVICE_CREATE_ARGS_DEFAULT, &sdis));
+
+ super_shape = create_super_shape();
+ solid = create_solid(sdis);
+ dummy = create_dummy(sdis);
+ interf = create_interface(sdis, solid, dummy);
+ scn = create_scene(sdis, super_shape, interf);
+
+ check(scn, super_shape);
+
+ OK(s3dut_mesh_ref_put(super_shape));
+ OK(sdis_device_ref_put(sdis));
+ OK(sdis_interface_ref_put(interf));
+ OK(sdis_medium_ref_put(solid));
+ OK(sdis_medium_ref_put(dummy));
+ OK(sdis_scene_ref_put(scn));
+
+ CHK(mem_allocated_size() == 0);
+ return 0;
+}
diff --git a/src/test_sdis_primkey_2d.c b/src/test_sdis_primkey_2d.c
@@ -0,0 +1,295 @@
+/* Copyright (C) 2016-2024 |Méso|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_mesh.h"
+
+#include <star/s2d.h>
+
+#include <rsys/float2.h>
+
+/*
+ * The system is a solid supershape whose boundary temperature is set to a
+ * constant. The temperature of the solid is therefore this same temperature.
+ * This simplistic test case is not used to verify a Monte Carlo estimate, but
+ * to ensure that the caller can recover the internal representation of the
+ * geometric primitives from his own data.
+ *
+ * /\
+ * ___/ \___
+ * \ /
+ * /_ __ _\
+ * \/ \/
+ *
+ */
+
+/*******************************************************************************
+ * Super shape
+ ******************************************************************************/
+static struct mesh
+create_super_shape(void)
+{
+ struct mesh sshape = MESH_NULL;
+
+ const unsigned nslices = 128;
+ const double a = 1.0;
+ const double b = 1.0;
+ const double n1 = 1.0;
+ const double n2 = 1.0;
+ const double n3 = 1.0;
+ const double m = 6.0;
+ size_t i = 0;
+
+ FOR_EACH(i, 0, nslices) {
+ const double theta = (double)i * (2.0*PI / (double)nslices);
+ const double tmp0 = pow(fabs(1.0/a * cos(m/4.0*theta)), n2);
+ const double tmp1 = pow(fabs(1.0/b * sin(m/4.0*theta)), n3);
+ const double tmp2 = pow(tmp0 + tmp1, 1.0/n1);
+ const double r = 1.0 / tmp2;
+ const double x = cos(theta) * r;
+ const double y = sin(theta) * r;
+ const size_t j = (i + 1) % nslices;
+
+ sa_push(sshape.positions, x);
+ sa_push(sshape.positions, y);
+ sa_push(sshape.indices, i);
+ sa_push(sshape.indices, j);
+ }
+
+ return sshape;
+}
+
+/*******************************************************************************
+ * Solid, i.e. medium of the super shape
+ ******************************************************************************/
+#define SOLID_PROP(Prop, Val) \
+ static double \
+ solid_get_##Prop \
+ (const struct sdis_rwalk_vertex* vtx, \
+ struct sdis_data* data) \
+ { \
+ (void)vtx, (void)data; /* Avoid the "unused variable" warning */ \
+ return Val; \
+ }
+SOLID_PROP(calorific_capacity, 1) /* [J/K/Kg] */
+SOLID_PROP(thermal_conductivity, 1) /* [W/m/K] */
+SOLID_PROP(volumic_mass, 1) /* [kg/m^3] */
+SOLID_PROP(temperature, SDIS_TEMPERATURE_NONE) /* [K] */
+SOLID_PROP(delta, 1.0/20.0) /* [m/fp_to_meter] */
+#undef SOLID_PROP
+
+static struct sdis_medium*
+create_solid(struct sdis_device* sdis)
+{
+ struct sdis_solid_shader shader = SDIS_SOLID_SHADER_NULL;
+ struct sdis_medium* solid = NULL;
+
+ shader.calorific_capacity = solid_get_calorific_capacity;
+ shader.thermal_conductivity = solid_get_thermal_conductivity;
+ shader.volumic_mass = solid_get_volumic_mass;
+ shader.delta = solid_get_delta;
+ shader.temperature = solid_get_temperature;
+ OK(sdis_solid_create(sdis, &shader, NULL, &solid));
+ return solid;
+}
+
+/*******************************************************************************
+ * Dummy environment, i.e. environment surrounding the super shape. It is
+ * defined only for Stardis compliance: in Stardis, an interface must divide 2
+ * media.
+ ******************************************************************************/
+static struct sdis_medium*
+create_dummy(struct sdis_device* sdis)
+{
+ struct sdis_fluid_shader shader = SDIS_FLUID_SHADER_NULL;
+ struct sdis_medium* dummy = NULL;
+
+ shader.calorific_capacity = dummy_medium_getter;
+ shader.volumic_mass = dummy_medium_getter;
+ shader.temperature = dummy_medium_getter;
+ OK(sdis_fluid_create(sdis, &shader, NULL, &dummy));
+ return dummy;
+}
+
+/*******************************************************************************
+ * Interface
+ ******************************************************************************/
+static double
+interface_get_temperature
+ (const struct sdis_interface_fragment* frag,
+ struct sdis_data* data)
+{
+ (void)frag, (void)data; /* Avoid the "unused variable" warning */
+ return 300; /* [K] */
+}
+
+static struct sdis_interface*
+create_interface
+ (struct sdis_device* sdis,
+ struct sdis_medium* front,
+ struct sdis_medium* back)
+{
+ struct sdis_interface* interf = NULL;
+ struct sdis_interface_shader shader = SDIS_INTERFACE_SHADER_NULL;
+
+ shader.front.temperature = interface_get_temperature;
+ shader.back.temperature = interface_get_temperature;
+ OK(sdis_interface_create(sdis, front, back, &shader, NULL, &interf));
+ return interf;
+}
+
+/*******************************************************************************
+ * Scene, i.e. the system to simulate
+ ******************************************************************************/
+struct scene_context {
+ const struct mesh* sshape;
+ struct sdis_interface* interf;
+};
+static const struct scene_context SCENE_CONTEXT_NULL = {NULL, NULL};
+
+static void
+scene_get_indices(const size_t iseg, size_t ids[2], void* ctx)
+{
+ struct scene_context* context = ctx;
+ CHK(ids && context);
+ /* Flip the indices to ensure that the normal points into the super shape */
+ ids[0] = (unsigned)context->sshape->indices[iseg*2+1];
+ ids[1] = (unsigned)context->sshape->indices[iseg*2+0];
+}
+
+static void
+scene_get_interface(const size_t iseg, struct sdis_interface** interf, void* ctx)
+{
+ struct scene_context* context = ctx;
+ (void)iseg;
+ CHK(interf && context);
+ *interf = context->interf;
+}
+
+static void
+scene_get_position(const size_t ivert, double pos[2], void* ctx)
+{
+ struct scene_context* context = ctx;
+ CHK(pos && context);
+ pos[0] = context->sshape->positions[ivert*2+0];
+ pos[1] = context->sshape->positions[ivert*2+1];
+}
+
+static struct sdis_scene*
+create_scene
+ (struct sdis_device* sdis,
+ const struct mesh* sshape,
+ struct sdis_interface* interf)
+{
+ struct sdis_scene* scn = NULL;
+ struct sdis_scene_create_args scn_args = SDIS_SCENE_CREATE_ARGS_DEFAULT;
+ struct scene_context context = SCENE_CONTEXT_NULL;
+
+ context.interf = interf;
+ context.sshape = sshape;
+
+ scn_args.get_indices = scene_get_indices;
+ scn_args.get_interface = scene_get_interface;
+ scn_args.get_position = scene_get_position;
+ scn_args.nprimitives = mesh_2d_nsegments(sshape);
+ scn_args.nvertices = mesh_2d_nvertices(sshape);
+ scn_args.context = &context;
+ OK(sdis_scene_2d_create(sdis, &scn_args, &scn));
+ return scn;
+}
+
+/*******************************************************************************
+ * Validation
+ ******************************************************************************/
+static void
+check(struct sdis_scene* scn, const struct mesh* mesh)
+{
+ struct s2d_primitive prim = S2D_PRIMITIVE_NULL;
+ struct sdis_primkey key = SDIS_PRIMKEY_NULL;
+ size_t iprim = 0;
+ size_t nprims = 0;
+
+ BA(sdis_scene_get_s2d_primitive(NULL, &key, &prim));
+ BA(sdis_scene_get_s2d_primitive(scn, NULL, &prim));
+ BA(sdis_scene_get_s2d_primitive(scn, &key, NULL));
+ BA(sdis_scene_get_s2d_primitive(scn, &key, &prim));
+
+ nprims = mesh_2d_nsegments(mesh);
+ FOR_EACH(iprim, 0, nprims) {
+ const double *v0, *v1;
+ float v0f[2], v1f[2];
+ struct s2d_attrib attr0, attr1;
+
+ /* Check that a primitive can be obtained from the key constructed on the
+ * user side */
+ v0 = mesh->positions + mesh->indices[iprim*2 + 0]*2;
+ v1 = mesh->positions + mesh->indices[iprim*2 + 1]*2;
+ sdis_primkey_2d_setup(&key, v0, v1);
+ OK(sdis_scene_get_s2d_primitive(scn, &key, &prim));
+
+ /* Check that the primitive on the solver side is the same as that on the
+ * user side. On the solver side, vertices are stored in simple precision in
+ * Star-3D view. We therefore need to take care of this conversion to check
+ * that the vertices are the same */
+ OK(s2d_segment_get_vertex_attrib(&prim, 0, S2D_POSITION, &attr0));
+ OK(s2d_segment_get_vertex_attrib(&prim, 1, S2D_POSITION, &attr1));
+ f2_set_d2(v0f, v0);
+ f2_set_d2(v1f, v1);
+
+ /* The vertices have been inverted on the user's side to reverse the normal
+ * orientation. Below it is taken into account */
+ CHK(f2_eq(v0f, attr1.value));
+ CHK(f2_eq(v1f, attr0.value));
+ }
+}
+
+/*******************************************************************************
+ * The test
+ ******************************************************************************/
+int
+main(int argc, char** argv)
+{
+ /* Stardis */
+ struct sdis_device* sdis = NULL;
+ struct sdis_interface* interf = NULL;
+ struct sdis_medium* solid = NULL;
+ struct sdis_medium* dummy = NULL; /* Medium surrounding the solid */
+ struct sdis_scene* scn = NULL;
+
+ /* Miscellaneous */
+ struct mesh sshape = MESH_NULL;
+ (void)argc, (void)argv;
+
+ OK(sdis_device_create(&SDIS_DEVICE_CREATE_ARGS_DEFAULT, &sdis));
+
+ sshape = create_super_shape();
+ solid = create_solid(sdis);
+ dummy = create_dummy(sdis);
+ interf = create_interface(sdis, solid, dummy);
+ scn = create_scene(sdis, &sshape, interf);
+
+ check(scn, &sshape);
+
+ mesh_release(&sshape);
+ OK(sdis_device_ref_put(sdis));
+ OK(sdis_interface_ref_put(interf));
+ OK(sdis_medium_ref_put(solid));
+ OK(sdis_medium_ref_put(dummy));
+ OK(sdis_scene_ref_put(scn));
+
+ CHK(mem_allocated_size() == 0);
+ return 0;
+}
diff --git a/src/test_sdis_scene.c b/src/test_sdis_scene.c
@@ -100,6 +100,8 @@ test_scene_3d
struct context ctx;
struct senc2d_scene* scn2d;
struct senc3d_scene* scn3d;
+ struct s2d_scene_view* view2d;
+ struct s3d_scene_view* view3d;
struct sdis_scene_create_args scn_args = SDIS_SCENE_CREATE_ARGS_DEFAULT;
struct sdis_scene_find_closest_point_args closest_pt_args =
SDIS_SCENE_FIND_CLOSEST_POINT_ARGS_NULL;
@@ -266,9 +268,13 @@ test_scene_3d
BA(sdis_scene_get_senc3d_scene(scn, NULL));
BA(sdis_scene_get_senc3d_scene(NULL, &scn3d));
OK(sdis_scene_get_senc3d_scene(scn, &scn3d));
- OK(senc3d_scene_ref_put(scn3d));
- /* No 2D available */
- BA(sdis_scene_get_senc2d_scene(scn, &scn2d));
+ BA(sdis_scene_get_senc2d_scene(scn, &scn2d)); /* No 2D available */
+
+ BA(sdis_scene_get_s3d_scene_view(NULL, NULL));
+ BA(sdis_scene_get_s3d_scene_view(NULL, &view3d));
+ BA(sdis_scene_get_s3d_scene_view(scn, NULL));
+ OK(sdis_scene_get_s3d_scene_view(scn, &view3d));
+ BA(sdis_scene_get_s2d_scene_view(scn, &view2d)); /* No 2D available */
BA(sdis_scene_get_radiative_env(NULL, &radenv));
BA(sdis_scene_get_radiative_env(scn, NULL));
@@ -311,6 +317,8 @@ test_scene_2d
struct context ctx;
struct senc2d_scene* scn2d;
struct senc3d_scene* scn3d;
+ struct s2d_scene_view* view2d;
+ struct s3d_scene_view* view3d;
size_t nsegs, npos;
size_t i;
size_t iprim;
@@ -502,9 +510,13 @@ test_scene_2d
BA(sdis_scene_get_senc2d_scene(scn, NULL));
BA(sdis_scene_get_senc2d_scene(NULL, &scn2d));
OK(sdis_scene_get_senc2d_scene(scn, &scn2d));
- OK(senc2d_scene_ref_put(scn2d));
- /* No 3D available */
- BA(sdis_scene_get_senc3d_scene(scn, &scn3d));
+ BA(sdis_scene_get_senc3d_scene(scn, &scn3d)); /* No 3D available */
+
+ BA(sdis_scene_get_s2d_scene_view(NULL, NULL));
+ BA(sdis_scene_get_s2d_scene_view(NULL, &view2d));
+ BA(sdis_scene_get_s2d_scene_view(scn, NULL));
+ OK(sdis_scene_get_s2d_scene_view(scn, &view2d));
+ BA(sdis_scene_get_s3d_scene_view(scn, &view3d)); /* No 3D available */
BA(sdis_scene_get_radiative_env(NULL, NULL));
BA(sdis_scene_get_radiative_env(scn, NULL));
diff --git a/src/test_sdis_solve_probe_boundary_list.c b/src/test_sdis_solve_probe_boundary_list.c
@@ -33,6 +33,7 @@
* estimate the temperature at the sphere boundary at several probe points. We
* should find by Monte Carlo the temperature of the trilinear profile at the
* position of the probe. It's the test.
+ *
* /\ <-- T(x,y,z)
* ___/ \___
* T(z) \ __ /
diff --git a/src/test_sdis_utils.h b/src/test_sdis_utils.h
@@ -195,6 +195,7 @@ static const struct sdis_solid_shader DUMMY_SOLID_SHADER = {
dummy_medium_getter, /* Delta */
dummy_medium_getter, /* Volumic power */
dummy_medium_getter, /* Temperature */
+ NULL, /* sample path */
0 /* Initial time */
};
