commit 743a809ffd4b2b60b877131342ce7930f6939fa8
parent 0af29fc052ff4d520686da8e7dd91dae4bddd07a
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Fri, 22 Dec 2023 09:15:02 +0100
Add the sdis_source API
A source is an external source whose contribution will be calculated by
Stardis at solid/fluid interfaces in the form of external net fluxes.
This data structure will replace the interface functors currently used
to manage external sources: managing them in this way is a shame, given
the complexity left to the caller to describe just one source.
At present, only spherical sources are proposed. But the reference
counter functions indicate that future source types will share the same
API. Sampling the spherical source and testing the intersection of a ray
with it are handled by two internal functions.
Note that although the code compiles, nothing has yet been tested.
Diffstat:
| M | Makefile | | | 1 | + |
| M | src/sdis.h | | | 76 | +++++++++++++++++++++++++++++++++++++++++++++++++++++++++------------------- |
| A | src/sdis_source.c | | | 274 | +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
| A | src/sdis_source_c.h | | | 55 | +++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
4 files changed, 387 insertions(+), 19 deletions(-)
diff --git a/Makefile b/Makefile
@@ -46,6 +46,7 @@ SRC =\
src/sdis_scene.c\
src/sdis_solve.c\
src/sdis_solve_camera.c\
+ src/sdis_source.c\
src/sdis_tile.c\
$($(DISTRIB_PARALLELISM)_SRC)
OBJ = $(SRC:.c=.o)
diff --git a/src/sdis.h b/src/sdis.h
@@ -69,6 +69,7 @@ struct sdis_green_function;
struct sdis_interface;
struct sdis_medium;
struct sdis_scene;
+struct sdis_source;
/* Forward declaration of non ref counted types */
struct sdis_green_path;
@@ -113,25 +114,6 @@ struct sdis_interface_fragment {
static const struct sdis_interface_fragment SDIS_INTERFACE_FRAGMENT_NULL =
SDIS_INTERFACE_FRAGMENT_NULL__;
-/*******************************************************************************
- * Estimation data types
- ******************************************************************************/
-enum sdis_estimator_type {
- SDIS_ESTIMATOR_TEMPERATURE, /* In Kelvin */
- SDIS_ESTIMATOR_FLUX, /* In Watt/m^2 */
- SDIS_ESTIMATOR_POWER, /* In Watt */
- SDIS_ESTIMATOR_TYPES_COUNT__
-};
-
-/* Monte-Carlo estimation */
-struct sdis_mc {
- double E; /* Expected value */
- double V; /* Variance */
- double SE; /* Standard error */
-};
-#define SDIS_MC_NULL__ {0, 0, 0}
-static const struct sdis_mc SDIS_MC_NULL = SDIS_MC_NULL__;
-
/* Input arguments of the sdis_device_create function */
struct sdis_device_create_args {
struct logger* logger; /* NULL <=> default logger */
@@ -158,6 +140,44 @@ struct sdis_info {
#define SDIS_INFO_NULL__ {0}
static const struct sdis_info SDIS_INFO_NULL = SDIS_INFO_NULL__;
+/* Type of functor used to retrieve the source's position relative to time. */
+typedef void
+(*sdis_get_position_T)
+ (const double time,
+ double pos[3],
+ struct sdis_data* data);
+
+/* Input arguments of the sdis_spherical_source_create function */
+struct sdis_spherical_source_create_args {
+ sdis_get_position_T position; /* [m] */
+ struct sdis_data* data; /* Data sent to the position functor */
+ double radius; /* [m] */
+ double power; /* Total power [W] */
+};
+#define SDIS_SPHERICAL_SOURCE_CREATE_ARGS_NULL__ {NULL, NULL, 0, 0}
+static const struct sdis_spherical_source_create_args
+SDIS_SPHERICAL_SOURCE_CREATE_ARGS_NULL =
+ SDIS_SPHERICAL_SOURCE_CREATE_ARGS_NULL__;
+
+/*******************************************************************************
+ * Estimation data types
+ ******************************************************************************/
+enum sdis_estimator_type {
+ SDIS_ESTIMATOR_TEMPERATURE, /* In Kelvin */
+ SDIS_ESTIMATOR_FLUX, /* In Watt/m^2 */
+ SDIS_ESTIMATOR_POWER, /* In Watt */
+ SDIS_ESTIMATOR_TYPES_COUNT__
+};
+
+/* Monte-Carlo estimation */
+struct sdis_mc {
+ double E; /* Expected value */
+ double V; /* Variance */
+ double SE; /* Standard error */
+};
+#define SDIS_MC_NULL__ {0, 0, 0}
+static const struct sdis_mc SDIS_MC_NULL = SDIS_MC_NULL__;
+
/*******************************************************************************
* Data type used to describe physical properties
******************************************************************************/
@@ -971,6 +991,24 @@ sdis_interface_get_id
(const struct sdis_interface* interf);
/*******************************************************************************
+ * External source API. When a scene has external sources, an external flux
+ * (in both its direct and diffuse parts) is imposed on the interfaces.
+ ******************************************************************************/
+SDIS_API res_T
+sdis_spherical_source_create
+ (struct sdis_device* dev,
+ struct sdis_spherical_source_create_args* args,
+ struct sdis_source** source);
+
+SDIS_API res_T
+sdis_source_ref_get
+ (struct sdis_source* source);
+
+SDIS_API res_T
+sdis_source_ref_put
+ (struct sdis_source* source);
+
+/*******************************************************************************
* A scene is a collection of primitives. Each primitive is the geometric
* support of the interface between 2 media.
******************************************************************************/
diff --git a/src/sdis_source.c b/src/sdis_source.c
@@ -0,0 +1,274 @@
+/* Copyright (C) 2016-2023 |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_device_c.h"
+#include "sdis_log.h"
+#include "sdis_source_c.h"
+
+#include <rsys/mem_allocator.h>
+#include <rsys/ref_count.h>
+
+#include <star/ssp.h>
+
+struct sdis_source {
+ struct sdis_spherical_source_create_args spherical;
+
+ struct sdis_device* dev;
+ ref_T ref;
+};
+
+/*******************************************************************************
+ * Helper functions
+ ******************************************************************************/
+static res_T
+check_spherical_source_create_args
+ (struct sdis_device* dev,
+ const char* func_name,
+ struct sdis_spherical_source_create_args* args)
+{
+ ASSERT(func_name);
+ if(!args) return RES_BAD_ARG;
+
+ if(!args->position) {
+ log_err(dev, "%s: the position functor is missing.\n", func_name);
+ return RES_BAD_ARG;
+ }
+
+ if(args->radius < 0) {
+ log_err(dev, "%s: invalid source radius '%g' m. It cannot be negative.\n",
+ func_name, args->radius);
+ return RES_BAD_ARG;
+ }
+
+ if(args->power < 0) {
+ log_err(dev, "%s: invalid source power '%g' W. It cannot be negative.\n",
+ func_name, args->power);
+ return RES_BAD_ARG;
+ }
+
+ return RES_OK;
+}
+
+static void
+release_source(ref_T* ref)
+{
+ struct sdis_device* dev = NULL;
+ struct sdis_source* src = CONTAINER_OF(ref, struct sdis_source, ref);
+ ASSERT(ref);
+ dev = src->dev;
+ SDIS(data_ref_put(src->spherical.data));
+ MEM_RM(dev->allocator, src);
+ SDIS(device_ref_put(dev));
+}
+
+/*******************************************************************************
+ * Exported symbols
+ ******************************************************************************/
+res_T
+sdis_spherical_source_create
+ (struct sdis_device* dev,
+ struct sdis_spherical_source_create_args* args,
+ struct sdis_source** out_src)
+{
+ struct sdis_source* src = NULL;
+ res_T res = RES_OK;
+
+ if(!dev || !out_src) { res = RES_BAD_ARG; goto error; }
+ res = check_spherical_source_create_args(dev, FUNC_NAME, args);
+ if(res != RES_OK) goto error;
+
+ src = MEM_CALLOC(dev->allocator, 1, sizeof(*src));
+ if(!src) {
+ log_err(dev, "%s: cannot allocate spherical source.\n", FUNC_NAME);
+ res = RES_OK;
+ goto error;
+ }
+ ref_init(&src->ref);
+ SDIS(device_ref_get(dev));
+ SDIS(data_ref_get(args->data));
+ src->spherical = *args;
+ src->dev = dev;
+
+exit:
+ if(out_src) *out_src = src;
+ return res;
+error:
+ if(src) { SDIS(source_ref_put(src)); src = NULL; }
+ goto exit;
+}
+
+res_T
+sdis_source_ref_get(struct sdis_source* src)
+{
+ if(!src) return RES_BAD_ARG;
+ ref_get(&src->ref);
+ return RES_OK;
+}
+
+res_T
+sdis_source_ref_put(struct sdis_source* src)
+{
+ if(!src) return RES_BAD_ARG;
+ ref_put(&src->ref, release_source);
+ return RES_OK;
+}
+
+/*******************************************************************************
+ * Local functions
+ ******************************************************************************/
+res_T
+source_sample
+ (const struct sdis_source* src,
+ struct ssp_rng* rng,
+ const double pos[3],
+ const double time,
+ struct source_sample* sample)
+{
+ double src_pos[3]; /* [m] */
+ double main_dir[3];
+ double half_angle; /* [radians] */
+ double cos_half_angle; /* [radians] */
+ double dst; /* [m] */
+ double radius; /* Source radius [m] */
+ double area; /* Source area [m^2] */
+ res_T res = RES_OK;
+ ASSERT(src && rng && pos && sample);
+
+ /* Retrieve current source position and radius */
+ src->spherical.position(time, src_pos, src->spherical.data);
+ radius = src->spherical.radius;
+
+ area = 4*PI*radius*radius; /* [m^2] */
+
+ /* compute the direction of `pos' toward the center of the source */
+ d3_sub(main_dir, src_pos, pos);
+
+ /* Normalize the direction and keep the distance from `pos' to the center of
+ * the source */
+ dst = d3_normalize(main_dir, main_dir);
+ if(dst <= radius) {
+ log_err(src->dev,
+ "%s: the position from which the external source is sampled "
+ "is included in the source:\n"
+ "\tsource position = %g, %g, %g\n"
+ "\tsource radius = %g\n"
+ "\tposition = %g, %g, %g\n"
+ "\ttime = %g\n"
+ "\tdistance from position to source = %g\n",
+ FUNC_NAME, SPLIT3(src_pos), radius, SPLIT3(pos), time, dst);
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ /* Point source */
+ if(area == 0) {
+ d3_set(sample->dir, main_dir);
+ sample->pdf = 1;
+ sample->dst = dst;
+ sample->radiance = src->spherical.power / (4*PI*dst*dst); /* [W/m^2/sr] */
+
+ /* Spherical source */
+ } else {
+ /* Sample the source according to its solid angle,
+ * i.e. 2*PI*(1 - cos(half_angle)) */
+ half_angle = asin(radius/dst);
+ cos_half_angle = cos(half_angle);
+ ssp_ran_sphere_cap_uniform /* pdf = 1/(2*PI*(1-cos(half_angle))) */
+ (rng, main_dir, cos_half_angle, sample->dir, &sample->pdf);
+
+
+ /* Set other sample variables */
+ sample->dst = dst - radius; /* From pos to source boundaries [m] */
+ sample->radiance = src->spherical.power / (PI*area); /* [W/m^2/sr] */
+ }
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
+res_T
+source_trace_to
+ (const struct sdis_source* src,
+ const double pos[3], /* Ray origin */
+ const double dir[3], /* Ray direction */
+ const double time, /* Time at which ray is traced */
+ struct source_sample* sample)
+{
+ double src_pos[3]; /* [m] */
+ double main_dir[3];
+ double radius; /* [m] */
+ double dst; /* Distance from pos to the source center [m] */
+ double half_angle; /* [radian] */
+ res_T res = RES_OK;
+ ASSERT(src && pos && dir && sample);
+ ASSERT(d3_is_normalized(dir));
+
+ radius = src->spherical.radius;
+
+ /* Point sources cannot be targeted */
+ if(radius == 0) {
+ *sample = SOURCE_SAMPLE_NULL;
+ goto exit;
+ }
+
+ /* Retrieve current source position and radius */
+ src->spherical.position(time, src_pos, src->spherical.data);
+
+ /* compute the direction of `pos' toward the center of the source */
+ d3_sub(main_dir, src_pos, pos);
+
+ /* Normalize the direction and keep the distance from `pos' to the center of
+ * the source */
+ dst = d3_normalize(main_dir, main_dir);
+ if(dst <= radius) {
+ log_err(src->dev,
+ "%s: the position from which the external source is targeted "
+ "is included in the source:\n"
+ "\tsource position = %g, %g, %g\n"
+ "\tsource radius = %g\n"
+ "\tposition = %g, %g, %g\n"
+ "\ttime = %g\n"
+ "\tdistance from position to source = %g\n",
+ FUNC_NAME, SPLIT3(src_pos), radius, SPLIT3(pos), time, dst);
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ /* Compute the the half angle of the source as seen from pos */
+ half_angle = asin(radius/dst);
+
+ /* The source is missed */
+ if(d3_dot(dir, main_dir) < cos(half_angle)) {
+ *sample = SOURCE_SAMPLE_NULL;
+
+ /* The source is intersected */
+ } else {
+ const double area = 4*PI*radius*radius; /* [m^2] */
+
+ d3_set(sample->dir, dir);
+ sample->pdf = 1;
+ sample->dst = dst - radius; /* From pos to source boundaries [m] */
+ sample->radiance = src->spherical.power / (PI*area); /* [W/m^2/sr] */
+ }
+
+exit:
+ return res;
+error:
+ *sample = SOURCE_SAMPLE_NULL;
+ goto exit;
+}
diff --git a/src/sdis_source_c.h b/src/sdis_source_c.h
@@ -0,0 +1,55 @@
+/* Copyright (C) 2016-2023 |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/>. */
+
+#ifndef SDIS_SOURCE_C_H
+#define SDIS_SOURCE_C_H
+
+#include <rsys/rsys.h>
+
+struct sdis_source;
+struct ssp_rng;
+
+struct source_sample {
+ double dir[3]; /* Direction _to_ the source */
+ double pdf; /* pdf of sampled direction */
+ double dst; /* Distance to the source [m] */
+ double radiance; /* [W/m^2/sr] */
+};
+#define SOURCE_SAMPLE_NULL__ {{0,0,0}, 0, 0, 0}
+static const struct source_sample SOURCE_SAMPLE_NULL = SOURCE_SAMPLE_NULL__;
+
+/* Helper macro used to define whether a sample is valid or not */
+#define SOURCE_SAMPLE_NONE(Sample) ((Sample)->pdf == 0)
+
+extern LOCAL_SYM res_T
+source_sample
+ (const struct sdis_source* source,
+ struct ssp_rng* rng,
+ const double pos[3], /* Position from which the source is sampled */
+ const double time, /* Time at which the source is sampled */
+ struct source_sample* sample);
+
+/* Trace a ray toward the source. The returned sample has a pdf of 1 or 0
+ * whether the source is intersected by the ray or not, respectively. You can
+ * use the SOURCE_SAMPLE_NONE macro to check this */
+extern LOCAL_SYM res_T
+source_trace_to
+ (const struct sdis_source* source,
+ const double pos[3], /* Ray origin */
+ const double dir[3], /* Ray direction */
+ const double time, /* Time at which ray is traced */
+ struct source_sample* sample); /* pdf == 0 if no source is reached */
+
+#endif /* SDIS_SOURCE_C_H */
