stardis-solver

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

sdis_source.c (10031B)


      1 /* Copyright (C) 2016-2025 |Méso|Star> (contact@meso-star.com)
      2  *
      3  * This program is free software: you can redistribute it and/or modify
      4  * it under the terms of the GNU General Public License as published by
      5  * the Free Software Foundation, either version 3 of the License, or
      6  * (at your option) any later version.
      7  *
      8  * This program is distributed in the hope that it will be useful,
      9  * but WITHOUT ANY WARRANTY; without even the implied warranty of
     10  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
     11  * GNU General Public License for more details.
     12  *
     13  * You should have received a copy of the GNU General Public License
     14  * along with this program. If not, see <http://www.gnu.org/licenses/>. */
     15 
     16 #include "sdis.h"
     17 #include "sdis_device_c.h"
     18 #include "sdis_log.h"
     19 #include "sdis_source_c.h"
     20 
     21 #include <rsys/free_list.h>
     22 #include <rsys/mem_allocator.h>
     23 #include <rsys/ref_count.h>
     24 
     25 #include <star/ssp.h>
     26 
     27 struct sdis_source {
     28   struct sdis_spherical_source_shader spherical;
     29   struct sdis_data* data;
     30 
     31   struct fid id; /* Unique identifier of the source */
     32   struct sdis_device* dev;
     33   ref_T ref;
     34 };
     35 
     36 /*******************************************************************************
     37  * Helper functions
     38  ******************************************************************************/
     39 static res_T
     40 check_spherical_source_shader
     41   (struct sdis_device* dev,
     42    const char* func_name,
     43    const struct sdis_spherical_source_shader* shader)
     44 {
     45   ASSERT(func_name);
     46   if(!shader) return RES_BAD_ARG;
     47 
     48   if(!shader->position) {
     49     log_err(dev, "%s: the position functor is missing.\n", func_name);
     50     return RES_BAD_ARG;
     51   }
     52 
     53   if(!shader->power) {
     54     log_err(dev, "%s: the power functor is missing.\n", func_name);
     55     return RES_BAD_ARG;
     56   }
     57 
     58   if(shader->radius < 0) {
     59     log_err(dev, "%s: invalid source radius '%g' m. It cannot be negative.\n",
     60       func_name, shader->radius);
     61     return RES_BAD_ARG;
     62   }
     63 
     64   return RES_OK;
     65 }
     66 
     67 static void
     68 release_source(ref_T* ref)
     69 {
     70   struct sdis_device* dev = NULL;
     71   struct sdis_source* src = CONTAINER_OF(ref, struct sdis_source, ref);
     72   ASSERT(ref);
     73   dev = src->dev;
     74   if(src->data) SDIS(data_ref_put(src->data));
     75   flist_name_del(&dev->source_names, src->id);
     76   MEM_RM(dev->allocator, src);
     77   SDIS(device_ref_put(dev));
     78 }
     79 
     80 /*******************************************************************************
     81  * Exported symbols
     82  ******************************************************************************/
     83 res_T
     84 sdis_spherical_source_create
     85   (struct sdis_device* dev,
     86    const struct sdis_spherical_source_shader* shader,
     87    struct sdis_data* data,
     88    struct sdis_source** out_src)
     89 {
     90   struct sdis_source* src = NULL;
     91   res_T res = RES_OK;
     92 
     93   if(!dev || !out_src) { res = RES_BAD_ARG; goto error; }
     94   res = check_spherical_source_shader(dev, FUNC_NAME, shader);
     95   if(res != RES_OK) goto error;
     96 
     97   src = MEM_CALLOC(dev->allocator, 1, sizeof(*src));
     98   if(!src) {
     99     log_err(dev, "%s: cannot allocate spherical source.\n", FUNC_NAME);
    100     res = RES_OK;
    101     goto error;
    102   }
    103   ref_init(&src->ref);
    104   SDIS(device_ref_get(dev));
    105   if(data) SDIS(data_ref_get(data));
    106   src->spherical = *shader;
    107   src->data = data;
    108   src->dev = dev;
    109   src->id = flist_name_add(&dev->source_names);
    110   flist_name_get(&dev->source_names, src->id)->mem = src;
    111 
    112 exit:
    113   if(out_src) *out_src = src;
    114   return res;
    115 error:
    116   if(src) { SDIS(source_ref_put(src)); src = NULL; }
    117   goto exit;
    118 }
    119 
    120 res_T
    121 sdis_spherical_source_get_shader
    122   (const struct sdis_source* source,
    123    struct sdis_spherical_source_shader* shader)
    124 {
    125   if(!source || !shader) return RES_BAD_ARG;
    126   *shader = source->spherical;
    127   return RES_OK;
    128 }
    129 
    130 res_T
    131 sdis_source_ref_get(struct sdis_source* src)
    132 {
    133   if(!src) return RES_BAD_ARG;
    134   ref_get(&src->ref);
    135   return RES_OK;
    136 }
    137 
    138 res_T
    139 sdis_source_ref_put(struct sdis_source* src)
    140 {
    141   if(!src) return RES_BAD_ARG;
    142   ref_put(&src->ref, release_source);
    143   return RES_OK;
    144 }
    145 
    146 struct sdis_data*
    147 sdis_source_get_data(struct sdis_source* src)
    148 {
    149   ASSERT(src);
    150   return src->data;
    151 }
    152 
    153 unsigned
    154 sdis_source_get_id(const struct sdis_source* source)
    155 {
    156   ASSERT(source);
    157   return source->id.index;
    158 }
    159 
    160 /*******************************************************************************
    161  * Local functions
    162  ******************************************************************************/
    163 res_T
    164 source_sample
    165   (const struct sdis_source* src,
    166    const struct source_props* props,
    167    struct ssp_rng* rng,
    168    const double pos[3],
    169    struct source_sample* sample)
    170 {
    171   double main_dir[3];
    172   double half_angle; /* [radians] */
    173   double cos_half_angle; /* [radians] */
    174   double dst; /* [m] */
    175   res_T res = RES_OK;
    176   ASSERT(src && rng && pos && sample);
    177 
    178   /* compute the direction of `pos' toward the center of the source */
    179   d3_sub(main_dir, props->pos, pos);
    180 
    181   /* Normalize the direction and keep the distance from `pos' to the center of
    182    * the source */
    183   dst = d3_normalize(main_dir, main_dir);
    184   if(dst <= props->radius) {
    185     log_err(src->dev,
    186       "%s: the position from which the external source is sampled "
    187       "is included in the source:\n"
    188       "\tsource position = %g, %g, %g\n"
    189       "\tsource radius = %g\n"
    190       "\tposition = %g, %g, %g\n"
    191       "\ttime = %g\n"
    192       "\tdistance from position to source = %g\n",
    193       FUNC_NAME, SPLIT3(props->pos), props->radius, SPLIT3(pos), props->time, dst);
    194     res = RES_BAD_ARG;
    195     goto error;
    196   }
    197 
    198   /* Point source */
    199   if(props->area == 0) {
    200     d3_set(sample->dir, main_dir);
    201     sample->pdf = 1;
    202     sample->dst = dst;
    203     sample->radiance_term = 1.0 / (4*PI*dst*dst); /* [W/m^2/sr] */
    204     sample->radiance = props->power * sample->radiance_term; /* [W/m^2/sr] */
    205 
    206   /* Spherical source */
    207   } else {
    208     /* Sample the source according to its solid angle,
    209      * i.e. 2*PI*(1 - cos(half_angle)) */
    210     half_angle = asin(props->radius/dst);
    211     cos_half_angle = cos(half_angle);
    212     ssp_ran_sphere_cap_uniform /* pdf = 1/(2*PI*(1-cos(half_angle))) */
    213       (rng, main_dir, cos_half_angle, sample->dir, &sample->pdf);
    214 
    215     /* Set other sample variables */
    216     sample->dst = dst - props->radius; /* From pos to source boundaries [m] */
    217     sample->radiance_term = 1.0 / (PI*props->area); /* [W/m^2/sr] */
    218     sample->radiance = props->power * sample->radiance_term; /* [W/m^2/sr] */
    219   }
    220 
    221 exit:
    222   return res;
    223 error:
    224   goto exit;
    225 }
    226 
    227 res_T
    228 source_trace_to
    229   (const struct sdis_source* src,
    230    const struct source_props* props,
    231    const double pos[3], /* Ray origin */
    232    const double dir[3], /* Ray direction */
    233    struct source_sample* sample)
    234 {
    235   double main_dir[3];
    236   double dst; /* Distance from pos to the source center [m] */
    237   double half_angle; /* [radian] */
    238   res_T res = RES_OK;
    239   ASSERT(src && props && pos && dir && sample);
    240   ASSERT(d3_is_normalized(dir));
    241 
    242   /* Point sources cannot be targeted */
    243   if(props->radius == 0) {
    244     *sample = SOURCE_SAMPLE_NULL;
    245     goto exit;
    246   }
    247 
    248   /* compute the direction of `pos' toward the center of the source */
    249   d3_sub(main_dir, props->pos, pos);
    250 
    251   /* Normalize the direction and keep the distance from `pos' to the center of
    252    * the source */
    253   dst = d3_normalize(main_dir, main_dir);
    254   if(dst <= props->radius) {
    255     log_err(src->dev,
    256       "%s: the position from which the external source is targeted "
    257       "is included in the source:\n"
    258       "\tsource position = %g, %g, %g\n"
    259       "\tsource radius = %g\n"
    260       "\tposition = %g, %g, %g\n"
    261       "\ttime = %g\n"
    262       "\tdistance from position to source = %g\n",
    263       FUNC_NAME, SPLIT3(props->pos), props->radius, SPLIT3(pos), props->time, dst);
    264     res = RES_BAD_ARG;
    265     goto error;
    266   }
    267 
    268   /* Compute the half angle of the source as seen from pos */
    269   half_angle = asin(props->radius/dst);
    270 
    271   /* The source is missed */
    272   if(d3_dot(dir, main_dir) < cos(half_angle)) {
    273     *sample = SOURCE_SAMPLE_NULL;
    274 
    275   /* The source is intersected */
    276   } else {
    277     d3_set(sample->dir, dir);
    278     sample->pdf = 1;
    279     sample->dst = dst - props->radius; /* From pos to source boundaries [m] */
    280     sample->radiance_term = 1.0 / (PI*props->area); /* [W/m^2/sr] */
    281     sample->radiance = props->power * sample->radiance_term; /* [W/m^2/sr] */
    282   }
    283 
    284 exit:
    285   return res;
    286 error:
    287   *sample = SOURCE_SAMPLE_NULL;
    288   goto exit;
    289 }
    290 
    291 res_T
    292 source_get_props
    293   (const struct sdis_source* src,
    294    const double time, /* [s] */
    295    struct source_props* props)
    296 {
    297   res_T res = RES_OK;
    298   ASSERT(src && props);
    299 
    300   /* Retrieve the source properties */
    301   src->spherical.position(time, props->pos, src->data);
    302   props->power = src->spherical.power(time, src->data);
    303   props->radius = src->spherical.radius;
    304 
    305   if(props->power < 0) {
    306     log_err(src->dev, "%s: invalid source power '%g' W. It cannot be negative.\n",
    307       FUNC_NAME, props->power);
    308     res = RES_BAD_ARG;
    309     goto error;
    310   }
    311 
    312   props->area = 4*PI*props->radius*props->radius; /* [m^2] */
    313   props->time = time; /* [s] */
    314 
    315 exit:
    316   return res;
    317 error:
    318   goto exit;
    319 }
    320 
    321 double /* [W] */
    322 source_get_power(const struct sdis_source* src, const double time /* [s] */)
    323 {
    324   ASSERT(src);
    325   return src->spherical.power(time, src->data);
    326 }
    327 
    328 double /* [W/perpendicular m^2/sr] */
    329 source_get_diffuse_radiance
    330   (const struct sdis_source* src,
    331    const double time /* [s] */,
    332    const double dir[3])
    333 {
    334   ASSERT(src);
    335   if(src->spherical.diffuse_radiance == NULL) {
    336     return 0;
    337   } else {
    338     return src->spherical.diffuse_radiance(time, dir, src->data);
    339   }
    340 }
    341 
    342 void
    343 source_compute_signature(const struct sdis_source* src, hash256_T hash)
    344 {
    345   struct sha256_ctx ctx;
    346   ASSERT(src && hash);
    347 
    348   /* Calculate the source signature. Currently, it is only the source radius.
    349    * But the Source API is designed to be independent of source type. In the
    350    * future, the source will not necessarily be spherical, so the data to be
    351    * hashed will depend on the type of source. This function anticipate this by
    352    * calculating a hash even if it is currently dispensable. */
    353   sha256_ctx_init(&ctx);
    354   sha256_ctx_update
    355     (&ctx, (const char*)&src->spherical.radius, sizeof(src->spherical.radius));
    356   sha256_ctx_finalize(&ctx, hash);
    357 }