regular1d.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__REGULAR1D_H
#define PLASK__REGULAR1D_H
 
#include "mesh.hpp"
#include "../math.hpp"
#include "../utils/iterators.hpp"
#include "../utils/interpolation.hpp"
#include "../utils/stl.hpp"
 
#include "axis1d.hpp"
 
namespace plask {
 
class PLASK_API RegularAxis: public MeshAxis {
 
    double lo, _step;
    std::size_t points_count;
 
  public:
 
    typedef double PointType;
 
    typedef IndexedIterator<const RegularAxis, PointType> native_const_iterator;
 
    native_const_iterator begin() const { return native_const_iterator(this, 0); }
 
    native_const_iterator end() const { return native_const_iterator(this, points_count); }
 
    RegularAxis():
        lo(0.), _step(0.), points_count(0) {}
 
    RegularAxis(const RegularAxis& src):
        lo(src.lo), _step(src._step), points_count(src.points_count) {}
 
    RegularAxis(double first, double last, std::size_t points_count):
        lo(first), _step( (last - first) / ((points_count>1)?double(points_count-1):1.) ),
        points_count(points_count) {}
 
    RegularAxis& operator=(const RegularAxis& src);
 
    void reset(double first, double last, std::size_t points_count);
 
    double first() const { return lo; }
 
    double last() const { return lo + _step * double(points_count-1); }
 
    double step() const { return _step; }
 
    std::size_t size() const override { return points_count; }
 
    bool operator==(const RegularAxis& to_compare) const {
        return this->lo == to_compare.lo && this->_step == to_compare._step && this->points_count == to_compare.points_count;
    }
 
    bool operator!=(const RegularAxis& to_compare) const {
        return !(*this == to_compare);
    }
 
    bool empty() const override { return points_count == 0; }
 
    double operator[](std::size_t index) const { return lo + double(index) * _step; }
 
    double at(std::size_t index) const override { return lo + double(index) * _step; }
 
    void clear() {
        if (empty()) return;
        points_count = 0;
        fireResized();
    }
 
    std::size_t findIndex(double to_find) const override {
        return clamp(std::ptrdiff_t(std::ceil((to_find - lo) / _step)), std::ptrdiff_t(0), std::ptrdiff_t(points_count));
    }
 
    native_const_iterator find(double to_find) const {
        return begin() + findIndex(to_find);
    }
 
    std::size_t findUpIndex(double to_find) const override;
 
    native_const_iterator findUp(double to_find) const {
        return begin() + findUpIndex(to_find);
    }
 
    native_const_iterator findNearest(double to_find) const {
        return find_nearest_using_lower_bound(begin(), end(), to_find, find(to_find));
    }
 
    std::size_t findNearestIndex(double to_find) const override { return findNearest(to_find) - begin(); }
 
    shared_ptr<MeshAxis> clone() const override { return plask::make_shared<RegularAxis>(*this); }
 
    void writeXML(XMLElement& object) const override;
 
    bool isIncreasing() const override;
 
    shared_ptr<MeshAxis> getMidpointAxis() const override;
 
    template <typename RandomAccessContainer>
    auto interpolateLinear(const RandomAccessContainer& data, double point) -> typename std::remove_reference<decltype(data[0])>::type;
 
protected:
    bool hasSameNodes(const MeshD<1>& to_compare) const override;
 
};
 
// RegularAxis method templates implementation
template <typename RandomAccessContainer>
auto RegularAxis::interpolateLinear(const RandomAccessContainer& data, double point) -> typename std::remove_reference<decltype(data[0])>::type {
    std::size_t index = findIndex(point);
    if (index == size()) return data[index - 1];     //TODO what should do if mesh is empty?
    if (index == 0 || this->operator[](index) == point) return data[index]; //hit exactly
    // here: points[index-1] < point < points[index]
    return interpolation::linear(this->operator[](index-1), data[index-1], this->operator[](index), data[index], point);
}
 
typedef RegularAxis RegularMesh1D;
 
PLASK_API shared_ptr<RegularMesh1D> readRegularMeshAxis(XMLReader& reader);
 
}   // namespace plask
 
#endif // PLASK__REGULAR1D_H