tools/mesh_converter/MeshElement.hpp

#ifndef ALLPIX_MESHELEMENT_H
#define ALLPIX_MESHELEMENT_H

#include <Eigen/Eigen>
#include <array>
#include <cmath>
#include <utility>

#include "core/utils/log.h"
#include "octree/Octree.hpp"

namespace mesh_converter {
    // Simple point class to store data
    class Point {
    public:
        Point() noexcept = default;
        Point(double px, double py, double pz) noexcept : x(px), y(py), z(pz), dim(3){};
        Point(double py, double pz) noexcept : y(py), z(pz), dim(2){};

        bool isFinite() const { return std::isfinite(x) && std::isfinite(y) && std::isfinite(z); }

        double x{0}, y{0}, z{0};
        unsigned int dim{0};

        friend std::ostream& operator<<(std::ostream& out, const Point& pt) {
            out << "(" << pt.x << "," << pt.y << "," << pt.z << ")";
            return out;
        }
    };

    class MeshElement {
    public:
        MeshElement() = delete;

        MeshElement(size_t dimension,
                    const std::array<Point, 4>& vertices_tetrahedron,
                    const std::array<Point, 4>& efield_vertices_tetrahedron)
            : dimension_(dimension), vertices_(vertices_tetrahedron), e_field_(efield_vertices_tetrahedron) {
            calculate_volume();
        }

        bool isValid(double volume_cut, Point& qp) const;

        Point getObservable(Point& qp) const;

        std::string print(Point& qp) const;

    private:
        double get_distance(size_t index, Point& qp) const;

        void calculate_volume();

        double get_sub_volume(size_t index, Point& p) const;

        size_t dimension_{3};
        std::array<Point, 4> vertices_{};
        std::array<Point, 4> e_field_{};

        double volume_{0};
    };

    class Combination {
        const std::vector<Point>* grid_;
        const std::vector<Point>* field_;
        Point reference_;
        Point result_;
        bool valid_{};
        double cut_;

        std::array<Point, 4> grid_elements;
        std::array<Point, 4> field_elements;

    public:
        explicit Combination(const std::vector<Point>* points,
                             const std::vector<Point>* field,
                             const Point& q,
                             const double volume_cut)
            : grid_(points), field_(field), reference_(q), cut_(volume_cut) {}

        template <class It> bool operator()(It begin, It end) {
            // Dimensionality is number of iterator elements minus one:
            size_t dimensions = static_cast<size_t>(end - begin) - 1;
            size_t idx = 0;

            LOG(TRACE) << "Constructing " << dimensions << "D element at " << reference_ << " with mesh points:";
            for(; begin < end; begin++) {
                LOG(TRACE) << "\t\t" << (*grid_)[*begin];
                grid_elements[idx] = (*grid_)[*begin];
                field_elements[idx++] = (*field_)[*begin];
            }

            MeshElement element(dimensions, grid_elements, field_elements);
            valid_ = element.isValid(cut_, reference_);
            if(valid_) {
                LOG(DEBUG) << element.print(reference_);
                result_ = element.getObservable(reference_);
                if(!result_.isFinite()) {
                    LOG(WARNING) << "Interpolated result not a finite number at " << reference_;
                    return false;
                }
            }

            return valid_;
        }

        bool valid() const { return valid_; }

        const Point& result() const { return result_; }
    };

} // namespace mesh_converter

#endif // ALLPIX_MESHELEMENT_H