test_sgf_cube.c (9220B)
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 <rsys/stretchy_array.h> 22 23 #include <star/ssp.h> 24 25 #include <omp.h> 26 27 #define NSTEPS 10000L 28 29 /* 30 * Analytic Gebhart Factor matrix 31 * 32 * 0.065871765 0.245053213 0.281090450 0.210375872 0.0838339939 0.113774706 33 * 0.183789910 0.100632183 0.291901621 0.218467251 0.0870583783 0.118150656 34 * 0.187393634 0.259468108 0.121154594 0.222750923 0.0887654054 0.120467336 35 * 0.180322176 0.249676859 0.286394044 0.082269756 0.0854157674 0.115921399 36 * 0.167667988 0.232155676 0.266296216 0.199303457 0.0267900995 0.107786564 37 * 0.170662059 0.236301313 0.271051506 0.202862448 0.0808399227 0.038282752 38 */ 39 40 static const float vertices[] = { 41 0.f, 0.f, 0.f, 42 1.f, 0.f, 0.f, 43 0.f, 1.f, 0.f, 44 1.f, 1.f, 0.f, 45 0.f, 0.f, 1.f, 46 1.f, 0.f, 1.f, 47 0.f, 1.f, 1.f, 48 1.f, 1.f, 1.f 49 }; 50 static const size_t nvertices = sizeof(vertices) / (3*sizeof(float)); 51 52 /* Front faces are CW. The normals point into the cube */ 53 static const unsigned indices[] = { 54 0, 2, 1, 1, 2, 3, /* Front */ 55 0, 4, 2, 2, 4, 6, /* Left */ 56 4, 5, 6, 6, 5, 7, /* Back */ 57 3, 7, 1, 1, 7, 5, /* Right */ 58 2, 6, 3, 3, 6, 7, /* Top */ 59 0, 1, 4, 4, 1, 5 /* Bottom */ 60 }; 61 static const size_t nprims = sizeof(indices) / (3*sizeof(unsigned)); 62 63 static const double emissivity[] = { 64 0.6, 0.6, /* Front */ 65 0.8, 0.8, /* Left */ 66 0.9, 0.9, /* Back */ 67 0.7, 0.7, /* Right */ 68 0.3, 0.3, /* Top */ 69 0.4, 0.4 /* Bottom */ 70 }; 71 72 /* Emissivity used to simulate 2 infinite planes */ 73 static const double emissivity_inf_bottom_top[] = { 74 0, 0, /* Front */ 75 0, 0, /* Left */ 76 0, 0, /* Back */ 77 0, 0, /* Right */ 78 1, 1, /* Top */ 79 1, 1 /* Bottom */ 80 }; 81 82 static const double specularity[] = { 83 0.0, 0.0, /* Front */ 84 0.0, 0.0, /* Left */ 85 0.0, 0.0, /* Back */ 86 0.0, 0.0, /* Right */ 87 0.0, 0.0, /* Top */ 88 0.0, 0.0 /* Bottom */ 89 }; 90 91 /* Specularity used to simulate 2 infinite planes */ 92 static const double specularity_inf_bottom_top[] = { 93 1.0, 1.0, /* Front */ 94 1.0, 1.0, /* Left */ 95 1.0, 1.0, /* Back */ 96 1.0, 1.0, /* Right */ 97 0.0, 0.0, /* Top */ 98 0.0, 0.0 /* Bottom */ 99 }; 100 101 /* Check the estimation of the bottom/top & bottom/medium Gebhart factors */ 102 static void 103 check_bottom_top_medium_gf 104 (struct sgf_scene* scn, 105 struct ssp_rng* rng, 106 const double gf_bottom_top, /* Ref of the bottom/top Gebhart factor */ 107 const double gf_bottom_medium) /* Ref of the bottom/medium Gebhart Factor */ 108 { 109 /* Indices of the top/bottom primitives */ 110 const size_t TOP0 = 8; 111 const size_t TOP1 = 9; 112 const size_t BOTTOM0 = 10; 113 const size_t BOTTOM1 = 11; 114 115 struct sgf_estimator* estimator; 116 struct sgf_status status[6]; 117 double E, SE; 118 119 CHK(sgf_scene_begin_integration(scn) == RES_OK); 120 121 /* Estimate the Gebhart factors for the 1st triangle of the bottom face */ 122 CHK(sgf_integrate(scn, BOTTOM0, rng, NSTEPS, &estimator) == RES_OK); 123 CHK(sgf_estimator_get_status(estimator, TOP0, 0, status + 0) == RES_OK); 124 CHK(sgf_estimator_get_status(estimator, TOP1, 0, status + 1) == RES_OK); 125 CHK(sgf_estimator_get_status_medium(estimator, 0, status + 2) == RES_OK); 126 CHK(sgf_estimator_get_status_medium(estimator, 0, status + 3) == RES_OK); 127 CHK(sgf_estimator_ref_put(estimator) == RES_OK); 128 129 /* Estimate the Gebhart factors for the 2nd triangle of the bottom face */ 130 CHK(sgf_integrate(scn, BOTTOM1, rng, NSTEPS, &estimator) == RES_OK); 131 CHK(sgf_estimator_get_status(estimator, TOP0, 0, status + 4) == RES_OK); 132 CHK(sgf_estimator_get_status(estimator, TOP1, 0, status + 5) == RES_OK); 133 CHK(sgf_estimator_ref_put(estimator) == RES_OK); 134 135 /* Check the Gebhart factor between the bottom and the top plane. */ 136 E = (status[0].E + status[1].E + status[4].E + status[5].E)/2; 137 SE = (status[0].SE + status[1].SE + status[4].SE + status[5].SE)/2; 138 CHK(eq_eps(E, gf_bottom_top, SE) == 1); 139 140 /* Check the Gebhart factor between the bottom plane and the medium */ 141 E = (status[2].E + status[3].E)/2; 142 SE = (status[2].SE + status[3].SE)/2; 143 CHK(eq_eps(E, gf_bottom_medium, SE*3) == 1); 144 145 CHK(sgf_scene_end_integration(scn) == RES_OK); 146 } 147 148 int 149 main(int argc, char** argv) 150 { 151 struct mem_allocator allocator; 152 struct shape_context shape; 153 struct sgf_scene* scn; 154 struct sgf_scene_desc desc = SGF_SCENE_DESC_NULL; 155 struct sgf_device* sgf = NULL; 156 struct sgf_status* status = NULL; 157 struct ssp_rng_proxy* proxy; 158 struct ssp_rng** rngs = NULL; 159 double ka[12]; 160 int iprim; 161 unsigned i; 162 unsigned nbuckets; 163 (void)argc, (void)argv; 164 165 mem_init_proxy_allocator(&allocator, &mem_default_allocator); 166 nbuckets = (unsigned)omp_get_num_procs(); 167 168 CHK(ssp_rng_proxy_create(&allocator, SSP_RNG_THREEFRY, nbuckets, &proxy) == RES_OK); 169 CHK(sgf_device_create(NULL, &allocator, 1, &sgf) == RES_OK); 170 CHK(sgf_scene_create(sgf, &scn) == RES_OK); 171 172 shape.vertices = vertices; 173 shape.nvertices = nvertices; 174 shape.indices = indices; 175 shape.nprimitives = nprims; 176 shape.emissivity = emissivity; 177 shape.specularity = specularity; 178 shape.dim = 3; 179 180 desc.get_position = get_position; 181 desc.get_indices = get_indices; 182 desc.get_emissivity = get_emissivity; 183 desc.get_specularity = get_specularity; 184 desc.get_reflectivity = get_reflectivity; 185 desc.nprims = (unsigned)nprims; 186 desc.nverts = (unsigned)nvertices; 187 desc.nbands = 1; 188 desc.context = &shape; 189 190 CHK(sgf_scene_setup_3d(scn, &desc) == RES_OK); 191 192 status = sa_add(status, nprims*nprims); 193 rngs = sa_add(rngs, nbuckets); 194 195 FOR_EACH(i, 0, nbuckets) 196 CHK(ssp_rng_proxy_create_rng(proxy, i, rngs + i) == RES_OK); 197 198 CHK(sgf_scene_begin_integration(scn) == RES_OK); 199 200 /* Integrate the Gebhart Factors */ 201 #pragma omp parallel for 202 for(iprim = 0; iprim < (int)nprims; ++iprim) { 203 size_t iprim2; 204 struct sgf_status* row = status + (size_t)iprim * nprims; 205 struct sgf_estimator* est = NULL; 206 struct sgf_status infinity; 207 const int ithread = omp_get_thread_num(); 208 209 CHK(sgf_integrate 210 (scn, (size_t)iprim, rngs[ithread], NSTEPS, &est) == RES_OK); 211 212 FOR_EACH(iprim2, 0, nprims) { 213 CHK(sgf_estimator_get_status(est, iprim2, 0, row + iprim2) == RES_OK); 214 CHK(row[iprim2].nsteps == NSTEPS); 215 } 216 217 CHK(sgf_estimator_get_status_infinity(est, 0, &infinity) == RES_OK); 218 CHK(eq_eps(infinity.E, 0, infinity.SE) == 1); 219 220 CHK(sgf_estimator_ref_put(est) == RES_OK); 221 } 222 223 /* Merge the radiative flux of coplanar primitives */ 224 for(iprim=0; iprim < (int)(nprims/2); ++iprim) { 225 const struct sgf_status* row_src0 = status + (size_t)iprim * 2 * nprims; 226 const struct sgf_status* row_src1 = row_src0 + nprims; 227 size_t icol; 228 double sum = 0; 229 230 FOR_EACH(icol, 0, nprims/2) { 231 const struct sgf_status* src0 = row_src0 + icol * 2; 232 const struct sgf_status* src1 = src0 + 1; 233 const struct sgf_status* src2 = row_src1 + icol * 2; 234 const struct sgf_status* src3 = src2 + 1; 235 double E = (src0->E + src1->E + src2->E + src3->E) / 2; 236 237 sum += E; 238 printf("%.6f ", E); 239 } 240 printf("\n"); 241 CHK(eq_eps(sum, 1.0, 1.e-3) == 1); /* Ensure the energy conservation */ 242 } 243 244 CHK(sgf_scene_end_integration(scn) == RES_OK); 245 246 /* 247 * Check medium attenuation with 2 parallel infinite planes. To simulate 248 * this configuration, the top and bottom faces of the cube are fully 249 * emissive. The other ones are fully specular and have no emissivity 250 */ 251 252 shape.absorption = ka; 253 shape.emissivity = emissivity_inf_bottom_top; 254 shape.specularity = specularity_inf_bottom_top; 255 desc.get_absorption = get_absorption; 256 257 FOR_EACH(i, 0, 12) ka[i] = 0; 258 CHK(sgf_scene_setup_3d(scn, &desc) == RES_OK); 259 check_bottom_top_medium_gf(scn, rngs[0], 1, 0); 260 261 FOR_EACH(i, 0, 12) ka[i] = 0.1; 262 CHK(sgf_scene_setup_3d(scn, &desc) == RES_OK); 263 check_bottom_top_medium_gf(scn, rngs[0], 0.832583, 0.167417); 264 265 FOR_EACH(i, 0, 12) ka[i] = 1; 266 CHK(sgf_scene_setup_3d(scn, &desc) == RES_OK); 267 check_bottom_top_medium_gf(scn, rngs[0], 0.219384, 0.780616); 268 269 FOR_EACH(i, 0, 12) ka[i] = 10; 270 CHK(sgf_scene_setup_3d(scn, &desc) == RES_OK); 271 check_bottom_top_medium_gf(scn, rngs[0], 7.0975e-6, 0.999992902); 272 273 CHK(ssp_rng_proxy_ref_put(proxy) == RES_OK); 274 CHK(sgf_device_ref_put(sgf) == RES_OK); 275 CHK(sgf_scene_ref_put(scn) == RES_OK); 276 FOR_EACH(i, 0, nbuckets) 277 CHK(ssp_rng_ref_put(rngs[i]) == RES_OK); 278 sa_release(rngs); 279 sa_release(status); 280 281 check_memory_allocator(&allocator); 282 mem_shutdown_proxy_allocator(&allocator); 283 CHK(mem_allocated_size() == 0); 284 return 0; 285 } 286