ordered1d.hpp

Go to the documentation of this file.

/* 
 * This file is part of PLaSK (https://plask.app) by Photonics Group at TUL
 * Copyright (c) 2022 Lodz University of Technology
 * 
 * This program is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, version 3.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
 * GNU General Public License for more details.
 */
#ifndef PLASK__RECTILINEAR1D_H
#define PLASK__RECTILINEAR1D_H
 
#include <vector>
#include <algorithm>
#include <initializer_list>
 
#include "mesh.hpp"
#include "../vec.hpp"
#include "../log/log.hpp"
#include "../utils/iterators.hpp"
#include "../utils/interpolation.hpp"
 
#include "axis1d.hpp"
 
namespace plask {
 
class PLASK_API OrderedAxis: public MeshAxis {
 
    std::vector<double> points;
 
    void sortPointsAndRemoveNonUnique(double min_dist);
 
public:
 
    struct WarningOff {
        OrderedAxis* axis;
        bool prev_state;
        WarningOff(OrderedAxis& axis): axis(&axis), prev_state(axis.warn_too_close) { axis.warn_too_close = false; }
        WarningOff(OrderedAxis* axis): axis(axis), prev_state(axis->warn_too_close) { axis->warn_too_close = false; }
        WarningOff(const shared_ptr<OrderedAxis>& axis): axis(axis.get()), prev_state(axis->warn_too_close) { axis->warn_too_close = false; }
        ~WarningOff() { axis->warn_too_close = prev_state; }
    };
 
    bool warn_too_close;
 
    constexpr static double MIN_DISTANCE = 1e-6; // 1 picometer
 
    typedef double PointType;
 
    typedef std::vector<double>::const_iterator native_const_iterator;
 
    native_const_iterator begin() const { return points.begin(); }
 
    native_const_iterator end() const { return points.end(); }
 
    const std::vector<double>& getPointsVector() const { return points; }
 
    native_const_iterator find(double to_find) const;
 
    std::size_t findIndex(double to_find) const override { return std::size_t(find(to_find) - begin()); }
 
    native_const_iterator findUp(double to_find) const;
 
    std::size_t findUpIndex(double to_find) const override { return std::size_t(findUp(to_find) - begin()); }
 
    native_const_iterator findNearest(double to_find) const;
 
    std::size_t findNearestIndex(double to_find) const override { return findNearest(to_find) - begin(); }
 
    OrderedAxis(): warn_too_close(true) {}
 
    OrderedAxis(const OrderedAxis& src): points(src.points), warn_too_close(true) {}
 
    OrderedAxis(OrderedAxis&& src): points(std::move(src.points)), warn_too_close(true) {}
 
    OrderedAxis(const MeshAxis& src): points(src.size()), warn_too_close(true) {
        if (src.isIncreasing())
            std::copy(src.begin(), src.end(), points.begin());
        else
            std::reverse_copy(src.begin(), src.end(), points.begin());
    }
 
    OrderedAxis(std::initializer_list<PointType> points, double min_dist=MIN_DISTANCE);
 
    OrderedAxis(const std::vector<PointType>& points, double min_dist=MIN_DISTANCE);
 
    OrderedAxis(std::vector<PointType>&& points, double min_dist=MIN_DISTANCE);
 
    OrderedAxis& operator=(const OrderedAxis& src);
 
    OrderedAxis& operator=(OrderedAxis&& src);
 
    OrderedAxis& operator=(const MeshAxis& src);
 
    bool operator==(const OrderedAxis& to_compare) const;
 
    void writeXML(XMLElement& object) const override;
 
    std::size_t size() const override { return points.size(); }
 
    double at(std::size_t index) const override { assert(index < points.size()); return points[index]; }
 
    // @return true only if there are no points in mesh
    //bool empty() const { return points.empty(); }
 
    bool addPoint(double new_node_cord, double min_dist);
 
    bool addPoint(double new_node_cord) {
        return addPoint(new_node_cord, MIN_DISTANCE);
    }
 
    void removePoint(std::size_t index);
 
    void removePoints(std::size_t start, std::size_t stop);
 
    void removePoints(std::size_t start, std::size_t stop, std::ptrdiff_t step);
 
    template <typename IteratorT>
    void addOrderedPoints(IteratorT begin, IteratorT end, std::size_t points_count_hint, double min_dist=MIN_DISTANCE);
 
    //TODO use iterator traits and write version for input iterator
    template <typename RandomAccessIteratorT>
    void addOrderedPoints(RandomAccessIteratorT begin, RandomAccessIteratorT end) { addOrderedPoints(begin, end, end - begin); }
 
    void addPointsLinear(double first, double last, std::size_t points_count);
 
    //const double& operator[](std::size_t index) const { return points[index]; }
 
    void clear();
 
    shared_ptr<MeshAxis> clone() const override;
 
    bool isIncreasing() const override { return true; }
 
    template <typename RandomAccessContainer>
    auto interpolateLinear(const RandomAccessContainer& data, double point) const -> typename std::remove_reference<decltype(data[0])>::type;
 
};
 
// OrderedAxis method templates implementation
template <typename RandomAccessContainer>
inline auto OrderedAxis::interpolateLinear(const RandomAccessContainer& data, double point) const -> typename std::remove_reference<decltype(data[0])>::type {
    std::size_t index = findIndex(point);
    if (index == size()) return data[index-1];     //TODO what should it do if mesh is empty?
    if (index == 0 || points[index] == point) return data[index]; // hit exactly
    // here: points[index-1] < point < points[index]
    return interpolation::linear(points[index-1], data[index-1], points[index], data[index], point);
}
 
template <typename IteratorT>
inline void OrderedAxis::addOrderedPoints(IteratorT begin, IteratorT end, std::size_t points_count_hint, double min_dist) {
    std::vector<double> result;
    result.reserve(this->size() + points_count_hint);
    std::set_union(this->points.begin(), this->points.end(), begin, end, std::back_inserter(result));
    this->points = std::move(result);
    // Remove points too close to each other
    auto almost_equal = [min_dist, this](const double& x, const double& y) -> bool {
        bool remove = std::abs(x-y) < min_dist;
        if (warn_too_close && remove) writelog(LOG_WARNING, "Points in ordered mesh too close, skipping point at {0}", y);
        return remove;
    };
    this->points.erase(std::unique(this->points.begin(), this->points.end(), almost_equal), this->points.end());
    fireResized();
}
 
typedef OrderedAxis OrderedMesh1D;
 
shared_ptr<OrderedMesh1D> readRectilinearMeshAxis(XMLReader& reader);
 
}   // namespace plask
 
#endif // PLASK__RECTILINEAR1D_H