commit 9fe8666d10ed9142ecf6f823f5ebee9395e25cc6
parent 74fdffe2c6ad99adca85a8aa923961c5ae261524
Author: Christophe Coustet <christophe.coustet@meso-star.com>
Date: Wed, 22 Mar 2017 18:23:41 +0100
Add a new test that shows a bug in the solver (MC weight for quadrics).
Diffstat:
5 files changed, 298 insertions(+), 17 deletions(-)
diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt
@@ -169,6 +169,7 @@ if(NOT NO_TEST)
new_test(test_ssol_solver5)
new_test(test_ssol_solver6)
new_test(test_ssol_solver7)
+ new_test(test_ssol_solver8)
new_test(test_ssol_sun)
build_test(test_ssol_draw)
diff --git a/src/test_ssol_circ2D_geometry.h b/src/test_ssol_circ2D_geometry.h
@@ -0,0 +1,54 @@
+/* Copyright (C) CNRS 2016-2017
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. */
+
+#include "test_ssol_geometries.h"
+#include <rsys/math.h>
+
+
+/*******************************************************************************
+* Circle polygon
+******************************************************************************/
+#if !defined(RADIUS)
+#error "Missing the RADIUS macro defining the circle radius"
+#endif
+#if !defined(NVERTS)
+#define NVERTS 36
+#endif
+#if !defined(CIRCLE_NAME)
+#error "Missing the CIRCLE_NAME macro defining the circle name"
+#endif
+
+#define EDGES__ CONCAT(CIRCLE_NAME, _EDGES__)
+#define INIT_FUNC__ CONCAT(init_, CIRCLE_NAME)
+
+/* should be const but scpr expects non-const data */
+static double EDGES__ [2*NVERTS];
+
+static void INIT_FUNC__() {
+ int n;
+ /* radius that give the same area than a circle */
+ double r = sqrt(2 * RADIUS * RADIUS * PI / (sin(2 * PI / NVERTS) * NVERTS));
+ for (n = 0; n < NVERTS; n++) {
+ EDGES__[2 * n] = r * cos((double)-n * 2 * PI / (double)NVERTS);
+ EDGES__[2 * n + 1] = r * sin((double)-n * 2 * PI / (double) NVERTS);
+ }
+}
+
+#undef EDGES__
+#undef INIT_FUNC__
+
+#undef RADIUS
+#undef NVERTS
+#undef CIRCLE_NAME
diff --git a/src/test_ssol_rect2D_geometry.h b/src/test_ssol_rect2D_geometry.h
@@ -18,25 +18,41 @@
/*******************************************************************************
* Rectangle polygon
******************************************************************************/
-#if !defined(HALF_X)
-#error "Missing the HALF_X macro defining the rectangle size"
+#if !defined(HALF_X) && !(defined(X_MIN) && defined(X_MAX))
+#error "Missing the HALF_X or X_MIN and X_MAX macros defining the rectangle size"
#endif
-#if !defined(HALF_Y)
-#error "Missing the HALF_Y macro defining the rectangle size"
+#if !defined(HALF_Y) && !(defined(Y_MIN) && defined(Y_MAX))
+#error "Missing the HALF_Y or Y_MIN and Y_MAX macros defining the rectangle size"
#endif
#if !defined(POLYGON_NAME)
-#error "Missing the POLYGON_NAME macro defining the polygon name"
+#error "Missing the POLYGON_NAME macro defining the rectangle name"
#endif
#define EDGES__ CONCAT(POLYGON_NAME, _EDGES__)
#define POLY_NVERTS__ CONCAT(POLYGON_NAME, _NVERTS__)
+#if !defined(X_MIN)
+#define X_MIN (float)(-(HALF_X))
+#endif
+
+#if !defined(X_MAX)
+#define X_MAX (float)(HALF_X)
+#endif
+
+#if !defined(Y_MIN)
+#define Y_MIN (float)(-(HALF_Y))
+#endif
+
+#if !defined(Y_MAX)
+#define Y_MAX (float)(HALF_Y)
+#endif
+
/* should be const but scpr expects non-const data */
static double EDGES__ [] = {
- -HALF_X, -HALF_Y,
- -HALF_X, HALF_Y,
- HALF_X, HALF_Y,
- HALF_X, -HALF_Y,
+ X_MIN, Y_MIN,
+ X_MIN, Y_MAX,
+ X_MAX, Y_MAX,
+ X_MAX, Y_MIN
};
const unsigned POLY_NVERTS__ = sizeof(EDGES__) / sizeof(double[2]);
@@ -46,4 +62,8 @@ const unsigned POLY_NVERTS__ = sizeof(EDGES__) / sizeof(double[2]);
#undef HALF_X
#undef HALF_Y
+#undef X_MIN
+#undef X_MAX
+#undef Y_MIN
+#undef Y_MAX
#undef POLYGON_NAME
diff --git a/src/test_ssol_rect_geometry.h b/src/test_ssol_rect_geometry.h
@@ -18,11 +18,11 @@
/*******************************************************************************
* Rectangle plane
******************************************************************************/
-#if !defined(HALF_X)
-#error "Missing the HALF_X macro defining the rectangle size"
+#if !defined(HALF_X) && !(defined(X_MIN) && defined(X_MAX))
+#error "Missing the HALF_X or X_MIN and X_MAX macros defining the rectangle size"
#endif
-#if !defined(HALF_Y)
-#error "Missing the HALF_Y macro defining the rectangle size"
+#if !defined(HALF_Y) && !(defined(Y_MIN) && defined(Y_MAX))
+#error "Missing the HALF_Y or Y_MIN and Y_MAX macros defining the rectangle size"
#endif
#if !defined(PLANE_NAME)
#error "Missing the DARRAY_NAME macro defining the rectangle name"
@@ -34,11 +34,27 @@
#define RECT_NVERTS__ CONCAT(PLANE_NAME, _NVERTS__)
#define RECT_NTRIS__ CONCAT(PLANE_NAME, _NTRIS__)
+#if !defined(X_MIN)
+#define X_MIN (float)(-(HALF_X))
+#endif
+
+#if !defined(X_MAX)
+#define X_MAX (float)(HALF_X)
+#endif
+
+#if !defined(Y_MIN)
+#define Y_MIN (float)(-(HALF_Y))
+#endif
+
+#if !defined(Y_MAX)
+#define Y_MAX (float)(HALF_Y)
+#endif
+
static const float EDGES__ [] = {
- (float) -HALF_X, (float) -HALF_Y, 0.f,
- (float) HALF_X, (float) -HALF_Y, 0.f,
- (float) HALF_X, (float) HALF_Y, 0.f,
- (float) -HALF_X, (float) HALF_Y, 0.f
+ X_MIN, Y_MIN, 0.f,
+ X_MAX, Y_MIN, 0.f,
+ X_MAX, Y_MAX, 0.f,
+ X_MIN, Y_MAX, 0.f
};
const unsigned RECT_NVERTS__ = sizeof(EDGES__) / sizeof(float[3]);
@@ -56,4 +72,8 @@ static const struct desc RECT_DESC__ = { EDGES__, TRG_IDS__ };
#undef HALF_X
#undef HALF_Y
+#undef X_MIN
+#undef X_MAX
+#undef Y_MIN
+#undef Y_MAX
#undef PLANE_NAME
diff --git a/src/test_ssol_solver8.c b/src/test_ssol_solver8.c
@@ -0,0 +1,186 @@
+/* Copyright (C) CNRS 2016-2017
+*
+* This program is free software: you can redistribute it and/or modify
+* it under the terms of the GNU General Public License as published by
+* the Free Software Foundation, either version 3 of the License, or
+* (at your option) any later version.
+*
+* This program is distributed in the hope that it will be useful,
+* but WITHOUT ANY WARRANTY; without even the implied warranty of
+* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+* GNU General Public License for more details.
+*
+* You should have received a copy of the GNU General Public License
+* along with this program. If not, see <http://www.gnu.org/licenses/>. */
+
+#include "ssol.h"
+#include "test_ssol_utils.h"
+#include "test_ssol_materials.h"
+
+#define PLANE_NAME TARGET
+#define HALF_X 10
+#define HALF_Y 10
+#include "test_ssol_rect_geometry.h"
+
+#define X_SZ 10
+#define Y_SZ 4
+#define POLYGON_NAME POLY
+#define HALF_X (X_SZ / 2)
+STATIC_ASSERT((HALF_X * 2 == X_SZ), ONLY_ENVEN_VALUES_FOR_X_SZ);
+#define Y_MIN 0
+#define Y_MAX Y_SZ
+#include "test_ssol_rect2D_geometry.h"
+
+#include <rsys/double33.h>
+
+#include <star/s3d.h>
+#include <star/ssp.h>
+
+static void
+get_wlen(const size_t i, double* wlen, double* data, void* ctx)
+{
+ double wavelengths[3] = { 1, 2, 3 };
+ double intensities[3] = { 1, 0.8, 1 };
+ CHECK(i < 3, 1);
+ (void) ctx;
+ *wlen = wavelengths[i];
+ *data = intensities[i];
+}
+
+int
+main(int argc, char** argv)
+{
+ struct mem_allocator allocator;
+ struct ssol_device* dev;
+ struct ssp_rng* rng;
+ struct ssol_scene* scene;
+ struct ssol_shape* square;
+ struct ssol_vertex_data attribs[1] = { SSOL_VERTEX_DATA_NULL__ };
+ struct ssol_shape* quad_square;
+ struct ssol_carving carving = SSOL_CARVING_NULL;
+ struct ssol_quadric quadric = SSOL_QUADRIC_DEFAULT;
+ struct ssol_punched_surface punched = SSOL_PUNCHED_SURFACE_NULL;
+ struct ssol_material* m_mtl;
+ struct ssol_material* v_mtl;
+ struct ssol_mirror_shader shader = SSOL_MIRROR_SHADER_NULL;
+ struct ssol_object* m_object;
+ struct ssol_object* t_object;
+ struct ssol_instance* heliostat;
+ struct ssol_instance* target;
+ struct ssol_sun* sun;
+ struct ssol_spectrum* spectrum;
+ struct ssol_estimator* estimator;
+ struct ssol_mc_global mc_global;
+ struct ssol_mc_receiver mc_rcv;
+ double dir[3];
+ double transform[12]; /* 3x4 column major matrix */
+ size_t count;
+ FILE* tmp;
+ uint32_t r_id;
+
+ (void) argc, (void) argv;
+ d3_splat(transform + 9, 0);
+ d33_rotation_pitch(transform, PI); /* flip faces: invert normal */
+ transform[11] = 4; /* set it just above the parabolic cylinder */
+
+ mem_init_proxy_allocator(&allocator, &mem_default_allocator);
+
+ CHECK(ssol_device_create
+ (NULL, &allocator, SSOL_NTHREADS_DEFAULT, 0, &dev), RES_OK);
+
+#define DNI 1000
+ CHECK(ssp_rng_create(&allocator, &ssp_rng_threefry, &rng), RES_OK);
+ CHECK(ssol_spectrum_create(dev, &spectrum), RES_OK);
+ CHECK(ssol_spectrum_setup(spectrum, get_wlen, 3, NULL), RES_OK);
+ CHECK(ssol_sun_create_directional(dev, &sun), RES_OK);
+ CHECK(ssol_sun_set_direction(sun, d3(dir, 0, 1, -1)), RES_OK);
+ CHECK(ssol_sun_set_spectrum(sun, spectrum), RES_OK);
+ CHECK(ssol_sun_set_dni(sun, DNI), RES_OK);
+ CHECK(ssol_scene_create(dev, &scene), RES_OK);
+ CHECK(ssol_scene_attach_sun(scene, sun), RES_OK);
+
+ /* Create scene content */
+
+ CHECK(ssol_shape_create_mesh(dev, &square), RES_OK);
+ attribs[0].usage = SSOL_POSITION;
+ attribs[0].get = get_position;
+ CHECK(ssol_mesh_setup(square, TARGET_NTRIS__, get_ids,
+ TARGET_NVERTS__, attribs, 1, (void*) &TARGET_DESC__), RES_OK);
+
+ CHECK(ssol_shape_create_punched_surface(dev, &quad_square), RES_OK);
+ carving.get = get_polygon_vertices;
+ carving.operation = SSOL_AND;
+ carving.nb_vertices = POLY_NVERTS__;
+ carving.context = &POLY_EDGES__;
+ quadric.type = SSOL_QUADRIC_PARABOLIC_CYLINDER;
+ quadric.data.parabol.focal = 1;
+ punched.nb_carvings = 1;
+ punched.quadric = &quadric;
+ punched.carvings = &carving;
+ CHECK(ssol_punched_surface_setup(quad_square, &punched), RES_OK);
+
+ CHECK(ssol_material_create_mirror(dev, &m_mtl), RES_OK);
+ shader.normal = get_shader_normal;
+ shader.reflectivity = get_shader_reflectivity;
+ shader.roughness = get_shader_roughness;
+ CHECK(ssol_mirror_set_shader(m_mtl, &shader), RES_OK);
+ CHECK(ssol_material_create_virtual(dev, &v_mtl), RES_OK);
+
+ CHECK(ssol_object_create(dev, &m_object), RES_OK);
+ CHECK(ssol_object_add_shaded_shape(m_object, quad_square, m_mtl, m_mtl), RES_OK);
+ CHECK(ssol_object_instantiate(m_object, &heliostat), RES_OK);
+ CHECK(ssol_scene_attach_instance(scene, heliostat), RES_OK);
+
+ CHECK(ssol_object_create(dev, &t_object), RES_OK);
+ CHECK(ssol_object_add_shaded_shape(t_object, square, v_mtl, v_mtl), RES_OK);
+ CHECK(ssol_object_instantiate(t_object, &target), RES_OK);
+ CHECK(ssol_instance_set_transform(target, transform), RES_OK);
+ CHECK(ssol_instance_set_receiver(target, SSOL_FRONT, 0), RES_OK);
+ CHECK(ssol_instance_sample(target, 0), RES_OK);
+ CHECK(ssol_scene_attach_instance(scene, target), RES_OK);
+
+ NCHECK(tmp = tmpfile(), 0);
+#define N__ 100000
+#define GET_MC_RCV ssol_estimator_get_mc_receiver
+ CHECK(ssol_solve(scene, rng, N__, 0, tmp, &estimator), RES_OK);
+ CHECK(ssol_instance_get_id(target, &r_id), RES_OK);
+ CHECK(ssol_estimator_get_count(estimator, &count), RES_OK);
+ CHECK(count, N__);
+ CHECK(fclose(tmp), 0);
+ CHECK(ssol_estimator_get_failed_count(estimator, &count), RES_OK);
+ CHECK(count, 0);
+#define S (sqrt(2) * X_SZ * Y_SZ)
+ CHECK(ssol_estimator_get_mc_global(estimator, &mc_global), RES_OK);
+ printf("Cos = %g +/- %g\n", mc_global.cos_loss.E, mc_global.cos_loss.SE);
+ printf("Shadows = %g +/- %g\n", mc_global.shadowed.E, mc_global.shadowed.SE);
+ printf("Missing = %g +/- %g\n", mc_global.missing.E, mc_global.missing.SE);
+ CHECK(eq_eps(mc_global.cos_loss.E, 0, 1e-4), 1);
+ CHECK(eq_eps(mc_global.shadowed.E, 0, 1e-4), 1);
+ CHECK(eq_eps(mc_global.missing.E, 0, 1e-4), 1);
+ CHECK(GET_MC_RCV(estimator, target, SSOL_FRONT, &mc_rcv), RES_OK);
+ printf("Ir(target1) = %g +/- %g\n",
+ mc_rcv.integrated_irradiance.E, mc_rcv.integrated_irradiance.SE);
+ CHECK(eq_eps(mc_rcv.integrated_irradiance.E, S * DNI, 2 * mc_rcv.integrated_irradiance.SE), 1);
+
+ /* Free data */
+ CHECK(ssol_instance_ref_put(heliostat), RES_OK);
+ CHECK(ssol_instance_ref_put(target), RES_OK);
+ CHECK(ssol_object_ref_put(m_object), RES_OK);
+ CHECK(ssol_object_ref_put(t_object), RES_OK);
+ CHECK(ssol_shape_ref_put(square), RES_OK);
+ CHECK(ssol_shape_ref_put(quad_square), RES_OK);
+ CHECK(ssol_material_ref_put(m_mtl), RES_OK);
+ CHECK(ssol_material_ref_put(v_mtl), RES_OK);
+ CHECK(ssol_estimator_ref_put(estimator), RES_OK);
+ CHECK(ssol_device_ref_put(dev), RES_OK);
+ CHECK(ssol_scene_ref_put(scene), RES_OK);
+ CHECK(ssp_rng_ref_put(rng), RES_OK);
+ CHECK(ssol_spectrum_ref_put(spectrum), RES_OK);
+ CHECK(ssol_sun_ref_put(sun), RES_OK);
+
+ check_memory_allocator(&allocator);
+ mem_shutdown_proxy_allocator(&allocator);
+ CHECK(mem_allocated_size(), 0);
+
+ return 0;
+}
