test_sgf_tetrahedron.c (4881B)
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 100000L 28 29 static const float vertices[] = { 30 0.57735026918962576450f, 0.f, 0.f, 31 -0.28867513459481288225f, 0.5f, 0.f, 32 -0.28867513459481288225f, -0.5f, 0.f, 33 0.f, 0.f, 0.81649658092772603273f 34 }; 35 static const size_t nvertices = sizeof(vertices) / (3*sizeof(float)); 36 37 /* Front faces are CW. The normals point into the cube */ 38 static const unsigned indices[] = { 39 0, 1, 3, 40 0, 2, 1, 41 1, 2, 3, 42 0, 3, 2 43 }; 44 static const size_t nprims = sizeof(indices) / (3*sizeof(unsigned)); 45 46 static const double emissivity[] = { 0.5, 1.0, 1.0, 1.0 }; 47 static const double specularity[] = { 0.0, 0.0, 0.0, 0.0 }; 48 49 int 50 main(int argc, char** argv) 51 { 52 struct mem_allocator allocator; 53 struct shape_context shape; 54 struct sgf_device* sgf = NULL; 55 struct sgf_scene_desc desc = SGF_SCENE_DESC_NULL; 56 struct sgf_scene* scn = NULL; 57 struct sgf_status* status = NULL; 58 struct ssp_rng_proxy* proxy = NULL; 59 struct ssp_rng** rngs = NULL; 60 size_t i; 61 size_t nbuckets; 62 int iprim; 63 (void)argc, (void)argv; 64 65 mem_init_proxy_allocator(&allocator, &mem_default_allocator); 66 nbuckets = (unsigned)omp_get_num_procs(); 67 CHK(ssp_rng_proxy_create(&allocator, SSP_RNG_THREEFRY, nbuckets, &proxy) == RES_OK); 68 CHK(sgf_device_create(NULL, &allocator, 1, &sgf) == RES_OK); 69 CHK(sgf_scene_create(sgf, &scn) == RES_OK); 70 71 shape.vertices = vertices; 72 shape.nvertices = nvertices; 73 shape.indices = indices; 74 shape.nprimitives = nprims; 75 shape.emissivity = emissivity; 76 shape.specularity = specularity; 77 shape.dim = 3; 78 79 desc.get_position = get_position; 80 desc.get_indices = get_indices; 81 desc.get_emissivity = get_emissivity; 82 desc.get_specularity = get_specularity; 83 desc.get_reflectivity = get_reflectivity; 84 desc.nprims = (unsigned)nprims; 85 desc.nverts = (unsigned)nvertices; 86 desc.nbands = 1; 87 desc.context = &shape; 88 89 CHK(sgf_scene_setup_3d(scn, &desc) == RES_OK); 90 91 status = sa_add(status, nprims*nprims); 92 rngs = sa_add(rngs, nbuckets); 93 94 FOR_EACH(i, 0, nbuckets) 95 CHK(ssp_rng_proxy_create_rng(proxy, i, rngs + i) == RES_OK); 96 97 CHK(sgf_scene_begin_integration(scn) == RES_OK); 98 99 /* Gebhart Factor integration */ 100 #pragma omp parallel for 101 for(iprim = 0; iprim < (int)nprims; ++iprim) { 102 struct sgf_status* row = status + (size_t)iprim*nprims; 103 struct sgf_estimator* estimator = NULL; 104 struct sgf_status infinity; 105 size_t iprim2; 106 const int ithread = omp_get_thread_num(); 107 108 CHK(sgf_integrate 109 (scn, (size_t)iprim, rngs[ithread], NSTEPS, &estimator) == RES_OK); 110 111 FOR_EACH(iprim2, 0, nprims) { 112 CHK(sgf_estimator_get_status(estimator, iprim2, 0, row + iprim2) == RES_OK); 113 CHK(row[iprim2].nsteps == NSTEPS); 114 } 115 116 CHK(sgf_estimator_get_status_infinity(estimator, 0, &infinity) == RES_OK); 117 CHK(eq_eps(infinity.E, 0, infinity.SE) == 1); 118 119 CHK(sgf_estimator_ref_put(estimator) == RES_OK); 120 } 121 CHK(sgf_scene_end_integration(scn) == RES_OK); 122 123 /* Expected value */ 124 for(iprim=0; iprim < (int)nprims; ++iprim) { 125 const struct sgf_status* row = status + (size_t)iprim * nprims; 126 unsigned icol; 127 double sum = 0.0; 128 FOR_EACH(icol, 0, nprims) { 129 sum += row[icol].E; 130 printf("%.6f ", row[icol].E); 131 } 132 printf("\n"); 133 /* Ensure the energy conservation property */ 134 CHK(eq_eps(sum, 1.0, 1.e-3) == 1); 135 } 136 printf("\n"); 137 138 /* Standard error */ 139 for(iprim=0; iprim < (int)nprims; ++iprim) { 140 const struct sgf_status* row = status + (size_t)iprim * nprims; 141 unsigned icol; 142 FOR_EACH(icol, 0, nprims) { 143 printf("%.6f ", row[icol].SE); 144 } 145 printf("\n"); 146 } 147 148 CHK(ssp_rng_proxy_ref_put(proxy) == RES_OK); 149 CHK(sgf_device_ref_put(sgf) == RES_OK); 150 CHK(sgf_scene_ref_put(scn) == RES_OK); 151 FOR_EACH(i, 0, nbuckets) 152 CHK(ssp_rng_ref_put(rngs[i]) == RES_OK); 153 154 sa_release(rngs); 155 sa_release(status); 156 157 check_memory_allocator(&allocator); 158 mem_shutdown_proxy_allocator(&allocator); 159 CHK(mem_allocated_size() == 0); 160 return 0; 161 } 162