star-gf

Compute Gebhart factors
git clone git://git.meso-star.fr/star-gf.git
Log | Files | Refs | README | LICENSE

test_sgf_square.c (6398B)


      1 /* Copyright (C) 2021, 2024 |Meso|Star> (contact@meso-star.com)
      2  * Copyright (C) 2015-2018 EDF S.A., France (syrthes-support@edf.fr)
      3  *
      4  * This program is free software: you can redistribute it and/or modify
      5  * it under the terms of the GNU General Public License as published by
      6  * the Free Software Foundation, either version 3 of the License, or
      7  * (at your option) any later version.
      8  *
      9  * This program is distributed in the hope that it will be useful,
     10  * but WITHOUT ANY WARRANTY; without even the implied warranty of
     11  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
     12  * GNU General Public License for more details.
     13  *
     14  * You should have received a copy of the GNU General Public License
     15  * along with this program. If not, see <http://www.gnu.org/licenses/>. */
     16 
     17 #include "sgf.h"
     18 #include "test_sgf_utils.h"
     19 
     20 #include <rsys/logger.h>
     21 #include <star/s2d.h>
     22 #include <star/ssp.h>
     23 
     24 #include <limits.h>
     25 
     26 #define NSTEPS 100000
     27 
     28 static const float vertices[] = {
     29   1.f, 0.f,
     30   0.f, 0.f,
     31   0.f, 1.f,
     32   1.f, 1.f
     33 };
     34 const unsigned nverts = sizeof(vertices)/(2*sizeof(float));
     35 
     36 const unsigned indices[] = {
     37   0, 1, /* Bottom */
     38   1, 2, /* Left */
     39   2, 3, /* Top */
     40   3, 0  /* Right */
     41 };
     42 const unsigned nsegs = sizeof(indices)/(2*sizeof(unsigned));
     43 
     44 static const double emissivity[] = {
     45   1.0, /* Bottom */
     46   1.0, /* Left */
     47   1.0, /* Top */
     48   1.0  /* Right */
     49 };
     50 
     51 /* Emissivity used to simulate 2 infinite planes */
     52 static const double emissivity_inf_bottom_top[] = {
     53   1, /* Bottom */
     54   0, /* Left */
     55   1, /* Top */
     56   0 /* Right */
     57 };
     58 
     59 static const double specularity[] = {
     60   0.0, /* Bottom */
     61   0.0, /* Left */
     62   0.0, /* Top */
     63   0.0  /* Right */
     64 };
     65 
     66 /* Specularity used to simulate 2 infinite segments */
     67 static const double specularity_inf_bottom_top[] = {
     68   0.0, /* Bottom */
     69   1.0, /* Left */
     70   0.0, /* Top */
     71   1.0  /* Right */
     72 };
     73 
     74 /* Check the estimation of the bottom/top & bottom/medium Gebhart factors */
     75 static void
     76 check_bottom_top_medium_gf
     77   (struct sgf_scene* scn,
     78    struct ssp_rng* rng,
     79    const double gf_bottom_top, /* Ref of the bottom/top Gebhart factor */
     80    const double gf_bottom_medium) /* Ref of the bottom/medium Gebhart Factor */
     81 {
     82   /* Indices of the top/bottom primitives */
     83   const size_t TOP = 2;
     84   const size_t BOTTOM = 0;
     85 
     86   struct sgf_estimator* estimator;
     87   struct sgf_status status[2];
     88 
     89   CHK(sgf_scene_begin_integration(scn) == RES_OK);
     90 
     91   /* Estimate the Gebhart factors for the bottom segment */
     92   CHK(sgf_integrate(scn, BOTTOM, rng, NSTEPS, &estimator) == RES_OK);
     93   CHK(sgf_estimator_get_status(estimator, TOP, 0, status + 0) == RES_OK);
     94   CHK(sgf_estimator_get_status_medium(estimator, 0, status + 1) == RES_OK);
     95   CHK(sgf_estimator_ref_put(estimator) == RES_OK);
     96 
     97   /* Check the Gebhart factor between the bottom and the top plane.  */
     98   CHK(eq_eps(status[0].E, gf_bottom_top, status[0].SE) == 1);
     99 
    100   /* Check the Gebhart factor between the bottom plane and the medium */
    101   CHK(eq_eps(status[1].E, gf_bottom_medium, status[1].SE) == 1);
    102 
    103   CHK(sgf_scene_end_integration(scn) == RES_OK);
    104 }
    105 
    106 int
    107 main(int argc, char** argv)
    108 {
    109   struct mem_allocator allocator;
    110   struct shape_context shape;
    111   struct sgf_device* sgf;
    112   struct sgf_scene* scn;
    113   struct sgf_scene_desc desc = SGF_SCENE_DESC_NULL;
    114   struct sgf_status status;
    115   struct sgf_estimator* estimator;
    116   struct ssp_rng* rng;
    117   double ka[4];
    118   size_t iprim, i;
    119   (void)argc, (void)argv;
    120 
    121   mem_init_proxy_allocator(&allocator, &mem_default_allocator);
    122 
    123   CHK(ssp_rng_create(&allocator, SSP_RNG_THREEFRY, &rng) == RES_OK);
    124   CHK(sgf_device_create(NULL, &allocator, 1, &sgf) == RES_OK);
    125   CHK(sgf_scene_create(sgf, &scn) == RES_OK);
    126 
    127   shape.vertices = vertices;
    128   shape.nvertices = nverts;
    129   shape.indices = indices;
    130   shape.nprimitives = nsegs;
    131   shape.emissivity = emissivity;
    132   shape.specularity = specularity;
    133   shape.dim = 2;
    134 
    135   desc.get_position = get_position;
    136   desc.get_indices = get_indices;
    137   desc.get_emissivity = get_emissivity;
    138   desc.get_specularity = get_specularity;
    139   desc.get_reflectivity = get_reflectivity;
    140   desc.nprims = (unsigned)nsegs;
    141   desc.nverts = (unsigned)nverts;
    142   desc.nbands = 1;
    143   desc.context = &shape;
    144 
    145   CHK(sgf_scene_setup_2d(scn, &desc) == RES_OK);
    146   CHK(sgf_integrate(scn, 0, rng, NSTEPS, &estimator) == RES_BAD_OP);
    147 
    148   CHK(sgf_scene_begin_integration(scn) == RES_OK);
    149 
    150   FOR_EACH(iprim, 0, 3) {
    151     CHK(sgf_integrate(scn, iprim, rng, NSTEPS, &estimator) == RES_OK);
    152     FOR_EACH(i, 0, 3) {
    153       CHK(sgf_estimator_get_status(estimator, i, 0, &status) == RES_OK);
    154 
    155       if(i == iprim) {
    156         CHK(eq_eps(status.E, 0, status.SE) == 1);
    157       } else if(i == iprim + 2 || i == iprim - 2) { /* parallel faces */
    158         CHK(eq_eps(status.E, sqrt(2) - 1, 2*status.SE) == 1);
    159       } else { /* Orthogonal faces */
    160         CHK(iprim == i+1 || iprim == i-1);
    161         CHK(eq_eps(status.E, 1 - sqrt(2)/2, 2*status.SE) == 1);
    162       }
    163     }
    164     CHK(sgf_estimator_get_status_infinity(estimator, 0, &status) == RES_OK);
    165     CHK(eq_eps(status.E, 0, status.SE) == 1);
    166 
    167     CHK(sgf_estimator_ref_put(estimator) == RES_OK);
    168   }
    169 
    170   CHK(sgf_integrate(scn, 4, rng,  NSTEPS, &estimator) == RES_BAD_ARG);
    171 
    172   CHK(sgf_scene_end_integration(scn) == RES_OK);
    173 
    174   /*
    175    * Check medium attenuation with 2 parallel infinite segments. To simulate
    176    * this configuration, the top and bottom segments of the square are fully
    177    * emissive.  The other ones are fully specular and have no emissivity
    178    */
    179 
    180   shape.absorption = ka;
    181   shape.emissivity = emissivity_inf_bottom_top;
    182   shape.specularity = specularity_inf_bottom_top;
    183   desc.get_absorption = get_absorption;
    184 
    185   FOR_EACH(i, 0, 4) ka[i] = 0;
    186   CHK(sgf_scene_setup_2d(scn, &desc) == RES_OK);
    187   check_bottom_top_medium_gf(scn, rng, 1, 0);
    188 
    189   FOR_EACH(i, 0, 4) ka[i] = 0.1;
    190   CHK(sgf_scene_setup_2d(scn, &desc) == RES_OK);
    191   check_bottom_top_medium_gf(scn, rng, 0.832583, 0.167417);
    192 
    193   FOR_EACH(i, 0, 4) ka[i] = 1;
    194   CHK(sgf_scene_setup_2d(scn, &desc) == RES_OK);
    195   check_bottom_top_medium_gf(scn, rng, 0.219384, 0.780616);
    196 
    197   FOR_EACH(i, 0, 4) ka[i] = 10;
    198   CHK(sgf_scene_setup_2d(scn, &desc) == RES_OK);
    199   check_bottom_top_medium_gf(scn, rng, 7.0975e-6, 0.999992902);
    200 
    201   CHK(ssp_rng_ref_put(rng) == RES_OK);
    202   CHK(sgf_device_ref_put(sgf) == RES_OK);
    203   CHK(sgf_scene_ref_put(scn) == RES_OK);
    204 
    205   check_memory_allocator(&allocator);
    206   mem_shutdown_proxy_allocator(&allocator);
    207   CHK(mem_allocated_size() == 0);
    208   return 0;
    209 }
    210