diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index 7d02d3cc..56b6cdb6 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -61,6 +61,43 @@ jobs: sudo apt-get update sudo apt-get install -y libfontconfig1-dev + - name: Install embree (Ubuntu) + if: matrix.os == 'ubuntu-latest' + run: | + sudo apt-get update + sudo apt-get install -y libembree-dev + + - name: Install embree (macOS) + if: matrix.os == 'macOS-latest' + run: | + brew install embree + + - name: Install embree (Windows) + if: matrix.os == 'windows-latest' + shell: pwsh + run: | + $embreeVersion = "4.4.0" + $embreeUrl = "https://github.com/embree/embree/releases/download/v$embreeVersion/embree-$embreeVersion.x64.windows.zip" + $embreeZip = "$env:TEMP\embree.zip" + $embreeDir = "C:\embree" + + Write-Host "Downloading Embree $embreeVersion..." + Invoke-WebRequest -Uri $embreeUrl -OutFile $embreeZip + + Write-Host "Extracting Embree..." + Expand-Archive -Path $embreeZip -DestinationPath $embreeDir -Force + + # $embreeRoot = Get-ChildItem -Path $embreeDir -Directory -Filter "embree*" | Select-Object -First 1 + $embreeCMake = Join-Path $embreeDir "lib\cmake\embree-$embreeVersion" + $embreeLib = Join-Path $embreeDir "bin" + + # Write-Host "Embree root: $($embreeRoot.FullName)" + Write-Host "Embree CMake: $embreeCMake" + Write-Host "Embree Lib: $embreeLib" + + echo "embree_DIR=$embreeCMake" >> $env:GITHUB_ENV + echo "$embreeLib" >> $env:GITHUB_PATH + - name: Set reusable strings # Turn repeated input strings (such as the build output directory) into step outputs. These step outputs can be used throughout the workflow file. id: strings @@ -77,6 +114,7 @@ jobs: -DCMAKE_C_COMPILER=${{ matrix.c_compiler }} -DCMAKE_BUILD_TYPE=${{ matrix.build_type }} -DSOLTRACE_BUILD_GUI=OFF + -DSOLTRACE_BUILD_EMBREE_SUPPORT=ON -S ${{ github.workspace }} - name: Build-Windows diff --git a/.github/workflows/coverage.yml b/.github/workflows/coverage.yml index 25d00466..a818b3e0 100644 --- a/.github/workflows/coverage.yml +++ b/.github/workflows/coverage.yml @@ -19,7 +19,7 @@ jobs: - name: Install dependencies run: | sudo apt-get update - sudo apt-get install -y lcov libfontconfig1-dev + sudo apt-get install -y lcov libfontconfig1-dev libembree-dev - name: Set reusable strings id: strings @@ -43,6 +43,7 @@ jobs: -DSOLTRACE_BUILD_CORETRACE=OFF -DSOLTRACE_BUILD_GUI=OFF -DENABLE_COVERAGE=ON + -DSOLTRACE_BUILD_EMBREE_SUPPORT=ON -S ${{ github.workspace }} - name: Build diff --git a/CMakeLists.txt b/CMakeLists.txt index fae21153..c28d99ad 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -6,8 +6,14 @@ set(CMAKE_CXX_STANDARD 17) set(CMAKE_CXX_STANDARD_REQUIRED True) option(SOLTRACE_BUILD_GUI "Build the GUI parts of Soltrace" OFF) -option(SOLTRACE_BUILD_CORETRACE "Build the coretrace library (ray tracing core)" ON) +option( + SOLTRACE_BUILD_CORETRACE + "Legacy ray tracing backend including the strace executable" + OFF) + +option(SOLTRACE_BUILD_EMBREE_SUPPORT "Build Embree support for ray tracing" OFF) option(SOLTRACE_BUILD_OPTIX_SUPPORT "Build the OptiX support for ray tracing" OFF) + option(SOLTRACE_BUILD_PERF_TEST "Build performance tests (not built for debug)" OFF) option(FULL_FIELD_VALIDATION_TEST "Include full field validation tests" OFF) diff --git a/README.md b/README.md index 9d7155f8..3f716574 100644 --- a/README.md +++ b/README.md @@ -21,7 +21,7 @@ The `pysoltrace` API is capable of running multi-threaded simulations, generatin The API requires the compiled coretrace library. Project files for building this library are generated using CMake as outlined in the steps below. It is possible to build only coretrace and not build the graphical interface by following the steps 1-7, but only building the `coretrace_api` project in step 7.vii. -## Steps for Building SolTrace +## Steps for Building SolTrace (Legacy) These are the general steps you need to follow to set up your computer for developing SolTrace: @@ -78,6 +78,62 @@ These are the general steps you need to follow to set up your computer for devel Note that output is NOT stored in the ```build-soltrace/``` directory! +## Steps for Building SolTrace (work in progress) + +SolTrace has been updated to use multiple ray tracing engines in addition to the prior implementation. Currently, there is no graphical user interface (it is underdevelopment). + +Building SolTrace (develop branch) requires a C++-17 capable compiler and cmake 3.19 or greater. Once these are available, building can be done in the normal pattern of configure and build: + +```sh +git clone https://github.com/NatLabRockies/SolTrace.git +cd SolTrace +mkdir build +cd build +cmake .. +cmake --build . -j4 +``` + +Note the `-j4` instructs cmake to use 4 processes to compile the source code. This can be adjusted if the number of processors available is greater (or less) than 4. + +### Building with Intel's Embree Ray Tracing Library + +Information about Embree can be found on the [Embree webpage](https://www.embree.org/) or on the [Embree github page](https://github.com/RenderKit/embree). + +Prior to building SolTrace, you will need to install Embree v4.x.x. On Linux this is best done with a package manager. E.g., + +```sh +sudo apt-get update +sudo apt-get install libembree-dev +``` +On Mac, you can use Homebrew + +```sh +brew install embree +``` + +On Windows (this works for Linux and MacOS as well), you need to download the binaries from the [github page](https://github.com/RenderKit/embree) (under the appropriate installation header). Follow the corresponding install instructions found there. + +Once Embree is installed, clone the SolTrace repo, configure with embree enabled, and build: + +```sh +git clone https://github.com/NatLabRockies/SolTrace.git +cd SolTrace +mkdir build +cd build +cmake .. -DSOLTRACE_BUILD_EMBREE_SUPPORT=ON +cmake --build . -j4 +``` + +If cmake is having trouble locating the Embree install, you specify Embree's install location passing the `embree_DIR` variable to cmake. In this case the configure command would be + +```sh +cmake .. -DSOLTRACE_BUILD_EMBREE_SUPPORT=ON -Dembree_DIR= +``` + +### Building with Nvidia's Optix Ray Tracing Library + +TODO. + ## Contributing If you would like to report an issue with SolTrace or make a feature request, please let us know by adding a new issue on the [issues page](https://github.com/NREL/SolTrace/issues). diff --git a/coretrace/simulation_data/CMakeLists.txt b/coretrace/simulation_data/CMakeLists.txt index 8380658a..e317b9ee 100644 --- a/coretrace/simulation_data/CMakeLists.txt +++ b/coretrace/simulation_data/CMakeLists.txt @@ -25,14 +25,14 @@ set(SIMDATA_SRC stage_element.cpp sun.cpp surface.cpp - virtual_element.cpp + utilities.cpp vector3d.cpp + virtual_element.cpp cst_templates/arclength.cpp cst_templates/heliostat.cpp cst_templates/linear_fresnel.cpp cst_templates/parabolic_dish.cpp cst_templates/parabolic_trough.cpp - cst_templates/utilities.cpp solar_position_calculators/solar_position_calculator.cpp solar_position_calculators/lib_irradproc.cpp solar_position_calculators/solpos.cpp @@ -60,14 +60,14 @@ set(SIMDATA_HDRS stage_element.hpp sun.hpp surface.hpp + utilities.hpp vector3d.hpp - # virtual_element.hpp + virtual_element.hpp cst_templates/arclength.hpp cst_templates/heliostat.hpp cst_templates/linear_fresnel.hpp cst_templates/parabolic_dish.hpp cst_templates/parabolic_trough.hpp - cst_templates/utilities.hpp solar_position_calculators/solar_position_calculator.hpp solar_position_calculators/lib_irradproc.h solar_position_calculators/solpos00.h diff --git a/coretrace/simulation_data/aperture.cpp b/coretrace/simulation_data/aperture.cpp index 689aca2d..fbdf78da 100644 --- a/coretrace/simulation_data/aperture.cpp +++ b/coretrace/simulation_data/aperture.cpp @@ -10,781 +10,899 @@ #include "constants.hpp" -namespace SolTrace::Data { - -aperture_ptr Aperture::make_aperture_from_type(ApertureType type, - const std::vector &args) +namespace SolTrace::Data { - // std::cout << "Type: " << type - // << " NArgs: " << args.size() - // << std::endl; - switch (type) - { - case ApertureType::ANNULUS: - if (args.size() < 3) - break; - return make_aperture(args[0], args[1], args[2]); - case ApertureType::CIRCLE: - if (args.size() < 1) - break; - return make_aperture(args[0]); - case ApertureType::HEXAGON: - if (args.size() < 1) - break; - return make_aperture(args[0]); - case ApertureType::RECTANGLE: - if (args.size() < 2) - break; - return make_aperture(args[0], args[1]); // This is assuming centered around the origin - case ApertureType::EQUILATERAL_TRIANGLE: - if (args.size() < 1) - break; - return make_aperture(args[0]); - case ApertureType::SINGLE_AXIS_CURVATURE_SECTION: - if (args.size() < 3) - break; - return make_aperture( - args[1] - args[0], args[2], - -0.5 * (args[1] - args[0]), -0.5 * args[2]); - // return make_aperture(args[1] - args[0], args[2]); - case ApertureType::IRREGULAR_TRIANGLE: - if (args.size() < 6) - break; - return make_aperture( - args[0], args[1], args[2], args[3], args[4], args[5]); - case ApertureType::IRREGULAR_QUADRILATERAL: - if (args.size() < 8) - break; - return make_aperture( - args[0], args[1], args[2], args[3], - args[4], args[5], args[6], args[7]); - default: + + aperture_ptr Aperture::make_aperture_from_type(ApertureType type, + const std::vector &args) + { + // std::cout << "Type: " << type + // << " NArgs: " << args.size() + // << std::endl; + switch (type) + { + case ApertureType::ANNULUS: + if (args.size() < 3) + break; + return make_aperture(args[0], args[1], args[2]); + case ApertureType::CIRCLE: + if (args.size() < 1) + break; + return make_aperture(args[0]); + case ApertureType::HEXAGON: + if (args.size() < 1) + break; + return make_aperture(args[0]); + case ApertureType::RECTANGLE: + if (args.size() < 2) + break; + return make_aperture(args[0], args[1]); // This is assuming centered around the origin + case ApertureType::EQUILATERAL_TRIANGLE: + if (args.size() < 1) + break; + return make_aperture(args[0]); + case ApertureType::SINGLE_AXIS_CURVATURE_SECTION: + if (args.size() < 3) + break; + return make_aperture( + args[1] - args[0], args[2], + -0.5 * (args[1] - args[0]), -0.5 * args[2]); + // return make_aperture(args[1] - args[0], args[2]); + case ApertureType::IRREGULAR_TRIANGLE: + if (args.size() < 6) + break; + return make_aperture( + args[0], args[1], args[2], args[3], args[4], args[5]); + case ApertureType::IRREGULAR_QUADRILATERAL: + if (args.size() < 8) + break; + return make_aperture( + args[0], args[1], args[2], args[3], + args[4], args[5], args[6], args[7]); + default: + // TODO handle error + // Unsupported case + return nullptr; + // break; + } + // TODO handle error - // Unsupported case + // Wrong number of arguments + return nullptr; - // break; + // return aperture_ptr(); } - // TODO handle error - // Wrong number of arguments - - return nullptr; - // return aperture_ptr(); -} - -aperture_ptr Aperture::make_aperture_from_json(const nlohmann::ordered_json& jnode) -{ - if (!jnode.contains("aperture_type")) - throw std::invalid_argument("Missing aperture_type"); - std::string type_str = jnode.at("aperture_type"); - ApertureType aperture_type = get_enum_from_string(type_str, ApertureTypeMap, ApertureType::APERTURE_UNKNOWN); - switch (aperture_type) - { - case ApertureType::ANNULUS: return make_aperture(jnode); - case ApertureType::CIRCLE: return make_aperture(jnode); - case ApertureType::HEXAGON: return make_aperture(jnode); - case ApertureType::RECTANGLE: return make_aperture(jnode); - case ApertureType::EQUILATERAL_TRIANGLE: return make_aperture(jnode); - case ApertureType::IRREGULAR_TRIANGLE: return make_aperture(jnode); - case ApertureType::IRREGULAR_QUADRILATERAL:return make_aperture(jnode); + aperture_ptr Aperture::make_aperture_from_json(const nlohmann::ordered_json &jnode) + { + if (!jnode.contains("aperture_type")) + throw std::invalid_argument("Missing aperture_type"); + std::string type_str = jnode.at("aperture_type"); + ApertureType aperture_type = get_enum_from_string(type_str, ApertureTypeMap, ApertureType::APERTURE_UNKNOWN); + switch (aperture_type) + { + case ApertureType::ANNULUS: + return make_aperture(jnode); + case ApertureType::CIRCLE: + return make_aperture(jnode); + case ApertureType::HEXAGON: + return make_aperture(jnode); + case ApertureType::RECTANGLE: + return make_aperture(jnode); + case ApertureType::EQUILATERAL_TRIANGLE: + return make_aperture(jnode); + case ApertureType::IRREGULAR_TRIANGLE: + return make_aperture(jnode); + case ApertureType::IRREGULAR_QUADRILATERAL: + return make_aperture(jnode); default: throw std::invalid_argument("Unsupported aperture_type: " + type_str); + } } -} -Annulus::Annulus(const nlohmann::ordered_json& jnode) - : Aperture(ApertureType::ANNULUS) -{ - this->inner_radius = jnode.at("inner_radius"); - this->outer_radius = jnode.at("outer_radius"); - this->arc_angle = jnode.at("arc_angle"); -} - -Aperture::Point Aperture::midpoint(const Aperture::Point& v0, - const Aperture::Point& v1) const { - return Aperture::Point((v0.x + v1.x) / 2, (v0.y + v1.y) / 2); -} -std::vector Aperture::subdivide(Aperture::Triangle tri, - int n) const { - Point v0 = tri.a; - Point v1 = tri.b; - Point v2 = tri.c; - Point m01 = midpoint(v0, v1); - Point m12 = midpoint(v1, v2); - Point m20 = midpoint(v2, v0); - std::vector result; - if (n - 1 == 0) { - result.push_back(Triangle(v0, m01, m20)); - result.push_back(Triangle(m01, v1, m12)); - result.push_back(Triangle(m12, m01, m20)); - result.push_back(Triangle(m12, m20, v2)); + Aperture::Point Aperture::midpoint(const Aperture::Point &v0, + const Aperture::Point &v1) const + { + return Aperture::Point((v0.x + v1.x) / 2, (v0.y + v1.y) / 2); + } + + std::vector Aperture::subdivide(Aperture::Triangle tri, + int n) const + { + Point v0 = tri.a; + Point v1 = tri.b; + Point v2 = tri.c; + Point m01 = midpoint(v0, v1); + Point m12 = midpoint(v1, v2); + Point m20 = midpoint(v2, v0); + std::vector result; + if (n - 1 == 0) + { + result.push_back(Triangle(v0, m01, m20)); + result.push_back(Triangle(m01, v1, m12)); + result.push_back(Triangle(m12, m01, m20)); + result.push_back(Triangle(m12, m20, v2)); + return result; + } + auto t1 = subdivide(Triangle(v0, m01, m20), n - 1); + auto t2 = subdivide(Triangle(v0, m01, m20), n - 1); + auto t3 = subdivide(Triangle(v0, m01, m20), n - 1); + auto t4 = subdivide(Triangle(v0, m01, m20), n - 1); + result.insert(result.end(), t1.begin(), t1.end()); + result.insert(result.end(), t2.begin(), t2.end()); + result.insert(result.end(), t3.begin(), t3.end()); + result.insert(result.end(), t4.begin(), t4.end()); return result; } - auto t1 = subdivide(Triangle(v0, m01, m20), n - 1); - auto t2 = subdivide(Triangle(v0, m01, m20), n - 1); - auto t3 = subdivide(Triangle(v0, m01, m20), n - 1); - auto t4 = subdivide(Triangle(v0, m01, m20), n - 1); - result.insert(result.end(), t1.begin(), t1.end()); - result.insert(result.end(), t2.begin(), t2.end()); - result.insert(result.end(), t3.begin(), t3.end()); - result.insert(result.end(), t4.begin(), t4.end()); - return result; -} -int Aperture::index_of(std::vector& v, - const Aperture::Point& p) const { - auto it = find(v.begin(), v.end(), p); - if (it == v.end()) { - v.push_back(p); - return v.size() - 1; - } - return std::distance(v.begin(), it); -} -std::tuple, std::vector> Aperture::indexed_triangles( - const std::vector& triangles) const { - std::vector indices; - std::vector points; - std::vector flattened; - for (const Triangle& tri : triangles) { - indices.push_back(index_of(points, tri.a)); - indices.push_back(index_of(points, tri.b)); - indices.push_back(index_of(points, tri.c)); - } - for (const Point& p : points) { - flattened.push_back(p.x); - flattened.push_back(p.y); - } - return std::make_pair(flattened, indices); -} - -double Annulus::aperture_area() const -{ - // TODO: input.cpp on line 219 uses the formula - // elm->ParameterC*(ACOSM1O180)*(elm->ParameterB - elm->ParameterA); - // = \theta * (r - R) - // This seems to be wrong... - double R = this->outer_radius; - double r = this->inner_radius; - // Convert to radians - double arc = this->arc_angle * D2R; - return 0.5 * arc * (R * R - r * r); -} - -double Annulus::diameter_circumscribed_circle() const -{ - return 2.0 * this->outer_radius; -} -bool Annulus::is_in(double x, double y) const -{ - double r = sqrt(x * x + y * y); - bool inside = false; - if (this->inner_radius <= r && - r <= this->outer_radius) - { - double theta = atan2(y, x); - // Arc is split across x-axis, hence the 0.5 - double arc = 0.5 * this->arc_angle * D2R; - inside = (-arc <= theta && theta <= arc); - } - return inside; -} - -std::tuple, std::vector> -Annulus::triangulation() const { - const int resolution = 32; - std::vector verts; - std::vector indices; - for (int i = 0; i <= resolution; i++) { - const double u = i / resolution * PI * 2; - verts.push_back(inner_radius * std::cos(u)); - verts.push_back(inner_radius * std::sin(u)); - verts.push_back(outer_radius * std::cos(u)); - verts.push_back(outer_radius * std::sin(u)); - } - for (int i = 0; i < resolution - 3; i += 2) { - const int a = i; - const int b = i + 1; - const int c = i + 2; - const int d = i + 3; - // Generate two triangles for each quad in the mesh - // Adjust order to be counter-clockwise - indices.push_back(a); - indices.push_back(d); - indices.push_back(b); - indices.push_back(b); - indices.push_back(d); - indices.push_back(c); - } - return std::make_tuple(verts, indices); -} - -aperture_ptr Annulus::make_copy() const -{ - // Invokes the implicit copy constructor - return make_aperture(*this); -} + int Aperture::index_of(std::vector &v, + const Aperture::Point &p) const + { + auto it = find(v.begin(), v.end(), p); + if (it == v.end()) + { + v.push_back(p); + return v.size() - 1; + } + return std::distance(v.begin(), it); + } -void Annulus::write_json(nlohmann::ordered_json& jnode) const -{ - ApertureType type = ApertureType::ANNULUS; - jnode["aperture_type"] = ApertureTypeMap.at(type); - jnode["inner_radius"] = this->inner_radius; - jnode["outer_radius"] = this->outer_radius; - jnode["arc_angle"] = this->arc_angle; -} - -Circle::Circle(const nlohmann::ordered_json& jnode) - : Aperture(ApertureType::CIRCLE) -{ - this->diameter = jnode.at("diameter"); -} + std::tuple, std::vector> Aperture::indexed_triangles( + const std::vector &triangles) const + { + std::vector indices; + std::vector points; + std::vector flattened; + for (const Triangle &tri : triangles) + { + indices.push_back(index_of(points, tri.a)); + indices.push_back(index_of(points, tri.b)); + indices.push_back(index_of(points, tri.c)); + } + for (const Point &p : points) + { + flattened.push_back(p.x); + flattened.push_back(p.y); + } + return std::make_pair(flattened, indices); + } -double Circle::aperture_area() const -{ - return 0.25 * PI * this->diameter * this->diameter; -} + Annulus::Annulus(const nlohmann::ordered_json &jnode) + : Aperture(ApertureType::ANNULUS) + { + this->inner_radius = jnode.at("inner_radius"); + this->outer_radius = jnode.at("outer_radius"); + this->arc_angle = jnode.at("arc_angle"); + } -double Circle::diameter_circumscribed_circle() const -{ - return this->diameter; -} + double Annulus::aperture_area() const + { + // TODO: input.cpp on line 219 uses the formula + // elm->ParameterC*(ACOSM1O180)*(elm->ParameterB - elm->ParameterA); + // = \theta * (r - R) + // This seems to be wrong... + double R = this->outer_radius; + double r = this->inner_radius; + // Convert to radians + double arc = this->arc_angle * D2R; + return 0.5 * arc * (R * R - r * r); + } -bool Circle::is_in(double x, double y) const -{ - double r = sqrt(x * x + y * y); - return r <= this->radius_circumscribed_circle(); -} + double Annulus::diameter_circumscribed_circle() const + { + return 2.0 * this->outer_radius; + } -aperture_ptr Circle::make_copy() const -{ - // Invokes the implicit copy constructor - return make_aperture(*this); -} + void Annulus::bounding_box(double &xmin, + double &xmax, + double &ymin, + double &ymax) const + { + xmin = -this->outer_radius; + xmax = this->outer_radius; + ymin = -this->outer_radius; + ymax = this->outer_radius; + return; + } -void Circle::write_json(nlohmann::ordered_json& jnode) const -{ - ApertureType type = ApertureType::CIRCLE; - jnode["aperture_type"] = ApertureTypeMap.at(type); - jnode["diameter"] = this->diameter; -} + bool Annulus::is_in(double x, double y) const + { + double r = sqrt(x * x + y * y); + bool inside = false; + if (this->inner_radius <= r && + r <= this->outer_radius) + { + double theta = atan2(y, x); + // Arc is split across x-axis, hence the 0.5 + double arc = 0.5 * this->arc_angle * D2R; + inside = (-arc <= theta && theta <= arc); + } + return inside; + } -EqualateralTriangle::EqualateralTriangle(const nlohmann::ordered_json& jnode) - : Aperture(ApertureType::EQUILATERAL_TRIANGLE) -{ - this->circumscribe_diameter = jnode.at("circumscribe_diameter"); -} - -std::tuple, std::vector> -Circle::triangulation() const { - // Using a fixed Delaunay triangulation of the unit circle - std::vector verts = { - 1.00000000e+00, 0.00000000e+00, 9.51056516e-01, 3.09016994e-01, - 8.09016994e-01, 5.87785252e-01, 5.87785252e-01, 8.09016994e-01, - 3.09016994e-01, 9.51056516e-01, 6.12323400e-17, 1.00000000e+00, - -3.09016994e-01, 9.51056516e-01, -5.87785252e-01, 8.09016994e-01, - -8.09016994e-01, 5.87785252e-01, -9.51056516e-01, 3.09016994e-01, - -1.00000000e+00, 1.22464680e-16, -9.51056516e-01, -3.09016994e-01, - -8.09016994e-01, -5.87785252e-01, -5.87785252e-01, -8.09016994e-01, - -3.09016994e-01, -9.51056516e-01, -1.83697020e-16, -1.00000000e+00, - 3.09016994e-01, -9.51056516e-01, 5.87785252e-01, -8.09016994e-01, - 8.09016994e-01, -5.87785252e-01, 9.51056516e-01, -3.09016994e-01, - 0.00000000e+00, 0.00000000e+00, -5.00000000e-01, -5.00000000e-01, - -5.00000000e-01, 0.00000000e+00, -5.00000000e-01, 5.00000000e-01, - 0.00000000e+00, -5.00000000e-01, 0.00000000e+00, 5.00000000e-01, - 5.00000000e-01, -5.00000000e-01, 5.00000000e-01, 0.00000000e+00, - 5.00000000e-01, 5.00000000e-01 - }; - std::vector indices = { - 22, 11, 21, 11, 22, 10, 17, 18, 26, 14, 24, 21, 13, 14, 21, 24, 14, 15, - 25, 6, 23, 6, 25, 5, 9, 22, 23, 8, 9, 23, 22, 9, 10, 27, 1, 28, - 1, 27, 0, 19, 27, 26, 18, 19, 26, 27, 19, 0, 4, 25, 28, 3, 4, 28, - 25, 4, 5, 12, 13, 21, 11, 12, 21, 24, 16, 26, 16, 17, 26, 16, 24, 15, - 7, 8, 23, 6, 7, 23, 2, 3, 28, 1, 2, 28, 24, 22, 21, 22, 24, 20, - 24, 27, 20, 27, 24, 26, 25, 22, 20, 22, 25, 23, 27, 25, 20, 25, 27, 28 - }; - // scale from the unit cirle to our circle - std::transform( - verts.begin(), verts.end(), verts.begin(), [this](double element) { - return element *= this->diameter / 2.0; - }); - return std::make_tuple(verts, indices); -} - -double EqualateralTriangle::aperture_area() const -{ - double r = 0.5 * this->circumscribe_diameter; - return 0.75 * sqrt(3) * r * r; -} + std::tuple, std::vector> + Annulus::triangulation() const + { + const int resolution = 32; + std::vector verts; + std::vector indices; + for (int i = 0; i <= resolution; i++) + { + const double u = i / resolution * PI * 2; + verts.push_back(inner_radius * std::cos(u)); + verts.push_back(inner_radius * std::sin(u)); + verts.push_back(outer_radius * std::cos(u)); + verts.push_back(outer_radius * std::sin(u)); + } + for (int i = 0; i < resolution - 3; i += 2) + { + const int a = i; + const int b = i + 1; + const int c = i + 2; + const int d = i + 3; + // Generate two triangles for each quad in the mesh + // Adjust order to be counter-clockwise + indices.push_back(a); + indices.push_back(d); + indices.push_back(b); + indices.push_back(b); + indices.push_back(d); + indices.push_back(c); + } + return std::make_tuple(verts, indices); + } -double EqualateralTriangle::diameter_circumscribed_circle() const -{ - return this->circumscribe_diameter; -} + aperture_ptr Annulus::make_copy() const + { + // Invokes the implicit copy constructor + return make_aperture(*this); + } -bool EqualateralTriangle::is_in(double x, double y) const -{ - double r = sqrt(x * x + y * y); - double ro = this->radius_circumscribed_circle(); - if (r > ro) - return false; + void Annulus::write_json(nlohmann::ordered_json &jnode) const + { + ApertureType type = ApertureType::ANNULUS; + jnode["aperture_type"] = ApertureTypeMap.at(type); + jnode["inner_radius"] = this->inner_radius; + jnode["outer_radius"] = this->outer_radius; + jnode["arc_angle"] = this->arc_angle; + } - double ri = 0.5 * ro; - if (r <= ri) - return true; + Circle::Circle(const nlohmann::ordered_json &jnode) + : Aperture(ApertureType::CIRCLE) + { + this->diameter = jnode.at("diameter"); + } - double y0; - // double a = ro / sqrt(3.0) = 2 * ri / sqrt(3.0); - if (0.0 <= x && x <= ro) + double Circle::aperture_area() const { - // y0 = -sqrt(3.0) * (x - a); - y0 = ro - sqrt(3.0) * x; - return (-ri <= y && y <= y0); + return 0.25 * PI * this->diameter * this->diameter; } - else if (-ro <= x && x < 0.0) + + double Circle::diameter_circumscribed_circle() const { - // y0 = sqrt(3.0) * (x + a); - y0 = sqrt(3.0) * x + ro; - return (-ri <= y && y <= y0); + return this->diameter; } - return false; -} + void Circle::bounding_box(double &xmin, + double &xmax, + double &ymin, + double &ymax) const + { + double r = this->radius_circumscribed_circle(); + xmin = -r; + xmax = r; + ymin = -r; + ymax = r; + return; + } -std::tuple, std::vector> -EqualateralTriangle::triangulation() const { - double r = circumscribe_diameter / 2.0; - Triangle tri(Point(0, r), - Point(r * cos(-PI / 6.0), r * sin(-PI / 6.0)), - Point(r * cos(7 * PI / 6.0), r * sin(7 * PI / 6.0))); - return indexed_triangles(subdivide(tri, 3)); -} + bool Circle::is_in(double x, double y) const + { + double r = sqrt(x * x + y * y); + return r <= this->radius_circumscribed_circle(); + } -aperture_ptr EqualateralTriangle::make_copy() const -{ - // Invokes the implicit copy constructor - return make_aperture(*this); -} + aperture_ptr Circle::make_copy() const + { + // Invokes the implicit copy constructor + return make_aperture(*this); + } -void EqualateralTriangle::write_json(nlohmann::ordered_json& jnode) const -{ - ApertureType type = ApertureType::EQUILATERAL_TRIANGLE; - jnode["aperture_type"] = ApertureTypeMap.at(type); - jnode["circumscribe_diameter"] = this->circumscribe_diameter; -} + void Circle::write_json(nlohmann::ordered_json &jnode) const + { + ApertureType type = ApertureType::CIRCLE; + jnode["aperture_type"] = ApertureTypeMap.at(type); + jnode["diameter"] = this->diameter; + } -Hexagon::Hexagon(const nlohmann::ordered_json& jnode) - : Aperture(ApertureType::HEXAGON) -{ - this->circumscribe_diameter = jnode.at("circumscribe_diameter"); -} + std::tuple, std::vector> + Circle::triangulation() const + { + // Using a fixed Delaunay triangulation of the unit circle + std::vector verts = { + 1.00000000e+00, 0.00000000e+00, 9.51056516e-01, 3.09016994e-01, + 8.09016994e-01, 5.87785252e-01, 5.87785252e-01, 8.09016994e-01, + 3.09016994e-01, 9.51056516e-01, 6.12323400e-17, 1.00000000e+00, + -3.09016994e-01, 9.51056516e-01, -5.87785252e-01, 8.09016994e-01, + -8.09016994e-01, 5.87785252e-01, -9.51056516e-01, 3.09016994e-01, + -1.00000000e+00, 1.22464680e-16, -9.51056516e-01, -3.09016994e-01, + -8.09016994e-01, -5.87785252e-01, -5.87785252e-01, -8.09016994e-01, + -3.09016994e-01, -9.51056516e-01, -1.83697020e-16, -1.00000000e+00, + 3.09016994e-01, -9.51056516e-01, 5.87785252e-01, -8.09016994e-01, + 8.09016994e-01, -5.87785252e-01, 9.51056516e-01, -3.09016994e-01, + 0.00000000e+00, 0.00000000e+00, -5.00000000e-01, -5.00000000e-01, + -5.00000000e-01, 0.00000000e+00, -5.00000000e-01, 5.00000000e-01, + 0.00000000e+00, -5.00000000e-01, 0.00000000e+00, 5.00000000e-01, + 5.00000000e-01, -5.00000000e-01, 5.00000000e-01, 0.00000000e+00, + 5.00000000e-01, 5.00000000e-01}; + std::vector indices = { + 22, 11, 21, 11, 22, 10, 17, 18, 26, 14, 24, 21, 13, 14, 21, 24, 14, 15, + 25, 6, 23, 6, 25, 5, 9, 22, 23, 8, 9, 23, 22, 9, 10, 27, 1, 28, + 1, 27, 0, 19, 27, 26, 18, 19, 26, 27, 19, 0, 4, 25, 28, 3, 4, 28, + 25, 4, 5, 12, 13, 21, 11, 12, 21, 24, 16, 26, 16, 17, 26, 16, 24, 15, + 7, 8, 23, 6, 7, 23, 2, 3, 28, 1, 2, 28, 24, 22, 21, 22, 24, 20, + 24, 27, 20, 27, 24, 26, 25, 22, 20, 22, 25, 23, 27, 25, 20, 25, 27, 28}; + // scale from the unit cirle to our circle + std::transform( + verts.begin(), verts.end(), verts.begin(), [this](double element) + { return element *= this->diameter / 2.0; }); + return std::make_tuple(verts, indices); + } -double Hexagon::aperture_area() const -{ - // TODO: input.cpp on line 210 uses the formula - // 5*sqr(elm->ParameterA/2.0)*cos(30.0*(ACOSM1O180))*sin(30.0*(ACOSM1O180)); - // = 5*(d/2)^2*cos(pi/6)*sin(pi/6) - // = 5*(d/2)^2*sqrt(3)/2*1/2 - // = 5*sqrt(3)/4 * (d/2)^2 - // = 1.25*sqrt(3) * (d/2)^2 - // This seems to be wrong... - double r = 0.5 * this->circumscribe_diameter; - return 1.5 * sqrt(3) * r * r; -} - -double Hexagon::diameter_circumscribed_circle() const -{ - return circumscribe_diameter; -} + EqualateralTriangle::EqualateralTriangle(const nlohmann::ordered_json &jnode) + : Aperture(ApertureType::EQUILATERAL_TRIANGLE) + { + this->circumscribe_diameter = jnode.at("circumscribe_diameter"); + } + + double EqualateralTriangle::aperture_area() const + { + double r = 0.5 * this->circumscribe_diameter; + return 0.75 * sqrt(3) * r * r; + } + + double EqualateralTriangle::diameter_circumscribed_circle() const + { + return this->circumscribe_diameter; + } + + void EqualateralTriangle::bounding_box(double &xmin, + double &xmax, + double &ymin, + double &ymax) const + { + double r = this->radius_circumscribed_circle(); + double side = sqrt(3.0) * r; + double h = 1.5 * r; + xmin = -0.5 * side; + xmax = 0.5 * side; + ymin = r - h; + ymax = r; + return; + } + + bool EqualateralTriangle::is_in(double x, double y) const + { + double r = sqrt(x * x + y * y); + double ro = this->radius_circumscribed_circle(); + if (r > ro) + return false; + + double ri = 0.5 * ro; + if (r <= ri) + return true; + + double y0; + // double a = ro / sqrt(3.0) = 2 * ri / sqrt(3.0); + if (0.0 <= x && x <= ro) + { + // y0 = -sqrt(3.0) * (x - a); + y0 = ro - sqrt(3.0) * x; + return (-ri <= y && y <= y0); + } + else if (-ro <= x && x < 0.0) + { + // y0 = sqrt(3.0) * (x + a); + y0 = sqrt(3.0) * x + ro; + return (-ri <= y && y <= y0); + } -bool Hexagon::is_in(double x, double y) const -{ - double r = sqrt(x * x + y * y); - double ro = this->radius_circumscribed_circle(); - if (r > ro) return false; + } - double ri = 0.5 * sqrt(3.0) * ro; - if (r <= ri) - return true; - - // NOTE: Old code used - // xl = sqrt(ro^2 - ri^2) - // where `ro` is the radius of the circumscribing circle and `ri` is - // the radius of the inscribing circle. But this is equivalent to - // xl = 0.5 * ro - - double xl = 0.5 * this->radius_circumscribed_circle(); - double y1, y2; - if (xl < x && x <= ro) - { - y1 = sqrt(3.0) * (x - ro); - y2 = -y1; - // if (y1 <= y && y <= y2) return true; - return (y1 <= y && y <= y2); - } - else if (-xl <= x && x <= xl) - { - return (-ri <= y && y <= ri); - } - else if (-ro <= x && x < -xl) - { - y1 = sqrt(3.0) * (x + ro); - y2 = -y1; - return (y2 <= y && y <= y1); - } - - return false; -} - -std::tuple, std::vector> -Hexagon::triangulation() const { - double r = circumscribe_diameter / 2.0; - std::vector t0 = - subdivide(Triangle(Point(0, 0), - Point(r * cos(PI / 3.0), r * sin(PI / 3.0)), - Point(r, 0)), - 2); - std::vector t1 = subdivide( - Triangle(Point(0, 0), - Point(r * cos(2 * PI / 3.0), r * sin(2 * PI / 3.0)), - Point(r * cos(PI / 3.0), r * sin(PI / 3.0))), - 2); - std::vector t2 = subdivide( - Triangle(Point(0, 0), - Point(-r, 0), - Point(r * cos(2 * PI / 3.0), r * sin(2 * PI / 3.0))), - 2); - std::vector t3 = subdivide( - Triangle(Point(0, 0), - Point(-r, 0), - Point(r * cos(4 * PI / 3.0), r * sin(4 * PI / 3.0))), - 2); - std::vector t4 = subdivide( - Triangle(Point(0, 0), - Point(r * cos(5 * PI / 3.0), r * sin(5 * PI / 3.0)), - Point(r * cos(4 * PI / 3.0), r * sin(4 * PI / 3.0))), - 2); - std::vector t5 = subdivide( - Triangle(Point(0, 0), - Point(r, 0), - Point(r * cos(5 * PI / 3.0), r * sin(5 * PI / 3.0))), - 2); - std::vector triangles; - triangles.insert(triangles.end(), t0.begin(), t0.end()); - triangles.insert(triangles.end(), t1.begin(), t1.end()); - triangles.insert(triangles.end(), t2.begin(), t2.end()); - triangles.insert(triangles.end(), t3.begin(), t3.end()); - triangles.insert(triangles.end(), t4.begin(), t4.end()); - triangles.insert(triangles.end(), t5.begin(), t5.end()); - return indexed_triangles(triangles); -} - -aperture_ptr Hexagon::make_copy() const -{ - // Invokes the implicit copy constructor - return make_aperture(*this); -} + std::tuple, std::vector> + EqualateralTriangle::triangulation() const + { + double r = circumscribe_diameter / 2.0; + Triangle tri(Point(0, r), + Point(r * cos(-PI / 6.0), r * sin(-PI / 6.0)), + Point(r * cos(7 * PI / 6.0), r * sin(7 * PI / 6.0))); + return indexed_triangles(subdivide(tri, 3)); + } -void Hexagon::write_json(nlohmann::ordered_json& jnode) const -{ - ApertureType type = ApertureType::HEXAGON; - jnode["aperture_type"] = ApertureTypeMap.at(type); - jnode["circumscribe_diameter"] = this->circumscribe_diameter; -} - -IrregularTriangle::IrregularTriangle(double x1, double y1, - double x2, double y2, - double x3, double y3) - : Aperture(ApertureType::IRREGULAR_TRIANGLE), - x1(x1), y1(y1), - x2(x2), y2(y2), - x3(x3), y3(y3) -{ -} + aperture_ptr EqualateralTriangle::make_copy() const + { + // Invokes the implicit copy constructor + return make_aperture(*this); + } -IrregularTriangle::IrregularTriangle(const nlohmann::ordered_json& jnode) - : Aperture(ApertureType::IRREGULAR_TRIANGLE) -{ - this->x1 = jnode.at("x1"); - this->y1 = jnode.at("y1"); - this->x2 = jnode.at("x2"); - this->y2 = jnode.at("y2"); - this->x3 = jnode.at("x3"); - this->y3 = jnode.at("y3"); -} - -std::tuple, std::vector> -IrregularTriangle::triangulation() const { - Triangle tri(Point(x1, y1), Point(x2, y2), Point(x3, y3)); - return indexed_triangles(subdivide(tri, 3)); -} - -double IrregularTriangle::aperture_area() const -{ - double v11 = this->x1 - this->x2; - double v12 = this->y1 - this->y2; - double v21 = this->x3 - this->x2; - double v22 = this->y3 - this->y2; + void EqualateralTriangle::write_json(nlohmann::ordered_json &jnode) const + { + ApertureType type = ApertureType::EQUILATERAL_TRIANGLE; + jnode["aperture_type"] = ApertureTypeMap.at(type); + jnode["circumscribe_diameter"] = this->circumscribe_diameter; + } - double v1m = sqrt(v11 * v11 + v12 * v12); - double v2m = sqrt(v21 * v21 + v22 * v22); + Hexagon::Hexagon(const nlohmann::ordered_json &jnode) + : Aperture(ApertureType::HEXAGON) + { + this->circumscribe_diameter = jnode.at("circumscribe_diameter"); + } - double theta = acos((v11 * v21 + v12 * v22) / (v1m * v2m)); - double area = 0.5 * v1m * v2m * sin(theta); + double Hexagon::aperture_area() const + { + // TODO: input.cpp on line 210 uses the formula + // 5*sqr(elm->ParameterA/2.0)*cos(30.0*(ACOSM1O180))*sin(30.0*(ACOSM1O180)); + // = 5*(d/2)^2*cos(pi/6)*sin(pi/6) + // = 5*(d/2)^2*sqrt(3)/2*1/2 + // = 5*sqrt(3)/4 * (d/2)^2 + // = 1.25*sqrt(3) * (d/2)^2 + // This seems to be wrong... + double r = 0.5 * this->circumscribe_diameter; + return 1.5 * sqrt(3) * r * r; + } - return area; -} + double Hexagon::diameter_circumscribed_circle() const + { + return circumscribe_diameter; + } -double IrregularTriangle::diameter_circumscribed_circle() const -{ - // TODO: Not sure this is exact. Is that a problem? - double xmax = std::max(std::max(x1, x2), x3); - double ymax = std::max(std::max(y1, y2), y3); - double xmin = std::min(std::min(x1, x2), x3); - double ymin = std::min(std::min(y1, y2), y3); - double dx = xmax - xmin; - double dy = ymax - ymin; - return sqrt(dx * dx + dy * dy); -} - -bool IrregularTriangle::is_in(double x, double y) const -{ - return intri(x1, y1, x2, y2, x3, y3, x, y); -} + void Hexagon::bounding_box(double &xmin, + double &xmax, + double &ymin, + double &ymax) const + { + double r = this->radius_circumscribed_circle(); + double apothem = 0.5 * r * sqrt(3.0); + xmin = -r; + xmax = r; + ymin = -apothem; + ymax = apothem; + return; + } -aperture_ptr IrregularTriangle::make_copy() const -{ - // Invokes the implicit copy constructor - return make_aperture(*this); -} + bool Hexagon::is_in(double x, double y) const + { + double r = sqrt(x * x + y * y); + double ro = this->radius_circumscribed_circle(); + if (r > ro) + return false; + + double ri = 0.5 * sqrt(3.0) * ro; + if (r <= ri) + return true; + + // NOTE: Old code used + // xl = sqrt(ro^2 - ri^2) + // where `ro` is the radius of the circumscribing circle and `ri` is + // the radius of the inscribing circle. But this is equivalent to + // xl = 0.5 * ro + + double xl = 0.5 * this->radius_circumscribed_circle(); + double y1, y2; + if (xl < x && x <= ro) + { + y1 = sqrt(3.0) * (x - ro); + y2 = -y1; + // if (y1 <= y && y <= y2) return true; + return (y1 <= y && y <= y2); + } + else if (-xl <= x && x <= xl) + { + return (-ri <= y && y <= ri); + } + else if (-ro <= x && x < -xl) + { + y1 = sqrt(3.0) * (x + ro); + y2 = -y1; + return (y2 <= y && y <= y1); + } -void IrregularTriangle::write_json(nlohmann::ordered_json& jnode) const -{ - ApertureType type = ApertureType::IRREGULAR_TRIANGLE; - jnode["aperture_type"] = ApertureTypeMap.at(type); - jnode["x1"] = this->x1; - jnode["y1"] = this->y1; - jnode["x2"] = this->x2; - jnode["y2"] = this->y2; - jnode["x3"] = this->x3; - jnode["y3"] = this->y3; -} - -IrregularQuadrilateral::IrregularQuadrilateral(double x1, double y1, - double x2, double y2, - double x3, double y3, - double x4, double y4) - : Aperture(ApertureType::IRREGULAR_QUADRILATERAL), - x1(x1), y1(y1), - x2(x2), y2(y2), - x3(x3), y3(y3), - x4(x4), y4(y4) -{ -} + return false; + } -IrregularQuadrilateral::IrregularQuadrilateral(const nlohmann::ordered_json& jnode) - : Aperture(ApertureType::IRREGULAR_QUADRILATERAL) -{ - this->x1 = jnode.at("x1"); - this->y1 = jnode.at("y1"); - this->x2 = jnode.at("x2"); - this->y2 = jnode.at("y2"); - this->x3 = jnode.at("x3"); - this->y3 = jnode.at("y3"); - this->x4 = jnode.at("x4"); - this->y4 = jnode.at("y4"); -} - -double IrregularQuadrilateral::aperture_area() const -{ - double v11 = this->x1 - this->x2; - double v12 = this->y1 - this->y2; - double v21 = this->x3 - this->x2; - double v22 = this->y3 - this->y2; - double v31 = this->x3 - this->x4; - double v32 = this->y3 - this->y4; - double v41 = this->x1 - this->x4; - double v42 = this->y1 - this->y4; - - double v1m = sqrt(v11 * v11 + v12 * v12); - double v2m = sqrt(v21 * v21 + v22 * v22); - double v3m = sqrt(v31 * v31 + v32 * v32); - double v4m = sqrt(v41 * v41 + v42 * v42); - - double theta1 = acos((v11 * v21 + v12 * v22) / (v1m * v2m)); - double theta2 = acos((v31 * v41 + v32 * v42) / (v3m * v4m)); - - double area = 0.5 * (v1m * v2m * sin(theta1) + v3m * v4m * sin(theta2)); - return area; -} - -std::tuple, std::vector> -IrregularQuadrilateral::triangulation() const { - std::vector t0 = - subdivide(Triangle(Point(x1, y1), Point(x3, y3), Point(x2, y2)), 2); - std::vector t1 = - subdivide(Triangle(Point(x1, y1), Point(x4, y4), Point(x3, y3)), 2); - std::vector triangles; - triangles.insert(triangles.end(), t0.begin(), t0.end()); - triangles.insert(triangles.end(), t1.begin(), t1.end()); - return indexed_triangles(triangles); -} - -double IrregularQuadrilateral::diameter_circumscribed_circle() const -{ - // TODO: Not sure this is exact. Is that a problem? - double xmax = std::max(std::max(x1, x2), std::max(x3, x4)); - double ymax = std::max(std::max(y1, y2), std::max(y3, y4)); - double xmin = std::min(std::min(x1, x2), std::min(x3, x4)); - double ymin = std::min(std::min(y1, y2), std::min(y3, y4)); - double dx = xmax - xmin; - double dy = ymax - ymin; - return sqrt(dx * dx + dy * dy); -} - -bool IrregularQuadrilateral::is_in(double x, double y) const -{ - return inquad(x1, y1, x2, y2, x3, y3, x4, y4, x, y); -} + std::tuple, std::vector> + Hexagon::triangulation() const + { + double r = circumscribe_diameter / 2.0; + std::vector t0 = + subdivide(Triangle(Point(0, 0), + Point(r * cos(PI / 3.0), r * sin(PI / 3.0)), + Point(r, 0)), + 2); + std::vector t1 = subdivide( + Triangle(Point(0, 0), + Point(r * cos(2 * PI / 3.0), r * sin(2 * PI / 3.0)), + Point(r * cos(PI / 3.0), r * sin(PI / 3.0))), + 2); + std::vector t2 = subdivide( + Triangle(Point(0, 0), + Point(-r, 0), + Point(r * cos(2 * PI / 3.0), r * sin(2 * PI / 3.0))), + 2); + std::vector t3 = subdivide( + Triangle(Point(0, 0), + Point(-r, 0), + Point(r * cos(4 * PI / 3.0), r * sin(4 * PI / 3.0))), + 2); + std::vector t4 = subdivide( + Triangle(Point(0, 0), + Point(r * cos(5 * PI / 3.0), r * sin(5 * PI / 3.0)), + Point(r * cos(4 * PI / 3.0), r * sin(4 * PI / 3.0))), + 2); + std::vector t5 = subdivide( + Triangle(Point(0, 0), + Point(r, 0), + Point(r * cos(5 * PI / 3.0), r * sin(5 * PI / 3.0))), + 2); + std::vector triangles; + triangles.insert(triangles.end(), t0.begin(), t0.end()); + triangles.insert(triangles.end(), t1.begin(), t1.end()); + triangles.insert(triangles.end(), t2.begin(), t2.end()); + triangles.insert(triangles.end(), t3.begin(), t3.end()); + triangles.insert(triangles.end(), t4.begin(), t4.end()); + triangles.insert(triangles.end(), t5.begin(), t5.end()); + return indexed_triangles(triangles); + } -aperture_ptr IrregularQuadrilateral::make_copy() const -{ - // Invokes the implicit copy constructor - return make_aperture(*this); -} + aperture_ptr Hexagon::make_copy() const + { + // Invokes the implicit copy constructor + return make_aperture(*this); + } -void IrregularQuadrilateral::write_json(nlohmann::ordered_json& jnode) const -{ - ApertureType type = ApertureType::IRREGULAR_QUADRILATERAL; - jnode["aperture_type"] = ApertureTypeMap.at(type); - jnode["x1"] = this->x1; - jnode["y1"] = this->y1; - jnode["x2"] = this->x2; - jnode["y2"] = this->y2; - jnode["x3"] = this->x3; - jnode["y3"] = this->y3; - jnode["x4"] = this->x4; - jnode["y4"] = this->y4; -} - -Rectangle::Rectangle(double xlen, double ylen) - : Aperture(ApertureType::RECTANGLE), - x_length(xlen), - y_length(ylen) -{ - // Default to rectangle centered at the origin. - this->x_coord = -0.5 * this->x_length; - this->y_coord = -0.5 * this->y_length; - return; -} - -Rectangle::Rectangle(const nlohmann::ordered_json& jnode) - : Aperture(ApertureType::RECTANGLE) -{ - this->x_length = jnode.at("x_length"); - this->y_length = jnode.at("y_length"); - this->x_coord = jnode.at("x_coord"); - this->y_coord = jnode.at("y_coord"); -} + void Hexagon::write_json(nlohmann::ordered_json &jnode) const + { + ApertureType type = ApertureType::HEXAGON; + jnode["aperture_type"] = ApertureTypeMap.at(type); + jnode["circumscribe_diameter"] = this->circumscribe_diameter; + } -double Rectangle::aperture_area() const -{ - return this->x_length * this->y_length; -} + IrregularTriangle::IrregularTriangle(double x1, double y1, + double x2, double y2, + double x3, double y3) + : Aperture(ApertureType::IRREGULAR_TRIANGLE), + x1(x1), y1(y1), + x2(x2), y2(y2), + x3(x3), y3(y3) + { + } -double Rectangle::diameter_circumscribed_circle() const -{ - return sqrt(x_length * x_length + y_length * y_length); -} + IrregularTriangle::IrregularTriangle(const nlohmann::ordered_json &jnode) + : Aperture(ApertureType::IRREGULAR_TRIANGLE) + { + this->x1 = jnode.at("x1"); + this->y1 = jnode.at("y1"); + this->x2 = jnode.at("x2"); + this->y2 = jnode.at("y2"); + this->x3 = jnode.at("x3"); + this->y3 = jnode.at("y3"); + } -aperture_ptr Rectangle::make_copy() const -{ - // Invokes the implicit copy constructor - return make_aperture(*this); -} + std::tuple, std::vector> + IrregularTriangle::triangulation() const + { + Triangle tri(Point(x1, y1), Point(x2, y2), Point(x3, y3)); + return indexed_triangles(subdivide(tri, 3)); + } -bool Rectangle::is_in(double x, double y) const -{ - double xl = this->x_coord; - double yl = this->y_coord; - double xu = xl + this->x_length; - double yu = yl + this->y_length; - return (xl <= x && x <= xu && yl <= y && y <= yu); -} - -Rectangle::Rectangle(double xlen, double ylen, double xl, double yl) - : Aperture(ApertureType::RECTANGLE), - x_length(xlen), - y_length(ylen), - x_coord(xl), - y_coord(yl) -{ -} + double IrregularTriangle::aperture_area() const + { + double v11 = this->x1 - this->x2; + double v12 = this->y1 - this->y2; + double v21 = this->x3 - this->x2; + double v22 = this->y3 - this->y2; -void Rectangle::write_json(nlohmann::ordered_json& jnode) const -{ - ApertureType type = ApertureType::RECTANGLE; - jnode["aperture_type"] = ApertureTypeMap.at(type); - jnode["x_length"] = this->x_length; - jnode["y_length"] = this->y_length; - jnode["x_coord"] = this->x_coord; - jnode["y_coord"] = this->y_coord; -} - - -std::tuple, std::vector> -Rectangle::triangulation() const { - const int segments = 5; - std::vector verts; - std::vector indices; - for (int i = 0; i <= segments; ++i) { - for (int j = 0; j <= segments; ++j) { - const double x = i * x_length / segments + x_coord; - const double y = j * y_length / segments + y_coord; - verts.push_back(x); - verts.push_back(y); - } + double v1m = sqrt(v11 * v11 + v12 * v12); + double v2m = sqrt(v21 * v21 + v22 * v22); + + double theta = acos((v11 * v21 + v12 * v22) / (v1m * v2m)); + double area = 0.5 * v1m * v2m * sin(theta); + + return area; } - for (int i = 0; i < segments; ++i) { - for (int j = 0; j < segments; ++j) { - const int a = (segments + 1) * i + j; - const int c = (segments + 1) * (i + 1) + j; - const int d = (segments + 1) * (i + 1) + j + 1; - const int b = (segments + 1) * i + j + 1; - // Generate two triangles for each quad in the mesh - // Adjust order to be counter-clockwise - indices.push_back(a); - indices.push_back(c); - indices.push_back(b); - indices.push_back(b); - indices.push_back(c); - indices.push_back(d); + + double IrregularTriangle::diameter_circumscribed_circle() const + { + // TODO: Not sure this is exact. Is that a problem? + double xmax = std::max(std::max(x1, x2), x3); + double ymax = std::max(std::max(y1, y2), y3); + double xmin = std::min(std::min(x1, x2), x3); + double ymin = std::min(std::min(y1, y2), y3); + double dx = xmax - xmin; + double dy = ymax - ymin; + return sqrt(dx * dx + dy * dy); + } + + void IrregularTriangle::bounding_box(double &xmin, + double &xmax, + double &ymin, + double &ymax) const + { + xmin = std::min(x1, std::min(x2, x3)); + xmax = std::max(x1, std::max(x2, x3)); + ymin = std::min(y1, std::min(y2, y3)); + ymax = std::max(y1, std::max(y2, y3)); + return; + } + + bool IrregularTriangle::is_in(double x, double y) const + { + return intri(x1, y1, x2, y2, x3, y3, x, y); + } + + aperture_ptr IrregularTriangle::make_copy() const + { + // Invokes the implicit copy constructor + return make_aperture(*this); + } + + void IrregularTriangle::write_json(nlohmann::ordered_json &jnode) const + { + ApertureType type = ApertureType::IRREGULAR_TRIANGLE; + jnode["aperture_type"] = ApertureTypeMap.at(type); + jnode["x1"] = this->x1; + jnode["y1"] = this->y1; + jnode["x2"] = this->x2; + jnode["y2"] = this->y2; + jnode["x3"] = this->x3; + jnode["y3"] = this->y3; + } + + IrregularQuadrilateral::IrregularQuadrilateral(double x1, double y1, + double x2, double y2, + double x3, double y3, + double x4, double y4) + : Aperture(ApertureType::IRREGULAR_QUADRILATERAL), + x1(x1), y1(y1), + x2(x2), y2(y2), + x3(x3), y3(y3), + x4(x4), y4(y4) + { + } + + IrregularQuadrilateral::IrregularQuadrilateral(const nlohmann::ordered_json &jnode) + : Aperture(ApertureType::IRREGULAR_QUADRILATERAL) + { + this->x1 = jnode.at("x1"); + this->y1 = jnode.at("y1"); + this->x2 = jnode.at("x2"); + this->y2 = jnode.at("y2"); + this->x3 = jnode.at("x3"); + this->y3 = jnode.at("y3"); + this->x4 = jnode.at("x4"); + this->y4 = jnode.at("y4"); + } + + double IrregularQuadrilateral::aperture_area() const + { + double v11 = this->x1 - this->x2; + double v12 = this->y1 - this->y2; + double v21 = this->x3 - this->x2; + double v22 = this->y3 - this->y2; + double v31 = this->x3 - this->x4; + double v32 = this->y3 - this->y4; + double v41 = this->x1 - this->x4; + double v42 = this->y1 - this->y4; + + double v1m = sqrt(v11 * v11 + v12 * v12); + double v2m = sqrt(v21 * v21 + v22 * v22); + double v3m = sqrt(v31 * v31 + v32 * v32); + double v4m = sqrt(v41 * v41 + v42 * v42); + + double theta1 = acos((v11 * v21 + v12 * v22) / (v1m * v2m)); + double theta2 = acos((v31 * v41 + v32 * v42) / (v3m * v4m)); + + double area = 0.5 * (v1m * v2m * sin(theta1) + v3m * v4m * sin(theta2)); + return area; + } + + std::tuple, std::vector> + IrregularQuadrilateral::triangulation() const + { + std::vector t0 = + subdivide(Triangle(Point(x1, y1), Point(x3, y3), Point(x2, y2)), 2); + std::vector t1 = + subdivide(Triangle(Point(x1, y1), Point(x4, y4), Point(x3, y3)), 2); + std::vector triangles; + triangles.insert(triangles.end(), t0.begin(), t0.end()); + triangles.insert(triangles.end(), t1.begin(), t1.end()); + return indexed_triangles(triangles); + } + + double IrregularQuadrilateral::diameter_circumscribed_circle() const + { + // TODO: Not sure this is exact. Is that a problem? + double xmax = std::max(std::max(x1, x2), std::max(x3, x4)); + double ymax = std::max(std::max(y1, y2), std::max(y3, y4)); + double xmin = std::min(std::min(x1, x2), std::min(x3, x4)); + double ymin = std::min(std::min(y1, y2), std::min(y3, y4)); + double dx = xmax - xmin; + double dy = ymax - ymin; + return sqrt(dx * dx + dy * dy); + } + + void IrregularQuadrilateral::bounding_box(double &xmin, + double &xmax, + double &ymin, + double &ymax) const + { + xmin = std::min(std::min(x1, x2), std::min(x3, x4)); + xmax = std::max(std::max(x1, x2), std::max(x3, x4)); + ymin = std::min(std::min(y1, y2), std::min(y3, y4)); + ymax = std::max(std::max(y1, y2), std::max(y3, y4)); + return; + } + + bool IrregularQuadrilateral::is_in(double x, double y) const + { + return inquad(x1, y1, x2, y2, x3, y3, x4, y4, x, y); + } + + aperture_ptr IrregularQuadrilateral::make_copy() const + { + // Invokes the implicit copy constructor + return make_aperture(*this); + } + + void IrregularQuadrilateral::write_json(nlohmann::ordered_json &jnode) const + { + ApertureType type = ApertureType::IRREGULAR_QUADRILATERAL; + jnode["aperture_type"] = ApertureTypeMap.at(type); + jnode["x1"] = this->x1; + jnode["y1"] = this->y1; + jnode["x2"] = this->x2; + jnode["y2"] = this->y2; + jnode["x3"] = this->x3; + jnode["y3"] = this->y3; + jnode["x4"] = this->x4; + jnode["y4"] = this->y4; + } + + Rectangle::Rectangle(double xlen, double ylen) + : Aperture(ApertureType::RECTANGLE), + x_length(xlen), + y_length(ylen) + { + // Default to rectangle centered at the origin. + this->x_coord = -0.5 * this->x_length; + this->y_coord = -0.5 * this->y_length; + return; + } + + Rectangle::Rectangle(double xlen, double ylen, double xl, double yl) + : Aperture(ApertureType::RECTANGLE), + x_length(xlen), + y_length(ylen), + x_coord(xl), + y_coord(yl) + { + } + + Rectangle::Rectangle(const nlohmann::ordered_json &jnode) + : Aperture(ApertureType::RECTANGLE) + { + this->x_length = jnode.at("x_length"); + this->y_length = jnode.at("y_length"); + this->x_coord = jnode.at("x_coord"); + this->y_coord = jnode.at("y_coord"); + } + + double Rectangle::aperture_area() const + { + return this->x_length * this->y_length; + } + + double Rectangle::diameter_circumscribed_circle() const + { + return sqrt(x_length * x_length + y_length * y_length); + } + + void Rectangle::bounding_box(double &xmin, + double &xmax, + double &ymin, + double &ymax) const + { + xmin = this->x_coord; + xmax = this->x_coord + this->x_length; + ymin = this->y_coord; + ymax = this->y_coord + this->y_length; + return; + } + + bool Rectangle::is_in(double x, double y) const + { + double xl = this->x_coord; + double yl = this->y_coord; + double xu = xl + this->x_length; + double yu = yl + this->y_length; + return (xl <= x && x <= xu && yl <= y && y <= yu); + } + + aperture_ptr Rectangle::make_copy() const + { + // Invokes the implicit copy constructor + return make_aperture(*this); + } + + void Rectangle::write_json(nlohmann::ordered_json &jnode) const + { + ApertureType type = ApertureType::RECTANGLE; + jnode["aperture_type"] = ApertureTypeMap.at(type); + jnode["x_length"] = this->x_length; + jnode["y_length"] = this->y_length; + jnode["x_coord"] = this->x_coord; + jnode["y_coord"] = this->y_coord; + } + + std::tuple, std::vector> + Rectangle::triangulation() const + { + const int segments = 5; + std::vector verts; + std::vector indices; + for (int i = 0; i <= segments; ++i) + { + for (int j = 0; j <= segments; ++j) + { + const double x = i * x_length / segments + x_coord; + const double y = j * y_length / segments + y_coord; + verts.push_back(x); + verts.push_back(y); + } + } + for (int i = 0; i < segments; ++i) + { + for (int j = 0; j < segments; ++j) + { + const int a = (segments + 1) * i + j; + const int c = (segments + 1) * (i + 1) + j; + const int d = (segments + 1) * (i + 1) + j + 1; + const int b = (segments + 1) * i + j + 1; + // Generate two triangles for each quad in the mesh + // Adjust order to be counter-clockwise + indices.push_back(a); + indices.push_back(c); + indices.push_back(b); + indices.push_back(b); + indices.push_back(c); + indices.push_back(d); + } } + return make_pair(verts, indices); } - return make_pair(verts, indices); -} -bool intri(double x1, double y1, - double x2, double y2, - double x3, double y3, - double xt, double yt) -{ - double a = (x1 - xt) * (y2 - yt) - (x2 - xt) * (y1 - yt); - double b = (x2 - xt) * (y3 - yt) - (x3 - xt) * (y2 - yt); - double c = (x3 - xt) * (y1 - yt) - (x1 - xt) * (y3 - yt); - return (std::signbit(a) == std::signbit(b) && - std::signbit(b) == std::signbit(c)); - // return (sign(a) == sign(b) && sign(b) == sign(c)); -} - -bool inquad(double x1, double y1, - double x2, double y2, - double x3, double y3, - double x4, double y4, - double xt, double yt) -{ - return (intri(x1, y1, x2, y2, x3, y3, xt, yt) || - intri(x1, y1, x3, y3, x4, y4, xt, yt)); -} + bool intri(double x1, double y1, + double x2, double y2, + double x3, double y3, + double xt, double yt) + { + double a = (x1 - xt) * (y2 - yt) - (x2 - xt) * (y1 - yt); + double b = (x2 - xt) * (y3 - yt) - (x3 - xt) * (y2 - yt); + double c = (x3 - xt) * (y1 - yt) - (x1 - xt) * (y3 - yt); + return (std::signbit(a) == std::signbit(b) && + std::signbit(b) == std::signbit(c)); + // return (sign(a) == sign(b) && sign(b) == sign(c)); + } + + bool inquad(double x1, double y1, + double x2, double y2, + double x3, double y3, + double x4, double y4, + double xt, double yt) + { + return (intri(x1, y1, x2, y2, x3, y3, xt, yt) || + intri(x1, y1, x3, y3, x4, y4, xt, yt)); + } } // namespace SolTrace::Data diff --git a/coretrace/simulation_data/aperture.hpp b/coretrace/simulation_data/aperture.hpp index d4a961be..625baade 100644 --- a/coretrace/simulation_data/aperture.hpp +++ b/coretrace/simulation_data/aperture.hpp @@ -42,17 +42,16 @@ namespace SolTrace::Data }; inline const std::map ApertureTypeMap = - { - {ApertureType::ANNULUS, "ANNULUS"}, - {ApertureType::CIRCLE, "CIRCLE"}, - {ApertureType::HEXAGON, "HEXAGON"}, - {ApertureType::RECTANGLE, "RECTANGLE"}, - {ApertureType::EQUILATERAL_TRIANGLE, "EQUILATERAL_TRIANGLE"}, - {ApertureType::SINGLE_AXIS_CURVATURE_SECTION, "SINGLE_AXIS_CURVATURE_SECTION"}, - {ApertureType::IRREGULAR_TRIANGLE, "IRREGULAR_TRIANGLE"}, - {ApertureType::IRREGULAR_QUADRILATERAL, "IRREGULAR_QUADRILATERAL"}, - {ApertureType::APERTURE_UNKNOWN, "APERTURE_UNKNOWN"} - }; + { + {ApertureType::ANNULUS, "ANNULUS"}, + {ApertureType::CIRCLE, "CIRCLE"}, + {ApertureType::HEXAGON, "HEXAGON"}, + {ApertureType::RECTANGLE, "RECTANGLE"}, + {ApertureType::EQUILATERAL_TRIANGLE, "EQUILATERAL_TRIANGLE"}, + {ApertureType::SINGLE_AXIS_CURVATURE_SECTION, "SINGLE_AXIS_CURVATURE_SECTION"}, + {ApertureType::IRREGULAR_TRIANGLE, "IRREGULAR_TRIANGLE"}, + {ApertureType::IRREGULAR_QUADRILATERAL, "IRREGULAR_QUADRILATERAL"}, + {ApertureType::APERTURE_UNKNOWN, "APERTURE_UNKNOWN"}}; struct Aperture; using aperture_ptr = std::shared_ptr; @@ -93,11 +92,11 @@ namespace SolTrace::Data const std::vector &args); /** - * @brief Factory method to create apertures from json - * @param jnode the json containing necessary parameters for each aperture type - * @return Shared pointer to the created aperture - */ - static aperture_ptr make_aperture_from_json(const nlohmann::ordered_json& jnode); + * @brief Factory method to create apertures from json + * @param jnode the json containing necessary parameters for each aperture type + * @return Shared pointer to the created aperture + */ + static aperture_ptr make_aperture_from_json(const nlohmann::ordered_json &jnode); /** * @brief Get the aperture type @@ -150,53 +149,70 @@ namespace SolTrace::Data */ virtual aperture_ptr make_copy() const = 0; + virtual void bounding_box(double &xmin, double &xmax, + double &ymin, double &ymax) const = 0; + /** * @brief Write aperture parameters to json * @param jnode JSON node */ - virtual void write_json(nlohmann::ordered_json& jnode) const = 0; + virtual void write_json(nlohmann::ordered_json &jnode) const = 0; /** * @brief Get the aperture type string * @return The aperture type string */ - inline std::string get_type_string() const { - switch (my_type) { - case ANNULUS: return "Annulus"; - case CIRCLE: return "Circle"; - case HEXAGON: return "Hexagon"; - case RECTANGLE: return "Rectangle"; - case EQUILATERAL_TRIANGLE: return "Regular Triangle"; - case SINGLE_AXIS_CURVATURE_SECTION: return "Single Axis Curvature"; - case IRREGULAR_TRIANGLE: return "Triangle"; - case IRREGULAR_QUADRILATERAL: return "Quad"; - case APERTURE_UNKNOWN: return "Unknown"; + inline std::string get_type_string() const + { + switch (my_type) + { + case ANNULUS: + return "Annulus"; + case CIRCLE: + return "Circle"; + case HEXAGON: + return "Hexagon"; + case RECTANGLE: + return "Rectangle"; + case EQUILATERAL_TRIANGLE: + return "Regular Triangle"; + case SINGLE_AXIS_CURVATURE_SECTION: + return "Single Axis Curvature"; + case IRREGULAR_TRIANGLE: + return "Triangle"; + case IRREGULAR_QUADRILATERAL: + return "Quad"; + case APERTURE_UNKNOWN: + return "Unknown"; } return "Unknown"; } protected: - struct Point { + struct Point + { public: double x; double y; - Point(double ix, double iy) : x(ix), y(iy) { } - bool operator==(const Point& p) const { + Point(double ix, double iy) : x(ix), y(iy) {} + bool operator==(const Point &p) const + { return x == p.x && y == p.y; } }; - struct Triangle { + struct Triangle + { public: Point a; Point b; Point c; - Triangle(Point ia, Point ib, Point ic) : a(ia), b(ib), c(ic) { } + Triangle(Point ia, Point ib, Point ic) : a(ia), b(ib), c(ic) {} }; /** * @brief Compute the mipoint between to points * @return The midpoint */ - Point midpoint(const Point& v0, const Point& v1) const; + Point midpoint(const Point &v0, const Point &v1) const; /** * @brief Recursively Subdivide a triangle by midpoints * @param tri The triangle to subdivide @@ -211,18 +227,18 @@ namespace SolTrace::Data * @param p The point of interest * @return index of p in v */ - int index_of(std::vector& v, const Point& p) const; + int index_of(std::vector &v, const Point &p) const; /** * @brief Convert a list of Triangles in indexed (flattened) faceset * @param triangles The list to convert * @return indexed faceset */ std::tuple, std::vector> - indexed_triangles(const std::vector& triangles) const; + indexed_triangles(const std::vector &triangles) const; }; - - struct Annulus : public Aperture { + struct Annulus : public Aperture + { double inner_radius; double outer_radius; double arc_angle; @@ -251,7 +267,7 @@ namespace SolTrace::Data * @brief Json-based constructor for annulus aperture * @param jnode contains ri, ro, and arc angle */ - Annulus(const nlohmann::ordered_json& jnode); + Annulus(const nlohmann::ordered_json &jnode); virtual ~Annulus() {} @@ -267,6 +283,11 @@ namespace SolTrace::Data */ virtual double diameter_circumscribed_circle() const override; + virtual void bounding_box(double &xmin, + double &xmax, + double &ymin, + double &ymax) const override; + /** * @brief Test if point is inside annulus aperture * @param x X coordinate @@ -285,7 +306,7 @@ namespace SolTrace::Data * @brief Write aperture parameters to json * @param jnode JSON node */ - virtual void write_json(nlohmann::ordered_json& jnode) const override; + virtual void write_json(nlohmann::ordered_json &jnode) const override; /** * @brief Triangulate the annulus shape @@ -311,7 +332,7 @@ namespace SolTrace::Data * @brief Json-based constructor for circular aperture * @param jnode contains diameter */ - Circle(const nlohmann::ordered_json& jnode); + Circle(const nlohmann::ordered_json &jnode); virtual ~Circle() {} @@ -327,6 +348,11 @@ namespace SolTrace::Data */ virtual double diameter_circumscribed_circle() const override; + virtual void bounding_box(double &xmin, + double &xmax, + double &ymin, + double &ymax) const override; + /** * @brief Test if point is inside circular aperture * @param x X coordinate @@ -345,7 +371,7 @@ namespace SolTrace::Data * @brief Write aperture parameters to json * @param jnode JSON node */ - virtual void write_json(nlohmann::ordered_json& jnode) const override; + virtual void write_json(nlohmann::ordered_json &jnode) const override; /** * @brief Triangulate the circle shape @@ -377,7 +403,7 @@ namespace SolTrace::Data * @brief Json-based constructor for equilateral triangle aperture * @param jnode contains diameter of circumscribed circle */ - EqualateralTriangle(const nlohmann::ordered_json& jnode); + EqualateralTriangle(const nlohmann::ordered_json &jnode); virtual ~EqualateralTriangle() {} @@ -393,6 +419,11 @@ namespace SolTrace::Data */ virtual double diameter_circumscribed_circle() const override; + virtual void bounding_box(double &xmin, + double &xmax, + double &ymin, + double &ymax) const override; + /** * @brief Test if point is inside triangular aperture * @param x X coordinate @@ -411,7 +442,7 @@ namespace SolTrace::Data * @brief Write aperture parameters to json * @param jnode JSON node */ - virtual void write_json(nlohmann::ordered_json& jnode) const override; + virtual void write_json(nlohmann::ordered_json &jnode) const override; /** * @brief Triangulate the triangle shape @@ -439,7 +470,7 @@ namespace SolTrace::Data * @brief Json-based constructor for hexagonal aperture * @param jnode contains diameter of circumscribed circle */ - Hexagon(const nlohmann::ordered_json& jnode); + Hexagon(const nlohmann::ordered_json &jnode); virtual ~Hexagon() {} @@ -455,6 +486,11 @@ namespace SolTrace::Data */ virtual double diameter_circumscribed_circle() const override; + virtual void bounding_box(double &xmin, + double &xmax, + double &ymin, + double &ymax) const override; + /** * @brief Test if point is inside hexagonal aperture * @param x X coordinate @@ -473,7 +509,7 @@ namespace SolTrace::Data * @brief Write aperture parameters to json * @param jnode JSON node */ - virtual void write_json(nlohmann::ordered_json& jnode) const override; + virtual void write_json(nlohmann::ordered_json &jnode) const override; /** * @brief Triangulate the hexagon shape @@ -512,7 +548,7 @@ namespace SolTrace::Data * @brief Json-based constructor for rectangular aperture * @param jnode contains xlen, ylen, xl, yl */ - Rectangle(const nlohmann::ordered_json& jnode); + Rectangle(const nlohmann::ordered_json &jnode); virtual ~Rectangle() {} @@ -528,6 +564,11 @@ namespace SolTrace::Data */ virtual double diameter_circumscribed_circle() const override; + virtual void bounding_box(double &xmin, + double &xmax, + double &ymin, + double &ymax) const override; + /** * @brief Test if point is inside rectangular aperture * @param x X coordinate @@ -546,7 +587,7 @@ namespace SolTrace::Data * @brief Write aperture parameters to json * @param jnode JSON node */ - virtual void write_json(nlohmann::ordered_json& jnode) const override; + virtual void write_json(nlohmann::ordered_json &jnode) const override; /** * @brief Triangulate the rectangle shape @@ -556,10 +597,10 @@ namespace SolTrace::Data triangulation() const override; }; - struct SingleAxisCurvatureSection : public Aperture - { - // TODO: Implement this? - }; + // struct SingleAxisCurvatureSection : public Aperture + // { + // // TODO: Implement this? + // }; struct IrregularTriangle : public Aperture { @@ -588,7 +629,7 @@ namespace SolTrace::Data * @brief Json-based constructor for irregular triangle aperture * @param jnode contains x1, y1, x2, y2, x3, y3 */ - IrregularTriangle(const nlohmann::ordered_json& jnode); + IrregularTriangle(const nlohmann::ordered_json &jnode); ~IrregularTriangle() {} @@ -604,6 +645,11 @@ namespace SolTrace::Data */ virtual double diameter_circumscribed_circle() const override; + virtual void bounding_box(double &xmin, + double &xmax, + double &ymin, + double &ymax) const override; + /** * @brief Test if point is inside irregular triangle * @param x X coordinate @@ -622,7 +668,7 @@ namespace SolTrace::Data * @brief Write aperture parameters to json * @param jnode JSON node */ - virtual void write_json(nlohmann::ordered_json& jnode) const override; + virtual void write_json(nlohmann::ordered_json &jnode) const override; /** * @brief Triangulate the triangle shape @@ -664,7 +710,7 @@ namespace SolTrace::Data * @brief Json-based constructor for irregular quadrilateral aperture * @param jnode contains x1, y1, x2, y2, x3, y3, x4, y4 */ - IrregularQuadrilateral(const nlohmann::ordered_json& jnode); + IrregularQuadrilateral(const nlohmann::ordered_json &jnode); ~IrregularQuadrilateral() {} @@ -680,6 +726,11 @@ namespace SolTrace::Data */ virtual double diameter_circumscribed_circle() const override; + virtual void bounding_box(double &xmin, + double &xmax, + double &ymin, + double &ymax) const override; + /** * @brief Test if point is inside irregular quadrilateral * @param x X coordinate @@ -698,7 +749,7 @@ namespace SolTrace::Data * @brief Write aperture parameters to json * @param jnode JSON node */ - virtual void write_json(nlohmann::ordered_json& jnode) const override; + virtual void write_json(nlohmann::ordered_json &jnode) const override; /** * @brief Triangulate the quad shape diff --git a/coretrace/simulation_data/composite_element.cpp b/coretrace/simulation_data/composite_element.cpp index 084fe814..6adcab20 100644 --- a/coretrace/simulation_data/composite_element.cpp +++ b/coretrace/simulation_data/composite_element.cpp @@ -5,8 +5,6 @@ #include "single_element.hpp" - - namespace SolTrace::Data { @@ -17,16 +15,16 @@ namespace SolTrace::Data return; } - CompositeElement::CompositeElement(const nlohmann::ordered_json& jnode) : ElementBase(jnode), - number_of_elements(0), - my_elements() + CompositeElement::CompositeElement(const nlohmann::ordered_json &jnode) : ElementBase(jnode), + number_of_elements(0), + my_elements() { using json = nlohmann::ordered_json; // Common parameters are set in ElementBase constructor before here json jelements = jnode.at("elements"); - for (auto& [key, jelement] : jelements.items()) + for (auto &[key, jelement] : jelements.items()) { bool is_single = jelement.at("is_single"); if (is_single) @@ -127,6 +125,11 @@ namespace SolTrace::Data { this->number_of_elements += el->get_number_of_elements(); el->set_reference_element(this); + // Mark any added elements as virtual if needed + if (this->is_virtual()) + { + el->mark_virtual(); + } } return id; } @@ -187,7 +190,7 @@ namespace SolTrace::Data } // Stage and Composite should have the same function - void CompositeElement::write_json(nlohmann::ordered_json& jnode) const + void CompositeElement::write_json(nlohmann::ordered_json &jnode) const { using json = nlohmann::ordered_json; diff --git a/coretrace/simulation_data/cst_templates/linear_fresnel.cpp b/coretrace/simulation_data/cst_templates/linear_fresnel.cpp index 85636c79..aef1a126 100644 --- a/coretrace/simulation_data/cst_templates/linear_fresnel.cpp +++ b/coretrace/simulation_data/cst_templates/linear_fresnel.cpp @@ -7,8 +7,7 @@ #include "constants.hpp" #include "element.hpp" #include "surface.hpp" - -#include "cst_templates/utilities.hpp" +#include "utilities.hpp" namespace SolTrace::Data { diff --git a/coretrace/simulation_data/matvec.hpp b/coretrace/simulation_data/matvec.hpp index 1562a0ee..e6b7aca8 100644 --- a/coretrace/simulation_data/matvec.hpp +++ b/coretrace/simulation_data/matvec.hpp @@ -70,6 +70,25 @@ namespace SolTrace::Data void MatrixVectorMult(const double M[3][3], const double V[3], double MxV[3]); + + template + inline void MatrixVectorMult_generic(const R M[3][3], + const T V[3], + T MxV[3]) + { + MxV[0] = M[0][0] * V[0] + M[0][1] * V[1] + M[0][2] * V[2]; + MxV[1] = M[1][0] * V[0] + M[1][1] * V[1] + M[1][2] * V[2]; + MxV[2] = M[2][0] * V[0] + M[2][1] * V[1] + M[2][2] * V[2]; + return; + } + + // void MatrixVectorMult(const double M[3][3], + // const double V[3], + // double MxV[3]) + // { + // return MatrixVectorMult_generic(M, V, MxV); + // } + void MatrixTranspose(const double InputMatrix[3][3], int NumRowsCols, double OutputMatrix[3][3]); diff --git a/coretrace/simulation_data/simdata_io.cpp b/coretrace/simulation_data/simdata_io.cpp index da3a1bf5..21baa72e 100644 --- a/coretrace/simulation_data/simdata_io.cpp +++ b/coretrace/simulation_data/simdata_io.cpp @@ -583,6 +583,14 @@ bool process_stages( stage->set_aim_vector(AX, AY, AZ); stage->set_zrot(ZRot); stage->compute_coordinate_rotations(); + if (virt) + { + stage->mark_virtual(); + } + else + { + stage->unmark_virtual(); + } // Loop through elements for (int i_element = 0; i_element < count_element; i_element++) diff --git a/coretrace/simulation_data/simulation_data_api.hpp b/coretrace/simulation_data/simulation_data_api.hpp index 1130e746..de8a9dc1 100644 --- a/coretrace/simulation_data/simulation_data_api.hpp +++ b/coretrace/simulation_data/simulation_data_api.hpp @@ -11,6 +11,7 @@ #include "stage_element.hpp" #include "sun.hpp" #include "surface.hpp" +#include "utilities.hpp" #include "vector3d.hpp" #include "virtual_element.hpp" diff --git a/coretrace/simulation_data/simulation_data_export.hpp b/coretrace/simulation_data/simulation_data_export.hpp index a6bd03db..fdd319d4 100644 --- a/coretrace/simulation_data/simulation_data_export.hpp +++ b/coretrace/simulation_data/simulation_data_export.hpp @@ -9,6 +9,7 @@ using SolTrace::Data::Aperture; using SolTrace::Data::ApertureType; using SolTrace::Data::Circle; using SolTrace::Data::CompositeElement; +using SolTrace::Data::Cone; using SolTrace::Data::Cylinder; using SolTrace::Data::DistributionType; using SolTrace::Data::EqualateralTriangle; @@ -35,7 +36,7 @@ using SolTrace::Data::Vector3d; using SolTrace::Data::VirtualElement; using SolTrace::Data::VirtualPlane; -// Template Types +// CST Template Types using SolTrace::Data::Heliostat; using SolTrace::Data::LinearFresnel; using SolTrace::Data::ParabolicDish; @@ -50,7 +51,7 @@ using SolTrace::Data::ray_source_ptr; using SolTrace::Data::stage_ptr; using SolTrace::Data::surface_ptr; -// Functions +// Construction Functions using SolTrace::Data::make_aperture; using SolTrace::Data::make_element; using SolTrace::Data::make_ray_source; @@ -78,6 +79,15 @@ using SolTrace::Data::CalculateTransformMatrices; using SolTrace::Data::TransformToLocal; using SolTrace::Data::TransformToReference; +// Utility Functions +using SolTrace::Data::abs_max; +using SolTrace::Data::abs_min; +using SolTrace::Data::is_approx; +using SolTrace::Data::project_onto_plane; +using SolTrace::Data::project_onto_vector; +using SolTrace::Data::rotate_vector_degrees; +using SolTrace::Data::rotate_vector_radians; + // Status Constants using SolTrace::Data::ELEMENT_ERROR; using SolTrace::Data::ELEMENT_ID_UNASSIGNED; diff --git a/coretrace/simulation_data/surface.cpp b/coretrace/simulation_data/surface.cpp index dca20fd6..d47e5879 100644 --- a/coretrace/simulation_data/surface.cpp +++ b/coretrace/simulation_data/surface.cpp @@ -4,151 +4,267 @@ #include -namespace SolTrace::Data { +#include -surface_ptr make_surface_from_type(SurfaceType type, const std::vector &args) +namespace SolTrace::Data { - surface_ptr retval = nullptr; - unsigned nargs = args.size(); - - switch (type) - { - case SurfaceType::CONE: - retval = nargs < 1 ? nullptr : make_surface(args[0]); - break; - case SurfaceType::CYLINDER: - retval = nargs < 1 ? nullptr : make_surface(1.0 / args[0]); - break; - case SurfaceType::FLAT: - retval = make_surface(); - break; - case SurfaceType::PARABOLA: - { - if (nargs < 2) - retval = nullptr; - else + + surface_ptr make_surface_from_type(SurfaceType type, const std::vector &args) + { + surface_ptr retval = nullptr; + unsigned nargs = args.size(); + + switch (type) + { + case SurfaceType::CONE: + retval = nargs < 1 ? nullptr : make_surface(args[0]); + break; + case SurfaceType::CYLINDER: + retval = nargs < 1 ? nullptr : make_surface(1.0 / args[0]); + break; + case SurfaceType::FLAT: + retval = make_surface(); + break; + case SurfaceType::PARABOLA: { - double cx = args[0]; - double cy = args[1]; - double fx = 1.0 / (2.0 * cx); - double fy = 1.0 / (2.0 * cy); - retval = make_surface(fx, fy); + if (nargs < 2) + retval = nullptr; + else + { + double cx = args[0]; + double cy = args[1]; + double fx = 1.0 / (2.0 * cx); + double fy = 1.0 / (2.0 * cy); + retval = make_surface(fx, fy); + } + break; } - break; - } - case SurfaceType::SPHERE: - retval = nargs < 1 ? nullptr : make_surface(args[0]); - break; - case SurfaceType::HYPER: - case SurfaceType::GENERAL_SPENCER_MURTY: - case SurfaceType::TORUS: - default: - retval = nullptr; // Not implemented yet - break; - } - - return retval; -} - -surface_ptr make_surface_from_json(const nlohmann::ordered_json& jnode) -{ - if (!jnode.contains("surface_type")) - throw std::invalid_argument("Missing surface_type"); - std::string type_str = jnode.at("surface_type"); - SurfaceType surface_type = get_enum_from_string(type_str, SurfaceTypeMap, SurfaceType::SURFACE_UNKNOWN); - if (surface_type == SurfaceType::SURFACE_UNKNOWN) - throw std::invalid_argument("Unknown surface"); - switch (surface_type) - { - case SurfaceType::CONE: return make_surface(jnode); - case SurfaceType::CYLINDER: return make_surface(jnode); - case SurfaceType::FLAT: return make_surface(jnode); - case SurfaceType::PARABOLA: return make_surface(jnode); - case SurfaceType::SPHERE: return make_surface(jnode); + case SurfaceType::SPHERE: + retval = nargs < 1 ? nullptr : make_surface(args[0]); + break; + case SurfaceType::HYPER: + case SurfaceType::GENERAL_SPENCER_MURTY: + case SurfaceType::TORUS: + default: + retval = nullptr; // Not implemented yet + break; + } + + return retval; + } + + surface_ptr make_surface_from_json(const nlohmann::ordered_json &jnode) + { + if (!jnode.contains("surface_type")) + throw std::invalid_argument("Missing surface_type"); + std::string type_str = jnode.at("surface_type"); + SurfaceType surface_type = get_enum_from_string(type_str, SurfaceTypeMap, SurfaceType::SURFACE_UNKNOWN); + if (surface_type == SurfaceType::SURFACE_UNKNOWN) + throw std::invalid_argument("Unknown surface"); + switch (surface_type) + { + case SurfaceType::CONE: + return make_surface(jnode); + case SurfaceType::CYLINDER: + return make_surface(jnode); + case SurfaceType::FLAT: + return make_surface(jnode); + case SurfaceType::PARABOLA: + return make_surface(jnode); + case SurfaceType::SPHERE: + return make_surface(jnode); default: throw std::invalid_argument("Unsupported surface_type: " + type_str); + } } -} -Cone::Cone(const nlohmann::ordered_json& jnode) - : Surface(SurfaceType::CONE) -{ - this->half_angle = jnode.at("half_angle"); -} + Cone::Cone(const nlohmann::ordered_json &jnode) + : Surface(SurfaceType::CONE) + { + this->half_angle = jnode.at("half_angle"); + } -void Cone::write_json(nlohmann::ordered_json& jnode) const -{ - SurfaceType type = SurfaceType::CONE; - jnode["surface_type"] = SurfaceTypeMap.at(type); - jnode["half_angle"] = this->half_angle; -} + void Cone::bounding_box(const double x_minmax[2], + const double y_minmax[2], + double &z_min, + double &z_max) const + { + double theta = this->half_angle * D2R; + double x_abs = abs_max(x_minmax, 2); + double y_abs = abs_max(y_minmax, 2); + z_min = 0.0; + z_max = sqrt(x_abs * x_abs + y_abs * y_abs) / tan(theta); + return; + } -Cylinder::Cylinder(const nlohmann::ordered_json& jnode) - : Surface(SurfaceType::CYLINDER) -{ - this->radius = jnode.at("radius"); -} + surface_ptr Cone::make_copy() const + { + return make_surface(*this); + } -void Cylinder::write_json(nlohmann::ordered_json& jnode) const -{ - SurfaceType type = SurfaceType::CYLINDER; - jnode["surface_type"] = SurfaceTypeMap.at(type); - jnode["radius"] = this->radius; -} + void Cone::write_json(nlohmann::ordered_json &jnode) const + { + SurfaceType type = SurfaceType::CONE; + jnode["surface_type"] = SurfaceTypeMap.at(type); + jnode["half_angle"] = this->half_angle; + } -void Flat::write_json(nlohmann::ordered_json& jnode) const -{ - SurfaceType type = SurfaceType::FLAT; - jnode["surface_type"] = SurfaceTypeMap.at(type); -} + Cylinder::Cylinder(const nlohmann::ordered_json &jnode) + : Surface(SurfaceType::CYLINDER) + { + this->radius = jnode.at("radius"); + } -Parabola::Parabola(const nlohmann::ordered_json& jnode) - : Surface(SurfaceType::PARABOLA) -{ - this->focal_length_x = jnode.at("focal_length_x"); - this->focal_length_y = jnode.at("focal_length_y"); -} + void Cylinder::bounding_box(const double x_minmax[2], + const double y_minmax[2], + double &z_min, + double &z_max) const + { + double r = this->radius; + // Debug check only. This should be caught upstream before + // any bounding box calculations are done by a SimulationRunner. + // See, e.g., native_runner/cylinder_calculator.cpp. + assert(is_approx(x_minmax[0], -r, 1e-6)); + assert(is_approx(x_minmax[1], r, 1e-6)); + z_min = 0.0; + z_max = 2.0 * r; + return; + } -void Parabola::write_json(nlohmann::ordered_json& jnode) const -{ - SurfaceType type = SurfaceType::PARABOLA; - jnode["surface_type"] = SurfaceTypeMap.at(type); - jnode["focal_length_x"] = this->focal_length_x; - jnode["focal_length_y"] = this->focal_length_y; -} - -Sphere::Sphere(const nlohmann::ordered_json& jnode) - : Surface(SurfaceType::SPHERE) -{ - this->vertex_curv = jnode.at("vertex_curv"); -} + surface_ptr Cylinder::make_copy() const + { + return make_surface(*this); + } -void Sphere::write_json(nlohmann::ordered_json& jnode) const -{ - SurfaceType type = SurfaceType::SPHERE; - jnode["surface_type"] = SurfaceTypeMap.at(type); - jnode["vertex_curv"] = this->vertex_curv; -} - -double Cone::z(double x, double y) const { - return sqrt(x * x + y * y) / tan(half_angle); -} -double Cylinder::z(double x, double) const { - // TODO: Fix ? This is really only the top half of the cylinder - // Clyinder breaks the model since it is a mulit-valued fuction: each - // x values produces two z values Returning only the positive root - return radius + sqrt(x * x + radius * radius); -} -double Parabola::z(double x, double y) const { - // z(x,y) = (cx * x^2 + cy * y^2) / 2 - return x * x / focal_length_x + y * y / focal_length_y; -} -double Sphere::z(double x, double y) const { - // TODO: Double-check: This provides only the bottom-half of the sphere, - // aligning with the image in the documentation - // TODO: Verify how vertex_curv relates to radius, i.e., vertex_curve = 1/r? - return vertex_curv * (x * x + y * y) / - (1 + sqrt(1 - vertex_curv * vertex_curv * (x * x + y * y))); -} + void Cylinder::write_json(nlohmann::ordered_json &jnode) const + { + SurfaceType type = SurfaceType::CYLINDER; + jnode["surface_type"] = SurfaceTypeMap.at(type); + jnode["radius"] = this->radius; + } + + void Flat::bounding_box(const double x_minmax[2], + const double y_minmax[2], + double &z_min, + double &z_max) const + { + z_min = -1e-4; + z_max = 1e-4; + return; + } + + surface_ptr Flat::make_copy() const + { + return make_surface(*this); + } + + void Flat::write_json(nlohmann::ordered_json &jnode) const + { + SurfaceType type = SurfaceType::FLAT; + jnode["surface_type"] = SurfaceTypeMap.at(type); + } + + Parabola::Parabola(const nlohmann::ordered_json &jnode) + : Surface(SurfaceType::PARABOLA) + { + this->focal_length_x = jnode.at("focal_length_x"); + this->focal_length_y = jnode.at("focal_length_y"); + } + + void Parabola::bounding_box(const double x_minmax[2], + const double y_minmax[2], + double &z_min, + double &z_max) const + { + double cx = 0.5 / this->focal_length_x; + double cy = 0.5 / this->focal_length_y; + double x_max = abs_max(x_minmax, 2); + double y_max = abs_max(y_minmax, 2); + z_max = 0.5 * (cx * x_max * x_max + cy * y_max * y_max); + + if (x_minmax[0] <= 0.0 && 0.0 <= x_minmax[1] && + y_minmax[0] <= 0.0 && 0.0 <= y_minmax[1]) + { + z_min = 0.0; + } + else + { + double x_min = abs_min(x_minmax, 2); + double y_min = abs_min(y_minmax, 2); + z_min = 0.5 * (cx * x_min * x_min + cy * y_min * y_min); + } + + return; + } + + surface_ptr Parabola::make_copy() const + { + return make_surface(*this); + } + + void Parabola::write_json(nlohmann::ordered_json &jnode) const + { + SurfaceType type = SurfaceType::PARABOLA; + jnode["surface_type"] = SurfaceTypeMap.at(type); + jnode["focal_length_x"] = this->focal_length_x; + jnode["focal_length_y"] = this->focal_length_y; + } + + Sphere::Sphere(const nlohmann::ordered_json &jnode) + : Surface(SurfaceType::SPHERE) + { + this->vertex_curv = jnode.at("vertex_curv"); + } + + void Sphere::bounding_box(const double x_minmax[2], + const double y_minmax[2], + double &z_min, + double &z_max) const + { + z_min = 0.0; + double R = 1.0 / this->vertex_curv; + double x_max = abs_max(x_minmax, 2); + double y_max = abs_max(y_minmax, 2); + double rsq = x_max * x_max + y_max * y_max; + z_max = R > sqrt(rsq) ? R - sqrt(R * R - rsq) : R; + return; + } + + surface_ptr Sphere::make_copy() const + { + return make_surface(*this); + } + + void Sphere::write_json(nlohmann::ordered_json &jnode) const + { + SurfaceType type = SurfaceType::SPHERE; + jnode["surface_type"] = SurfaceTypeMap.at(type); + jnode["vertex_curv"] = this->vertex_curv; + } + + double Cone::z(double x, double y) const + { + return sqrt(x * x + y * y) / tan(half_angle); + } + + double Cylinder::z(double x, double) const + { + // TODO: Fix ? This is really only the top half of the cylinder. + // Cylinder breaks the model since it is a multi-valued function: each + // x value produces two z values. Returning only the positive root. + return radius + sqrt(x * x + radius * radius); + } + + double Parabola::z(double x, double y) const + { + // z(x,y) = (cx * x^2 + cy * y^2) / 2 + return x * x / focal_length_x + y * y / focal_length_y; + } + + double Sphere::z(double x, double y) const + { + return vertex_curv * (x * x + y * y) / + (1 + sqrt(1 - vertex_curv * vertex_curv * (x * x + y * y))); + } } // namespace SolTrace::Data diff --git a/coretrace/simulation_data/surface.hpp b/coretrace/simulation_data/surface.hpp index 27068a80..04cec56f 100644 --- a/coretrace/simulation_data/surface.hpp +++ b/coretrace/simulation_data/surface.hpp @@ -18,165 +18,209 @@ #include #include -namespace SolTrace::Data { - -enum SurfaceType -{ - CONE, - CYLINDER, - FLAT, - PARABOLA, - SPHERE, - - HYPER, - GENERAL_SPENCER_MURTY, - TORUS, - - SURFACE_UNKNOWN -}; - -inline const std::map SurfaceTypeMap = +namespace SolTrace::Data { - {SurfaceType::CONE, "CONE"}, - {SurfaceType::CYLINDER, "CYLINDER"}, - {SurfaceType::FLAT, "FLAT"}, - {SurfaceType::PARABOLA, "PARABOLA"}, - {SurfaceType::SPHERE, "SPHERE"}, - {SurfaceType::HYPER, "HYPER"}, - {SurfaceType::GENERAL_SPENCER_MURTY, "GENERAL_SPENCER_MURTY"}, - {SurfaceType::TORUS, "TORUS"}, - {SurfaceType::SURFACE_UNKNOWN, "SURFACE_UNKNOWN"} -}; - -struct Surface -{ -public: - SurfaceType my_type; - - Surface(SurfaceType st) : my_type(st) {} - virtual ~Surface() {} - - SurfaceType get_type() { return my_type; } - - virtual void write_json(nlohmann::ordered_json& jnode) const = 0; - virtual double z(double x, double y) const { return 0; } - - inline std::string get_type_string() const { - switch (my_type) { - case CONE: return "Cone"; - case CYLINDER: return "Cylinder"; - case FLAT: return "Flat"; - case PARABOLA: return "Parabola"; - case SPHERE: return "Sphere"; - case HYPER: return "Hyper"; - case GENERAL_SPENCER_MURTY: return "General Spencer Murty"; - case TORUS: return "Torus"; - case SURFACE_UNKNOWN: return "Unknown"; - } - return "Unknown"; - } -}; - -struct Cone : public Surface -{ - // z(x,y) = sqrt(x^2 + y^2) / tan(theta) - // where theta = half_angle - double half_angle; - Cone(double ha) : Surface(SurfaceType::CONE), half_angle(ha) {} - Cone(const nlohmann::ordered_json& jnode); - virtual ~Cone() {} - virtual void write_json(nlohmann::ordered_json& jnode) const override; - - virtual double z(double x, double y) const override; -}; - -struct Cylinder : public Surface -{ - // x^2 + (z - r)^2 = r^2 - // where r = radius - double radius; - Cylinder(double r) : Surface(SurfaceType::CYLINDER), radius(r) + enum SurfaceType { - } - Cylinder(const nlohmann::ordered_json& jnode); - virtual ~Cylinder() {} - virtual void write_json(nlohmann::ordered_json& jnode) const override; - - virtual double z(double x, double y) const override; -}; - -struct Flat : public Surface -{ - Flat() : Surface(SurfaceType::FLAT) {} - Flat(const nlohmann::ordered_json& jnode) : Surface(SurfaceType::FLAT) {}; - virtual ~Flat() {} - virtual void write_json(nlohmann::ordered_json& jnode) const override; -}; - -struct Parabola : public Surface -{ - // z(x,y) = (cx * x^2 + cy * y^2) / 2 - // TODO: Assuming that vertex_x_curv gives cx and - // that vertex_y_curv gives cy - double focal_length_x; - double focal_length_y; - - Parabola(double focal_x, double focal_y) : Surface(SurfaceType::PARABOLA), - focal_length_x(focal_x), - focal_length_y(focal_y) + CONE, + CYLINDER, + FLAT, + PARABOLA, + SPHERE, + + HYPER, + GENERAL_SPENCER_MURTY, + TORUS, + + SURFACE_UNKNOWN + }; + + inline const std::map SurfaceTypeMap = + { + {SurfaceType::CONE, "CONE"}, + {SurfaceType::CYLINDER, "CYLINDER"}, + {SurfaceType::FLAT, "FLAT"}, + {SurfaceType::PARABOLA, "PARABOLA"}, + {SurfaceType::SPHERE, "SPHERE"}, + {SurfaceType::HYPER, "HYPER"}, + {SurfaceType::GENERAL_SPENCER_MURTY, "GENERAL_SPENCER_MURTY"}, + {SurfaceType::TORUS, "TORUS"}, + {SurfaceType::SURFACE_UNKNOWN, "SURFACE_UNKNOWN"}}; + + struct Surface; + using surface_ptr = std::shared_ptr; + + template + inline auto make_surface(Args &&...args) { + return std::make_shared(std::forward(args)...); } - Parabola(const nlohmann::ordered_json& jnode); - virtual ~Parabola() {} - virtual void write_json(nlohmann::ordered_json& jnode) const override; - virtual double z(double x, double y) const override; -}; - -struct Sphere : public Surface -{ - // z(x,y) = c(x^2 + y^2) / [1 + sqrt(1 - c^2{x^2 + y^2})] - // where c = 1/R. - // TODO: This form seems to be unnecessarily complicated. - // Could easily just use one of the equations - // z(x,y) = (1 - sqrt(1 - c^2 (x^2 + y^2))) / c - // = R - sqrt(R^2 - (x^2 + y^2)) - // Need to check on this. - double vertex_curv; - - Sphere(double curv) : Surface(SurfaceType::SPHERE), - vertex_curv(curv) + struct Surface { - } - Sphere(const nlohmann::ordered_json& jnode); - virtual ~Sphere() {} - virtual void write_json(nlohmann::ordered_json& jnode) const override; - - virtual double z(double x, double y) const override; -}; - -// TODO: Add other surface types. Documentation has the following: -// 1. Hyperboloid/Ellipsoid -// 2. Zernike Series -// 3. VSHOT data -// 4. Finite Element data -// 5. General Spencer & Murty Equation -// 6. Polynomial Series (rotationally symmetric) -// 7. Cubic Spline Interpolation (rotationally symmetric) - -using surface_ptr = std::shared_ptr; - -template -inline auto make_surface(Args &&...args) -{ - return std::make_shared(std::forward(args)...); -} - -surface_ptr make_surface_from_type(SurfaceType type, - const std::vector &args); + public: + SurfaceType my_type; + + Surface(SurfaceType st) : my_type(st) {} + virtual ~Surface() {} + + SurfaceType get_type() { return my_type; } + + virtual void bounding_box(const double x_minmax[2], + const double y_minmax[2], + double &z_min, + double &z_max) const = 0; + + virtual surface_ptr make_copy() const = 0; + + virtual void write_json(nlohmann::ordered_json &jnode) const = 0; + + virtual double z(double x, double y) const { return 0; } + + inline std::string get_type_string() const + { + switch (my_type) + { + case CONE: + return "Cone"; + case CYLINDER: + return "Cylinder"; + case FLAT: + return "Flat"; + case PARABOLA: + return "Parabola"; + case SPHERE: + return "Sphere"; + case HYPER: + return "Hyper"; + case GENERAL_SPENCER_MURTY: + return "General Spencer Murty"; + case TORUS: + return "Torus"; + case SURFACE_UNKNOWN: + return "Unknown"; + } + return "Unknown"; + } + }; -surface_ptr make_surface_from_json(const nlohmann::ordered_json& jnode); + struct Cone : public Surface + { + // z(x,y) = sqrt(x^2 + y^2) / tan(theta) + // where theta = half_angle + double half_angle; + Cone(double ha) : Surface(SurfaceType::CONE), half_angle(ha) {} + Cone(const nlohmann::ordered_json &jnode); + virtual ~Cone() {} + virtual void bounding_box(const double x_minmax[2], + const double y_minmax[2], + double &z_min, + double &z_max) const override; + virtual surface_ptr make_copy() const override; + virtual void write_json(nlohmann::ordered_json &jnode) const override; + + virtual double z(double x, double y) const override; + }; + + struct Cylinder : public Surface + { + // x^2 + (z - r)^2 = r^2 + // where r = radius + double radius; + Cylinder(double r) : Surface(SurfaceType::CYLINDER), radius(r) + { + } + Cylinder(const nlohmann::ordered_json &jnode); + virtual ~Cylinder() {} + virtual void bounding_box(const double x_minmax[2], + const double y_minmax[2], + double &z_min, + double &z_max) const override; + virtual surface_ptr make_copy() const override; + virtual void write_json(nlohmann::ordered_json &jnode) const override; + + virtual double z(double x, double y) const override; + }; + + struct Flat : public Surface + { + Flat() : Surface(SurfaceType::FLAT) {} + Flat(const nlohmann::ordered_json &jnode) : Surface(SurfaceType::FLAT) {}; + virtual ~Flat() {} + virtual void bounding_box(const double x_minmax[2], + const double y_minmax[2], + double &z_min, + double &z_max) const override; + virtual surface_ptr make_copy() const override; + virtual void write_json(nlohmann::ordered_json &jnode) const override; + }; + + struct Parabola : public Surface + { + // z(x,y) = (cx * x^2 + cy * y^2) / 2 + // TODO: Assuming that vertex_x_curv gives cx and + // that vertex_y_curv gives cy + double focal_length_x; + double focal_length_y; + + Parabola(double focal_x, double focal_y) : Surface(SurfaceType::PARABOLA), + focal_length_x(focal_x), + focal_length_y(focal_y) + { + } + Parabola(const nlohmann::ordered_json &jnode); + virtual ~Parabola() {} + virtual void bounding_box(const double x_minmax[2], + const double y_minmax[2], + double &z_min, + double &z_max) const override; + virtual surface_ptr make_copy() const override; + virtual void write_json(nlohmann::ordered_json &jnode) const override; + + virtual double z(double x, double y) const override; + }; + + struct Sphere : public Surface + { + // z(x,y) = c(x^2 + y^2) / [1 + sqrt(1 - c^2{x^2 + y^2})] + // where c = 1/R. + // TODO: This form seems to be unnecessarily complicated. + // Could easily just use one of the equations + // z(x,y) = (1 - sqrt(1 - c^2 (x^2 + y^2))) / c + // = R - sqrt(R^2 - (x^2 + y^2)) + // Need to check on this. + double vertex_curv; + + Sphere(double curv) : Surface(SurfaceType::SPHERE), + vertex_curv(curv) + { + } + Sphere(const nlohmann::ordered_json &jnode); + virtual ~Sphere() {} + virtual void bounding_box(const double x_minmax[2], + const double y_minmax[2], + double &z_min, + double &z_max) const override; + virtual surface_ptr make_copy() const override; + virtual void write_json(nlohmann::ordered_json &jnode) const override; + + virtual double z(double x, double y) const override; + }; + + // TODO: Add other surface types. Documentation has the following: + // 1. Hyperboloid/Ellipsoid + // 2. Zernike Series + // 3. VSHOT data + // 4. Finite Element data + // 5. General Spencer & Murty Equation + // 6. Polynomial Series (rotationally symmetric) + // 7. Cubic Spline Interpolation (rotationally symmetric) + + surface_ptr make_surface_from_type(SurfaceType type, + const std::vector &args); + + surface_ptr make_surface_from_json(const nlohmann::ordered_json &jnode); } // namespace SolTrace::Data diff --git a/coretrace/simulation_data/cst_templates/utilities.cpp b/coretrace/simulation_data/utilities.cpp similarity index 100% rename from coretrace/simulation_data/cst_templates/utilities.cpp rename to coretrace/simulation_data/utilities.cpp diff --git a/coretrace/simulation_data/cst_templates/utilities.hpp b/coretrace/simulation_data/utilities.hpp similarity index 66% rename from coretrace/simulation_data/cst_templates/utilities.hpp rename to coretrace/simulation_data/utilities.hpp index 78bad606..9dea49ac 100644 --- a/coretrace/simulation_data/cst_templates/utilities.hpp +++ b/coretrace/simulation_data/utilities.hpp @@ -2,10 +2,41 @@ #ifndef SOLTRACE_UTILITIES_H #define SOLTRACE_UTILITIES_H +#include +#include + #include "vector3d.hpp" namespace SolTrace::Data { +template +R abs_min(const R values[], uint_fast64_t size) +{ + R amin = std::numeric_limits::max(); + for (uint_fast64_t k=0; k < size; ++k) + { + amin = std::min(amin, std::abs(values[k])); + } + return amin; +} + +template +R abs_max(const R values[], uint_fast64_t size) +{ + R amax = 0.0; + for (uint_fast64_t k=0; k < size; ++k) + { + amax = std::max(amax, std::abs(values[k])); + } + return amax; +} + +template +bool is_approx(const R &x, const R &y, const R &atol=1e-6) +{ + return (fabs(x - y) < atol); +} + void project_onto_plane(const Vector3d &n, const Vector3d &u, Vector3d &uproj); diff --git a/coretrace/simulation_runner/CMakeLists.txt b/coretrace/simulation_runner/CMakeLists.txt index 959f6490..29b50821 100644 --- a/coretrace/simulation_runner/CMakeLists.txt +++ b/coretrace/simulation_runner/CMakeLists.txt @@ -1,5 +1,9 @@ add_subdirectory(native_runner) +if (SOLTRACE_BUILD_EMBREE_SUPPORT) + add_subdirectory(embree_runner) +endif() + # find cuda and optix libraries, if found then add optix_runner if (SOLTRACE_BUILD_OPTIX_SUPPORT) @@ -42,6 +46,4 @@ if (SOLTRACE_BUILD_OPTIX_SUPPORT) message(STATUS "OptiX or CUDA not found, skipping optix_runner") endif() - - -endif() \ No newline at end of file +endif() diff --git a/coretrace/simulation_runner/embree_runner/CMakeLists.txt b/coretrace/simulation_runner/embree_runner/CMakeLists.txt new file mode 100644 index 00000000..b1aaabc8 --- /dev/null +++ b/coretrace/simulation_runner/embree_runner/CMakeLists.txt @@ -0,0 +1,55 @@ +Project(embree_runner) + +# include(FetchContent) +# FetchContent_Declare( +# embree +# GIT_REPOSITORY https://github.com/RenderKit/embree.git +# GIT_TAG v4.4.0 +# ) + +include_directories( + . + ${CMAKE_SOURCE_DIR}/coretrace/simulation_data + ${CMAKE_SOURCE_DIR}/coretrace/simulation_runner + ${CMAKE_SOURCE_DIR}/coretrace/simulation_runner/native_runner + ${CMAKE_SOURCE_DIR}/coretrace/simulation_results +) + +set(EMBREE_RUNNER_SRC + bbox_calculator.cpp + find_element_hit_embree.cpp + embree_helper.cpp + embree_runner.cpp + trace_embree.cpp +) + +set(POSITION_INDEPENDENT_CODE ON) + +FIND_PACKAGE(embree 4 REQUIRED) + +if(WIN32) + # Hack to build on Windows + add_library(embree_runner STATIC ${EMBREE_RUNNER_SRC}) +else(WIN32) + add_library(embree_runner SHARED ${EMBREE_RUNNER_SRC}) +endif() + +set_target_properties( + embree_runner + PROPERTIES + DEBUG_POSTFIX "d" + PREFIX "" +) + +target_link_libraries( + embree_runner + PUBLIC + simdata + simresult + native_runner + embree +) + +if (NOT APPLE AND ENABLE_COVERAGE) + SET(CMAKE_CXX_FLAGS "-O0 -coverage") +endif() diff --git a/coretrace/simulation_runner/embree_runner/bbox_calculator.cpp b/coretrace/simulation_runner/embree_runner/bbox_calculator.cpp new file mode 100644 index 00000000..450725e5 --- /dev/null +++ b/coretrace/simulation_runner/embree_runner/bbox_calculator.cpp @@ -0,0 +1,378 @@ +#include "bbox_calculator.hpp" + +#include + +#include +#include + +namespace SolTrace::EmbreeRunner +{ + using SolTrace::Data::MatrixVectorMult_generic; + + using SolTrace::NativeRunner::TElement; + using SolTrace::NativeRunner::telement_ptr; + using SolTrace::NativeRunner::TStage; + using SolTrace::NativeRunner::tstage_ptr; + + // float get_absolute_minmax(const float values[], int size, bool is_max) + // { + // float minmax_abs = std::abs(values[0]); + + // if (is_max) + // { + // for (int i = 1; i < size; i++) + // { + // minmax_abs = std::max(std::abs(values[i]), + // minmax_abs); + // } + // } + // else + // { + // for (int i = 1; i < size; i++) + // { + // minmax_abs = std::min(std::abs(values[i]), + // minmax_abs); + // } + // } + + // return minmax_abs; + // } + + // void process_zernike_bounds(TElement* st_element, float x_minmax[2], float y_minmax[2], + // float(&z_minmax)[2]) + // { + // double z, x, y; + + // float& z_min = z_minmax[0]; + // float& z_max = z_minmax[1]; + + // z_min = std::numeric_limits::infinity(); + // z_max = -std::numeric_limits::infinity(); + + // float x_range[3] = { x_minmax[0], x_minmax[1], 0.f }; + // float y_range[3] = { y_minmax[0], y_minmax[1], 0.f }; + + // for (int xi = 0; xi < 3; ++xi) + // { + // for (int yi = 0; yi < 3; ++yi) + // { + // x = x_range[xi]; + // y = y_range[yi]; + // EvalMono(x, y, st_element->BCoefficients, st_element->FitOrder, 0.0, 0.0, &z); + + // if (z < z_min) z_min = z; + // if (z > z_max) z_max = z; + // } + // } + // } + + // void process_poly_bounds(TElement* st_element, float x_minmax[2], float y_minmax[2], + // float(&z_minmax)[2]) + // { + // double z, x, y; + + // float& z_min = z_minmax[0]; + // float& z_max = z_minmax[1]; + + // z_min = std::numeric_limits::infinity(); + // z_max = -std::numeric_limits::infinity(); + + // float x_range[3] = { x_minmax[0], x_minmax[1], 0.f }; + // float y_range[3] = { y_minmax[0], y_minmax[1], 0.f }; + + // for (int xi = 0; xi < 3; ++xi) + // { + // for (int yi = 0; yi < 3; ++yi) + // { + // x = x_range[xi]; + // y = y_range[yi]; + // EvalPoly(x, y, st_element->PolyCoeffs, st_element->FitOrder, &z); + + // if (z < z_min) z_min = z; + // if (z > z_max) z_max = z; + // } + // } + // } + + // bool find_spline_extrema(std::vector& xa, + // std::vector& ya, + // std::vector& y2a, + // double xMin, double xMax, + // double& yMin, double& yMax) + // { + // if (xa.size() < 2 || xa.size() != ya.size() || xa.size() != y2a.size()) + // return false; + + // int n = xa.size(); + // yMin = INFINITY; + // yMax = -INFINITY; + + // // Check all interval critical points + // for (int i = 0; i < n - 1; i++) { + // double h = xa[i + 1] - xa[i]; + // if (h == 0) continue; + + // // The spline derivative equal to zero gives us a quadratic equation + // // Coefficients derived from the derivative expression in splint function + // double A = h * (y2a[i + 1] - y2a[i]) / 2.0; + // double B = h * y2a[i] / 2.0 - (ya[i + 1] - ya[i]) / h; + // double C = -h * y2a[i] / 6.0; + + // // Solve quadratic equation: A*t² + B*t + C = 0 where t = (x-xa[i])/h + // double discriminant = B * B - 4 * A * C; + + // if (std::abs(A) < 1e-10) { + // // Linear case + // if (std::abs(B) > 1e-10) { + // double t = -C / B; + // double x = xa[i] + t * h; + // if (x >= xMin && x <= xMax && x >= xa[i] && x <= xa[i + 1]) { + // double y, dydx; + // if (splint(xa, ya, y2a, n, x, &y, &dydx)) { + // yMin = std::min(yMin, y); + // yMax = std::max(yMax, y); + // } + // } + // } + // } + // else if (discriminant >= 0) { + // // Two possible roots + // double t1 = (-B + sqrt(discriminant)) / (2 * A); + // double t2 = (-B - sqrt(discriminant)) / (2 * A); + + // double x1 = xa[i] + t1 * h; + // double x2 = xa[i] + t2 * h; + + // // Check if critical points are in this interval and the overall range + // if (x1 >= xa[i] && x1 <= xa[i + 1] && x1 >= xMin && x1 <= xMax) { + // double y, dydx; + // if (splint(xa, ya, y2a, n, x1, &y, &dydx)) { + // yMin = std::min(yMin, y); + // yMax = std::max(yMax, y); + // } + // } + + // if (x2 >= xa[i] && x2 <= xa[i + 1] && x2 >= xMin && x2 <= xMax) { + // double y, dydx; + // if (splint(xa, ya, y2a, n, x2, &y, &dydx)) { + // yMin = std::min(yMin, y); + // yMax = std::max(yMax, y); + // } + // } + // } + // } + + // // Also check endpoints and knot points within range + // for (int i = 0; i < n; i++) { + // if (xa[i] >= xMin && xa[i] <= xMax) { + // double y, dydx; + // if (splint(xa, ya, y2a, n, xa[i], &y, &dydx)) { + // yMin = std::min(yMin, y); + // yMax = std::max(yMax, y); + // } + // } + // } + + // // Check range endpoints if they're not knot points + // double y, dydx; + // if (splint(xa, ya, y2a, n, xMin, &y, &dydx)) { + // yMin = std::min(yMin, y); + // yMax = std::max(yMax, y); + // } + + // if (splint(xa, ya, y2a, n, xMax, &y, &dydx)) { + // yMin = std::min(yMin, y); + // yMax = std::max(yMax, y); + // } + + // // Check x = 0; + // if (splint(xa, ya, y2a, n, 0, &y, &dydx)) { + // yMin = std::min(yMin, y); + // yMax = std::max(yMax, y); + // } + + // return yMin != INFINITY && yMax != -INFINITY; + // } + + // void process_cubic_spline_bounds(TElement* st_element, float x_minmax[2], float y_minmax[2], + // float(&z_minmax)[2]) + // { + // double z, x, y; + + // float& z_min = z_minmax[0]; + // float& z_max = z_minmax[1]; + + // z_min = std::numeric_limits::infinity(); + // z_max = -std::numeric_limits::infinity(); + + // float x_min_abs = get_absolute_minmax(x_minmax, 2, false); + // float x_max_abs = get_absolute_minmax(x_minmax, 2, true); + // float y_min_abs = get_absolute_minmax(y_minmax, 2, false); + // float y_max_abs = get_absolute_minmax(y_minmax, 2, true); + + // double Rho_min = sqrt(x_min_abs * x_min_abs + y_min_abs * y_min_abs); + // double Rho_max = sqrt(x_max_abs * x_max_abs + y_max_abs * y_max_abs); + // double z_min_test; + // double z_max_test; + + // find_spline_extrema(st_element->CubicSplineXData, + // st_element->CubicSplineYData, + // st_element->CubicSplineY2Data, + // Rho_min, Rho_max, z_min_test, z_max_test); + + // /*for (int xi = 0; xi < 3; ++xi) + // { + // for (int yi = 0; yi < 3; ++yi) + // { + // x = x_range[xi]; + // y = y_range[yi]; + // double Rho = sqrt(x * x + y * y); + // double dummy; + // splint(st_element->CubicSplineXData, + // st_element->CubicSplineYData, + // st_element->CubicSplineY2Data, + // st_element->CubicSplineXData.size(), + // Rho, &z, &dummy); + + // if (z < z_min) z_min = z; + // if (z > z_max) z_max = z; + // } + // }*/ + + // z_min = z_min_test; + // z_max = z_max_test; + // } + + // void process_FE_bounds(TElement* st_element, float x_minmax[2], float y_minmax[2], + // float(&z_minmax)[2]) + // { + // float& z_min = z_minmax[0]; + // float& z_max = z_minmax[1]; + + // z_min = std::numeric_limits::infinity(); + // z_max = -std::numeric_limits::infinity(); + + // const MatDoub& xyz_nodes = st_element->FEData.nodes; + + // for (const std::vector& fe_node : xyz_nodes) + // { + // if (fe_node[2] > z_max) + // z_max = fe_node[2]; + // if (fe_node[2] < z_min) + // z_min = fe_node[2]; + // } + // } + + void transform_to_global(const float coord_element[3], + const tstage_ptr st_stage, + const TElement *st_element, + float (&coord_global)[3]) + { + float PosDumStage[3]; + float coord_stage[3]; + MatrixVectorMult_generic(st_element->RLocToRef, coord_element, PosDumStage); + for (int i = 0; i < 3; i++) + coord_stage[i] = PosDumStage[i] + st_element->Origin[i]; + + float PosDumGlob[3]; + MatrixVectorMult_generic(st_stage->RLocToRef, coord_stage, PosDumGlob); + for (int i = 0; i < 3; i++) + coord_global[i] = PosDumGlob[i] + st_stage->Origin[i]; + return; + } + + void transform_bounds(const float min_coord_element[3], + const float max_coord_element[3], + const tstage_ptr st_stage, + const TElement *st_element, + float (&min_coord_global)[3], + float (&max_coord_global)[3]) + { + // Transform min and max bounding box from element coordinates to global + float corners_element[8][3] = + { + {min_coord_element[0], min_coord_element[1], min_coord_element[2]}, + {min_coord_element[0], min_coord_element[1], max_coord_element[2]}, + {min_coord_element[0], max_coord_element[1], min_coord_element[2]}, + {min_coord_element[0], max_coord_element[1], max_coord_element[2]}, + {max_coord_element[0], min_coord_element[1], min_coord_element[2]}, + {max_coord_element[0], min_coord_element[1], max_coord_element[2]}, + {max_coord_element[0], max_coord_element[1], min_coord_element[2]}, + {max_coord_element[0], max_coord_element[1], max_coord_element[2]}}; + + // Convert corners to global coordinates + float corners_global[8][3]; + for (int i = 0; i < 8; i++) + transform_to_global(corners_element[i], st_stage, st_element, corners_global[i]); + + // Find min and max xyz + min_coord_global[0] = corners_global[0][0]; + min_coord_global[1] = corners_global[0][1]; + min_coord_global[2] = corners_global[0][2]; + max_coord_global[0] = corners_global[0][0]; + max_coord_global[1] = corners_global[0][1]; + max_coord_global[2] = corners_global[0][2]; + for (int i = 1; i < 8; i++) + { + for (int j = 0; j < 3; j++) + { + float val = corners_global[i][j]; + if (val < min_coord_global[j]) + min_coord_global[j] = val; + if (val > max_coord_global[j]) + max_coord_global[j] = val; + } + } + } + + bool get_bounds(const TElement *st_element, + float (&min_coord_global)[3], + float (&max_coord_global)[3]) + { + // Get stage + tstage_ptr st_stage = st_element->parent_stage; + + // Define element coord bounds + float min_coord_element[3] = {0.f, 0.f, 0.f}; + float max_coord_element[3] = {0.f, 0.f, 0.f}; + + double x_minmax[2] = {0.0, 0.0}; + double y_minmax[2] = {0.0, 0.0}; + double z_minmax[2] = {0.0, 0.0}; + + st_element->aperture->bounding_box(x_minmax[0], + x_minmax[1], + y_minmax[0], + y_minmax[1]); + + st_element->surface->bounding_box(x_minmax, + y_minmax, + z_minmax[0], + z_minmax[1]); + + // Expand bounding boxes slightly to account for float precision + const float expand = 1e-3f; + x_minmax[0] -= expand; + x_minmax[1] += expand; + y_minmax[0] -= expand; + y_minmax[1] += expand; + z_minmax[0] -= expand; + z_minmax[1] += expand; + + // Assign points to min/max coordinate element arrays + min_coord_element[0] = x_minmax[0]; + min_coord_element[1] = y_minmax[0]; + min_coord_element[2] = z_minmax[0]; + max_coord_element[0] = x_minmax[1]; + max_coord_element[1] = y_minmax[1]; + max_coord_element[2] = z_minmax[1]; + + // Convert local element bounds, to global xyz + transform_bounds(min_coord_element, max_coord_element, + st_stage, st_element, + min_coord_global, max_coord_global); + + return true; + } + +} // namespace SolTrace::EmbreeRunner diff --git a/coretrace/simulation_runner/embree_runner/bbox_calculator.hpp b/coretrace/simulation_runner/embree_runner/bbox_calculator.hpp new file mode 100644 index 00000000..98e319a8 --- /dev/null +++ b/coretrace/simulation_runner/embree_runner/bbox_calculator.hpp @@ -0,0 +1,16 @@ +#ifndef SOLTRACE_BBOX_CALCULATOR_H +#define SOLTRACE_BBOX_CALCULATOR_H + +// #include "types.h" + +#include + +namespace SolTrace::EmbreeRunner +{ + bool get_bounds(const SolTrace::NativeRunner::TElement *st_element, + float (&min_coord_global)[3], + float (&max_coord_global)[3]); + +} // namespace SolTrace::EmbreeRunner + +#endif diff --git a/coretrace/simulation_runner/embree_runner/embree_helper.cpp b/coretrace/simulation_runner/embree_runner/embree_helper.cpp new file mode 100644 index 00000000..643f175b --- /dev/null +++ b/coretrace/simulation_runner/embree_runner/embree_helper.cpp @@ -0,0 +1,288 @@ +#include "embree_helper.hpp" + +#include "bbox_calculator.hpp" + +// SimulationData headers +#include + +// NativeRunner headers +#include +#include + +#include + +#include +#include +#include +#include +#include + +namespace SolTrace::EmbreeRunner +{ + + using SolTrace::Data::CopyVec3; + using SolTrace::Data::TransformToLocal; + using SolTrace::Data::TransformToReference; + + using SolTrace::NativeRunner::TElement; + using SolTrace::NativeRunner::telement_ptr; + using SolTrace::NativeRunner::TStage; + using SolTrace::NativeRunner::tstage_ptr; + + template + bool compare_Vec3(T vec1[3], T vec2[3], double tol_diff) + { + for (int i = 0; i < 3; i++) + { + if ((std::abs(vec1[i] / vec2[i] - 1) > tol_diff) && + (std::abs(vec1[i] - vec2[i]) > 1e-5)) + return false; + } + return true; + } + + void error_function(void *userPtr, + RTCError error, + const char *str) + { + printf("error %d: %s\n", error, str); + } + + int bounds_error(const RTCBoundsFunctionArguments *args) + { + // Set bounds to zero and report error + RTCBounds *bounds_o = args->bounds_o; + bounds_o->lower_x = std::numeric_limits::quiet_NaN(); + bounds_o->upper_x = std::numeric_limits::quiet_NaN(); + bounds_o->lower_y = std::numeric_limits::quiet_NaN(); + bounds_o->upper_y = std::numeric_limits::quiet_NaN(); + bounds_o->lower_z = std::numeric_limits::quiet_NaN(); + bounds_o->upper_z = std::numeric_limits::quiet_NaN(); + + error_function(args->geometryUserPtr, + RTC_ERROR_INVALID_OPERATION, + "Bounding box computation failed for some element"); + + throw std::runtime_error("Embree scene commit failed"); + + return 0; + } + + void bounds_function(const RTCBoundsFunctionArguments *args) + { + // Get element + TElement *st_element = (TElement *)args->geometryUserPtr; + + // Get bounds + float min_coord_global[3]; + float max_coord_global[3]; + bool success = get_bounds(st_element, + min_coord_global, + max_coord_global); + + // Check error + if (!success) + { + bounds_error(args); + } + + // Assign bounds + RTCBounds *bounds_o = args->bounds_o; + bounds_o->lower_x = min_coord_global[0]; + bounds_o->upper_x = max_coord_global[0]; + bounds_o->lower_y = min_coord_global[1]; + bounds_o->upper_y = max_coord_global[1]; + bounds_o->lower_z = min_coord_global[2]; + bounds_o->upper_z = max_coord_global[2]; + } + + void intersect_function(const RTCIntersectFunctionNArguments *args) + { + // Get payload object + RayIntersectPayload *payload = (RayIntersectPayload *)args->context; + + double PosRayGlob[3], CosRayGlob[3]; + CopyVec3(PosRayGlob, payload->PosRayGlobIn); + CopyVec3(CosRayGlob, payload->CosRayGlobIn); + + // Get Element data + TElement *st_element = (TElement *)args->geometryUserPtr; + tstage_ptr st_stage = st_element->parent_stage; + + // First, convert ray coordinates to element + // Global -> stage -> element + // transform the global incoming ray to local stage coordinates + double PosRayStage[3], CosRayStage[3]; + TransformToLocal(PosRayGlob, CosRayGlob, + st_stage->Origin, st_stage->RRefToLoc, + PosRayStage, CosRayStage); + + // {Transform ray to element[j] coord system of Stage[i]} + double PosRayElement[3], CosRayElement[3]; + TransformToLocal(PosRayStage, CosRayStage, + st_element->Origin, st_element->RRefToLoc, + PosRayElement, CosRayElement); + + // Increment position by tiny amount to get off the element if + // tracing to the same element. + PosRayElement[0] = PosRayElement[0] + 1.0e-4 * CosRayElement[0]; + PosRayElement[1] = PosRayElement[1] + 1.0e-4 * CosRayElement[1]; + PosRayElement[2] = PosRayElement[2] + 1.0e-4 * CosRayElement[2]; + + // Call DeterminElementIntersectionNew + double PosRaySurfElement[3] = {0.0, 0.0, 0.0}; + double CosRaySurfElement[3] = {0.0, 0.0, 0.0}; + double DFXYZ[3] = {0.0, 0.0, 0.0}; + double PathLength = 0; + + int InterceptFlag = 0; + int HitBackSide = 0; + SolTrace::NativeRunner::DetermineElementIntersectionNew( + st_element, PosRayElement, CosRayElement, + PosRaySurfElement, CosRaySurfElement, DFXYZ, + &PathLength, &payload->ErrorFlag, &InterceptFlag, &HitBackSide); + + // Update rayhit info (if hit) + if (InterceptFlag != 0) + { + // Get rayhit data + RTCRayHit *rayhit_out = (RTCRayHit *)args->rayhit; + + // Check if hit is closer than other hits + // Using payload pathlength for double precision + // If pathlength == lastpathlength, use which element number is lower + // to match results with original code + if ((PathLength < payload->LastPathLength) || + ((PathLength == payload->LastPathLength) && + (st_element->element_number < payload->element_number))) + { + + if (PosRaySurfElement[2] <= st_element->ZAperture) + { + + // Transform ray back to stage coordinate system + double PosRaySurfStage[3] = {0.0, 0.0, 0.0}; + double CosRaySurfStage[3] = {0.0, 0.0, 0.0}; + TransformToReference(PosRaySurfElement, CosRaySurfElement, + st_element->Origin, st_element->RLocToRef, + PosRaySurfStage, CosRaySurfStage); + + rayhit_out->ray.tfar = (float)PathLength; // Update intersection distance + rayhit_out->hit.geomID = args->geomID; + rayhit_out->hit.primID = 0; // Single primitive + + // Define Ng (will not be used) + rayhit_out->hit.Ng_x = std::numeric_limits::quiet_NaN(); + rayhit_out->hit.Ng_y = std::numeric_limits::quiet_NaN(); + rayhit_out->hit.Ng_z = std::numeric_limits::quiet_NaN(); + + // Assign custom outputs + payload->LastHitBackSide = HitBackSide; + CopyVec3(payload->LastDFXYZ, DFXYZ); + // CopyVec3(payload->LastPosRaySurfGlob, PosRaySurfGlob); + // CopyVec3(payload->LastCosRaySurfGlob, CosRaySurfGlob); + CopyVec3(payload->LastPosRaySurfStage, PosRaySurfStage); + CopyVec3(payload->LastCosRaySurfStage, CosRaySurfStage); + CopyVec3(payload->LastPosRaySurfElement, PosRaySurfElement); + CopyVec3(payload->LastCosRaySurfElement, CosRaySurfElement); + payload->element_number = st_element->element_number; + payload->LastPathLength = PathLength; + } + } + } + + return; + } + + RTCScene make_scene(RTCDevice &device, + SolTrace::NativeRunner::TSystem &system) + { + // Make scene + RTCScene scene = rtcNewScene(device); + + // Loop through stages + unsigned int stage_index = 1; + for (tstage_ptr stage : system.StageList) + { + // Loop through elements in each stage + unsigned int stage_mask = 1u << stage_index; + for (telement_ptr st_element : stage->ElementList) + { + // Make embree geometry for each element + RTCGeometry embree_geom = rtcNewGeometry(device, RTC_GEOMETRY_TYPE_USER); + rtcSetGeometryUserPrimitiveCount(embree_geom, 1); + + // Attach st element to embree geometry + rtcSetGeometryUserData(embree_geom, st_element.get()); + + // Assign bounds function + rtcSetGeometryBoundsFunction(embree_geom, bounds_function, NULL); + + // Assign intersect function + rtcSetGeometryIntersectFunction(embree_geom, intersect_function); + + // Set mask (stage number) + rtcSetGeometryMask(embree_geom, stage_mask); + + // Commit geometry + rtcCommitGeometry(embree_geom); + + // Attach geometry + unsigned int geomID = rtcAttachGeometry(scene, embree_geom); + + // Release geometry (it is owned by the scene now) + rtcReleaseGeometry(embree_geom); + } + + stage_index++; + } + + return scene; + } + + bool validate_intersect(double (&LastPosRaySurfElement1)[3], + double (&LastCosRaySurfElement1)[3], + double (&LastDFXYZ1)[3], + uint_fast64_t &LastElementNumber1, + uint_fast64_t &LastRayNumber1, + double (&LastPosRaySurfStage1)[3], + double (&LastCosRaySurfStage1)[3], + int &ErrorFlag1, + int &LastHitBackSide1, + bool &StageHit1, + double (&LastPosRaySurfElement2)[3], + double (&LastCosRaySurfElement2)[3], + double (&LastDFXYZ2)[3], + uint_fast64_t &LastElementNumber2, + uint_fast64_t &LastRayNumber2, + double (&LastPosRaySurfStage2)[3], + double (&LastCosRaySurfStage2)[3], + int &ErrorFlag2, + int &LastHitBackSide2, + bool &StageHit2) + { + if (StageHit1 == false && StageHit2 == false) + return true; + + double tol_diff = 1e-4; + if (StageHit1 != StageHit2) + return false; + if (LastElementNumber1 != LastElementNumber2) + return false; + if (LastHitBackSide1 != LastHitBackSide2) + return false; + if (!compare_Vec3(LastPosRaySurfElement1, LastPosRaySurfElement2, tol_diff)) + return false; + if (!compare_Vec3(LastCosRaySurfElement1, LastCosRaySurfElement2, tol_diff)) + return false; + if (!compare_Vec3(LastDFXYZ1, LastDFXYZ2, tol_diff)) + return false; + if (!compare_Vec3(LastPosRaySurfStage1, LastPosRaySurfStage2, tol_diff)) + return false; + if (!compare_Vec3(LastCosRaySurfStage1, LastCosRaySurfStage2, tol_diff)) + return false; + + return true; + } + +} // namespace SolTrace::EmbreeRunner \ No newline at end of file diff --git a/coretrace/simulation_runner/embree_runner/embree_helper.hpp b/coretrace/simulation_runner/embree_runner/embree_helper.hpp new file mode 100644 index 00000000..ceaad9a4 --- /dev/null +++ b/coretrace/simulation_runner/embree_runner/embree_helper.hpp @@ -0,0 +1,75 @@ +#ifndef SOLTRACE_EMBREE_HELPER_H +#define SOLTRACE_EMBREE_HELPER_H + +#include + +#include + +#include +#include +#include + +namespace SolTrace::EmbreeRunner +{ + + struct RayIntersectPayload + { + RTCRayQueryContext context; // Embree built-in context (MUST come first) + + int LastHitBackSide = -1; + double LastDFXYZ[3] = {0.0, 0.0, 0.0}; + + // double LastPosRaySurfGlob[3] = { 0.0, 0.0, 0.0 }; + // double LastCosRaySurfGlob[3] = { 0.0, 0.0, 0.0 }; + + double LastPosRaySurfStage[3] = {0.0, 0.0, 0.0}; + double LastCosRaySurfStage[3] = {0.0, 0.0, 0.0}; + + double LastPosRaySurfElement[3] = {0.0, 0.0, 0.0}; + double LastCosRaySurfElement[3] = {0.0, 0.0, 0.0}; + + int_fast64_t element_number = -1; + int ErrorFlag = 0; + + double PosRayGlobIn[3] = {0.0, 0.0, 0.0}; + double CosRayGlobIn[3] = {0.0, 0.0, 0.0}; + + double LastPathLength = std::numeric_limits::infinity(); + }; + + inline unsigned int embree_mask(int_fast64_t stage_index) + { + return 1u << stage_index; + } + + void error_function(void *userPtr, + RTCError error, + const char *str); + + RTCScene make_scene(RTCDevice &device, + SolTrace::NativeRunner::TSystem &system); + + bool validate_intersect(double (&LastPosRaySurfElement1)[3], + double (&LastCosRaySurfElement1)[3], + double (&LastDFXYZ1)[3], + uint_fast64_t &LastElementNumber1, + uint_fast64_t &LastRayNumber1, + double (&LastPosRaySurfStage1)[3], + double (&LastCosRaySurfStage1)[3], + int &ErrorFlag1, + int &LastHitBackSide1, + bool &StageHit1, + double (&LastPosRaySurfElement2)[3], + double (&LastCosRaySurfElement2)[3], + double (&LastDFXYZ2)[3], + uint_fast64_t &LastElementNumber2, + uint_fast64_t &LastRayNumber2, + double (&LastPosRaySurfStage2)[3], + double (&LastCosRaySurfStage2)[3], + int &ErrorFlag2, + int &LastHitBackSide2, + bool &StageHit2); + +} // namespace SolTrace::EmbreeRunner + +#endif diff --git a/coretrace/simulation_runner/embree_runner/embree_runner.cpp b/coretrace/simulation_runner/embree_runner/embree_runner.cpp new file mode 100644 index 00000000..52484840 --- /dev/null +++ b/coretrace/simulation_runner/embree_runner/embree_runner.cpp @@ -0,0 +1,94 @@ + +#include "embree_runner.hpp" + +#include "trace_embree.hpp" + +#include +#include +#include +#include + +namespace SolTrace::EmbreeRunner +{ + using SolTrace::Runner::RunnerStatus; + + using SolTrace::NativeRunner::telement_ptr; + using SolTrace::NativeRunner::TRayData; + using SolTrace::NativeRunner::tstage_ptr; + using SolTrace::NativeRunner::TSystem; + + using SolTrace::Result::SimulationResult; + + EmbreeRunner::EmbreeRunner() : NativeRunner(), + embree_device(nullptr), + embree_scene(nullptr) + { + return; + } + + EmbreeRunner::~EmbreeRunner() + { + this->clean_embree(); + return; + } + + RunnerStatus EmbreeRunner::setup_simulation(const SimulationData *data) + { + + RunnerStatus sts = NativeRunner::setup_simulation(data); + + make_embree_scene(this->my_logger, + &this->tsys, + this->embree_device, + this->embree_scene, + this->number_of_threads); + + return sts; + } + + RunnerStatus EmbreeRunner::update_simulation(const SimulationData *data) + { + // TODO: Do a more efficient implementation of this? + this->clean_embree(); + NativeRunner::update_simulation(data); + return RunnerStatus::SUCCESS; + } + + RunnerStatus EmbreeRunner::run_simulation() + { + this->set_seeds(); + + RunnerStatus sts = trace_embree( + this->my_manager, + this->my_logger, + &this->tsys, + // this->tsys.seed, + this->seeds, + this->number_of_threads, + this->tsys.sim_raycount, + this->tsys.sim_raymax, + this->tsys.sim_errors_sunshape, + this->tsys.sim_errors_optical, + this->embree_scene); + + return sts; + } + + void EmbreeRunner::clean_embree() + { + // Clean embree + if (this->embree_scene != nullptr) + { + rtcReleaseScene(embree_scene); + embree_scene = nullptr; + } + + if (this->embree_device != nullptr) + { + rtcReleaseDevice(embree_device); + embree_device = nullptr; + } + return; + } + +} // namespace SolTrace::EmbreeRunner diff --git a/coretrace/simulation_runner/embree_runner/embree_runner.hpp b/coretrace/simulation_runner/embree_runner/embree_runner.hpp new file mode 100644 index 00000000..0b398992 --- /dev/null +++ b/coretrace/simulation_runner/embree_runner/embree_runner.hpp @@ -0,0 +1,63 @@ +#ifndef SOLTRACE_EMBREE_RUNNER_H +#define SOLTRACE_EMBREE_RUNNER_H + +#include + +#include +#include + +#include +#include + +namespace SolTrace::EmbreeRunner +{ + using SolTrace::Runner::RunnerStatus; + using SolTrace::Runner::SimulationRunner; + + // NativeRunner types that are used here that we want to make visible + // through the EmbreeRunner namespace + using SolTrace::NativeRunner::TRayData; + using SolTrace::NativeRunner::TSystem; + using SolTrace::NativeRunner::tstage_ptr; + using SolTrace::NativeRunner::telement_ptr; + + class EmbreeRunner : public SolTrace::NativeRunner::NativeRunner + { + public: + EmbreeRunner(); + virtual ~EmbreeRunner(); + + EmbreeRunner(const EmbreeRunner &) = delete; + EmbreeRunner(EmbreeRunner &&) = delete; + + virtual RunnerStatus setup_simulation(const SolTrace::Data::SimulationData *data); + virtual RunnerStatus run_simulation(); + virtual RunnerStatus update_simulation(const SolTrace::Data::SimulationData *data); + + // TODO: Do we want loud errors when a user calls these? + void disable_power_tower() = delete; + void enable_power_tower() = delete; + void disable_point_focus() = delete; + void enable_point_focus() = delete; + + private: + // // Number of threads to use when tracing + // uint_fast64_t number_of_threads; + // std::vector seeds; + + // SolTrace::NativeRunner::ElementParameters eparams; + // SolTrace::NativeRunner::TSystem tsys; + + // bool set_aperture_planes(SolTrace::NativeRunner::TSystem *tsys); + // bool set_aperture_planes(SolTrace::NativeRunner::tstage_ptr stage); + // bool aperture_plane(SolTrace::NativeRunner::telement_ptr Element); + + RTCDevice embree_device; + RTCScene embree_scene; + + void clean_embree(); + }; + +} // namespace SolTrace::EmbreeRunner + +#endif diff --git a/coretrace/simulation_runner/embree_runner/find_element_hit_embree.cpp b/coretrace/simulation_runner/embree_runner/find_element_hit_embree.cpp new file mode 100644 index 00000000..3f94fd66 --- /dev/null +++ b/coretrace/simulation_runner/embree_runner/find_element_hit_embree.cpp @@ -0,0 +1,90 @@ +#include "find_element_hit_embree.hpp" + +#include + +// Embree header +#include + +#include + +#include "embree_helper.hpp" + +namespace SolTrace::EmbreeRunner +{ + void FindElementHit_embree( + // Embree args + const RTCScene &scene, + // Ray info + const int i, + const uint_fast64_t RayNumber, + const double (&PosRayGlob)[3], + const double (&CosRayGlob)[3], + // outputs + double (&LastPosRaySurfElement)[3], + double (&LastCosRaySurfElement)[3], + double (&LastDFXYZ)[3], + uint_fast64_t &LastElementNumber, + uint_fast64_t &LastRayNumber, + double (&LastPosRaySurfStage)[3], + double (&LastCosRaySurfStage)[3], + int &ErrorFlag, + int &LastHitBackSide, + bool &StageHit) + { + // Initialize outputs + StageHit = false; + + // Make payload object to store intersect outputs + RayIntersectPayload ray_payload; + CopyVec3(ray_payload.PosRayGlobIn, PosRayGlob); // Copy position (with full double precision) + CopyVec3(ray_payload.CosRayGlobIn, CosRayGlob); // Copy direction (with full double precision) + rtcInitRayQueryContext(&ray_payload.context); + RTCIntersectArguments args; + rtcInitIntersectArguments(&args); + args.context = &ray_payload.context; + ray_payload.LastPathLength = std::numeric_limits::infinity(); + + // Make rayhit object + RTCRayHit rayhit; + rayhit.ray.org_x = PosRayGlob[0]; + rayhit.ray.org_y = PosRayGlob[1]; + rayhit.ray.org_z = PosRayGlob[2]; + rayhit.ray.dir_x = CosRayGlob[0]; + rayhit.ray.dir_y = CosRayGlob[1]; + rayhit.ray.dir_z = CosRayGlob[2]; + + // Define rayhit outputs + rayhit.ray.tnear = 0; + rayhit.ray.tfar = std::numeric_limits::infinity(); + rayhit.ray.mask = 1u << (i + 1); + rayhit.ray.flags = 0; + rayhit.hit.geomID = RTC_INVALID_GEOMETRY_ID; + rayhit.hit.instID[0] = RTC_INVALID_GEOMETRY_ID; + + // Find intersection + rtcIntersect1(scene, &rayhit, &args); + + // Check if ray hit any elements + if (rayhit.hit.geomID != RTC_INVALID_GEOMETRY_ID) + { + // Collect intersection outputs + StageHit = true; + CopyVec3(LastPosRaySurfElement, ray_payload.LastPosRaySurfElement); + CopyVec3(LastCosRaySurfElement, ray_payload.LastCosRaySurfElement); + CopyVec3(LastDFXYZ, ray_payload.LastDFXYZ); + LastElementNumber = ray_payload.element_number; + LastRayNumber = RayNumber; + CopyVec3(LastPosRaySurfStage, ray_payload.LastPosRaySurfStage); + CopyVec3(LastCosRaySurfStage, ray_payload.LastCosRaySurfStage); + ErrorFlag = ray_payload.ErrorFlag; + LastHitBackSide = ray_payload.LastHitBackSide; + } + + // No hit + else + { + StageHit = false; + } + } + +} // namespace SolTrace::EmbreeRunner diff --git a/coretrace/simulation_runner/embree_runner/find_element_hit_embree.hpp b/coretrace/simulation_runner/embree_runner/find_element_hit_embree.hpp new file mode 100644 index 00000000..9f3842d3 --- /dev/null +++ b/coretrace/simulation_runner/embree_runner/find_element_hit_embree.hpp @@ -0,0 +1,31 @@ +#ifndef SOLTRACE_FIND_ELEMENT_HIT_H +#define SOLTRACE_FIND_ELEMENT_HIT_H + +#include + +#include + +namespace SolTrace::EmbreeRunner +{ + void FindElementHit_embree( + // Embree args + const RTCScene &scene, + // Ray info + const int i, + const uint_fast64_t RayNumber, + const double (&PosRayGlob)[3], + const double (&CosRayGlob)[3], + // outputs + double (&LastPosRaySurfElement)[3], + double (&LastCosRaySurfElement)[3], + double (&LastDFXYZ)[3], + uint_fast64_t &LastElementNumber, + uint_fast64_t &LastRayNumber, + double (&LastPosRaySurfStage)[3], + double (&LastCosRaySurfStage)[3], + int &ErrorFlag, + int &LastHitBackSide, + bool &StageHit); +} // namespace SolTrace::EmbreeRunner + +#endif diff --git a/coretrace/simulation_runner/embree_runner/trace_embree.cpp b/coretrace/simulation_runner/embree_runner/trace_embree.cpp new file mode 100644 index 00000000..c87bd23c --- /dev/null +++ b/coretrace/simulation_runner/embree_runner/trace_embree.cpp @@ -0,0 +1,601 @@ +#include "trace_embree.hpp" + +#include +#include +#include +#include +#include + +// Embree header +#include + +// SimulationData header +#include + +// SimulationRunner header +#include + +// NativeRunner header(s) +#include +#include +#include +#include +#include +#include +#include +#include + +#include "embree_helper.hpp" +#include "find_element_hit_embree.hpp" + +namespace SolTrace::EmbreeRunner +{ + using SolTrace::NativeRunner::GlobalRay_refactored; + using SolTrace::NativeRunner::MTRand; + using SolTrace::NativeRunner::TElement; + using SolTrace::NativeRunner::telement_ptr; + using SolTrace::NativeRunner::thread_manager_ptr; + using SolTrace::NativeRunner::ThreadManager; + using SolTrace::NativeRunner::trace_logger_ptr; + using SolTrace::NativeRunner::tstage_ptr; + using SolTrace::NativeRunner::TSystem; + + using SolTrace::Runner::RunnerStatus; + + using SolTrace::Result::RayEvent; + + RunnerStatus make_embree_scene(trace_logger_ptr logger, + TSystem *System, + RTCDevice &embree_device, + RTCScene &embree_scene, + unsigned nthreads) + { + RunnerStatus sts = RunnerStatus::SUCCESS; + // // Initialize Embree vars + // RTCDevice embree_device = nullptr; + // RTCScene embree_scene = nullptr; + // bool use_shared_embree = false; + + // Make device + // std::cout << "Making embree device..." << std::endl; + std::stringstream ss; + ss << "threads=" << nthreads; + embree_device = rtcNewDevice(ss.str().c_str()); + + // std::cout << "Setting error function..." << std::endl; + rtcSetDeviceErrorFunction(embree_device, error_function, NULL); + + // Convert st stages into scene + // std::cout << "Making scene..." << std::endl; + embree_scene = make_scene(embree_device, *System); + + // std::cout << "Committing scene..." << std::endl; + rtcCommitScene(embree_scene); + + // Validate bounds + RTCError err = rtcGetDeviceError(embree_device); + if (err != RTC_ERROR_NONE) + { + // int asdg = 0; + // return RunnerStatus::ERROR; + sts = RunnerStatus::ERROR; + logger->error_log("Error setting up Embree scene"); + } + + return sts; + } + + RunnerStatus trace_embree( + thread_manager_ptr manager, + trace_logger_ptr logger, + TSystem *System, + const std::vector &seeds, + unsigned nthreads, + uint_fast64_t NumberOfRays, + uint_fast64_t MaxNumberOfRays, + bool IncludeSunShape, + bool IncludeErrors, + const RTCScene &embree_scene) + { + + System->RayData.SetUp(nthreads, NumberOfRays); + System->SunRayCount = 0; + + // Initialize Sun + Vector3d PosSunStage; + bool status = SolTrace::NativeRunner::SunToPrimaryStage( + logger, System, System->StageList[0].get(), + &System->Sun, PosSunStage.data); + + if (!status) + return RunnerStatus::ERROR; + + uint_fast64_t rem = NumberOfRays % nthreads; + uint_fast64_t nrays_per_thread = NumberOfRays / nthreads; + uint_fast64_t nrays; + + for (unsigned k = 0; k < nthreads; ++k) + { + nrays = k < rem ? nrays_per_thread + 1 : nrays_per_thread; + ThreadManager::future my_future = std::async( + std::launch::async, + trace_embree_single_thread, + k, + manager, + logger, + System, + seeds[k], + nrays, + MaxNumberOfRays / nthreads + 1, + IncludeSunShape, + IncludeErrors, + PosSunStage, + embree_scene); + manager->manage(k, std::move(my_future)); + } + + return manager->monitor_until_completion(); + } + + RunnerStatus trace_embree_single_thread( + unsigned thread_id, + thread_manager_ptr manager, + trace_logger_ptr logger, + TSystem *System, + unsigned seed, + uint_fast64_t NumberOfRays, + uint_fast64_t MaxNumberOfRays, + bool IncludeSunShape, + bool IncludeErrors, + const SolTrace::Data::Vector3d &PosSunStage, + const RTCScene &embree_scene) + { + // std::cout << "Thread " << thread_id << " with seed " << seed + // << std::endl; + // Initialize Internal State Variables + MTRand myrng(seed); + + uint_fast64_t update_rate = std::min( + std::max(static_cast(1), NumberOfRays / 10), + static_cast(1000)); + uint_fast64_t update_count = 0; + double total_work = System->StageList.size() * NumberOfRays; + + uint_fast64_t RayNumber = 1; // Ray Number of current ray + bool PreviousStageHasRays = false; + uint_fast64_t LastRayNumberInPreviousStage = NumberOfRays; + + // Define IncomingRays + std::vector IncomingRays; + IncomingRays.resize(NumberOfRays); + + // Initialize stage variables + uint_fast64_t StageDataArrayIndex = 0; + uint_fast64_t PreviousStageDataArrayIndex = 0; + uint_fast64_t n_rays_active = NumberOfRays; + uint_fast64_t sun_ray_count_local = 0; + + // Loop through stages + for (uint_fast64_t i = 0; i < System->StageList.size(); i++) + { + // Check if previous stage has rays + bool StageHasRays = true; + if (i > 0 && PreviousStageHasRays == false) + { + StageHasRays = false; + } + + // Get Current Stage + tstage_ptr Stage = System->StageList[i]; + + // Initialize stage variables + StageDataArrayIndex = 0; + PreviousStageDataArrayIndex = 0; + + // Loop through rays + while (StageHasRays) + { + // Initialize Global Coordinates + double PosRayGlob[3] = {0.0, 0.0, 0.0}; + double CosRayGlob[3] = {0.0, 0.0, 0.0}; + + // Initialize Stage Coordinates + double PosRayStage[3] = {0.0, 0.0, 0.0}; + double CosRayStage[3] = {0.0, 0.0, 0.0}; + + // Get Ray + if (i == 0) + { + // Make ray (if first stage) + double PosRaySun[3]; + SolTrace::NativeRunner::GenerateRay( + myrng, PosSunStage.data, Stage->Origin, + Stage->RLocToRef, &System->Sun, + PosRayGlob, CosRayGlob, PosRaySun); + System->SunRayCount++; + } + else + { + // Get ray from previous stage + RayNumber = IncomingRays[StageDataArrayIndex].Num; + CopyVec3(PosRayGlob, IncomingRays[StageDataArrayIndex].Pos); + CopyVec3(CosRayGlob, IncomingRays[StageDataArrayIndex].Cos); + StageDataArrayIndex++; + } + + // transform the global incoming ray to local stage coordinates + TransformToLocal(PosRayGlob, CosRayGlob, + Stage->Origin, Stage->RRefToLoc, + PosRayStage, CosRayStage); + + // Initialize internal variables for ray intersection tracing + bool RayInStage = true; + bool in_multi_hit_loop = false; + double LastPosRaySurfElement[3] = {0.0, 0.0, 0.0}; + double LastCosRaySurfElement[3] = {0.0, 0.0, 0.0}; + double LastPosRaySurfStage[3] = {0.0, 0.0, 0.0}; + double LastCosRaySurfStage[3] = {0.0, 0.0, 0.0}; + double LastDFXYZ[3] = {0.0, 0.0, 0.0}; + uint_fast64_t LastElementNumber = 0; + uint_fast64_t LastRayNumber = 0; + int ErrorFlag; + int LastHitBackSide; + bool StageHit; + int MultipleHitCount = 0; + double PosRayOutElement[3] = {0.0, 0.0, 0.0}; + double CosRayOutElement[3] = {0.0, 0.0, 0.0}; + + // Start Loop to trace ray until it leaves stage + bool RayIsAbsorbed = false; + while (RayInStage) + { + + FindElementHit_embree(embree_scene, i, RayNumber, + PosRayGlob, CosRayGlob, + LastPosRaySurfElement, LastCosRaySurfElement, + LastDFXYZ, LastElementNumber, LastRayNumber, + LastPosRaySurfStage, LastCosRaySurfStage, + ErrorFlag, LastHitBackSide, StageHit); + + // Breakout if ray left stage + if (!StageHit) + { + RayInStage = false; + break; + } + + // Increment MultipleHitCount? + MultipleHitCount++; + + if (i == 0 && MultipleHitCount == 1) + { + // Add ray to Stage RayData + auto r = System->RayData.Append(thread_id, + PosRayGlob, + CosRayGlob, + ELEMENT_NULL, + i + 1, + LastRayNumber, + RayEvent::CREATE); + + if (r == nullptr) + { + std::stringstream ss; + ss << "Failed to record ray data.\n"; + logger->error_log(ss.str()); + } + } + + // Get optics and check for absorption + const OpticalProperties *optics = 0; + RayEvent rev = RayEvent::VIRTUAL; + if (Stage->Virtual) + { + // If stage is virtual, there is no interaction + CopyVec3(PosRayOutElement, LastPosRaySurfElement); + CopyVec3(CosRayOutElement, LastCosRaySurfElement); + } + else + { + telement_ptr optelm = + Stage->ElementList[LastElementNumber - 1]; + if (LastHitBackSide) + optics = &optelm->Optics.Back; + else + optics = &optelm->Optics.Front; + + bool good = + SolTrace::NativeRunner::determine_interaction_type( + logger, + i, + thread_id, + myrng, + optics, + LastDFXYZ, + LastCosRaySurfElement, + rev); + + if (!good) + { + return RunnerStatus::ERROR; + } + + if (rev == RayEvent::ABSORB) + { + RayIsAbsorbed = true; + break; + } + } + + // Process Interaction + int k = LastElementNumber - 1; + SolTrace::NativeRunner::ProcessInteraction( + System, + myrng, + IncludeSunShape, + optics, + IncludeErrors, + i, + Stage, + MultipleHitCount, + LastDFXYZ, + LastCosRaySurfElement, + ErrorFlag, + CosRayOutElement, + LastPosRaySurfElement, + PosRayOutElement); + + // Transform ray back to stage coordinate system + TransformToReference(PosRayOutElement, + CosRayOutElement, + Stage->ElementList[k]->Origin, + Stage->ElementList[k]->RLocToRef, + PosRayStage, + CosRayStage); + TransformToReference(PosRayStage, + CosRayStage, + Stage->Origin, + Stage->RLocToRef, + PosRayGlob, + CosRayGlob); + + System->RayData.Append(thread_id, + PosRayGlob, + CosRayGlob, + LastElementNumber, + i + 1, + LastRayNumber, + rev); + + // Break out if multiple hits are not allowed + if (!Stage->MultiHitsPerRay) + { + StageHit = false; + break; + } + else + { + in_multi_hit_loop = true; + } + } + + if (MultipleHitCount > 0) + ++update_count; + + if (update_count % update_rate == 0) + { + double progress = update_count / total_work; + manager->progress_update(thread_id, progress); + if (manager->terminate(thread_id)) + return RunnerStatus::CANCEL; + } + + // Handle if Ray was absorbed + if (RayIsAbsorbed) + { + TransformToReference(LastPosRaySurfStage, + LastCosRaySurfStage, + Stage->Origin, + Stage->RLocToRef, + PosRayGlob, + CosRayGlob); + + System->RayData.Append(thread_id, + PosRayGlob, + CosRayGlob, + LastElementNumber, + i + 1, + LastRayNumber, + RayEvent::ABSORB); + + n_rays_active--; + + // ray was fully absorbed + if (RayNumber == LastRayNumberInPreviousStage) + { + PreviousStageHasRays = false; + if (PreviousStageDataArrayIndex > 0) + { + PreviousStageDataArrayIndex--; + PreviousStageHasRays = true; + } + // goto Label_EndStageLoop; + break; + } + else + { + if (i == 0) + { + if (RayNumber == NumberOfRays) + // goto Label_EndStageLoop; + break; + else + RayNumber++; + } + + // Next ray in loop + continue; + } + } + + // Ray has left the stage + bool FlagMiss = false; + if (i == 0) + { + if (MultipleHitCount == 0) + { + // Ray in first stage missed stage entirely + // Generate new ray + continue; + } + else + { + // Ray hit an element, so save it for next stage + CopyVec3(IncomingRays[PreviousStageDataArrayIndex].Pos, PosRayGlob); + CopyVec3(IncomingRays[PreviousStageDataArrayIndex].Cos, CosRayGlob); + IncomingRays[PreviousStageDataArrayIndex].Num = RayNumber; + + // Is Ray the last in the stage? + if (RayNumber == NumberOfRays) + { + StageHasRays = false; + break; + } + + PreviousStageDataArrayIndex++; + PreviousStageHasRays = true; + + // Move on to next ray + RayNumber++; + continue; + } + } + else + { + // After the first stage + // Ray hit element OR is traced through stage + if (Stage->TraceThrough || MultipleHitCount > 0) + { + // Ray is saved for the next stage + CopyVec3(IncomingRays[PreviousStageDataArrayIndex].Pos, PosRayGlob); + CopyVec3(IncomingRays[PreviousStageDataArrayIndex].Cos, CosRayGlob); + IncomingRays[PreviousStageDataArrayIndex].Num = RayNumber; + + // Check if ray is last in stage + if (RayNumber == LastRayNumberInPreviousStage) + { + StageHasRays = false; + break; + } + + PreviousStageDataArrayIndex++; + PreviousStageHasRays = true; + + if (MultipleHitCount == 0) + { + FlagMiss = true; + } + + // Go to next ray + continue; + } + // Ray missed stage entirely and is not traced + else + { + FlagMiss = true; + } + + // Handle FlagMiss condition ( + if (FlagMiss == true) + { + LastRayNumber = RayNumber; + + System->RayData.Append(thread_id, + PosRayGlob, + CosRayGlob, + ELEMENT_NULL, + i + 1, + LastRayNumber, + RayEvent::EXIT); + + n_rays_active--; + + if (RayNumber == LastRayNumberInPreviousStage) + { + if (!Stage->TraceThrough) + { + PreviousStageHasRays = false; + if (PreviousStageDataArrayIndex > 0) + { + PreviousStageHasRays = true; + PreviousStageDataArrayIndex--; // last ray was previous one + } + } + + // Exit stage + StageHasRays = false; + break; + } + else + { + if (i == 0) + RayNumber++; // generate new sun ray + + // Start new ray + continue; + } + } + } + } + + // EndStage section... + + // skipping save_st_data logic + + if (!PreviousStageHasRays) + { + LastRayNumberInPreviousStage = 0; + continue; // No rays to carry forward + } + + if (PreviousStageDataArrayIndex < IncomingRays.size()) + { + LastRayNumberInPreviousStage = IncomingRays[PreviousStageDataArrayIndex].Num; + if (LastRayNumberInPreviousStage == 0) + { + // size_t pp = IncomingRays[PreviousStageDataArrayIndex - 1].Num; + // System->errlog("LastRayNumberInPreviousStage=0, stage %d, PrevIdx=%d, CurIdx=%d, pp=%d", i + 1, + // PreviousStageDataArrayIndex, StageDataArrayIndex, pp); + return RunnerStatus::ERROR; + } + } + else + { + // System->errlog("Invalid PreviousStageDataArrayIndex: %u, @ stage %d", + // PreviousStageDataArrayIndex, i + 1); + return RunnerStatus::ERROR; + } + } + + // Close out any remaining rays as misses + unsigned idx = System->StageList.size() - 1; + tstage_ptr Stage = System->StageList[idx]; + for (uint_fast64_t k = 0; k < n_rays_active; ++k) + { + GlobalRay_refactored ray = IncomingRays[k]; + System->RayData.Append(thread_id, + ray.Pos, + ray.Cos, + ELEMENT_NULL, + idx + 1, + ray.Num, + RayEvent::EXIT); + } + + // System->SunRayCount is atomic so this is thread safe + System->SunRayCount += sun_ray_count_local; + + return RunnerStatus::SUCCESS; + } + +} // namespace SolTrace::EmbreeRunner diff --git a/coretrace/simulation_runner/embree_runner/trace_embree.hpp b/coretrace/simulation_runner/embree_runner/trace_embree.hpp new file mode 100644 index 00000000..568b1d33 --- /dev/null +++ b/coretrace/simulation_runner/embree_runner/trace_embree.hpp @@ -0,0 +1,50 @@ +#ifndef SOLTRACE_TRACE_EMBREE_H +#define SOLTRACE_TRACE_EMBREE_H + +#include +#include + +#include +#include +#include +#include +#include + +namespace SolTrace::EmbreeRunner +{ + + SolTrace::Runner::RunnerStatus make_embree_scene( + SolTrace::NativeRunner::trace_logger_ptr logger, + SolTrace::NativeRunner::TSystem *System, + RTCDevice &embree_device, + RTCScene &embree_scene, + unsigned nthreads); + + SolTrace::Runner::RunnerStatus trace_embree( + SolTrace::NativeRunner::thread_manager_ptr manager, + SolTrace::NativeRunner::trace_logger_ptr logger, + SolTrace::NativeRunner::TSystem *System, + const std::vector &seeds, + unsigned nthreads, + uint_fast64_t NumberOfRays, + uint_fast64_t MaxNumberOfRays, + bool IncludeSunShape, + bool IncludeErrors, + const RTCScene &embree_scene); + + SolTrace::Runner::RunnerStatus trace_embree_single_thread( + unsigned thread_id, + SolTrace::NativeRunner::thread_manager_ptr manager, + SolTrace::NativeRunner::trace_logger_ptr logger, + SolTrace::NativeRunner::TSystem *System, + unsigned seed, + uint_fast64_t NumberOfRays, + uint_fast64_t MaxNumberOfRays, + bool IncludeSunShape, + bool IncludeErrors, + const SolTrace::Data::Vector3d &PosSunStage, + const RTCScene &embree_scene); + +} // namespace SolTrace::EmbreeRunner + +#endif diff --git a/coretrace/simulation_runner/native_runner/CMakeLists.txt b/coretrace/simulation_runner/native_runner/CMakeLists.txt index 0187e7a3..c893cfd1 100644 --- a/coretrace/simulation_runner/native_runner/CMakeLists.txt +++ b/coretrace/simulation_runner/native_runner/CMakeLists.txt @@ -11,6 +11,7 @@ set(NATIVE_RUNNER_SRC calculator_factory.cpp cylinder_calculator.cpp determine_element_intersection_new.cpp + determine_interaction_type.cpp find_element_hit.cpp flat_calculator.cpp generate_ray.cpp @@ -25,6 +26,7 @@ set(NATIVE_RUNNER_SRC thread_manager.cpp trace.cpp tracing_errors.cpp + trace_logger.cpp treemesh.cpp ) @@ -51,6 +53,5 @@ target_link_libraries( ) if (NOT APPLE AND ENABLE_COVERAGE) - # SET(CMAKE_CXX_FLAGS "-O0 -coverage -fkeep-inline-functions") SET(CMAKE_CXX_FLAGS "-O0 -coverage") endif() diff --git a/coretrace/simulation_runner/native_runner/calculator_factory.cpp b/coretrace/simulation_runner/native_runner/calculator_factory.cpp index 23787369..ad028082 100644 --- a/coretrace/simulation_runner/native_runner/calculator_factory.cpp +++ b/coretrace/simulation_runner/native_runner/calculator_factory.cpp @@ -22,6 +22,7 @@ calculator_ptr CalculatorFactory::make_calculator( surface_ptr surf, const ElementParameters &eparams) { + // TODO: Make logger accessible to intersection calculators... // Input validation if (surf == nullptr) { diff --git a/coretrace/simulation_runner/native_runner/determine_interaction_type.cpp b/coretrace/simulation_runner/native_runner/determine_interaction_type.cpp new file mode 100644 index 00000000..0170a357 --- /dev/null +++ b/coretrace/simulation_runner/native_runner/determine_interaction_type.cpp @@ -0,0 +1,123 @@ +#include "determine_interaction_type.hpp" + +#include + +#include + +#include "mtrand.hpp" +#include "native_runner_types.hpp" +#include "trace_logger.hpp" + +namespace SolTrace::NativeRunner +{ + using SolTrace::Result::RayEvent; + + bool determine_interaction_type( + trace_logger_ptr logger, + int_fast64_t stage, + unsigned thread_id, + MTRand &myrng, + const OpticalProperties *optics, + const double (&LastDFXYZ)[3], + const double (&LastCosRaySurfElement)[3], + // bool LastHitBackSide, + RayEvent &rev) + { + bool good = true; + rev = RayEvent::VIRTUAL; + + double TestValue; + double UnitLastDFXYZ[3] = {0.0, 0.0, 0.0}; + double IncidentAngle = 0; + // TODO: Implement tables... + switch (optics->my_type) + { + case InteractionType::REFRACTION: + // if (optics->UseTransmissivityTable) + // { + // int npoints = optics->TransmissivityTable.size(); + // int m = 0; + + // UnitLastDFXYZ[0] = -LastDFXYZ[0] / sqrt(DOT(LastDFXYZ, LastDFXYZ)); + // UnitLastDFXYZ[1] = -LastDFXYZ[1] / sqrt(DOT(LastDFXYZ, LastDFXYZ)); + // UnitLastDFXYZ[2] = -LastDFXYZ[2] / sqrt(DOT(LastDFXYZ, LastDFXYZ)); + // IncidentAngle = acos(DOT(LastCosRaySurfElement, UnitLastDFXYZ)) * 1000.; //[mrad] + // if (IncidentAngle >= optics->TransmissivityTable[npoints - 1].angle) + // { + // TestValue = optics->TransmissivityTable[npoints - 1].trans; + // } + // else + // { + // while (optics->TransmissivityTable[m].angle < IncidentAngle) + // m++; + + // if (m == 0) + // TestValue = optics->TransmissivityTable[m].trans; + // else + // TestValue = (optics->TransmissivityTable[m].trans + optics->TransmissivityTable[m - 1].trans) / 2.0; + // } + // } + // else + // { + // TestValue = optics->transmitivity; + // rev = RayEvent::TRANSMIT; + // } + TestValue = optics->transmitivity; + rev = RayEvent::TRANSMIT; + break; + case InteractionType::REFLECTION: + // if (optics->UseReflectivityTable) + // { + // int npoints = optics->ReflectivityTable.size(); + // int m = 0; + // UnitLastDFXYZ[0] = -LastDFXYZ[0] / sqrt(DOT(LastDFXYZ, LastDFXYZ)); + // UnitLastDFXYZ[1] = -LastDFXYZ[1] / sqrt(DOT(LastDFXYZ, LastDFXYZ)); + // UnitLastDFXYZ[2] = -LastDFXYZ[2] / sqrt(DOT(LastDFXYZ, LastDFXYZ)); + // IncidentAngle = acos(DOT(LastCosRaySurfElement, UnitLastDFXYZ)) * 1000.; //[mrad] + // if (IncidentAngle >= optics->ReflectivityTable[npoints - 1].angle) + // { + // TestValue = optics->ReflectivityTable[npoints - 1].refl; + // } + // else + // { + // while (optics->ReflectivityTable[m].angle < IncidentAngle) + // m++; + + // if (m == 0) + // TestValue = optics->ReflectivityTable[m].refl; + // else + // TestValue = (optics->ReflectivityTable[m].refl + optics->ReflectivityTable[m - 1].refl) / 2.0; + // } + // } + // else + // { + // TestValue = optics->reflectivity; + // rev = RayEvent::REFLECT; + // } + TestValue = optics->reflectivity; + rev = RayEvent::REFLECT; + break; + default: + good = false; + std::stringstream ss; + ss << "Bad optical interaction." + << " Type: " << static_cast(optics->my_type) + << " Stage: " << stage + << " Thread: " << thread_id + << "\n"; + logger->error_log(ss.str()); + break; + } + + // Apply MonteCarlo probability of absorption. Limited + // for now, but can make more complex later on if desired + if (TestValue <= myrng()) + { + // ray was fully absorbed + rev = RayEvent::ABSORB; + } + + return good; + } + +} // namespace SolTrace::NativeRunner diff --git a/coretrace/simulation_runner/native_runner/determine_interaction_type.hpp b/coretrace/simulation_runner/native_runner/determine_interaction_type.hpp new file mode 100644 index 00000000..d05def09 --- /dev/null +++ b/coretrace/simulation_runner/native_runner/determine_interaction_type.hpp @@ -0,0 +1,27 @@ +#ifndef SOLTRACE_DETERMINE_INTERACTION_TYPE_H +#define SOLTRACE_DETERMINE_INTERACTION_TYPE_H + +#include +#include + +#include "mtrand.hpp" +#include "native_runner_types.hpp" +#include "trace_logger.hpp" + +namespace SolTrace::NativeRunner +{ + + bool determine_interaction_type( + trace_logger_ptr logger, + int_fast64_t stage, + unsigned thread_id, + MTRand &myrng, + const SolTrace::Data::OpticalProperties *optics, + const double (&LastDFXYZ)[3], + const double (&LastCosRaySurfElement)[3], + // bool LastHitBackSide, + SolTrace::Result::RayEvent &rev); + +} // namespace SolTrace::NativeRunner + +#endif diff --git a/coretrace/simulation_runner/native_runner/find_element_hit.cpp b/coretrace/simulation_runner/native_runner/find_element_hit.cpp index 88a49335..38179a2f 100644 --- a/coretrace/simulation_runner/native_runner/find_element_hit.cpp +++ b/coretrace/simulation_runner/native_runner/find_element_hit.cpp @@ -127,7 +127,6 @@ namespace SolTrace::NativeRunner CopyVec3(LastPosRaySurfElement, PosRaySurfElement); CopyVec3(LastCosRaySurfElement, CosRaySurfElement); CopyVec3(LastDFXYZ, DFXYZ); - // LastElementNumber = ((i == 0 && !PT_override) ? Element->element_number : j + 1); // mjw change from j index to element id LastElementNumber = Element->element_number; LastRayNumber = RayNumber; TransformToReference(PosRaySurfElement, CosRaySurfElement, diff --git a/coretrace/simulation_runner/native_runner/native_runner.cpp b/coretrace/simulation_runner/native_runner/native_runner.cpp index 2328035d..3a0b1aca 100644 --- a/coretrace/simulation_runner/native_runner/native_runner.cpp +++ b/coretrace/simulation_runner/native_runner/native_runner.cpp @@ -17,6 +17,7 @@ // NativeRunner headers #include "native_runner_types.hpp" #include "trace.hpp" +#include "trace_logger.hpp" namespace SolTrace::NativeRunner { @@ -25,11 +26,16 @@ namespace SolTrace::NativeRunner as_power_tower(false), number_of_threads(1) { - this->my_manager = make_thread_manager(); + this->my_logger = make_trace_logger(); + this->my_manager = make_thread_manager(this->my_logger); + return; } NativeRunner::~NativeRunner() { + this->my_manager = nullptr; + this->my_logger = nullptr; + return; } RunnerStatus NativeRunner::initialize() @@ -155,7 +161,7 @@ namespace SolTrace::NativeRunner auto my_map = std::map(); // int_fast64_t current_stage_id = -1; tstage_ptr current_stage = nullptr; - int_fast64_t element_number = 1; + // int_fast64_t element_number = 1; bool element_found_before_stage = false; if (data->get_number_of_elements() <= 0) @@ -184,17 +190,17 @@ namespace SolTrace::NativeRunner { // TODO: Duplicate stage numbers. Need to make an error // message. + throw std::runtime_error("Duplicate stage numbers found."); sts = RunnerStatus::ERROR; } current_stage = stage; - element_number = 1; + // element_number = 1; } else if (el->is_enabled() && el->is_single()) { if (current_stage == nullptr) { - // throw std::runtime_error("No stage to add element to"); element_found_before_stage = true; continue; } @@ -205,10 +211,11 @@ namespace SolTrace::NativeRunner } telement_ptr elem = make_telement(iter->second, - element_number, + current_stage, this->eparams); - ++element_number; - current_stage->ElementList.push_back(elem); + // ++element_number; + // current_stage->ElementList.push_back(elem); + current_stage->add_element(elem); } } @@ -225,7 +232,7 @@ namespace SolTrace::NativeRunner // set to correspond to global coordinates. This is necessary // so that the element coordinate setup in make_element are // correct. - int_fast64_t element_number = 1; + // int_fast64_t element_number = 1; auto stage = make_tstage(this->eparams); stage->ElementList.reserve(data->get_number_of_elements()); for (auto iter = data->get_const_iterator(); @@ -236,10 +243,11 @@ namespace SolTrace::NativeRunner if (el->is_enabled() && el->is_single()) { telement_ptr tel = make_telement(el, - element_number, + stage, this->eparams); - stage->ElementList.push_back(tel); - ++element_number; + // stage->ElementList.push_back(tel); + // ++element_number; + stage->add_element(tel); } } my_map.insert(std::make_pair(0, stage)); @@ -280,22 +288,11 @@ namespace SolTrace::NativeRunner RunnerStatus NativeRunner::run_simulation() { - if (this->seeds.empty() || - this->seeds.size() != this->number_of_threads) - { - this->seeds.clear(); - for (unsigned k = 0; k < this->number_of_threads; ++k) - { - this->seeds.push_back(this->tsys.seed + 123 * k); - } - } - else - { - ; // Intentional no-op - } + this->set_seeds(); RunnerStatus sts = trace_native( this->my_manager, + this->my_logger, &this->tsys, this->seeds, this->number_of_threads, @@ -450,4 +447,22 @@ namespace SolTrace::NativeRunner return true; } + void NativeRunner::set_seeds() + { + if (this->seeds.empty() || + this->seeds.size() != this->number_of_threads) + { + this->seeds.clear(); + for (unsigned k = 0; k < this->number_of_threads; ++k) + { + this->seeds.push_back(this->tsys.seed + 123 * k); + } + } + else + { + ; // Intentional no-op + } + return; + } + } // namespace SolTrace::NativeRunner diff --git a/coretrace/simulation_runner/native_runner/native_runner.hpp b/coretrace/simulation_runner/native_runner/native_runner.hpp index e7107515..abdf4583 100644 --- a/coretrace/simulation_runner/native_runner/native_runner.hpp +++ b/coretrace/simulation_runner/native_runner/native_runner.hpp @@ -80,7 +80,7 @@ namespace SolTrace::NativeRunner void print_log(std::ostream &os) { - this->my_manager->print_log(os); + this->my_logger->print_log(os); return; } @@ -109,7 +109,7 @@ namespace SolTrace::NativeRunner RunnerStatus setup_sun(const SolTrace::Data::SimulationData *data); RunnerStatus setup_elements(const SolTrace::Data::SimulationData *data); - private: + protected: // Use power tower speed ups bool as_power_tower; @@ -120,12 +120,17 @@ namespace SolTrace::NativeRunner ElementParameters eparams; + trace_logger_ptr my_logger; thread_manager_ptr my_manager; TSystem tsys; bool set_aperture_planes(TSystem *tsys); bool set_aperture_planes(tstage_ptr stage); bool aperture_plane(telement_ptr Element); + + void set_seeds(); + + private: }; } // namespace SolTrace::NativeRunner diff --git a/coretrace/simulation_runner/native_runner/native_runner_types.cpp b/coretrace/simulation_runner/native_runner/native_runner_types.cpp index 38acacce..5424ced5 100644 --- a/coretrace/simulation_runner/native_runner/native_runner_types.cpp +++ b/coretrace/simulation_runner/native_runner/native_runner_types.cpp @@ -407,7 +407,7 @@ namespace SolTrace::NativeRunner if (Query(i, pos, cos, &elm, &stage, &ray, &rev)) { printf(" [%llu] = { [%lg,%lg,%lg][%lg,%lg,%lg] %d %d %llu %s\(%d) }\n", - i, + static_cast(i), pos[0], pos[1], pos[2], cos[0], cos[1], cos[2], elm, @@ -469,7 +469,7 @@ namespace SolTrace::NativeRunner } telement_ptr make_telement(element_ptr el, - int_fast64_t el_num, + tstage_ptr my_stage, const ElementParameters &eparams) { // std::cout << "Name: " << el->get_name() @@ -487,6 +487,7 @@ namespace SolTrace::NativeRunner matrix_copy(telem->RLocToRef, el->get_local_to_stage()); telem->aperture = el->get_aperture()->make_copy(); + telem->surface = el->get_surface()->make_copy(); telem->icalc = CalculatorFactory::get()->make_calculator(telem->aperture, el->get_surface(), @@ -496,9 +497,10 @@ namespace SolTrace::NativeRunner telem->Optics.Front = *el->get_front_optical_properties(); telem->Optics.Back = *el->get_back_optical_properties(); - // telem->element_number = el->get_id(); telem->sim_data_id = el->get_id(); - telem->element_number = el_num; + telem->element_number = my_stage->next_element_number(); + + telem->parent_stage = my_stage; return telem; } diff --git a/coretrace/simulation_runner/native_runner/native_runner_types.hpp b/coretrace/simulation_runner/native_runner/native_runner_types.hpp index efd3fd62..bbc1da7a 100644 --- a/coretrace/simulation_runner/native_runner/native_runner_types.hpp +++ b/coretrace/simulation_runner/native_runner/native_runner_types.hpp @@ -69,6 +69,22 @@ namespace SolTrace::NativeRunner { + class GlobalRay_refactored + { + public: + GlobalRay_refactored() // : active(true) + { + Num = 0; + for (int i = 0; i < 3; i++) + Pos[i] = Cos[i] = 0.0; + } + + double Pos[3]; + double Cos[3]; + uint_fast64_t Num; + // bool active; + }; + // #define ACOSM1O180 0.017453292519943295 // acos(-1)/180.0 // #endif @@ -132,9 +148,8 @@ namespace SolTrace::NativeRunner uint_fast64_t newton_max_iters; }; - class TOpticalPropertySet + struct TOpticalPropertySet { - public: std::string Name; // TOpticalProperties Front; // TOpticalProperties Back; @@ -142,6 +157,10 @@ namespace SolTrace::NativeRunner SolTrace::Data::OpticalProperties Back; }; + // Forward declaration so TElement can have a stage pointer + struct TStage; + using tstage_ptr = typename std::shared_ptr; + struct TElement { TElement(); @@ -162,6 +181,7 @@ namespace SolTrace::NativeRunner SolTrace::Data::aperture_ptr aperture; /////////// SURFACE PARAMETERS /////////////// + SolTrace::Data::surface_ptr surface; calculator_ptr icalc; // double Kappa; @@ -181,11 +201,12 @@ namespace SolTrace::NativeRunner // of addition to element list int_fast64_t element_number; SolTrace::Data::element_id sim_data_id; + tstage_ptr parent_stage; }; using telement_ptr = typename std::shared_ptr; telement_ptr make_telement(SolTrace::Data::element_ptr el, - int_fast64_t el_num, + tstage_ptr my_stage, const ElementParameters &eparams); struct TSun @@ -306,6 +327,18 @@ namespace SolTrace::NativeRunner TStage(); ~TStage(); + void add_element(telement_ptr telem) + { + this->ElementList.push_back(telem); + return; + } + + int_fast64_t next_element_number() const + { + // element number is 1-based + return this->ElementList.size() + 1; + } + bool MultiHitsPerRay; bool Virtual; bool TraceThrough; @@ -328,7 +361,6 @@ namespace SolTrace::NativeRunner int_fast64_t stage_id; }; - using tstage_ptr = typename std::shared_ptr; tstage_ptr make_tstage(const ElementParameters &eparams); tstage_ptr make_tstage(SolTrace::Data::element_ptr el, const ElementParameters &eparams); diff --git a/coretrace/simulation_runner/native_runner/sun_to_primary_stage.cpp b/coretrace/simulation_runner/native_runner/sun_to_primary_stage.cpp index 190a1799..a28f3540 100644 --- a/coretrace/simulation_runner/native_runner/sun_to_primary_stage.cpp +++ b/coretrace/simulation_runner/native_runner/sun_to_primary_stage.cpp @@ -3,13 +3,13 @@ #include "native_runner_types.hpp" #include "simulation_data_export.hpp" -#include "thread_manager.hpp" +#include "trace_logger.hpp" namespace SolTrace::NativeRunner { bool SunToPrimaryStage( - thread_manager_ptr manager, + trace_logger_ptr logger, TSystem *System, TStage *Stage, TSun *Sun, @@ -61,7 +61,7 @@ namespace SolTrace::NativeRunner if (dtot == 0.0) { - manager->error_log("error calculating sun position in primary stage, dtot = 0.0\n"); + logger->error_log("error calculating sun position in primary stage, dtot = 0.0\n"); return false; } @@ -193,7 +193,7 @@ namespace SolTrace::NativeRunner } if (nelements == 0) - manager->error_log("error calculating sun position in primary stage because there are no elements"); + logger->error_log("error calculating sun position in primary stage because there are no elements"); return (nelements > 0); } diff --git a/coretrace/simulation_runner/native_runner/sun_to_primary_stage.hpp b/coretrace/simulation_runner/native_runner/sun_to_primary_stage.hpp index 2d23f6cc..aeff706c 100644 --- a/coretrace/simulation_runner/native_runner/sun_to_primary_stage.hpp +++ b/coretrace/simulation_runner/native_runner/sun_to_primary_stage.hpp @@ -2,13 +2,13 @@ #define SOLTRACE_SUN_TO_PRIMARY_STAGE_H #include "native_runner_types.hpp" -#include "thread_manager.hpp" +#include "trace_logger.hpp" namespace SolTrace::NativeRunner { bool SunToPrimaryStage( - thread_manager_ptr manager, + trace_logger_ptr logger, TSystem *System, TStage *Stage, TSun *Sun, diff --git a/coretrace/simulation_runner/native_runner/thread_manager.cpp b/coretrace/simulation_runner/native_runner/thread_manager.cpp index 55d3af38..99d27b3b 100644 --- a/coretrace/simulation_runner/native_runner/thread_manager.cpp +++ b/coretrace/simulation_runner/native_runner/thread_manager.cpp @@ -8,7 +8,7 @@ namespace SolTrace::NativeRunner { - ThreadManager::ThreadManager() + ThreadManager::ThreadManager(trace_logger_ptr log) : logger(log) { this->initialize(); return; @@ -26,6 +26,7 @@ namespace SolTrace::NativeRunner assert(this->threads.find(id) == this->threads.end()); this->threads[id] = std::move(f); + std::lock_guard lk(this->progress_mutex); this->progress[id] = 0.0; return id; } @@ -70,7 +71,7 @@ namespace SolTrace::NativeRunner << " returned status " << status_string(thread_status) << ". This is an error."; - this->error_log(ss.str()); + this->logger->error_log(ss.str()); // Unexpected behavior so we terminate this->cancel(); @@ -93,13 +94,6 @@ namespace SolTrace::NativeRunner return sts; } - void ThreadManager::error_log(const std::string &msg) - { - std::lock_guard lk(this->message_mutex); - this->messages.push_back(msg); - return; - } - void ThreadManager::progress_update(unsigned int id, double prog) { // std::cout << "Update from thread " << id << " reporting " @@ -152,18 +146,6 @@ namespace SolTrace::NativeRunner return; } - void ThreadManager::print_log(std::ostream &os) const - { - std::lock_guard lk(this->message_mutex); - for (auto iter = this->messages.cbegin(); - iter != this->messages.cend(); - ++iter) - { - os << *iter; - } - return; - } - void ThreadManager::initialize() { // Must call prior to manage and not while running!! @@ -177,10 +159,6 @@ namespace SolTrace::NativeRunner this->progress.clear(); this->threads.clear(); } - { - std::lock_guard lk(this->message_mutex); - this->messages.clear(); - } return; } diff --git a/coretrace/simulation_runner/native_runner/thread_manager.hpp b/coretrace/simulation_runner/native_runner/thread_manager.hpp index 19efe1a6..167bfe71 100644 --- a/coretrace/simulation_runner/native_runner/thread_manager.hpp +++ b/coretrace/simulation_runner/native_runner/thread_manager.hpp @@ -10,6 +10,8 @@ #include +#include "trace_logger.hpp" + namespace SolTrace::NativeRunner { @@ -18,7 +20,8 @@ namespace SolTrace::NativeRunner public: using ThreadStatus = SolTrace::Runner::RunnerStatus; using future = std::future; - ThreadManager(); + + ThreadManager(trace_logger_ptr logger); ~ThreadManager(); // Functions for whoever "hired" the manager. Pattern for use is @@ -27,6 +30,7 @@ namespace SolTrace::NativeRunner // the returned future to manage along with a unique id // 3. Call monitor_until_completion() to wait for everything // to finish. + // NOTE: Tasks should return type of ThreadManager::ThreadStatus // Must call prior to manage and not while running!! void initialize(); @@ -35,14 +39,12 @@ namespace SolTrace::NativeRunner ThreadStatus monitor_until_completion(); // For the worker threads that are being managed - void error_log(const std::string &msg); void progress_update(unsigned int id, double progress); bool terminate(unsigned int id); // For general public use -- thread safe and can be called whenever ThreadStatus status(double *progress = nullptr) const; void cancel() const; - void print_log(std::ostream &os) const; private: // unsigned int next_id; @@ -55,8 +57,7 @@ namespace SolTrace::NativeRunner std::map threads; - mutable std::mutex message_mutex; - std::vector messages; + trace_logger_ptr logger; }; using thread_manager_ptr = std::shared_ptr; diff --git a/coretrace/simulation_runner/native_runner/trace.cpp b/coretrace/simulation_runner/native_runner/trace.cpp index 7bbef786..61b715bb 100644 --- a/coretrace/simulation_runner/native_runner/trace.cpp +++ b/coretrace/simulation_runner/native_runner/trace.cpp @@ -65,6 +65,7 @@ #include // NativeRunner headers +#include "determine_interaction_type.hpp" #include "find_element_hit.hpp" #include "generate_ray.hpp" #include "native_runner_types.hpp" @@ -72,6 +73,7 @@ #include "pt_optimizations.hpp" #include "sun_to_primary_stage.hpp" #include "thread_manager.hpp" +#include "trace_logger.hpp" #include "treemesh.hpp" namespace SolTrace::NativeRunner @@ -83,9 +85,10 @@ namespace SolTrace::NativeRunner // Trace method RunnerStatus trace_native( thread_manager_ptr manager, + trace_logger_ptr logger, TSystem *System, const std::vector &seeds, - uint_fast64_t nthreads, + unsigned nthreads, uint_fast64_t NumberOfRays, uint_fast64_t MaxNumberOfRays, bool IncludeSunShape, @@ -94,7 +97,7 @@ namespace SolTrace::NativeRunner { // Initialize Sun Vector3d PosSunStage; - if (!SunToPrimaryStage(manager, + if (!SunToPrimaryStage(logger, System, System->StageList[0].get(), &System->Sun, @@ -124,6 +127,7 @@ namespace SolTrace::NativeRunner // having trouble with all the arguments... ThreadInfo my_info; my_info.manager = manager; + my_info.logger = logger; my_info.System = System; // my_info.NumberOfRays = NumberOfRays / nthreads; uint_fast64_t rem = NumberOfRays % nthreads; @@ -163,6 +167,7 @@ namespace SolTrace::NativeRunner RunnerStatus trace_single_thread( unsigned thread_id, thread_manager_ptr manager, + trace_logger_ptr logger, TSystem *System, unsigned int seed, uint_fast64_t NumberOfRays, @@ -208,11 +213,6 @@ namespace SolTrace::NativeRunner std::vector IncomingRays; // Vector of rays from previous stage, going into next stage IncomingRays.resize(NumberOfRays); - // Start the clock - // clock_t startTime = clock(); - // int rays_per_callback_estimate = 50; - // uint_fast64_t RaysTracedTotal = 0; - // Initialize stage variables uint_fast64_t StageDataArrayIndex = 0; uint_fast64_t PreviousStageDataArrayIndex = 0; @@ -363,7 +363,7 @@ namespace SolTrace::NativeRunner std::stringstream ss; ss << "Thread " << thread_id << " failed to record ray data.\n"; - manager->error_log(ss.str()); + logger->error_log(ss.str()); } } @@ -379,98 +379,31 @@ namespace SolTrace::NativeRunner else { // trace through the interaction - telement_ptr optelm = Stage->ElementList[LastElementNumber - 1]; + telement_ptr optelm = + Stage->ElementList[LastElementNumber - 1]; if (LastHitBackSide) optics = &optelm->Optics.Back; else optics = &optelm->Optics.Front; - double TestValue; - double UnitLastDFXYZ[3] = {0.0, 0.0, 0.0}; - double IncidentAngle = 0; - // switch (optelm->InteractionType) - switch (optics->my_type) + bool good = determine_interaction_type( + logger, + i, + 0, + myrng, + optics, + LastDFXYZ, + LastCosRaySurfElement, + rev); + + if (!good) { - case InteractionType::REFRACTION: // refraction - // TODO: Implement transmissivity table? - // if (optics->UseTransmissivityTable) - // { - // int npoints = optics->TransmissivityTable.size(); - // int m = 0; - - // UnitLastDFXYZ[0] = -LastDFXYZ[0] / sqrt(DOT(LastDFXYZ, LastDFXYZ)); - // UnitLastDFXYZ[1] = -LastDFXYZ[1] / sqrt(DOT(LastDFXYZ, LastDFXYZ)); - // UnitLastDFXYZ[2] = -LastDFXYZ[2] / sqrt(DOT(LastDFXYZ, LastDFXYZ)); - // IncidentAngle = acos(DOT(LastCosRaySurfElement, UnitLastDFXYZ)) * 1000.; //[mrad] - // if (IncidentAngle >= optics->TransmissivityTable[npoints - 1].angle) - // { - // TestValue = optics->TransmissivityTable[npoints - 1].trans; - // } - // else - // { - // while (optics->TransmissivityTable[m].angle < IncidentAngle) - // m++; - - // if (m == 0) - // TestValue = optics->TransmissivityTable[m].trans; - // else - // TestValue = (optics->TransmissivityTable[m].trans + optics->TransmissivityTable[m - 1].trans) / 2.0; - // } - // } - // else - // TestValue = optics->Transmissivity; - TestValue = optics->transmitivity; - rev = RayEvent::TRANSMIT; - break; - case InteractionType::REFLECTION: // reflection - // TODO: Implement reflectivity table? - // if (optics->UseReflectivityTable) - // { - // int npoints = optics->ReflectivityTable.size(); - // int m = 0; - // UnitLastDFXYZ[0] = -LastDFXYZ[0] / sqrt(DOT(LastDFXYZ, LastDFXYZ)); - // UnitLastDFXYZ[1] = -LastDFXYZ[1] / sqrt(DOT(LastDFXYZ, LastDFXYZ)); - // UnitLastDFXYZ[2] = -LastDFXYZ[2] / sqrt(DOT(LastDFXYZ, LastDFXYZ)); - // IncidentAngle = acos(DOT(LastCosRaySurfElement, UnitLastDFXYZ)) * 1000.; //[mrad] - // if (IncidentAngle >= optics->ReflectivityTable[npoints - 1].angle) - // { - // TestValue = optics->ReflectivityTable[npoints - 1].refl; - // } - // else - // { - // while (optics->ReflectivityTable[m].angle < IncidentAngle) - // m++; - - // if (m == 0) - // TestValue = optics->ReflectivityTable[m].refl; - // else - // TestValue = (optics->ReflectivityTable[m].refl + optics->ReflectivityTable[m - 1].refl) / 2.0; - // } - // } - // else - // TestValue = optics->Reflectivity; - TestValue = optics->reflectivity; - rev = RayEvent::REFLECT; - break; - default: - std::stringstream ss; - ss << "Bad optical interaction." - << " Type: " << static_cast(optics->my_type) - << " Stage: " << i - << " Thread: " << thread_id - << "\n"; - manager->error_log(ss.str()); return RunnerStatus::ERROR; } - // Apply MonteCarlo probability of absorption. Limited - // for now, but can make more complex later on if desired - // if (TestValue <= myrng()) - double flip = myrng(); - if (TestValue <= flip) + if (rev == RayEvent::ABSORB) { - // ray was fully absorbed RayIsAbsorbed = true; break; } @@ -483,7 +416,8 @@ namespace SolTrace::NativeRunner IncludeSunShape, optics, IncludeErrors, - i, Stage, // k, + i, + Stage, MultipleHitCount, LastDFXYZ, LastCosRaySurfElement, @@ -526,7 +460,9 @@ namespace SolTrace::NativeRunner } } - ++update_count; + if (MultipleHitCount > 0) + ++update_count; + if (update_count % update_rate == 0) { double progress = update_count / total_work; diff --git a/coretrace/simulation_runner/native_runner/trace.hpp b/coretrace/simulation_runner/native_runner/trace.hpp index 99c10598..31b79780 100644 --- a/coretrace/simulation_runner/native_runner/trace.hpp +++ b/coretrace/simulation_runner/native_runner/trace.hpp @@ -10,133 +10,17 @@ #include "mtrand.hpp" #include "native_runner_types.hpp" #include "thread_manager.hpp" +#include "trace_logger.hpp" #include "treemesh.hpp" namespace SolTrace::NativeRunner { - - class GlobalRay_refactored - { - public: - GlobalRay_refactored() // : active(true) - { - Num = 0; - for (int i = 0; i < 3; i++) - Pos[i] = Cos[i] = 0.0; - } - - double Pos[3]; - double Cos[3]; - uint_fast64_t Num; - // bool active; - }; - - // void SpencerandMurtySurfaceClosedForm(TElement *Element, - // double PosLoc[3], - // double CosLoc[3], - // double PosXYZ[3], - // double DFXYZ[3], - // double *PathLength, - // int *ErrorFlag); - - // bool LoadExistingStage0Ray(int index, - // std::vector> *raydat, - // double PosRayGlobal[3], - // double CosRayGlobal[3], - // st_uint_t &ElementNum, - // st_uint_t &RayNum); - - // bool LoadExistingStage1Ray(int index, - // std::vector> *raydat, - // double PosRayGlobal[3], - // double CosRayGlobal[3], - // int &raynum); - - // bool AperturePlane(TElement *Element); - - // void NewZStartforCubicSplineSurf(double CRadius, - // double PosLoc[3], - // double CosLoc[3], - // char AperShapeIndex, - // double *NewZStart, - // double *PLength, - // int *EFlag); - - // // the 0.0's are values for DeltaX and DeltaY; - // **[need to look at this further]** - // void EvalPoly(double ax, - // double ay, - // std::vector &Coeffs, - // int POrder, - // double *az); - - // void PolySlope(std::vector &Coeffs, - // int POrder, - // double ax, - // double ay, - // double *dzdx, - // double *dzdy); - - // void spline(std::vector &x, - // std::vector &y, - // int n, - // double yp1, double ypn, - // std::vector &y2); - - // void EvalMono(double ax, - // double ay, - // HPM2D &B, - // int order, - // double DeltaX, - // double DeltaY, - // double *az); - - // void FEInterpolate(double Xray, double Yray, double Delta, double Density, - // HPM2D &FEData, int NumFEPoints, - // double *z, double *zx, double *zy); - - // void MonoSlope(HPM2D &B, - // int order, - // double sxp, - // double syp, - // double *dzdx, - // double *dzdy); - - // void VSHOTInterpolateModShepard(double Xray, - // double Yray, double Density, - // HPM2D &VSHOTData, - // int NumVSHOTPoints, - // double *zx, - // double *zy, - // int *ErrorFlag); - - // void FEInterpNew(double Xray, - // double Yray, - // double Density, - // HPM2D &FEData, - // int NumFEPoints, - // double *zr); - // void FEInterpGM(double Xray, double Yray, GaussMarkov *gm, double *zr); - // void FEInterpKD(double Xray, - // double Yray, - // FEDataObj *kd, - // double step, - // double *zr, - // double *dzrdx, - // double *dzrdy); - - // bool InitGeometries(TSystem *sys); - // bool TranslateSurfaceParams( TElement *elm, double params[8]); - // bool ReadSurfaceFile( const char *file, TElement *elm ); - - // bool TranslateSurfaceParams(TSystem *sys, TElement *elm, double params[8]); - // bool ReadSurfaceFile(const char *file, TElement *elm, TSystem *sys); - SolTrace::Runner::RunnerStatus trace_native( thread_manager_ptr manager, + trace_logger_ptr logger, TSystem *System, const std::vector &seeds, - uint_fast64_t nthreads, + unsigned nthreads, uint_fast64_t NumberOfRays, uint_fast64_t MaxNumberOfRays, bool IncludeSunShape, @@ -146,6 +30,7 @@ namespace SolTrace::NativeRunner SolTrace::Runner::RunnerStatus trace_single_thread( unsigned thread_id, thread_manager_ptr manager, + trace_logger_ptr logger, TSystem *System, unsigned seed, uint_fast64_t NumberOfRays, @@ -161,6 +46,7 @@ namespace SolTrace::NativeRunner struct ThreadInfo { thread_manager_ptr manager; + trace_logger_ptr logger; TSystem *System; // unsigned int seed; uint_fast64_t NumberOfRays; @@ -183,6 +69,7 @@ namespace SolTrace::NativeRunner return trace_single_thread( thread_id, info.manager, + info.logger, info.System, seed, info.NumberOfRays, diff --git a/coretrace/simulation_runner/native_runner/trace_logger.cpp b/coretrace/simulation_runner/native_runner/trace_logger.cpp new file mode 100644 index 00000000..0817ab7f --- /dev/null +++ b/coretrace/simulation_runner/native_runner/trace_logger.cpp @@ -0,0 +1,37 @@ +#include "trace_logger.hpp" + +#include + +namespace SolTrace::NativeRunner +{ + TraceLogger::TraceLogger() + { + return; + } + + TraceLogger::~TraceLogger() + { + this->messages.clear(); + return; + } + + void TraceLogger::error_log(const std::string &msg) + { + std::lock_guard lk(this->message_mutex); + this->messages.push_back(msg); + return; + } + + void TraceLogger::print_log(std::ostream &os) const + { + std::lock_guard lk(this->message_mutex); + for (auto iter = this->messages.cbegin(); + iter != this->messages.cend(); + ++iter) + { + os << *iter; + } + return; + } + +} // namespace SolTrace::NativeRunner diff --git a/coretrace/simulation_runner/native_runner/trace_logger.hpp b/coretrace/simulation_runner/native_runner/trace_logger.hpp new file mode 100644 index 00000000..3d353aa3 --- /dev/null +++ b/coretrace/simulation_runner/native_runner/trace_logger.hpp @@ -0,0 +1,37 @@ +#ifndef SOLTRACE_TRACE_LOGGER_H +#define SOLTRACE_TRACE_LOGGER_H + +#include +#include +#include +#include +#include + +namespace SolTrace::NativeRunner +{ + + class TraceLogger + { + public: + TraceLogger(); + ~TraceLogger(); + + void error_log(const std::string &msg); + void print_log(std::ostream &os) const; + + private: + mutable std::mutex message_mutex; + std::vector messages; + }; + + using trace_logger_ptr = std::shared_ptr; + + template + inline auto make_trace_logger(Args &&...args) + { + return std::make_shared(std::forward(args)...); + } + +} // namespace SolTrace::NativeRunner + +#endif diff --git a/google-tests/regression-tests/CMakeLists.txt b/google-tests/regression-tests/CMakeLists.txt index 3801a996..8909cbff 100644 --- a/google-tests/regression-tests/CMakeLists.txt +++ b/google-tests/regression-tests/CMakeLists.txt @@ -5,11 +5,13 @@ include_directories( . # ../app/deploy/samples ../test_tools + ../unit-tests/common ../../coretrace/simulation_data ../../coretrace/simulation_data/cst_templates ../../coretrace/simulation_results ../../coretrace/simulation_runner ../../coretrace/simulation_runner/native_runner + ../../coretrace/simulation_runner/embree_runner ) # file(GLOB UNIT_TEST_SRC *.cpp) @@ -29,6 +31,13 @@ if(CMAKE_BUILD_TYPE STREQUAL "Release" OR SOLTRACE_BUILD_PERF_TEST) ) endif() +if(SOLTRACE_BUILD_EMBREE_SUPPORT) + list(APPEND + REGRESSION_TEST_SRC + embree_runner_validation_test.cpp + ) +endif() + add_executable(RegressionTests ${CSV_TOOL} ${REGRESSION_TEST_SRC} @@ -43,6 +52,13 @@ target_link_libraries( GTest::gtest_main ) +if(SOLTRACE_BUILD_EMBREE_SUPPORT) + target_link_libraries( + RegressionTests + embree_runner + ) +endif() + if (ENABLE_COVERAGE AND CMAKE_BUILD_TYPE STREQUAL "Debug") target_link_options(RegressionTests PRIVATE --coverage) target_link_libraries(RegressionTests gcov) diff --git a/google-tests/regression-tests/embree_runner_validation_test.cpp b/google-tests/regression-tests/embree_runner_validation_test.cpp new file mode 100644 index 00000000..a0000957 --- /dev/null +++ b/google-tests/regression-tests/embree_runner_validation_test.cpp @@ -0,0 +1,687 @@ +#include + +#include +#include + +#include +#include +#include + +#include "common.hpp" +#include "split_csv.h" + +using SolTrace::Runner::RunnerStatus; + +using SolTrace::EmbreeRunner::EmbreeRunner; +using SolTrace::EmbreeRunner::TRayData; +using SolTrace::EmbreeRunner::tstage_ptr; +using SolTrace::EmbreeRunner::TSystem; + +using SolTrace::Result::RayEvent; + +bool get_runner_element_and_stage(const EmbreeRunner *runner, + element_id id, + int_fast64_t &re, + int_fast64_t &rs) +{ + // Assumes that `id` is a unique identifier and that an element + // can belong only to a single stage + bool found = false; + const TSystem *sys = runner->get_system(); + for (unsigned stage_idx = 0; + stage_idx < sys->StageList.size(); + ++stage_idx) + { + tstage_ptr stage = sys->StageList[stage_idx]; + for (unsigned element_idx = 0; + element_idx < stage->ElementList.size(); + ++element_idx) + { + SolTrace::EmbreeRunner::telement_ptr tel = + stage->ElementList[element_idx]; + if (tel->sim_data_id == id) + { + // EmbreeRunner stage and element ids are 1 based + // so we add 1 to the vector index here. + found = true; + re = element_idx + 1; + rs = stage_idx + 1; + } + } + + if (found) + break; + } + + return found; +} + +// int_fast64_t count_element_event(const SimulationResult &res, element_id el, RayEvent rev) +// { +// int_fast64_t count = 0; + +// for (auto ray_idx = 0; +// ray_idx < res.get_number_of_records(); +// ++ray_idx) +// { +// auto rr = res[ray_idx]; +// for (auto event_idx = 0; +// event_idx < rr->get_number_of_interactions(); +// ++event_idx) +// { +// if (rr->get_event(event_idx) == rev && +// rr->get_element(event_idx) == el) +// { +// ++count; +// } +// } +// } + +// return count; +// } + +TEST(EmbreeRunner, ValidationTest1) +{ + // Pulling in path variable from CMake and creating path to .stinput sample file + std::string sample_path = std::string(PROJECT_DIR) + std::string("/High Flux Solar Furnace.stinput"); + + // Path to .csv exported from Soltrace as ground truth + std::string ground_csv_path = PROJECT_DIR + std::string("/hfsf_example_raydata.csv"); + + std::ifstream csv_file(ground_csv_path); + std::vector> ground_raydata = split_csv(ground_csv_path); + + // Create Simuluation Data + SimulationData sd; + + // Constants + const uint_fast64_t NRAYS = 10000; + const double TOL = 1e-4; + + // Read Input File + bool success = sd.import_from_file(sample_path); + EXPECT_TRUE(success); + EXPECT_TRUE(sd.get_number_of_elements() > 0); + EXPECT_TRUE(sd.get_number_of_ray_sources() > 0); + + // Parameters + SimulationParameters ¶ms = sd.get_simulation_parameters(); + params.include_optical_errors = false; + params.include_sun_shape_errors = false; + params.max_number_of_rays = NRAYS * 100; + params.number_of_rays = NRAYS; + params.seed = 1; + + // Run Ray Trace + EmbreeRunner runner; + runner.set_number_of_threads(1); + RunnerStatus sts; + sts = runner.initialize(); + EXPECT_EQ(sts, RunnerStatus::SUCCESS); + sts = runner.setup_simulation(&sd); + EXPECT_EQ(sts, RunnerStatus::SUCCESS); + + auto stage_list = runner.get_system()->StageList; + tstage_ptr tstage = nullptr; + for (auto iter = stage_list.begin(); + iter != stage_list.end(); + ++iter) + { + tstage = *iter; + tstage->MultiHitsPerRay = false; + tstage->TraceThrough = false; + } + + sts = runner.run_simulation(); + EXPECT_EQ(sts, RunnerStatus::SUCCESS); + + // const TSystem *sys = runner.get_system(); + // const TRayData *ray_data = &(runner.get_system()->RayData); + // size_t nrdata = ray_data->Count(); + + // ray_data->Print(); + + SimulationResult result; + sts = runner.report_simulation(&result, 0); + EXPECT_EQ(sts, RunnerStatus::SUCCESS); + EXPECT_EQ(result.get_number_of_records(), NRAYS); + + Vector3d point, cosines; + Vector3d pos_stage, dir_stage; + Vector3d temp; + int_fast64_t element; + int_fast64_t stage; + uint_fast64_t rayidx; + uint_fast64_t iidx; + element_ptr el = nullptr; + int_fast64_t run_element, run_stage; + + // size_t total_lines = ground_raydata[0].size(); + size_t total_lines = std::numeric_limits::max(); + size_t i; + for (i = 0; i < 9; ++i) + { + total_lines = std::min(total_lines, ground_raydata[i].size()); + } + + // Compare saved CSV values to runner values + for (i = 1; i < total_lines; ++i) + { + element = stoi(ground_raydata[6][i]); + stage = stoi(ground_raydata[7][i]); + // Legacy SolTrace and CSV file had 1-based ray IDs. SimulationResult + // has 0-based ray ID's so subtract 1 here. + rayidx = stoul(ground_raydata[8][i]) - 1; + + const ray_record_ptr rr = result[rayidx]; + if (element > 0) + { + bool found = false; + for (uint_fast64_t idx = 0; idx < rr->get_number_of_interactions(); ++idx) + { + run_element = -1; + run_stage = -1; + bool sts = get_runner_element_and_stage(&runner, + rr->get_element(idx), + run_element, + run_stage); + if (sts && element == run_element && stage == run_stage) + { + iidx = idx; + found = true; + break; + } + } + + // assert(found); + EXPECT_TRUE(found); + + if (found) + { + EXPECT_NE(rr->get_event(iidx), RayEvent::CREATE); + EXPECT_NE(rr->get_event(iidx), RayEvent::ABSORB); + EXPECT_NE(rr->get_event(iidx), RayEvent::EXIT); + } + else + { + std::cout << "CSV Line: " << i + 1 + << "Ray Record: " << *rr + << std::endl; + break; + } + + rr->get_position(iidx, point); + // Legacy SolTrace stored the incoming ray direction whereas + // EmbreeRunner/SimulationResult stores the exit direction so + // we take the direction for the previous ray event. + rr->get_direction(iidx - 1, cosines); + } + else + { + iidx = rr->get_number_of_interactions() - 1; + if (element == 0) + { + // Ray miss -- check only that it was noted + EXPECT_EQ(rr->get_event(iidx), RayEvent::EXIT); + continue; + } + else + { + EXPECT_EQ(rr->get_event(iidx), RayEvent::ABSORB); + get_runner_element_and_stage(&runner, + rr->get_element(iidx), + run_element, + run_stage); + EXPECT_EQ(run_element, abs(element)); + EXPECT_EQ(run_stage, stage); + + rr->get_position(iidx, point); + // Legacy SolTrace stored the incoming ray direction whereas + // EmbreeRunner/SimulationResult stores the exit direction so + // we take the direction for the previous ray event. + rr->get_direction(iidx - 1, cosines); + } + } + + el = sd.get_element(rr->get_element(iidx)); + EXPECT_NE(el, nullptr); + + if (el == nullptr) + { + std::cout << "CSV Line: " << i + 1 + << "\nRay Number: " << rayidx + 1 + << "\nElement: " << rr->get_element(iidx) + << " CSV Element: " << element + << " Runner Element: " << run_element + << "\nCSV Stage: " << stage + << " Runner Stage: " << run_stage + << std::endl; + break; + } + + // Runner and SimulationResult store everything in global + // coordinate whereas the CSV file is in stage coordinates + // as per legacy SolTrace + el->convert_global_to_reference(pos_stage, point); + + // See previous comment about coordinates + el->convert_vector_global_to_reference(dir_stage, cosines); + + EXPECT_NEAR(pos_stage[0], stod(ground_raydata[0][i]), TOL); + EXPECT_NEAR(pos_stage[1], stod(ground_raydata[1][i]), TOL); + EXPECT_NEAR(pos_stage[2], stod(ground_raydata[2][i]), TOL); + + EXPECT_NEAR(dir_stage[0], stod(ground_raydata[3][i]), TOL); + EXPECT_NEAR(dir_stage[1], stod(ground_raydata[4][i]), TOL); + EXPECT_NEAR(dir_stage[2], stod(ground_raydata[5][i]), TOL); + + // if (fabs(pos_stage[0] - stod(ground_raydata[0][i]) > TOL)) + // { + // Vector3d pos_csv(stod(ground_raydata[0][i]), + // stod(ground_raydata[1][i]), + // stod(ground_raydata[2][i])); + // Vector3d dir_csv(stod(ground_raydata[3][i]), + // stod(ground_raydata[4][i]), + // stod(ground_raydata[5][i])); + + // Vector3d pos_loc; + // Vector3d dir_loc; + + // Vector3d csv_glob; + // el->convert_stage_to_local(temp, pos_csv); + // el->convert_local_to_global(csv_glob, temp); + + // el->convert_global_to_local(pos_loc, point); + // el->convert_vector_global_to_local(dir_loc, cosines); + + // std::cout << "CSV Line: " << i + 1 + // << "\nRay Number: " << rayidx + 1 + // << "\nElement: " << rr->get_element(iidx) + // << " CSV Element: " << element + // << " Runner Element: " << run_element + // << "\nCSV Stage: " << stage + // << " Runner Stage: " << run_stage + // << std::endl; + + // std::cout << "CSV Pos: " << pos_csv + // << "\nCSV Global: " << csv_glob + // << "\nPos Global: " << point + // << "\nPos Stage: " << pos_stage + // << "\nPos Local: " << pos_loc + // << "\nCSV Dir: " << dir_csv + // << "\nDir Global: " << cosines + // << "\nDir Stage: " << dir_stage + // << "\nDir Local: " << dir_loc + // << std::endl; + + // std::cout << "Ray Record: " + // << *rr + // << std::endl; + + // std::cout << "Element Global Origin: " << el->get_origin_global() + // << "\nElement Local to Global: " << el->get_local_to_global() + // << "\nElement Ref to Local: " << el->get_local_to_reference() + // << std::endl; + + // break; + // } + } +} + +TEST(EmbreeRunner, ValidationTest2) +{ + // Pulling in path variable from CMake and creating path to .stinput sample file + std::string sample_path = std::string(PROJECT_DIR) + std::string("/Power-tower-surround_singlefacet.stinput"); + + // Path to .csv exported from Soltrace as ground truth + std::string ground_csv_path = PROJECT_DIR + std::string("/powertower_example_raydata.csv"); + + std::ifstream csv_file(ground_csv_path); + std::vector> ground_raydata = split_csv(ground_csv_path); + + // Create Simuluation Data + SimulationData sd; + + // Constants + const uint_fast64_t NRAYS = 50000; + const double TOL = 1e-4; + + // Read Input File + bool success = sd.import_from_file(sample_path); + EXPECT_TRUE(success); + EXPECT_TRUE(sd.get_number_of_elements() > 0); + EXPECT_TRUE(sd.get_number_of_ray_sources() > 0); + + std::cout << "Num Elements: " << sd.get_number_of_elements() << std::endl; + + // Parameters + SimulationParameters ¶ms = sd.get_simulation_parameters(); + params.include_optical_errors = false; + params.include_sun_shape_errors = false; + params.max_number_of_rays = NRAYS * 100; + params.number_of_rays = NRAYS; + params.seed = 1; + + // for (auto cit = sd.get_const_iterator(); + // !sd.is_at_end(cit); + // ++cit) + // { + // auto elem = cit->second; + // if (elem->is_stage()) + // { + // stage_ptr st = dynamic_pointer_cast(elem); + // assert(st != nullptr); + // std::cout << "Stage: " << st->get_stage() + // << "\nNumber of elements: " + // << st->get_number_of_elements() + // << std::endl; + // if (st->get_stage() == 1) + // { + // auto el = st->get_const_iterator()->second; + + // auto surf = dynamic_pointer_cast(el->get_surface()); + // assert(surf != nullptr); + // auto ap = dynamic_pointer_cast(el->get_aperture()); + // assert(ap != nullptr); + + // std::cout << "------------\n" + // << "Element ID: " << el->get_id() + // << "\nElement name: " << el->get_name() + // << "\nIs Stage: " << el->is_stage() + // << "\nIs Composite: " << el->is_composite() + // << "\nIs Single: " << el->is_single() + // << "\nOrigin (ref): " << el->get_origin_ref() + // << "\nOrigin (stage): " << el->get_origin_stage() + // << "\nOrigin (global): " << el->get_origin_global() + // << "\nAim (ref): " << el->get_aim_vector_ref() + // << "\nAim (stage): " << el->get_aim_vector_stage() + // << "\nAim (global): " << el->get_aim_vector_global() + // << "\nRefToLoc: " << el->get_reference_to_local() + // << "\nCylinder Radius: " << surf->radius + // << "\nAperture Dimension: " << ap->x_length << " x " << ap->y_length + // << " at (" << ap->x_coord << ", " << ap->y_coord << ")" + // << "\n** Front Optical Properties **\n" + // << *(el->get_front_optical_properties()) + // << "\n** Back Optical Properties **\n" + // << *(el->get_back_optical_properties()) + // << "\n------------" + // << std::endl; + // } + // } + // } + + // return; + + // Run Ray Trace + EmbreeRunner runner; + runner.set_number_of_threads(1); + RunnerStatus sts; + sts = runner.initialize(); + EXPECT_EQ(sts, RunnerStatus::SUCCESS); + sts = runner.setup_simulation(&sd); + EXPECT_EQ(sts, RunnerStatus::SUCCESS); + + // auto stage_list = runner.get_system()->StageList; + // stage_list[1]->MultiHitsPerRay = false; + // tstage_ptr tstage = nullptr; + // for (auto iter = stage_list.begin(); + // iter != stage_list.end(); + // ++iter) + // { + // tstage = *iter; + // tstage->MultiHitsPerRay = false; + // tstage->TraceThrough = false; + // } + + auto t0 = std::chrono::high_resolution_clock::now(); + sts = runner.run_simulation(); + auto t1 = std::chrono::high_resolution_clock::now(); + EXPECT_EQ(sts, RunnerStatus::SUCCESS); + + std::chrono::duration dur = t1 - t0; + EXPECT_TRUE(dur.count() < 75000.0); + + std::cout << "Time: " << dur.count() << " ms" << std::endl; + + // const TSystem *sys = runner.get_system(); + // const TRayData *ray_data = &(runner.get_system()->RayData); + // size_t nrdata = ray_data->Count(); + + // ray_data->Print(); + + SimulationResult result; + sts = runner.report_simulation(&result, 0); + EXPECT_EQ(sts, RunnerStatus::SUCCESS); + EXPECT_EQ(result.get_number_of_records(), NRAYS); + + element_id absorber_id = 6285; + int_fast64_t nabsorbed = count_element_event(result, absorber_id, RayEvent::ABSORB); + int_fast64_t nreflect = count_element_event(result, absorber_id, RayEvent::REFLECT); + int_fast64_t nevents = nabsorbed + nreflect; + + std::cout << "Total: " << nevents + << "\nAbsorb: " << nabsorbed << " (" + << static_cast(nabsorbed) / nevents << ")" + << "\nReflect: " << nreflect << " (" + << static_cast(nreflect) / nevents << ")" + << std::endl; + + // result.write_csv_file("native_runner_result_dump.csv"); + + Vector3d point, cosines; + Vector3d pos_stage, dir_stage; + Vector3d temp; + int_fast64_t element; + int_fast64_t stage; + uint_fast64_t rayidx; + uint_fast64_t iidx; + element_ptr el = nullptr; + int_fast64_t run_element, run_stage; + double RTOL = 0.0; + + // size_t total_lines = ground_raydata[0].size(); + size_t total_lines = std::numeric_limits::max(); + size_t i; + for (i = 0; i < 9; ++i) + { + total_lines = std::min(total_lines, ground_raydata[i].size()); + } + + // Compare saved CSV values to runner values + for (i = 1; i < total_lines; ++i) + { + element = stoi(ground_raydata[6][i]); + stage = stoi(ground_raydata[7][i]); + // Legacy SolTrace and CSV file had 1-based ray IDs. SimulationResult + // has 0-based ray ID's so subtract 1 here. + rayidx = stoul(ground_raydata[8][i]) - 1; + + const ray_record_ptr rr = result[rayidx]; + if (element > 0) + { + bool found = false; + for (uint_fast64_t idx = 0; idx < rr->get_number_of_interactions(); ++idx) + { + run_element = -1; + run_stage = -1; + bool sts = get_runner_element_and_stage(&runner, + rr->get_element(idx), + run_element, + run_stage); + if (sts && element == run_element && stage == run_stage) + { + iidx = idx; + found = true; + break; + } + } + + // assert(found); + EXPECT_TRUE(found); + + if (found) + { + EXPECT_NE(rr->get_event(iidx), RayEvent::CREATE); + EXPECT_NE(rr->get_event(iidx), RayEvent::ABSORB); + EXPECT_NE(rr->get_event(iidx), RayEvent::EXIT); + } + + if (!found || + rr->get_event(iidx) == RayEvent::CREATE || + rr->get_event(iidx) == RayEvent::ABSORB || + rr->get_event(iidx) == RayEvent::EXIT) + { + std::cout << "CSV Line: " << i + 1 + << "\nRay Record: " << *rr + << std::endl; + break; + } + + rr->get_position(iidx, point); + // Legacy SolTrace stored the incoming ray direction whereas + // EmbreeRunner/SimulationResult stores the exit direction so + // we take the direction for the previous ray event. + rr->get_direction(iidx - 1, cosines); + } + else + { + iidx = rr->get_number_of_interactions() - 1; + if (element == 0) + { + // Ray miss -- check only that it was noted + EXPECT_EQ(rr->get_event(iidx), RayEvent::EXIT); + continue; + } + else + { + EXPECT_EQ(rr->get_event(iidx), RayEvent::ABSORB); + + get_runner_element_and_stage(&runner, + rr->get_element(iidx), + run_element, + run_stage); + EXPECT_EQ(run_element, abs(element)); + EXPECT_EQ(run_stage, stage); + + if (run_element != abs(element) || run_stage != stage) + { + std::cout << "CSV Line: " << i + 1 + << "\nRay Number: " << rayidx + 1 + << "\nElement: " << rr->get_element(iidx) + << " CSV Element: " << element + << " Runner Element: " << run_element + << "\nCSV Stage: " << stage + << " Runner Stage: " << run_stage + << "\nRay Record: " << *rr + << std::endl; + break; + } + + rr->get_position(iidx, point); + // Legacy SolTrace stored the incoming ray direction whereas + // EmbreeRunner/SimulationResult stores the exit direction so + // we take the direction for the previous ray event. + rr->get_direction(iidx - 1, cosines); + } + } + + el = sd.get_element(rr->get_element(iidx)); + EXPECT_NE(el, nullptr); + + if (el == nullptr) + { + std::cout << "CSV Line: " << i + 1 + << "\nRay Number: " << rayidx + 1 + << "\nElement: " << rr->get_element(iidx) + << " CSV Element: " << element + << " Runner Element: " << run_element + << "\nCSV Stage: " << stage + << " Runner Stage: " << run_stage + << "\nRay Record: " << *rr + << std::endl; + break; + } + + // Runner and SimulationResult store everything in global + // coordinate whereas the CSV file is in stage coordinates + // as per legacy SolTrace + el->convert_global_to_reference(pos_stage, point); + + // See previous comment about coordinates + el->convert_vector_global_to_reference(dir_stage, cosines); + + RTOL = vector_norm(pos_stage) * TOL; + EXPECT_NEAR(pos_stage[0], stod(ground_raydata[0][i]), RTOL); + EXPECT_NEAR(pos_stage[1], stod(ground_raydata[1][i]), RTOL); + EXPECT_NEAR(pos_stage[2], stod(ground_raydata[2][i]), RTOL); + + EXPECT_NEAR(dir_stage[0], stod(ground_raydata[3][i]), TOL); + EXPECT_NEAR(dir_stage[1], stod(ground_raydata[4][i]), TOL); + EXPECT_NEAR(dir_stage[2], stod(ground_raydata[5][i]), TOL); + + if (fabs(pos_stage[0] - stod(ground_raydata[0][i])) > RTOL) + { + std::cout << "CSV Line: " << i + 1 + << "\nRay Number: " << rayidx + 1 + << "\nElement: " << rr->get_element(iidx) + << " CSV Element: " << element + << " Runner Element: " << run_element + << "\nCSV Stage: " << stage + << " Runner Stage: " << run_stage + << "\nRay Record: " << *rr + << std::endl; + break; + } + + // if (fabs(pos_stage[0] - stod(ground_raydata[0][i]) > TOL)) + // { + // Vector3d pos_csv(stod(ground_raydata[0][i]), + // stod(ground_raydata[1][i]), + // stod(ground_raydata[2][i])); + // Vector3d dir_csv(stod(ground_raydata[3][i]), + // stod(ground_raydata[4][i]), + // stod(ground_raydata[5][i])); + + // Vector3d pos_loc; + // Vector3d dir_loc; + + // Vector3d csv_glob; + // el->convert_stage_to_local(temp, pos_csv); + // el->convert_local_to_global(csv_glob, temp); + + // el->convert_global_to_local(pos_loc, point); + // el->convert_vector_global_to_local(dir_loc, cosines); + + // std::cout << "CSV Line: " << i + 1 + // << "\nRay Number: " << rayidx + 1 + // << "\nElement: " << rr->get_element(iidx) + // << " CSV Element: " << element + // << " Runner Element: " << run_element + // << "\nCSV Stage: " << stage + // << " Runner Stage: " << run_stage + // << std::endl; + + // std::cout << "CSV Pos: " << pos_csv + // << "\nCSV Global: " << csv_glob + // << "\nPos Global: " << point + // << "\nPos Stage: " << pos_stage + // << "\nPos Local: " << pos_loc + // << "\nCSV Dir: " << dir_csv + // << "\nDir Global: " << cosines + // << "\nDir Stage: " << dir_stage + // << "\nDir Local: " << dir_loc + // << std::endl; + + // std::cout << "Ray Record: " + // << *rr + // << std::endl; + + // std::cout << "Element Global Origin: " << el->get_origin_global() + // << "\nElement Local to Global: " << el->get_local_to_global() + // << "\nElement Ref to Local: " << el->get_local_to_reference() + // << std::endl; + + // break; + // } + } +} diff --git a/google-tests/regression-tests/native_runner_multithreading_performance_test.cpp b/google-tests/regression-tests/native_runner_multithreading_performance_test.cpp index c1f8c463..03a0f1bf 100644 --- a/google-tests/regression-tests/native_runner_multithreading_performance_test.cpp +++ b/google-tests/regression-tests/native_runner_multithreading_performance_test.cpp @@ -8,6 +8,8 @@ #include #include +#include "common.hpp" + TEST(NativeRunner, Multithreading) { const uint_fast64_t NRAYS = 50000; @@ -59,15 +61,23 @@ TEST(NativeRunner, Multithreading) EXPECT_TRUE(dur.count() < 20000.0); - // element_id absorber_id = 6285; - // int_fast64_t nabsorbed = count_element_event(result, absorber_id, RayEvent::ABSORB); - // int_fast64_t nreflect = count_element_event(result, absorber_id, RayEvent::REFLECT); - // int_fast64_t nevents = nabsorbed + nreflect; - - // std::cout << "Total: " << nevents - // << "\nAbsorb: " << nabsorbed << " (" - // << static_cast(nabsorbed) / nevents << ")" - // << "\nReflect: " << nreflect << " (" - // << static_cast(nreflect) / nevents << ")" - // << std::endl; + SimulationResult result; + sts = nr.report_simulation(&result, 0); + EXPECT_EQ(sts, SolTrace::Runner::RunnerStatus::SUCCESS); + + element_id absorber_id = 6285; + int_fast64_t nabsorbed = count_element_event(result, absorber_id, RayEvent::ABSORB); + int_fast64_t nreflect = count_element_event(result, absorber_id, RayEvent::REFLECT); + int_fast64_t nevents = nabsorbed + nreflect; + + std::cout << "Total: " << nevents + << "\nAbsorb: " << nabsorbed << " (" + << static_cast(nabsorbed) / nevents << ")" + << "\nReflect: " << nreflect << " (" + << static_cast(nreflect) / nevents << ")" + << std::endl; + + EXPECT_GT(nevents, 0); + EXPECT_GT(nabsorbed, 0); + EXPECT_GT(nreflect, 0); } diff --git a/google-tests/regression-tests/native_runner_validation_test.cpp b/google-tests/regression-tests/native_runner_validation_test.cpp index 1e62f269..c5680eee 100644 --- a/google-tests/regression-tests/native_runner_validation_test.cpp +++ b/google-tests/regression-tests/native_runner_validation_test.cpp @@ -7,9 +7,9 @@ #include #include +#include "common.hpp" #include "split_csv.h" - using SolTrace::Runner::RunnerStatus; using SolTrace::NativeRunner::NativeRunner; @@ -174,14 +174,6 @@ TEST(NativeRunner, ValidationTest1) // has 0-based ray ID's so subtract 1 here. rayidx = stoul(ground_raydata[8][i]) - 1; - if (stage == 3) - { - // TODO: Stage 3 is virtual in the input file but - // that has not been implemented yet. Remove after - // implementing this. - continue; - } - const ray_record_ptr rr = result[rayidx]; if (element > 0) { @@ -210,6 +202,11 @@ TEST(NativeRunner, ValidationTest1) EXPECT_NE(rr->get_event(iidx), RayEvent::CREATE); EXPECT_NE(rr->get_event(iidx), RayEvent::ABSORB); EXPECT_NE(rr->get_event(iidx), RayEvent::EXIT); + + if (rr->get_event(iidx) == RayEvent::ABSORB) + { + break; + } } else { @@ -643,6 +640,7 @@ TEST(NativeRunner, ValidationTest2) break; } + // -------- Test Debugging Code -------- // // if (fabs(pos_stage[0] - stod(ground_raydata[0][i]) > TOL)) // { // Vector3d pos_csv(stod(ground_raydata[0][i]), diff --git a/google-tests/unit-tests/common/common.cpp b/google-tests/unit-tests/common/common.cpp index de0b32b9..08078540 100644 --- a/google-tests/unit-tests/common/common.cpp +++ b/google-tests/unit-tests/common/common.cpp @@ -3,11 +3,8 @@ #include -#include -#include -#include -#include -#include +#include +#include bool is_identical(const Vector3d &x, const Vector3d &y) { @@ -45,3 +42,27 @@ element_ptr make_configured_element() el->set_surface(SolTrace::Data::make_surface()); return el; } + +int_fast64_t count_element_event(const SimulationResult &res, element_id el, RayEvent rev) +{ + int_fast64_t count = 0; + + for (auto ray_idx = 0; + ray_idx < res.get_number_of_records(); + ++ray_idx) + { + auto rr = res[ray_idx]; + for (auto event_idx = 0; + event_idx < rr->get_number_of_interactions(); + ++event_idx) + { + if (rr->get_event(event_idx) == rev && + rr->get_element(event_idx) == el) + { + ++count; + } + } + } + + return count; +} diff --git a/google-tests/unit-tests/common/common.hpp b/google-tests/unit-tests/common/common.hpp index 7223e8e3..4624f204 100644 --- a/google-tests/unit-tests/common/common.hpp +++ b/google-tests/unit-tests/common/common.hpp @@ -1,9 +1,8 @@ #ifndef SOLTRACE_UNITTESTS_COMMON_H #define SOLTRACE_UNITTESTS_COMMON_H -#include #include -#include +#include // Vectors exactly match component-wise bool is_identical(const SolTrace::Data::Vector3d &x, @@ -65,4 +64,6 @@ inline std::shared_ptr create_rectangle_aperture(double x_length = 2. return rect; } +int_fast64_t count_element_event(const SimulationResult &res, element_id el, RayEvent rev); + #endif diff --git a/google-tests/unit-tests/simulation_data/CMakeLists.txt b/google-tests/unit-tests/simulation_data/CMakeLists.txt index 4f7b5973..2b972c37 100644 --- a/google-tests/unit-tests/simulation_data/CMakeLists.txt +++ b/google-tests/unit-tests/simulation_data/CMakeLists.txt @@ -18,6 +18,7 @@ add_subdirectory(cst-templates) set(SIMULATION_DATA_TEST_SRC ../common/common.cpp linear_algebra_test.cpp + utilities_test.cpp aperture_test.cpp surface_test.cpp container_test.cpp diff --git a/google-tests/unit-tests/simulation_data/aperture_test.cpp b/google-tests/unit-tests/simulation_data/aperture_test.cpp index 4203f47a..95c65267 100644 --- a/google-tests/unit-tests/simulation_data/aperture_test.cpp +++ b/google-tests/unit-tests/simulation_data/aperture_test.cpp @@ -6,6 +6,79 @@ #include "common.hpp" +// Bounding box test helper functions +std::vector generate_grid(double a, double b, uint_fast64_t npoints) +{ + double h = (b - a) / (npoints - 1); + std::vector grid; + grid.resize(npoints); + + for (uint_fast64_t k = 0; k < npoints; ++k) + { + grid[k] = k * h + a; + } + + return grid; +} + +bool is_in_box(double x, double y, const double xbox[2], const double ybox[2]) +{ + return xbox[0] <= x && x <= xbox[1] && + ybox[0] <= y && y <= ybox[1]; +} + +void bounding_box_test(aperture_ptr ap, + double ax, double bx, uint_fast64_t nx, + double ay, double by, uint_fast64_t ny, + uint_fast64_t &nhit_outside, + uint_fast64_t &nhit_inside) +{ + std::vector xgrid = generate_grid(ax, bx, nx); + std::vector ygrid = generate_grid(ay, by, ny); + + double xbox[2] = {0.0, 0.0}; + double ybox[2] = {0.0, 0.0}; + ap->bounding_box(xbox[0], xbox[1], ybox[0], ybox[1]); + + // std::cout << "Grid bounds: [ax, bx] X [ay, by] = " + // << "[" << ax << ", " << bx << "] X [" + // << ay << ", " << by << "]" + // << std::endl; + // std::cout << "Bounding box: [ax, bx] X [ay, by] = " + // << "[" << xbox[0] << ", " << xbox[1] << "] X [" + // << ybox[0] << ", " << ybox[1] << "]" + // << std::endl; + + nhit_outside = 0; + nhit_inside = 0; + + bool hit = false; + bool in_box = false; + for (auto xi : xgrid) + { + for (auto yj : ygrid) + { + // std::cout << "Checking grid point (x,y) = (" << xi << ", " << yj << ")" << std::endl; + hit = ap->is_in(xi, yj); + if (hit) + { + // std::cout << "Grid point within aperture" << std::endl; + in_box = is_in_box(xi, yj, xbox, ybox); + if (in_box) + { + ++nhit_inside; + } + else + { + ++nhit_outside; + } + } + } + } + + return; +} + TEST(Aperture, ApertureBase) { struct TestAperture : public Aperture @@ -15,18 +88,28 @@ TEST(Aperture, ApertureBase) { } virtual ~TestAperture() {} - virtual double aperture_area() const { return my_value; } - virtual double diameter_circumscribed_circle() const { return 2.0; } - virtual bool is_in(double x, double y) const { return false; } - virtual aperture_ptr make_copy() const + virtual double aperture_area() const override { return my_value; } + virtual double diameter_circumscribed_circle() const override { return 2.0; } + virtual void bounding_box(double &xmin, + double &xmax, + double &ymin, + double &ymax) const override + { + return; + } + virtual bool is_in(double x, double y) const override { return false; } + virtual aperture_ptr make_copy() const override { return make_aperture(*this); } - virtual void write_json(nlohmann::ordered_json& jnode) const + virtual void write_json(nlohmann::ordered_json &jnode) const override { jnode["my_value"] = my_value; } - virtual std::tuple, std::vector> triangulation() const { return std::make_tuple(std::vector{}, std::vector{});} + virtual std::tuple, std::vector> triangulation() const override + { + return std::make_tuple(std::vector{}, std::vector{}); + } }; TestAperture ta1(1.2, ApertureType::CIRCLE); @@ -45,7 +128,7 @@ TEST(Aperture, ApertureBase) TEST(Aperture, Annulus) { - const double RO = 5.0; + constexpr double RO = 5.0; const double RI = 1.0; const double ARC1 = 90.0; const double ARC2 = 360.0; @@ -106,11 +189,33 @@ TEST(Aperture, Annulus) EXPECT_FALSE(a1->is_in(X3, Y3)); EXPECT_FALSE(a1->is_in(X4, Y4)); EXPECT_FALSE(a1->is_in(X5, Y5)); + + const uint_fast64_t NPOINTS = 1001; + constexpr double AX = -RO - 1.0; + constexpr double BX = RO + 1.0; + constexpr double AY = -RO - 1.0; + constexpr double BY = RO + 1.0; + constexpr double GRID_AREA = (BX - AX) * (BY - AY); + uint_fast64_t nhit_out = 0; + uint_fast64_t nhit_in = 0; + bounding_box_test(ann1, + AX, BX, NPOINTS, + AY, BY, NPOINTS, + nhit_out, nhit_in); + + double frac = ann1->aperture_area() / GRID_AREA; + + EXPECT_EQ(nhit_out, 0); + EXPECT_NEAR((double)(nhit_in) / (NPOINTS * NPOINTS), frac, 1e-3); + + std::cout << "Number of hits outside of bounding box: " << nhit_out << std::endl; + std::cout << "Number of hits inside bounding box: " << nhit_in << std::endl; + std::cout << "Expected fraction of hits: " << frac << std::endl; } TEST(Aperture, Circle) { - const double D = 2.0; + constexpr double D = 2.0; const double X1 = 0.5; const double Y1 = -0.5; const double X2 = 1.0; @@ -132,13 +237,35 @@ TEST(Aperture, Circle) EXPECT_EQ(ap->aperture_area(), cir->aperture_area()); EXPECT_TRUE(ap->is_in(X1, Y1)); EXPECT_FALSE(ap->is_in(X2, Y2)); + + const uint_fast64_t NPOINTS = 1001; + constexpr double AX = -0.5 * D - 1.0; + constexpr double BX = 0.5 * D + 1.0; + constexpr double AY = -0.5 * D - 1.0; + constexpr double BY = 0.5 * D + 1.0; + constexpr double GRID_AREA = (BX - AX) * (BY - AY); + uint_fast64_t nhit_out = 0; + uint_fast64_t nhit_in = 0; + bounding_box_test(cir, + AX, BX, NPOINTS, + AY, BY, NPOINTS, + nhit_out, nhit_in); + + double frac = cir->aperture_area() / GRID_AREA; + + EXPECT_EQ(nhit_out, 0); + EXPECT_NEAR((double)(nhit_in) / (NPOINTS * NPOINTS), frac, 1e-3); + + std::cout << "Number of hits outside of bounding box: " << nhit_out << std::endl; + std::cout << "Number of hits inside bounding box: " << nhit_in << std::endl; + std::cout << "Expected fraction of hits: " << frac << std::endl; } TEST(Aperture, EqualateralTriangle) { const double TOL = 1e-12; - const double D = 2.0; - const double R = 0.5 * D; + constexpr double D = 2.0; + constexpr double R = 0.5 * D; const double S = sqrt(3.0) * R; // Side length of triangle const double AREA = sqrt(27.0) * R * R / 4.0; @@ -178,6 +305,28 @@ TEST(Aperture, EqualateralTriangle) EXPECT_EQ(ap->aperture_area(), et->aperture_area()); EXPECT_TRUE(ap->is_in(X1, Y1)); EXPECT_FALSE(ap->is_in(X3, Y3)); + + const uint_fast64_t NPOINTS = 1001; + constexpr double AX = -R - 1.0; + constexpr double BX = R + 1.0; + constexpr double AY = -R - 1.0; + constexpr double BY = R + 1.0; + constexpr double GRID_AREA = (BX - AX) * (BY - AY); + uint_fast64_t nhit_out = 0; + uint_fast64_t nhit_in = 0; + bounding_box_test(et, + AX, BX, NPOINTS, + AY, BY, NPOINTS, + nhit_out, nhit_in); + + double frac = et->aperture_area() / GRID_AREA; + + EXPECT_EQ(nhit_out, 0); + EXPECT_NEAR((double)(nhit_in) / (NPOINTS * NPOINTS), frac, 1e-3); + + std::cout << "Number of hits outside of bounding box: " << nhit_out << std::endl; + std::cout << "Number of hits inside bounding box: " << nhit_in << std::endl; + std::cout << "Expected fraction of hits: " << frac << std::endl; } TEST(Aperture, Hexagon) @@ -231,6 +380,28 @@ TEST(Aperture, Hexagon) EXPECT_EQ(ap->aperture_area(), hex->aperture_area()); EXPECT_TRUE(ap->is_in(X2, Y2)); EXPECT_FALSE(ap->is_in(X4, Y4)); + + const uint_fast64_t NPOINTS = 1001; + const double AX = -R - 1.0; + const double BX = R + 1.0; + const double AY = -R - 1.0; + const double BY = R + 1.0; + const double GRID_AREA = (BX - AX) * (BY - AY); + uint_fast64_t nhit_out = 0; + uint_fast64_t nhit_in = 0; + bounding_box_test(hex, + AX, BX, NPOINTS, + AY, BY, NPOINTS, + nhit_out, nhit_in); + + double frac = hex->aperture_area() / GRID_AREA; + + EXPECT_EQ(nhit_out, 0); + EXPECT_NEAR((double)(nhit_in) / (NPOINTS * NPOINTS), frac, 1e-3); + + std::cout << "Number of hits outside of bounding box: " << nhit_out << std::endl; + std::cout << "Number of hits inside bounding box: " << nhit_in << std::endl; + std::cout << "Expected fraction of hits: " << frac << std::endl; } TEST(Aperture, Rectangle) @@ -320,6 +491,28 @@ TEST(Aperture, Rectangle) EXPECT_EQ(ap_shift->aperture_area(), rect_shift->aperture_area()); EXPECT_TRUE(ap_shift->is_in(X6, Y6)); EXPECT_FALSE(ap_shift->is_in(X8, Y8)); + + const uint_fast64_t NPOINTS = 1001; + const double AX = -LX - 1.0; + const double BX = LX + 1.0; + const double AY = -LY - 1.0; + const double BY = LY + 1.0; + const double GRID_AREA = (BX - AX) * (BY - AY); + uint_fast64_t nhit_out = 0; + uint_fast64_t nhit_in = 0; + bounding_box_test(rect, + AX, BX, NPOINTS, + AY, BY, NPOINTS, + nhit_out, nhit_in); + + double frac = rect->aperture_area() / GRID_AREA; + + EXPECT_EQ(nhit_out, 0); + EXPECT_NEAR((double)(nhit_in) / (NPOINTS * NPOINTS), frac, 1e-3); + + std::cout << "Number of hits outside of bounding box: " << nhit_out << std::endl; + std::cout << "Number of hits inside bounding box: " << nhit_in << std::endl; + std::cout << "Expected fraction of hits: " << frac << std::endl; } TEST(Aperture, IrregularTriangle) @@ -338,6 +531,28 @@ TEST(Aperture, IrregularTriangle) EXPECT_NEAR(ap->aperture_area(), 0.5 * y2 * (x3 - x1), TOL); EXPECT_TRUE(ap->is_in(1.0, 1.0)); EXPECT_FALSE(ap->is_in(1.5, 2.0)); + + const uint_fast64_t NPOINTS = 1001; + const double AX = x1 - 1.0; + const double BX = x3 + 1.0; + const double AY = y1 - 1.0; + const double BY = y2 + 1.0; + const double GRID_AREA = (BX - AX) * (BY - AY); + uint_fast64_t nhit_out = 0; + uint_fast64_t nhit_in = 0; + bounding_box_test(tri, + AX, BX, NPOINTS, + AY, BY, NPOINTS, + nhit_out, nhit_in); + + double frac = tri->aperture_area() / GRID_AREA; + + EXPECT_EQ(nhit_out, 0); + EXPECT_NEAR((double)(nhit_in) / (NPOINTS * NPOINTS), frac, 1e-3); + + std::cout << "Number of hits outside of bounding box: " << nhit_out << std::endl; + std::cout << "Number of hits inside bounding box: " << nhit_in << std::endl; + std::cout << "Expected fraction of hits: " << frac << std::endl; } TEST(Aperture, IrregularQuadrilateral) @@ -361,12 +576,34 @@ TEST(Aperture, IrregularQuadrilateral) EXPECT_TRUE(ap->is_in(3.0, 1.0)); EXPECT_TRUE(ap->is_in(1.0, 1.5)); EXPECT_FALSE(ap->is_in(4.0, 1.0)); + + const uint_fast64_t NPOINTS = 1001; + const double AX = x1 - 1.0; + const double BX = x2 + 1.0; + const double AY = y1 - 1.0; + const double BY = y3 + 1.0; + const double GRID_AREA = (BX - AX) * (BY - AY); + uint_fast64_t nhit_out = 0; + uint_fast64_t nhit_in = 0; + bounding_box_test(quad, + AX, BX, NPOINTS, + AY, BY, NPOINTS, + nhit_out, nhit_in); + + double frac = quad->aperture_area() / GRID_AREA; + + EXPECT_EQ(nhit_out, 0); + EXPECT_NEAR((double)(nhit_in) / (NPOINTS * NPOINTS), frac, 1e-3); + + std::cout << "Number of hits outside of bounding box: " << nhit_out << std::endl; + std::cout << "Number of hits inside bounding box: " << nhit_in << std::endl; + std::cout << "Expected fraction of hits: " << frac << std::endl; } TEST(Aperture, MakeApertureFromType) { const double TOL = 1e-12; - + // Test Circle creation std::vector circle_args = {2.0}; // diameter auto circle_ap = Aperture::make_aperture_from_type(ApertureType::CIRCLE, circle_args); @@ -374,7 +611,7 @@ TEST(Aperture, MakeApertureFromType) EXPECT_EQ(circle_ap->get_type(), ApertureType::CIRCLE); EXPECT_EQ(circle_ap->diameter_circumscribed_circle(), 2.0); EXPECT_NEAR(circle_ap->aperture_area(), PI, TOL); - + // Test Rectangle creation std::vector rect_args = {3.0, 2.0}; // width, height auto rect_ap = Aperture::make_aperture_from_type(ApertureType::RECTANGLE, rect_args); @@ -382,7 +619,7 @@ TEST(Aperture, MakeApertureFromType) EXPECT_EQ(rect_ap->get_type(), ApertureType::RECTANGLE); EXPECT_NEAR(rect_ap->aperture_area(), 6.0, TOL); EXPECT_TRUE(rect_ap->is_in(0.0, 0.0)); // Center should be inside - + // Test Annulus creation std::vector annulus_args = {1.0, 3.0, 180.0}; // inner_radius, outer_radius, arc_angle auto annulus_ap = Aperture::make_aperture_from_type(ApertureType::ANNULUS, annulus_args); @@ -390,39 +627,39 @@ TEST(Aperture, MakeApertureFromType) EXPECT_EQ(annulus_ap->get_type(), ApertureType::ANNULUS); EXPECT_EQ(annulus_ap->diameter_circumscribed_circle(), 6.0); EXPECT_NEAR(annulus_ap->aperture_area(), 0.5 * PI * (9.0 - 1.0), TOL); - + // Test Hexagon creation std::vector hex_args = {4.0}; // circumscribe_diameter auto hex_ap = Aperture::make_aperture_from_type(ApertureType::HEXAGON, hex_args); ASSERT_TRUE(hex_ap != nullptr); EXPECT_EQ(hex_ap->get_type(), ApertureType::HEXAGON); EXPECT_EQ(hex_ap->diameter_circumscribed_circle(), 4.0); - + // Test Equilateral Triangle creation std::vector tri_args = {2.0}; // circumscribe_diameter auto tri_ap = Aperture::make_aperture_from_type(ApertureType::EQUILATERAL_TRIANGLE, tri_args); ASSERT_TRUE(tri_ap != nullptr); EXPECT_EQ(tri_ap->get_type(), ApertureType::EQUILATERAL_TRIANGLE); EXPECT_EQ(tri_ap->diameter_circumscribed_circle(), 2.0); - + // Test Irregular Triangle creation std::vector irregular_tri_args = {0.0, 0.0, 1.0, 2.0, 2.0, 0.0}; // x1,y1, x2,y2, x3,y3 auto irregular_tri_ap = Aperture::make_aperture_from_type(ApertureType::IRREGULAR_TRIANGLE, irregular_tri_args); ASSERT_TRUE(irregular_tri_ap != nullptr); EXPECT_EQ(irregular_tri_ap->get_type(), ApertureType::IRREGULAR_TRIANGLE); EXPECT_NEAR(irregular_tri_ap->aperture_area(), 2.0, TOL); - + // Test Irregular Quadrilateral creation std::vector quad_args = {0.0, 0.0, 3.0, 0.0, 4.0, 2.0, 1.0, 2.0}; // x1,y1, x2,y2, x3,y3, x4,y4 auto quad_ap = Aperture::make_aperture_from_type(ApertureType::IRREGULAR_QUADRILATERAL, quad_args); ASSERT_TRUE(quad_ap != nullptr); EXPECT_EQ(quad_ap->get_type(), ApertureType::IRREGULAR_QUADRILATERAL); - + // Test insufficient arguments - should return null pointer std::vector insufficient_args = {1.0}; // Not enough args for annulus auto null_ap = Aperture::make_aperture_from_type(ApertureType::ANNULUS, insufficient_args); EXPECT_TRUE(null_ap == nullptr); - + // Test empty arguments std::vector empty_args; auto null_ap2 = Aperture::make_aperture_from_type(ApertureType::CIRCLE, empty_args); diff --git a/google-tests/unit-tests/simulation_data/cst-templates/CMakeLists.txt b/google-tests/unit-tests/simulation_data/cst-templates/CMakeLists.txt index 69759a20..9c05a1a1 100644 --- a/google-tests/unit-tests/simulation_data/cst-templates/CMakeLists.txt +++ b/google-tests/unit-tests/simulation_data/cst-templates/CMakeLists.txt @@ -14,7 +14,7 @@ include_directories( # CST Template test files set(CST_TEST_SRC - utilities_test.cpp + # utilities_test.cpp linear_fresnel_test.cpp parabolic_trough_test.cpp parabolic_dish_test.cpp diff --git a/google-tests/unit-tests/simulation_data/cst-templates/heliostat_field_test.cpp b/google-tests/unit-tests/simulation_data/cst-templates/heliostat_field_test.cpp index ff5a284d..16dbe789 100644 --- a/google-tests/unit-tests/simulation_data/cst-templates/heliostat_field_test.cpp +++ b/google-tests/unit-tests/simulation_data/cst-templates/heliostat_field_test.cpp @@ -13,11 +13,11 @@ #include #include #include +#include #include <../../hpvm.h> #include -#include #include "common.hpp" #include "count_absorbed_native.h" @@ -995,4 +995,4 @@ TEST_F(HeliostatFieldSimulation, multiFacet_BandFocused_BandCanted) set_scatter_aimpoints(); simulate_check_outputs("3b", "2", "8"); simulate_check_outputs("3b", "2", "12"); -} \ No newline at end of file +} diff --git a/google-tests/unit-tests/simulation_data/cst-templates/heliostat_test.cpp b/google-tests/unit-tests/simulation_data/cst-templates/heliostat_test.cpp index 2bc4d8c0..bc57e910 100644 --- a/google-tests/unit-tests/simulation_data/cst-templates/heliostat_test.cpp +++ b/google-tests/unit-tests/simulation_data/cst-templates/heliostat_test.cpp @@ -6,9 +6,9 @@ #include #include #include +#include #include -#include #include "common.hpp" #include "count_absorbed_native.h" @@ -660,4 +660,4 @@ TEST(Heliostat, UpdateGeometry) EXPECT_TRUE(n >= NRAYS); EXPECT_TRUE(num_absorbed > N_ABSORBED_THRESH); -} \ No newline at end of file +} diff --git a/google-tests/unit-tests/simulation_data/cst-templates/linear_fresnel_test.cpp b/google-tests/unit-tests/simulation_data/cst-templates/linear_fresnel_test.cpp index 81141cc3..7301080e 100644 --- a/google-tests/unit-tests/simulation_data/cst-templates/linear_fresnel_test.cpp +++ b/google-tests/unit-tests/simulation_data/cst-templates/linear_fresnel_test.cpp @@ -5,10 +5,10 @@ #include #include #include +#include #include #include -#include #include "common.hpp" #include "count_absorbed_native.h" diff --git a/google-tests/unit-tests/simulation_data/cst-templates/parabolic_dish_test.cpp b/google-tests/unit-tests/simulation_data/cst-templates/parabolic_dish_test.cpp index 52a917c3..1036e078 100644 --- a/google-tests/unit-tests/simulation_data/cst-templates/parabolic_dish_test.cpp +++ b/google-tests/unit-tests/simulation_data/cst-templates/parabolic_dish_test.cpp @@ -4,10 +4,10 @@ #include #include #include +#include #include #include -#include #include "common.hpp" #include "count_absorbed_native.h" diff --git a/google-tests/unit-tests/simulation_data/cst-templates/parabolic_trough_test.cpp b/google-tests/unit-tests/simulation_data/cst-templates/parabolic_trough_test.cpp index 56dacec0..6a79aebc 100644 --- a/google-tests/unit-tests/simulation_data/cst-templates/parabolic_trough_test.cpp +++ b/google-tests/unit-tests/simulation_data/cst-templates/parabolic_trough_test.cpp @@ -4,10 +4,10 @@ #include #include #include +#include #include #include -#include #include "common.hpp" #include "count_absorbed_native.h" diff --git a/google-tests/unit-tests/simulation_data/cst-templates/single_heliostat_test.cpp b/google-tests/unit-tests/simulation_data/cst-templates/single_heliostat_test.cpp index ea1abaaf..84fe6d89 100644 --- a/google-tests/unit-tests/simulation_data/cst-templates/single_heliostat_test.cpp +++ b/google-tests/unit-tests/simulation_data/cst-templates/single_heliostat_test.cpp @@ -8,11 +8,11 @@ #include #include #include +#include #include <../../hpvm.h> #include -#include #include "common.hpp" #include "count_absorbed_native.h" @@ -792,4 +792,4 @@ TEST_F(SingleHeliostatSimulation, MultiFacetFocused_SlantCanting_Southeast) EXPECT_NEAR(sun_height, 11.5183, 1.e-3); } -// TODO: add off-axis cases \ No newline at end of file +// TODO: add off-axis cases diff --git a/google-tests/unit-tests/simulation_data/element_test.cpp b/google-tests/unit-tests/simulation_data/element_test.cpp index 96b78a50..0b94cf3a 100644 --- a/google-tests/unit-tests/simulation_data/element_test.cpp +++ b/google-tests/unit-tests/simulation_data/element_test.cpp @@ -255,6 +255,14 @@ TEST(Element, CompositeElementAccessors) EXPECT_FALSE(cmp->is_virtual()); EXPECT_FALSE(elem2->is_virtual()); + cmp->mark_virtual(); + auto elem3 = make_configured_element(); + EXPECT_FALSE(elem3->is_virtual()); + cmp->add_element(elem3); + EXPECT_TRUE(elem3->is_virtual()); + cmp->unmark_virtual(); + EXPECT_FALSE(elem3->is_virtual()); + // Check that pass through functions are hooked up correctly auto iter = cmp->get_iterator(); while (!cmp->is_at_end(iter)) diff --git a/google-tests/unit-tests/simulation_data/surface_test.cpp b/google-tests/unit-tests/simulation_data/surface_test.cpp index 3d7feca4..ac99beaf 100644 --- a/google-tests/unit-tests/simulation_data/surface_test.cpp +++ b/google-tests/unit-tests/simulation_data/surface_test.cpp @@ -1,5 +1,6 @@ #include +#include #include TEST(Surface, Typing) @@ -7,8 +8,6 @@ TEST(Surface, Typing) // Test that each constructor properly sets the type auto cone = SolTrace::Data::make_surface(50.0); EXPECT_EQ(cone->get_type(), SolTrace::Data::CONE); - // cone = SolTrace::Data::make_surface(); - // EXPECT_EQ(cone->get_type(), CONE); auto cylinder = SolTrace::Data::make_surface(1.0); EXPECT_EQ(cylinder->get_type(), SolTrace::Data::CYLINDER); @@ -18,13 +17,9 @@ TEST(Surface, Typing) auto para = SolTrace::Data::make_surface(1.0, 1.0); EXPECT_EQ(para->get_type(), SolTrace::Data::PARABOLA); - // para = SolTrace::Data::make_surface(); - // EXPECT_EQ(para->get_type(), PARABOLA); auto sph = SolTrace::Data::make_surface(10.0); EXPECT_EQ(sph->get_type(), SolTrace::Data::SPHERE); - // sph = SolTrace::Data::make_surface(); - // EXPECT_EQ(sph->get_type(), SPHERE); } TEST(Surface, MakeSurfaceFromType) @@ -170,3 +165,182 @@ TEST(Surface, MakeSurfaceFromType) EXPECT_DOUBLE_EQ(sphere_cast->vertex_curv, 5.0); } } + +TEST(Cone, MakeCopy) +{ + const double xbox[2] = {-1.0, 2.0}; + const double ybox[2] = {0.0, 1.0}; + double z0, z1; + double z0copy, z1copy; + + auto cone = make_surface(35.0); + auto copy = cone->make_copy(); + + cone->bounding_box(xbox, ybox, z0, z1); + copy->bounding_box(xbox, ybox, z0copy, z1copy); + EXPECT_EQ(z0, z0copy); + EXPECT_EQ(z1, z1copy); + EXPECT_EQ(cone->z(xbox[0], ybox[1]), copy->z(xbox[0], ybox[1])); +} + +TEST(Cone, BoundingBox) +{ + const double TOL = 1e-12; + const double THETA = 30.0; + const double xbox[2] = {-1.0, 2.0}; + const double ybox[2] = {0.0, 1.0}; + double z0, z1; + + auto cone = make_surface(THETA); + cone->bounding_box(xbox, ybox, z0, z1); + EXPECT_NEAR(z0, 0.0, TOL); + EXPECT_NEAR(z1, sqrt(4.0 + 1.0) / tan(THETA * D2R), TOL); +} + +TEST(Cylinder, MakeCopy) +{ + const double xbox[2] = {-1.0, 1.0}; + const double ybox[2] = {0.0, 5.0}; + double z0, z1; + double z0copy, z1copy; + + auto cylin = make_surface(1.0); + auto copy = cylin->make_copy(); + + cylin->bounding_box(xbox, ybox, z0, z1); + copy->bounding_box(xbox, ybox, z0copy, z1copy); + EXPECT_EQ(z0, z0copy); + EXPECT_EQ(z1, z1copy); + EXPECT_EQ(cylin->z(xbox[0], ybox[1]), copy->z(xbox[0], ybox[1])); +} + +TEST(Cylinder, BoundingBox) +{ + const double TOL = 1e-12; + const double R = 1.5; + const double xbox[2] = {-R, R}; + const double ybox[2] = {0.0, 5.0}; + double z0, z1; + + auto cylin = make_surface(R); + auto copy = cylin->make_copy(); + + cylin->bounding_box(xbox, ybox, z0, z1); + EXPECT_NEAR(z0, 0.0, TOL); + EXPECT_NEAR(z1, 2.0 * R, TOL); +} + +TEST(Flat, MakeCopy) +{ + const double xbox[2] = {-1.0, 2.0}; + const double ybox[2] = {0.0, 1.0}; + double z0, z1; + double z0copy, z1copy; + + auto flat = make_surface(); + auto copy = flat->make_copy(); + + flat->bounding_box(xbox, ybox, z0, z1); + copy->bounding_box(xbox, ybox, z0copy, z1copy); + EXPECT_EQ(z0, z0copy); + EXPECT_EQ(z1, z1copy); + EXPECT_EQ(flat->z(xbox[0], ybox[1]), copy->z(xbox[0], ybox[1])); +} + +TEST(Flat, BoundingBox) +{ + const double TOL = 1e-12; + const double xbox[2] = {-PI, 2.0 * PI}; + const double ybox[2] = {0.0, 5.0}; + double z0, z1; + + auto flat = make_surface(); + auto copy = flat->make_copy(); + + flat->bounding_box(xbox, ybox, z0, z1); + EXPECT_NEAR(z0, -1e-4, TOL); + EXPECT_NEAR(z1, 1e-4, TOL); +} + +TEST(Parabola, MakeCopy) +{ + const double xbox[2] = {-1.0, 2.0}; + const double ybox[2] = {0.0, 1.0}; + double z0, z1; + double z0copy, z1copy; + + auto para = make_surface(1.0, 2.0); + auto copy = para->make_copy(); + + para->bounding_box(xbox, ybox, z0, z1); + copy->bounding_box(xbox, ybox, z0copy, z1copy); + EXPECT_EQ(z0, z0copy); + EXPECT_EQ(z1, z1copy); + EXPECT_EQ(para->z(xbox[0], ybox[1]), copy->z(xbox[0], ybox[1])); +} + +TEST(Parabola, BoundingBox) +{ + const double TOL = 1e-12; + const double xbox[2] = {-1.0, 2.0}; + const double ybox[2] = {0.0, 1.0}; + const double FX = 1.0; + const double FY = 2.0; + const double CX = 0.5 / FX; + const double CY = 0.5 / FY; + double z0, z1; + + auto para = make_surface(FX, FY); + + para->bounding_box(xbox, ybox, z0, z1); + EXPECT_NEAR(z0, 0.0, TOL); + EXPECT_NEAR(z1, 1.125, TOL); + + const double xbox2[2] = {1.0, 4.0}; + const double ybox2[2] = {0.5, 1.0}; + + para->bounding_box(xbox2, ybox2, z0, z1); + EXPECT_NEAR(z0, 0.5 * (CX * 1.0 + CY * 0.25), TOL); + EXPECT_NEAR(z1, 0.5 * (CX * 16.0 + CY * 1.0), TOL); +} + +TEST(Sphere, MakeCopy) +{ + const double xbox[2] = {-1.0, 2.0}; + const double ybox[2] = {0.0, 1.0}; + double z0, z1; + double z0copy, z1copy; + + auto sph = make_surface(0.5); + auto copy = sph->make_copy(); + + sph->bounding_box(xbox, ybox, z0, z1); + copy->bounding_box(xbox, ybox, z0copy, z1copy); + EXPECT_EQ(z0, z0copy); + EXPECT_EQ(z1, z1copy); + EXPECT_EQ(sph->z(xbox[0], ybox[1]), copy->z(xbox[0], ybox[1])); +} + +TEST(Sphere, BoundingBox) +{ + const double TOL = 1e-12; + const double CURV = 0.25; + const double R = 1.0 / CURV; + const double xbox[2] = {-1.0, 4.0}; + const double ybox[2] = {0.0, 3.0}; + double z0, z1; + + auto sph = make_surface(CURV); + + sph->bounding_box(xbox, ybox, z0, z1); + EXPECT_NEAR(z0, 0.0, TOL); + EXPECT_NEAR(z1, R, TOL); + + const double xbox2[2] = {0.0, 1.0}; + const double ybox2[2] = {0.0, 1.0}; + const double rsq = 2.0; + + sph->bounding_box(xbox2, ybox2, z0, z1); + EXPECT_NEAR(z0, 0.0, TOL); + EXPECT_NEAR(z1, R - sqrt(R*R - rsq), TOL); +} diff --git a/google-tests/unit-tests/simulation_data/cst-templates/utilities_test.cpp b/google-tests/unit-tests/simulation_data/utilities_test.cpp similarity index 64% rename from google-tests/unit-tests/simulation_data/cst-templates/utilities_test.cpp rename to google-tests/unit-tests/simulation_data/utilities_test.cpp index 87e4a602..3757e850 100644 --- a/google-tests/unit-tests/simulation_data/cst-templates/utilities_test.cpp +++ b/google-tests/unit-tests/simulation_data/utilities_test.cpp @@ -2,11 +2,30 @@ #include +#include #include -#include #include "common.hpp" +TEST(Utilities, IsApprox) +{ + const double TOL1 = 1e-12; + const double TOL2 = 1e-8; + const double x = 1.0; + const double y = x + 0.25 * TOL2; + + EXPECT_TRUE(is_approx(x, y, TOL2)); + EXPECT_FALSE(is_approx(x, y, TOL1)); +} + +TEST(Utilities, AbsMaxMin) +{ + const unsigned N = 5; + const double u[N] = {1.0, -2.0, -3.0, 4.0, -5.0}; + EXPECT_EQ(abs_max(u, N), 5.0); + EXPECT_EQ(abs_min(u, N), 1.0); +} + TEST(Utilities, Projection) { Vector3d v(1.0, 1.0, 1.0); diff --git a/google-tests/unit-tests/simulation_runner/CMakeLists.txt b/google-tests/unit-tests/simulation_runner/CMakeLists.txt index d599f4c9..8e0368d5 100644 --- a/google-tests/unit-tests/simulation_runner/CMakeLists.txt +++ b/google-tests/unit-tests/simulation_runner/CMakeLists.txt @@ -3,6 +3,11 @@ Project(SimulationRunnerUnitTests) # Add native runner tests (always built) add_subdirectory(native_runner) +# Add embree runner tests only if Embree is enabled +if(SOLTRACE_BUILD_EMBREE_SUPPORT) + add_subdirectory(embree_runner) +endif() + # Add optix runner tests only if OptiX is enabled if(SOLTRACE_BUILD_OPTIX_SUPPORT) add_subdirectory(optix_runner) diff --git a/google-tests/unit-tests/simulation_runner/embree_runner/CMakeLists.txt b/google-tests/unit-tests/simulation_runner/embree_runner/CMakeLists.txt new file mode 100644 index 00000000..7177de40 --- /dev/null +++ b/google-tests/unit-tests/simulation_runner/embree_runner/CMakeLists.txt @@ -0,0 +1,51 @@ +Project(EmbreeRunnerUnitTests) + +include_directories( + . + ../../common + ../../../test_tools + ../../../app/deploy/samples + ../../../../coretrace/simulation_data + ../../../../coretrace/simulation_results + ../../../../coretrace/simulation_runner + ../../../../coretrace/simulation_runner/native_runner + ../../../../coretrace/simulation_runner/embree_runner +) + +# Add test files for Embree_runner +set(EMBREE_RUNNER_TEST_SRC + ../../common/common.cpp + embree_runner_test.cpp +) + +if (NOT ENABLE_COVERAGE) + list(APPEND + EMBREE_RUNNER_TEST_SRC + embree_runner_multithreading_test.cpp + ) +endif() + +add_executable(EmbreeRunnerUnitTests + ${EMBREE_RUNNER_TEST_SRC} +) + +target_link_libraries( + EmbreeRunnerUnitTests + simdata + simresult + # native_runner + embree_runner + testtools + GTest::gtest_main +) + +if (ENABLE_COVERAGE AND CMAKE_BUILD_TYPE STREQUAL "Debug") + target_link_options(EmbreeRunnerUnitTests PRIVATE --coverage) + target_link_libraries(EmbreeRunnerUnitTests gcov) +endif() + +if(UNIX) + target_link_libraries(EmbreeRunnerUnitTests pthread) +endif() + +add_test(NAME EmbreeRunnerUnitTests COMMAND EmbreeRunnerUnitTests) diff --git a/google-tests/unit-tests/simulation_runner/embree_runner/embree_runner_multithreading_test.cpp b/google-tests/unit-tests/simulation_runner/embree_runner/embree_runner_multithreading_test.cpp new file mode 100644 index 00000000..7d9eedb6 --- /dev/null +++ b/google-tests/unit-tests/simulation_runner/embree_runner/embree_runner_multithreading_test.cpp @@ -0,0 +1,147 @@ +#include + +#include +#include +#include + +#include +#include +#include +#include + +using SolTrace::Runner::RunnerStatus; +using SolTrace::EmbreeRunner::EmbreeRunner; + +TEST(EmbreeRunner, StatusAndCancelMultiThread) +{ + std::string sample_path = std::string(PROJECT_DIR) + + std::string("/Power-tower-surround_singlefacet.stinput"); + + SimulationData sd; + EXPECT_TRUE(sd.import_from_file(sample_path)); + sd.set_number_of_rays(1000000); + + EmbreeRunner runner; + runner.set_number_of_threads(2); + RunnerStatus sts; + sts = runner.setup_simulation(&sd); + EXPECT_EQ(sts, RunnerStatus::SUCCESS); + + auto t0 = std::chrono::high_resolution_clock::now(); + + auto fsts = std::async(&EmbreeRunner::run_simulation, &runner); + + std::this_thread::sleep_for(std::chrono::milliseconds(200)); + sts = runner.status_simulation(); + EXPECT_EQ(sts, RunnerStatus::RUNNING); + + double prog = -1.0; + std::this_thread::sleep_for(std::chrono::milliseconds(200)); + sts = runner.status_simulation(&prog); + EXPECT_EQ(sts, RunnerStatus::RUNNING); + EXPECT_LT(prog, 1.0); + EXPECT_GT(prog, 0.0); + + runner.cancel_simulation(); + fsts.wait(); + + auto t1 = std::chrono::high_resolution_clock::now(); + std::chrono::duration dur = t1 - t0; + + EXPECT_EQ(fsts.get(), RunnerStatus::CANCEL); + EXPECT_LT(dur.count(), 10000.0); + + std::cout << "Time for run: " << dur.count() << std::endl; + std::cout << "Progress before cancel: " << prog << std::endl; +} + +TEST(EmbreeRunner, CancelMultithread) +{ + std::string sample_path = std::string(PROJECT_DIR) + + std::string("/Power-tower-surround_singlefacet.stinput"); + + SimulationData sd; + EXPECT_TRUE(sd.import_from_file(sample_path)); + sd.set_number_of_rays(1000000); + + EmbreeRunner runner; + runner.set_number_of_threads(4); + RunnerStatus sts; + sts = runner.setup_simulation(&sd); + ASSERT_EQ(sts, RunnerStatus::SUCCESS); + + auto fsts = std::async(&EmbreeRunner::run_simulation, &runner); + + // Give time to start processing + std::this_thread::sleep_for(std::chrono::milliseconds(500)); + sts = runner.status_simulation(); + EXPECT_EQ(sts, RunnerStatus::RUNNING); + + // Shut everything down to make sure it doesn't hang + auto t0 = std::chrono::high_resolution_clock::now(); + runner.cancel_simulation(); + ASSERT_EQ(fsts.wait_for(std::chrono::seconds(10)), std::future_status::ready); + auto t1 = std::chrono::high_resolution_clock::now(); + std::chrono::duration dur = t1 - t0; + + EXPECT_EQ(fsts.get(), RunnerStatus::CANCEL); + + std::cout << "Time to cancel: " << dur.count() << std::endl; +} + +TEST(EmbreeRunner, RayIdAssignmentMultiThread) +{ + const uint_fast64_t NRAYS = 50000; + const unsigned NTHREADS = 12; + + std::string project_path = std::string(PROJECT_DIR); + std::string sample_path = project_path + + std::string("/simple_test_case.stinput"); + + // Load simulation data from file + SimulationData sd; + ASSERT_TRUE(sd.import_from_file(sample_path)); + + sd.set_number_of_rays(NRAYS); + + // Create and run the native runner + EmbreeRunner runner; + runner.set_number_of_threads(NTHREADS); + RunnerStatus sts = runner.initialize(); + ASSERT_EQ(sts, RunnerStatus::SUCCESS); + sts = runner.setup_simulation(&sd); + ASSERT_EQ(sts, RunnerStatus::SUCCESS); + sts = runner.run_simulation(); + ASSERT_EQ(sts, RunnerStatus::SUCCESS); + + double prog = -1.0; + sts = runner.status_simulation(&prog); + EXPECT_LE(prog, 1.0); + EXPECT_GT(prog, 0.0); + std::cout << "Progress: " << prog << std::endl; + + SimulationResult result; + sts = runner.report_simulation(&result, 0); + ASSERT_EQ(sts, RunnerStatus::SUCCESS); + + // std::cout << result << std::endl; + + uint_fast64_t nrec = result.get_number_of_records(); + + std::cout << "Number of rays: " << NRAYS + << "\nNumber of records: " << nrec + << std::endl; + + if (nrec != NRAYS) + { + runner.print_log(std::cout); + } + + ASSERT_EQ(nrec, NRAYS); + + for (uint_fast64_t k = 0; k < nrec; ++k) + { + auto ray_rec = result[k]; + ASSERT_EQ(ray_rec->id, k + 1); + } +} diff --git a/google-tests/unit-tests/simulation_runner/embree_runner/embree_runner_test.cpp b/google-tests/unit-tests/simulation_runner/embree_runner/embree_runner_test.cpp new file mode 100644 index 00000000..8b1aaae2 --- /dev/null +++ b/google-tests/unit-tests/simulation_runner/embree_runner/embree_runner_test.cpp @@ -0,0 +1,328 @@ +#include + +#include +#include +#include + +#include "common.hpp" +#include "count_absorbed_native.h" + +using SolTrace::EmbreeRunner::EmbreeRunner; +using SolTrace::EmbreeRunner::RunnerStatus; +using SolTrace::EmbreeRunner::TSystem; +using SolTrace::EmbreeRunner::TRayData; + +using SolTrace::Result::ray_record_ptr; +using SolTrace::Result::RayEvent; +using SolTrace::Result::SimulationResult; + +TEST(EmbreeRunner, SingleRayValidationTest) +{ + const double TOL = 1e-8; + + constexpr double c = 0.09; + constexpr double R = 1.0 / c; + constexpr double x = -3.0621423346154577; + constexpr double y = 5.9286205128611948; + constexpr double z0 = 15.0; + + double zref = R - sqrt(R * R - (x * x + y * y)); + double z = z0 - zref; + + Vector3d u(0.0, 0.0, -1.0); + Vector3d v(2.0 * x, 2.0 * y, -2.0 * (z0 - z - R)); + v.make_unit(); + Vector3d w; + double alpha = dot_product(u, v); + vector_add(1.0, u, -2.0 * alpha, v, w); + + SimulationData sd; + + // Set parameters + const uint_fast64_t NRAYS = 1; + SimulationParameters ¶ms = sd.get_simulation_parameters(); + params.number_of_rays = NRAYS; + params.max_number_of_rays = params.number_of_rays * 100; + params.include_optical_errors = false; + params.include_sun_shape_errors = false; + params.seed = 1; + + // Sun + auto sun = SolTrace::Data::make_ray_source(); + sun->set_position(0.0, 0.0, 100.0); + sun->set_shape(SolTrace::Data::SunShape::PILLBOX, -1.0, 1.0, 0.0); + sd.add_ray_source(sun); + + auto sph = SolTrace::Data::make_element(); + Vector3d origin(0.0, 0.0, z0); + Vector3d aim(0.0, 0.0, -1.0); + double zrot = 0.0; + sph->set_reference_frame_geometry(origin, aim, zrot); + sph->set_aperture(SolTrace::Data::make_aperture(20.0)); + sph->set_surface(SolTrace::Data::make_surface(c)); + sph->get_front_optical_properties()->set_ideal_reflection(); + sph->get_back_optical_properties()->set_ideal_reflection(); + sph->set_name("Sphere"); + sd.add_element(sph); + + auto para = SolTrace::Data::make_element(); + origin.set_values(0.0, 0.0, -1.0); + aim.set_values(0.0, 0.0, 0.0); + zrot = 0.0; + para->set_reference_frame_geometry(origin, aim, zrot); + para->set_aperture(SolTrace::Data::make_aperture(31.0, 31.0)); + para->set_surface(SolTrace::Data::make_surface(0.5 / 0.03, 0.5 / 0.03)); + para->set_name("Parabola"); + sd.add_element(para); + + // std::cout << "Constructing..." << std::endl; + EmbreeRunner runner; + runner.set_number_of_threads(1); + // std::cout << "Initializing..." << std::endl; + RunnerStatus sts = runner.initialize(); + EXPECT_EQ(sts, RunnerStatus::SUCCESS); + // std::cout << "Setting up simulation..." << std::endl; + sts = runner.setup_simulation(&sd); + EXPECT_EQ(sts, RunnerStatus::SUCCESS); + // std::cout << "Running simulation..." << std::endl; + sts = runner.run_simulation(); + EXPECT_EQ(sts, RunnerStatus::SUCCESS); + + const TSystem *sys = runner.get_system(); + sys->RayData.Print(); + const TRayData *ray_data = &(sys->RayData); + size_t n = ray_data->Count(); + + // sys->RayData.Print(); + + EXPECT_EQ(n, 3); + + Vector3d ipoint, idir; + int element, stage; + uint_fast64_t raynum; + SolTrace::Result::RayEvent rev; + sys->RayData.Query(0, ipoint.data, idir.data, + &element, &stage, &raynum, &rev); + + EXPECT_EQ(raynum, 1); + EXPECT_EQ(rev, SolTrace::Result::RayEvent::CREATE); + + EXPECT_NEAR(ipoint[0], x, TOL); + EXPECT_NEAR(ipoint[1], y, TOL); + EXPECT_NEAR(ipoint[2], 10000.0, TOL); + + EXPECT_NEAR(idir[0], u[0], TOL); + EXPECT_NEAR(idir[1], u[1], TOL); + EXPECT_NEAR(idir[2], u[2], TOL); + + sys->RayData.Query(1, ipoint.data, idir.data, + &element, &stage, &raynum, &rev); + + EXPECT_EQ(raynum, 1); + EXPECT_EQ(rev, SolTrace::Result::RayEvent::REFLECT); + + EXPECT_NEAR(ipoint[0], x, TOL); + EXPECT_NEAR(ipoint[1], y, TOL); + EXPECT_NEAR(ipoint[2], z, TOL); + + EXPECT_NEAR(idir[0], w[0], TOL); + EXPECT_NEAR(idir[1], w[1], TOL); + EXPECT_NEAR(idir[2], w[2], TOL); +} + +TEST(EmbreeRunner, PowerTowerSmokeTest) +{ + SimulationData sd; + + // Sun + auto sun = SolTrace::Data::make_ray_source(); + sun->set_position(0.0, 0.0, 100.0); + sd.add_ray_source(sun); + + // Absorber -- Flat + auto absorber = SolTrace::Data::make_element(); + absorber->set_origin(0.0, 0.0, 10.0); + absorber->set_aim_vector(0.0, 0.0, 5.0); + absorber->set_zrot(0.0); + absorber->compute_coordinate_rotations(); + absorber->set_surface(SolTrace::Data::make_surface()); // surface(nullptr) + absorber->set_aperture(SolTrace::Data::make_aperture(2.0, 2.0)); + OpticalProperties *foptics = absorber->get_front_optical_properties(); + foptics->my_type = InteractionType::REFLECTION; + foptics->reflectivity = 0.0; + + // Make stage 1 -- second stage -- these can be added to SimulationData + // in any order but should be numbered in the desired order + auto st1 = SolTrace::Data::make_stage(1); + // Origin is initialized to zero but set it explicitly + st1->set_origin(0.0, 0.0, 0.0); + // Set aim vector so stage and global coordinates are identical + st1->set_aim_vector(0.0, 0.0, 1.0); + // Set z rotation so stage and global coordinates are identical + st1->set_zrot(0.0); + // Compute coordinate rotations + st1->compute_coordinate_rotations(); + st1->add_element(absorber); + // Optional -- to help the user identify things + + // Make stage 0 -- this will be the first stage if the runner uses stages + auto st0 = SolTrace::Data::make_stage(0); + st0->set_origin(0.0, 0.0, 0.0); + st0->set_aim_vector(0.0, 0.0, 1.0); + + Vector3d rvec, svec, avec; + Vector3d aim, pos; + + const int NUM_ELEMENTS = 10; + for (int k = 0; k < NUM_ELEMENTS; ++k) + { + auto el = SolTrace::Data::make_element(); + foptics = el->get_front_optical_properties(); + foptics->reflectivity = 1.0; + + pos.set_values(5 * sin(k * PI * 2.0 / NUM_ELEMENTS), + 5 * cos(k * PI * 2.0 / NUM_ELEMENTS), + 0.0); + vector_add(1.0, absorber->get_origin_global(), + -1.0, pos, + rvec); + make_unit_vector(rvec); + svec = sun->get_position(); + make_unit_vector(svec); + vector_add(0.5, rvec, 0.5, svec, avec); + vector_add(1.0, pos, 100.0, avec, aim); + + el->set_reference_frame_geometry(pos, aim, 0.0); + + el->set_surface(SolTrace::Data::make_surface()); + el->set_aperture(SolTrace::Data::make_aperture(2.0)); + + st0->add_element(el); + } + EXPECT_EQ(st0->get_number_of_elements(), NUM_ELEMENTS); + sd.add_stage(st0); + EXPECT_EQ(sd.get_number_of_elements(), NUM_ELEMENTS); + sd.add_stage(st1); + EXPECT_EQ(sd.get_number_of_elements(), NUM_ELEMENTS + 1); + + // Set parameters + const uint_fast64_t NRAYS = 10000; + SimulationParameters ¶ms = sd.get_simulation_parameters(); + params.number_of_rays = NRAYS; + params.max_number_of_rays = params.number_of_rays * 100; + params.include_optical_errors = false; + params.include_sun_shape_errors = false; + params.seed = 12345; + + EmbreeRunner runner; + runner.set_number_of_threads(1); + RunnerStatus sts = runner.initialize(); + EXPECT_EQ(sts, RunnerStatus::SUCCESS); + sts = runner.setup_simulation(&sd); + EXPECT_EQ(sts, RunnerStatus::SUCCESS); + sts = runner.run_simulation(); + EXPECT_EQ(sts, RunnerStatus::SUCCESS); + + const TSystem *sys = runner.get_system(); + // sys->RayData.Print(); + const TRayData *ray_data = &(sys->RayData); + uint_fast64_t num_absorbed = count_absorbed_native(ray_data); + uint_fast64_t ncreate = count_event_native(ray_data, RayEvent::CREATE); + uint_fast64_t nexit = count_event_native(ray_data, RayEvent::EXIT); + size_t n = ray_data->Count(); + + std::cout << "Number Created: " << ncreate << std::endl; + std::cout << "Number Absorbed: " << num_absorbed << std::endl; + std::cout << "Number Interactions: " << n << std::endl; + std::cout << "Number Exit: " << nexit << std::endl; + + scan_events_native(ray_data); + + EXPECT_EQ(ncreate, NRAYS); + EXPECT_EQ(num_absorbed + nexit, ncreate); + EXPECT_EQ(num_absorbed + nexit, NRAYS); + + // ray_data->Print(); + + EXPECT_TRUE(n >= NRAYS); + EXPECT_TRUE(num_absorbed > 0); + + SimulationResult result; + sts = runner.report_simulation(&result, 0); + EXPECT_EQ(sts, RunnerStatus::SUCCESS); + + EXPECT_EQ(result.get_number_of_records(), NRAYS); + for (int_fast64_t idx = 0; idx < result.get_number_of_records(); ++idx) + { + const ray_record_ptr rr = result[idx]; + EXPECT_TRUE(rr->get_number_of_interactions() > 0); + } + + std::cout << "Number of ray records: " + << result.get_number_of_records() + << std::endl; +} + +TEST(EmbreeRunner, PowerTowerTest) +{ + // Pulling in path variable from CMake and creating path to .stinput sample file + std::string sample_path = std::string(PROJECT_DIR) + std::string("/Power-tower-surround_singlefacet.stinput"); + + // Create Simuluation Data + SimulationData sd; + + // Constants + const uint_fast64_t NRAYS = 10000; + const double TOL = 1e-4; + + // Read Input File + bool success = sd.import_from_file(sample_path); + EXPECT_TRUE(success); + EXPECT_TRUE(sd.get_number_of_elements() > 0); + EXPECT_TRUE(sd.get_number_of_ray_sources() > 0); + + std::cout << "Num Elements: " << sd.get_number_of_elements() << std::endl; + + // Parameters + SimulationParameters ¶ms = sd.get_simulation_parameters(); + params.include_optical_errors = true; + params.include_sun_shape_errors = true; + params.max_number_of_rays = NRAYS * 100; + params.number_of_rays = NRAYS; + params.seed = 1; + + // Run Ray Trace + EmbreeRunner runner; + runner.set_number_of_threads(1); + RunnerStatus sts; + sts = runner.initialize(); + EXPECT_EQ(sts, RunnerStatus::SUCCESS); + sts = runner.setup_simulation(&sd); + EXPECT_EQ(sts, RunnerStatus::SUCCESS); + + auto t0 = std::chrono::high_resolution_clock::now(); + sts = runner.run_simulation(); + auto t1 = std::chrono::high_resolution_clock::now(); + EXPECT_EQ(sts, RunnerStatus::SUCCESS); + + std::chrono::duration dur = t1 - t0; + + std::cout << "Time: " << dur.count() << " ms" << std::endl; + + SimulationResult result; + sts = runner.report_simulation(&result, 0); + EXPECT_EQ(sts, RunnerStatus::SUCCESS); + EXPECT_EQ(result.get_number_of_records(), NRAYS); + + element_id absorber_id = 6285; + int_fast64_t nabsorbed = count_element_event(result, absorber_id, RayEvent::ABSORB); + int_fast64_t nreflect = count_element_event(result, absorber_id, RayEvent::REFLECT); + int_fast64_t nevents = nabsorbed + nreflect; + + std::cout << "Total: " << nevents + << "\nAbsorb: " << nabsorbed << " (" + << static_cast(nabsorbed) / nevents << ")" + << "\nReflect: " << nreflect << " (" + << static_cast(nreflect) / nevents << ")" + << std::endl; +} diff --git a/google-tests/unit-tests/simulation_runner/native_runner/CMakeLists.txt b/google-tests/unit-tests/simulation_runner/native_runner/CMakeLists.txt index 3b29bff1..aa1ff572 100644 --- a/google-tests/unit-tests/simulation_runner/native_runner/CMakeLists.txt +++ b/google-tests/unit-tests/simulation_runner/native_runner/CMakeLists.txt @@ -20,6 +20,7 @@ set(NATIVE_RUNNER_TEST_SRC sphere_calculator_test.cpp cylinder_calculator_test.cpp calculator_factory_test.cpp + thread_manager_test.cpp native_runner_test.cpp tower_demo.cpp ) diff --git a/google-tests/unit-tests/simulation_runner/native_runner/calculator_factory_test.cpp b/google-tests/unit-tests/simulation_runner/native_runner/calculator_factory_test.cpp index f0529142..33ec71c0 100644 --- a/google-tests/unit-tests/simulation_runner/native_runner/calculator_factory_test.cpp +++ b/google-tests/unit-tests/simulation_runner/native_runner/calculator_factory_test.cpp @@ -39,9 +39,19 @@ class CalculatorFactoryTest : public ::testing::Test struct UnknownSurface : public Surface { - UnknownSurface() : Surface(SurfaceType::SURFACE_UNKNOWN){} + UnknownSurface() : Surface(SurfaceType::SURFACE_UNKNOWN) {} ~UnknownSurface() {} - virtual void write_json(nlohmann::ordered_json& json) const override {} + virtual void bounding_box(const double x_minmax[2], + const double y_minmax[2], + double &z_min, + double &z_max) const override + { + z_max = 1.0; + z_min = 0.0; + return; + } + virtual surface_ptr make_copy() const override { return make_surface(*this); } + virtual void write_json(nlohmann::ordered_json &json) const override {} }; TEST_F(CalculatorFactoryTest, SingletonBehavior) diff --git a/google-tests/unit-tests/simulation_runner/native_runner/native_runner_multithreading_test.cpp b/google-tests/unit-tests/simulation_runner/native_runner/native_runner_multithreading_test.cpp index b13fcdf4..15a570de 100644 --- a/google-tests/unit-tests/simulation_runner/native_runner/native_runner_multithreading_test.cpp +++ b/google-tests/unit-tests/simulation_runner/native_runner/native_runner_multithreading_test.cpp @@ -37,7 +37,7 @@ TEST(NativeRunner, StatusAndCancelMultiThread) sts = runner.status_simulation(); EXPECT_EQ(sts, RunnerStatus::RUNNING); - double prog; + double prog = -1.0; std::this_thread::sleep_for(std::chrono::milliseconds(500)); sts = runner.status_simulation(&prog); EXPECT_EQ(sts, RunnerStatus::RUNNING); diff --git a/google-tests/unit-tests/simulation_runner/native_runner/native_runner_test.cpp b/google-tests/unit-tests/simulation_runner/native_runner/native_runner_test.cpp index f89a283c..b28e95e2 100644 --- a/google-tests/unit-tests/simulation_runner/native_runner/native_runner_test.cpp +++ b/google-tests/unit-tests/simulation_runner/native_runner/native_runner_test.cpp @@ -15,6 +15,7 @@ #include #include #include +#include #include "common.hpp" #include "count_absorbed_native.h" @@ -26,30 +27,6 @@ using SolTrace::NativeRunner::NativeRunner; using SolTrace::NativeRunner::TRayData; using SolTrace::NativeRunner::TSystem; -int_fast64_t count_element_event(const SimulationResult &res, element_id el, RayEvent rev) -{ - int_fast64_t count = 0; - - for (auto ray_idx = 0; - ray_idx < res.get_number_of_records(); - ++ray_idx) - { - auto rr = res[ray_idx]; - for (auto event_idx = 0; - event_idx < rr->get_number_of_interactions(); - ++event_idx) - { - if (rr->get_event(event_idx) == rev && - rr->get_element(event_idx) == el) - { - ++count; - } - } - } - - return count; -} - TEST(RandomNumberGenerator, SingleNumberMersenneTwister) { MTRand myrng(1); @@ -75,13 +52,13 @@ TEST(NativeRunnerTypes, TSun) EXPECT_EQ(sys->Sun.ShapeIndex, SunShape::PILLBOX); } -TEST(NativeRunnerTypes, MakeElement) -{ -} +// TEST(NativeRunnerTypes, MakeElement) +// { +// } -TEST(NativeRunnerTypes, MakeStage) -{ -} +// TEST(NativeRunnerTypes, MakeStage) +// { +// } // TEST(NativeRunnerTypes, TElement) // { @@ -130,12 +107,12 @@ TEST(NativeRunnerTypes, MakeStage) // // TODO: Implement test // } -TEST(ThreadManager, Logging) +TEST(TraceLogger, Logging) { - SolTrace::NativeRunner::ThreadManager manager; - manager.error_log("This is a test message to test logging"); + SolTrace::NativeRunner::TraceLogger logger; + logger.error_log("This is a test message to test logging"); std::stringstream ss; - manager.print_log(ss); + logger.print_log(ss); EXPECT_GT(ss.str().size(), 0); } diff --git a/google-tests/unit-tests/simulation_runner/native_runner/thread_manager_test.cpp b/google-tests/unit-tests/simulation_runner/native_runner/thread_manager_test.cpp new file mode 100644 index 00000000..0e420c40 --- /dev/null +++ b/google-tests/unit-tests/simulation_runner/native_runner/thread_manager_test.cpp @@ -0,0 +1,154 @@ +#include + +#include +#include + +#include +#include +#include + +using SolTrace::NativeRunner::make_thread_manager; +using SolTrace::NativeRunner::make_trace_logger; + +using SolTrace::NativeRunner::thread_manager_ptr; +using SolTrace::NativeRunner::ThreadManager; +using SolTrace::NativeRunner::trace_logger_ptr; + +// "Runs" nsteps step where each step takes at least tlen milliseconds +// checking the manager to see if it needs to terminate early after +// each step +ThreadManager::ThreadStatus timed_task(thread_manager_ptr manager, + trace_logger_ptr logger, + unsigned thread_id, + uint_fast64_t nsteps, + uint_fast64_t tlen) +{ + uint_fast64_t count = 0; + double thread_progress = 0.0; + + while (count < nsteps) + { + std::this_thread::sleep_for(std::chrono::milliseconds(tlen)); + ++count; + thread_progress = static_cast(count) / nsteps; + manager->progress_update(thread_id, count); + + if (manager->terminate(thread_id)) + { + return ThreadManager::ThreadStatus::CANCEL; + } + } + + return ThreadManager::ThreadStatus::SUCCESS; +} + +// Waits for tlen milliseconds and returns an error +ThreadManager::ThreadStatus error_task(thread_manager_ptr manager, + trace_logger_ptr logger, + unsigned thread_id, + uint_fast64_t tlen) +{ + std::this_thread::sleep_for(std::chrono::milliseconds(tlen)); + return ThreadManager::ThreadStatus::ERROR; +} + +TEST(ThreadManager, CleanExit) +{ + const unsigned NTHREADS = 2; + const uint_fast64_t NSTEPS = 5; + const uint_fast64_t T_MS = 100; + + auto logger = make_trace_logger(); + auto manager = make_thread_manager(logger); + + manager->initialize(); + for (auto thid = 0; thid < NTHREADS; ++thid) + { + auto f = std::async(std::launch::async, + timed_task, + manager, + logger, + thid, + NSTEPS, + T_MS); + manager->manage(thid, std::move(f)); + } + + std::this_thread::sleep_for(std::chrono::milliseconds(2 * T_MS)); + double progress = 0.0; + ThreadManager::ThreadStatus sts = manager->status(&progress); + EXPECT_EQ(sts, ThreadManager::ThreadStatus::RUNNING); + EXPECT_GT(progress, 0.0); + + sts = manager->monitor_until_completion(); + EXPECT_EQ(sts, ThreadManager::ThreadStatus::SUCCESS); +} + +TEST(ThreadManager, ErrorExit) +{ + const unsigned NTHREADS = 2; + const uint_fast64_t NSTEPS = 5; + const uint_fast64_t T_MS = 100; + + auto logger = make_trace_logger(); + auto manager = make_thread_manager(logger); + + manager->initialize(); + for (auto thid = 0; thid < NTHREADS; ++thid) + { + ThreadManager::future f; + if (thid < NTHREADS - 1) + { + f = std::async(std::launch::async, + timed_task, + manager, + logger, + thid, + NSTEPS, + T_MS); + manager->manage(thid, std::move(f)); + } + else + { + f = std::async(std::launch::async, + error_task, + manager, + logger, + thid, + T_MS); + manager->manage(thid, std::move(f)); + } + } + + ThreadManager::ThreadStatus sts = manager->monitor_until_completion(); + EXPECT_EQ(sts, ThreadManager::ThreadStatus::ERROR); +} + +TEST(ThreadManager, CancelExit) +{ + const unsigned NTHREADS = 2; + const uint_fast64_t NSTEPS = 5; + const uint_fast64_t T_MS = 100; + + auto logger = make_trace_logger(); + auto manager = make_thread_manager(logger); + + manager->initialize(); + for (auto thid = 0; thid < NTHREADS; ++thid) + { + auto f = std::async(std::launch::async, + timed_task, + manager, + logger, + thid, + NSTEPS, + T_MS); + manager->manage(thid, std::move(f)); + } + + std::this_thread::sleep_for(std::chrono::milliseconds(2 * T_MS)); + + manager->cancel(); + ThreadManager::ThreadStatus sts = manager->monitor_until_completion(); + EXPECT_EQ(sts, ThreadManager::ThreadStatus::CANCEL); +}