test_ssol_solver2b.c (9130B)
1 /* Copyright (C) 2018-2026 |Meso|Star> (contact@meso-star.com) 2 * Copyright (C) 2016, 2018 CNRS 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 "ssol.h" 18 #include "test_ssol_utils.h" 19 #include "test_ssol_materials.h" 20 21 #define PLANE_NAME RECT 22 #define HALF_X 0.5 23 #define HALF_Y 1 24 #include "test_ssol_rect_geometry.h" 25 26 #define POLYGON_NAME RECT_POLY 27 #define HALF_X 0.5 28 #define HALF_Y 1 29 #include "test_ssol_rect2D_geometry.h" 30 31 #define POLYGON_NAME SQUARE_POLY 32 #define HALF_X 1 33 #define HALF_Y 1 34 #include "test_ssol_rect2D_geometry.h" 35 36 #include <rsys/double33.h> 37 38 #include <star/s3d.h> 39 #include <star/ssp.h> 40 41 static void 42 get_wlen(const size_t i, double* wlen, double* data, void* ctx) 43 { 44 double wavelengths[3] = { 1, 2, 3 }; 45 double intensities[3] = { 1, 0.8, 1 }; 46 CHK(i < 3); 47 (void)ctx; 48 *wlen = wavelengths[i]; 49 *data = intensities[i]; 50 } 51 52 int 53 main(int argc, char** argv) 54 { 55 struct mem_allocator allocator; 56 struct ssol_device* dev; 57 struct ssp_rng* rng; 58 struct ssol_scene* scene; 59 struct ssol_shape* rect; 60 struct ssol_vertex_data attribs[1] = { SSOL_VERTEX_DATA_NULL__ }; 61 struct ssol_shape* quad_square; 62 struct ssol_shape* quad_rect; 63 struct ssol_carving carving = SSOL_CARVING_NULL; 64 struct ssol_quadric quadric = SSOL_QUADRIC_DEFAULT; 65 struct ssol_punched_surface punched = SSOL_PUNCHED_SURFACE_NULL; 66 struct ssol_material* m_mtl; 67 struct ssol_material* v_mtl; 68 struct ssol_mirror_shader shader = SSOL_MIRROR_SHADER_NULL; 69 struct ssol_object* m_object; 70 struct ssol_object* s_object; 71 struct ssol_object* t_object; 72 struct ssol_instance* heliostat1; 73 struct ssol_instance* heliostat2; 74 struct ssol_instance* secondary; 75 struct ssol_instance* target; 76 struct ssol_sun* sun; 77 struct ssol_spectrum* spectrum; 78 struct ssol_estimator* estimator; 79 struct ssol_mc_global mc_global; 80 struct ssol_mc_receiver mc_rcv; 81 double dir[3]; 82 double transform1[12]; /* 3x4 column major matrix */ 83 double transform2[12]; /* 3x4 column major matrix */ 84 double transform3[12]; /* 3x4 column major matrix */ 85 size_t count; 86 double m, std; 87 (void) argc, (void) argv; 88 89 d33_splat(transform1, 0); 90 d3_splat(transform1 + 9, 0); 91 d33_rotation_pitch(transform1, PI); /* flip faces: invert normal */ 92 transform1[9] = 2; /* +2 offset along X axis */ 93 transform1[11] = 2; /* +2 offset along Z axis */ 94 95 d33_set_identity(transform2); 96 d3_splat(transform2 + 9, 0); 97 transform2[9] = 4; /* +4 offset along X axis */ 98 99 d33_set_identity(transform3); 100 d3_splat(transform3 + 9, 0); 101 102 mem_init_proxy_allocator(&allocator, &mem_default_allocator); 103 104 CHK(ssol_device_create 105 (NULL, &allocator, SSOL_NTHREADS_DEFAULT, 0, &dev) == RES_OK); 106 107 CHK(ssp_rng_create(&allocator, SSP_RNG_THREEFRY, &rng) == RES_OK); 108 CHK(ssol_spectrum_create(dev, &spectrum) == RES_OK); 109 CHK(ssol_spectrum_setup(spectrum, get_wlen, 3, NULL) == RES_OK); 110 CHK(ssol_sun_create_directional(dev, &sun) == RES_OK); 111 CHK(ssol_sun_set_direction(sun, d3(dir, 1, 0, -1)) == RES_OK); 112 CHK(ssol_sun_set_spectrum(sun, spectrum) == RES_OK); 113 CHK(ssol_sun_set_dni(sun, 1000) == RES_OK); 114 CHK(ssol_scene_create(dev, &scene) == RES_OK); 115 CHK(ssol_scene_attach_sun(scene, sun) == RES_OK); 116 117 /* Create scene content */ 118 119 CHK(ssol_shape_create_mesh(dev, &rect) == RES_OK); 120 attribs[0].usage = SSOL_POSITION; 121 attribs[0].get = get_position; 122 CHK(ssol_mesh_setup(rect, RECT_NTRIS__, get_ids, 123 RECT_NVERTS__, attribs, 1, (void*) &RECT_DESC__) == RES_OK); 124 125 CHK(ssol_shape_create_punched_surface(dev, &quad_rect) == RES_OK); 126 carving.get = get_polygon_vertices; 127 carving.operation = SSOL_AND; 128 carving.nb_vertices = RECT_POLY_NVERTS__; 129 carving.context = &RECT_POLY_EDGES__; 130 quadric.type = SSOL_QUADRIC_PLANE; 131 punched.nb_carvings = 1; 132 punched.quadric = &quadric; 133 punched.carvings = &carving; 134 CHK(ssol_punched_surface_setup(quad_rect, &punched) == RES_OK); 135 136 CHK(ssol_shape_create_punched_surface(dev, &quad_square) == RES_OK); 137 carving.get = get_polygon_vertices; 138 carving.operation = SSOL_AND; 139 carving.nb_vertices = SQUARE_POLY_NVERTS__; 140 carving.context = &SQUARE_POLY_EDGES__; 141 quadric.type = SSOL_QUADRIC_PLANE; 142 punched.nb_carvings = 1; 143 punched.quadric = &quadric; 144 punched.carvings = &carving; 145 CHK(ssol_punched_surface_setup(quad_square, &punched) == RES_OK); 146 147 CHK(ssol_material_create_mirror(dev, &m_mtl) == RES_OK); 148 shader.normal = get_shader_normal; 149 shader.reflectivity = get_shader_reflectivity; 150 shader.roughness = get_shader_roughness; 151 CHK(ssol_mirror_setup(m_mtl, &shader, SSOL_MICROFACET_BECKMANN) == RES_OK); 152 CHK(ssol_material_create_virtual(dev, &v_mtl) == RES_OK); 153 154 CHK(ssol_object_create(dev, &m_object) == RES_OK); 155 CHK(ssol_object_add_shaded_shape(m_object, quad_rect, m_mtl, m_mtl) == RES_OK); 156 CHK(ssol_object_instantiate(m_object, &heliostat1) == RES_OK); 157 CHK(ssol_object_instantiate(m_object, &heliostat2) == RES_OK); 158 CHK(ssol_instance_set_receiver(heliostat1, SSOL_FRONT, 0) == RES_OK); 159 CHK(ssol_instance_set_receiver(heliostat2, SSOL_FRONT, 0) == RES_OK); 160 transform3[9] = -0.5; /* -0.5 offset along X axis */ 161 CHK(ssol_instance_set_transform(heliostat1, transform3) == RES_OK); 162 transform3[9] = +0.5; /* +0.5 offset along X axis */ 163 CHK(ssol_instance_set_transform(heliostat2, transform3) == RES_OK); 164 CHK(ssol_scene_attach_instance(scene, heliostat1) == RES_OK); 165 CHK(ssol_scene_attach_instance(scene, heliostat2) == RES_OK); 166 167 CHK(ssol_object_create(dev, &s_object) == RES_OK); 168 CHK(ssol_object_add_shaded_shape(s_object, quad_square, m_mtl, m_mtl) == RES_OK); 169 CHK(ssol_object_instantiate(s_object, &secondary) == RES_OK); 170 CHK(ssol_instance_set_receiver(secondary, SSOL_FRONT, 0) == RES_OK); 171 CHK(ssol_instance_set_transform(secondary, transform1) == RES_OK); 172 CHK(ssol_instance_sample(secondary, 0) == RES_OK); 173 CHK(ssol_scene_attach_instance(scene, secondary) == RES_OK); 174 175 CHK(ssol_object_create(dev, &t_object) == RES_OK); 176 CHK(ssol_object_add_shaded_shape(t_object, rect, v_mtl, v_mtl) == RES_OK); 177 CHK(ssol_object_instantiate(t_object, &target) == RES_OK); 178 CHK(ssol_instance_set_transform(target, transform2) == RES_OK); 179 CHK(ssol_instance_set_receiver(target, SSOL_FRONT, 0) == RES_OK); 180 CHK(ssol_instance_sample(target, 0) == RES_OK); 181 CHK(ssol_scene_attach_instance(scene, target) == RES_OK); 182 183 #define N__ 50000 184 CHK(ssol_solve(scene, rng, N__, 0, NULL, &estimator) == RES_OK); 185 CHK(ssol_estimator_get_realisation_count(estimator, &count) == RES_OK); 186 CHK(count == N__); 187 #define COS cos(PI / 4) 188 #define DNI_cos (1000 * COS) 189 m = 2 * DNI_cos; 190 #define SQR(x) ((x)*(x)) 191 std = sqrt((SQR(4 * DNI_cos) / 2 - SQR(2 * DNI_cos)) / (double)count); 192 CHK(ssol_estimator_get_mc_global(estimator, &mc_global) == RES_OK); 193 printf("Shadows = %g +/- %g\n", mc_global.shadowed.E, mc_global.shadowed.SE); 194 printf("Missing = %g +/- %g\n", mc_global.missing.E, mc_global.missing.SE); 195 printf("Cos = %g +/- %g\n", mc_global.cos_factor.E, mc_global.cos_factor.SE); 196 CHK(eq_eps(mc_global.shadowed.E, 0, 1e-4) == 1); 197 CHK(eq_eps(mc_global.missing.E, 4 * DNI_cos, 1e-4) == 1); /* nothing absorbed */ 198 CHK(eq_eps(mc_global.cos_factor.E, COS, 1e-4) == 1); 199 CHK(ssol_estimator_get_mc_receiver 200 (estimator, target, SSOL_FRONT, &mc_rcv) == RES_OK); 201 printf("Ir(target) = %g +/- %g\n", 202 mc_rcv.incoming_flux.E, mc_rcv.incoming_flux.SE); 203 CHK(eq_eps(mc_rcv.incoming_flux.E, m, 2 * std) == 1); 204 CHK(eq_eps(mc_rcv.incoming_flux.SE, std, 1e-3) == 1); 205 CHK(ssol_estimator_get_failed_count(estimator, &count) == RES_OK); 206 CHK(count == 0); 207 208 /* Free data */ 209 CHK(ssol_instance_ref_put(heliostat1) == RES_OK); 210 CHK(ssol_instance_ref_put(heliostat2) == RES_OK); 211 CHK(ssol_instance_ref_put(secondary) == RES_OK); 212 CHK(ssol_instance_ref_put(target) == RES_OK); 213 CHK(ssol_object_ref_put(m_object) == RES_OK); 214 CHK(ssol_object_ref_put(s_object) == RES_OK); 215 CHK(ssol_object_ref_put(t_object) == RES_OK); 216 CHK(ssol_shape_ref_put(rect) == RES_OK); 217 CHK(ssol_shape_ref_put(quad_rect) == RES_OK); 218 CHK(ssol_shape_ref_put(quad_square) == RES_OK); 219 CHK(ssol_material_ref_put(m_mtl) == RES_OK); 220 CHK(ssol_material_ref_put(v_mtl) == RES_OK); 221 CHK(ssol_estimator_ref_put(estimator) == RES_OK); 222 CHK(ssol_device_ref_put(dev) == RES_OK); 223 CHK(ssol_scene_ref_put(scene) == RES_OK); 224 CHK(ssp_rng_ref_put(rng) == RES_OK); 225 CHK(ssol_spectrum_ref_put(spectrum) == RES_OK); 226 CHK(ssol_sun_ref_put(sun) == RES_OK); 227 228 check_memory_allocator(&allocator); 229 mem_shutdown_proxy_allocator(&allocator); 230 CHK(mem_allocated_size() == 0); 231 232 return 0; 233 }