rectangular2d.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__RECTANGULAR2D_H
#define PLASK__RECTANGULAR2D_H
 
#include <iterator>
 
#include "rectangular_common.hpp"
#include "../utils/interpolation.hpp"
#include "../geometry/object.hpp"
#include "../geometry/space.hpp"
#include "../math.hpp"
 
namespace plask {
 
class PLASK_API RectangularMesh2D: public RectangularMeshBase2D {
 
    typedef std::size_t index_ft(const RectangularMesh2D* mesh, std::size_t axis0_index, std::size_t axis1_index);
    typedef std::size_t index01_ft(const RectangularMesh2D* mesh, std::size_t mesh_index);
 
    // Our own virtual table, changeable at run-time:
    index_ft* index_f;
    index01_ft* index0_f;
    index01_ft* index1_f;
 
    const shared_ptr<MeshAxis>* minor_axis; 
    const shared_ptr<MeshAxis>* major_axis; 
 
    void onAxisChanged(Event& e);
 
    void setChangeSignal(const shared_ptr<MeshAxis>& axis) { if (axis) axis->changedConnectMethod(this, &RectangularMesh2D::onAxisChanged); }
    void unsetChangeSignal(const shared_ptr<MeshAxis>& axis) { if (axis) axis->changedDisconnectMethod(this, &RectangularMesh2D::onAxisChanged); }
 
  public:
 
    class PLASK_API Element {
        const RectangularMesh2D& mesh;
        std::size_t index0, index1; // probably this form allows us to do most operation the fastest in average, low indexes of element corner or just element indexes
 
        public:
 
        Element(const RectangularMesh2D& mesh, std::size_t index0, std::size_t index1): mesh(mesh), index0(index0), index1(index1) {}
 
        Element(const RectangularMesh2D& mesh, std::size_t elementIndex): mesh(mesh) {
            std::size_t v = mesh.getElementMeshLowIndex(elementIndex);
            index0 = mesh.index0(v);
            index1 = mesh.index1(v);
        }
 
        inline std::size_t getIndex0() const { return index0; }
 
        inline std::size_t getIndex1() const { return index1; }
 
        inline std::size_t getLowerIndex0() const { return index0; }
 
        inline std::size_t getLowerIndex1() const { return index1; }
 
        inline double getLower0() const { return mesh.axis[0]->at(index0); }
 
        inline double getLower1() const { return mesh.axis[1]->at(index1); }
 
        inline std::size_t getUpperIndex0() const { return index0+1; }
 
        inline std::size_t getUpperIndex1() const { return index1+1; }
 
        inline double getUpper0() const { return mesh.axis[0]->at(getUpperIndex0()); }
 
        inline double getUpper1() const { return mesh.axis[1]->at(getUpperIndex1()); }
 
        inline double getSize0() const { return getUpper0() - getLower0(); }
 
        inline double getSize1() const { return getUpper1() - getLower1(); }
 
        inline Vec<2, double> getSize() const { return getUpUp() - getLoLo(); }
 
        inline Vec<2, double> getMidpoint() const { return mesh.getElementMidpoint(index0, index1); }
 
        inline std::size_t getIndex() const { return mesh.getElementIndexFromLowIndexes(getLowerIndex0(), getLowerIndex1()); }
 
        inline Box2D toBox() const { return mesh.getElementBox(index0, index1); }
 
        bool contains(Vec<2, double> p) const { return toBox().contains(p); }
 
        inline double getVolume() const { return getSize0() * getSize1(); }
 
        inline double getArea() const { return getVolume(); }
 
        inline std::size_t getLoLoIndex() const { return mesh.index(getLowerIndex0(), getLowerIndex1()); }
 
        inline std::size_t getLoUpIndex() const { return mesh.index(getLowerIndex0(), getUpperIndex1()); }
 
        inline std::size_t getUpLoIndex() const { return mesh.index(getUpperIndex0(), getLowerIndex1()); }
 
        inline std::size_t getUpUpIndex() const { return mesh.index(getUpperIndex0(), getUpperIndex1()); }
 
        inline Vec<2, double> getLoLo() const { return mesh(getLowerIndex0(), getLowerIndex1()); }
 
        inline Vec<2, double> getLoUp() const { return mesh(getLowerIndex0(), getUpperIndex1()); }
 
        inline Vec<2, double> getUpLo() const { return mesh(getUpperIndex0(), getLowerIndex1()); }
 
        inline Vec<2, double> getUpUp() const { return mesh(getUpperIndex0(), getUpperIndex1()); }
 
    };
 
    class PLASK_API Elements {
 
        static inline Element deref(const RectangularMesh2D& mesh, std::size_t index) { return mesh.getElement(index); }
    public:
        typedef IndexedIterator<const RectangularMesh2D, Element, deref> const_iterator;
        typedef const_iterator iterator;
 
        const RectangularMesh2D* mesh;
 
        Elements(const RectangularMesh2D* mesh): mesh(mesh) {}
 
        Element operator[](std::size_t i) const { return Element(*mesh, i); }
 
        Element operator()(std::size_t i0, std::size_t i1) const { return Element(*mesh, i0, i1); }
 
        std::size_t size() const { return mesh->getElementsCount(); }
 
        const_iterator begin() const { return const_iterator(mesh, 0); }
 
        const_iterator end() const { return const_iterator(mesh, size()); }
 
    };
 
    class ElementMesh;
 
    const shared_ptr<MeshAxis> axis[2];
 
    Elements elements() const { return Elements(this); }
    Elements getElements() const { return elements(); }
 
    Element element(std::size_t i0, std::size_t i1) const { return Element(*this, i0, i1); }
    Element getElement(std::size_t i0, std::size_t i1) const { return element(i0, i1); }
 
    Element element(std::size_t i) const { return Element(*this, i); }
 
    Element getElement(std::size_t i) const { return element(i); }
 
    enum IterationOrder { ORDER_10, ORDER_01 };
 
    void setIterationOrder(IterationOrder order);
 
    IterationOrder getIterationOrder() const;
 
    const char* getIterationOrderAsArray() const {
        return &"\1\0\0\1"[getIterationOrder() * 2];
        //       |   |   |
    };
 
    const char* getAxisToIterationOrder() const {
        return &"\1\0\0\1"[getIterationOrder() * 2];
        //       |   |   |
    };
 
    void setOptimalIterationOrder() {
        setIterationOrder(axis[0]->size() > axis[1]->size() ? ORDER_01 : ORDER_10);
    }
 
    explicit RectangularMesh2D(IterationOrder iterationOrder = ORDER_01);
 
    RectangularMesh2D(shared_ptr<MeshAxis> axis0, shared_ptr<MeshAxis> axis1, IterationOrder iterationOrder = ORDER_01);
 
    void reset(shared_ptr<MeshAxis> axis0, shared_ptr<MeshAxis> axis1, IterationOrder iterationOrder = ORDER_01);
 
    RectangularMesh2D(const RectangularMesh2D& src, bool clone_axes = false);
 
    void reset(const RectangularMesh2D& src, bool clone_axes = false);
 
    RectangularMesh2D& operator=(const RectangularMesh2D& src) { reset(src, true); return *this; }
 
    RectangularMesh2D& operator=(RectangularMesh2D&& src) { reset(src, false); return *this; }
 
    ~RectangularMesh2D();
 
    void setAxis(std::size_t axis_nr, shared_ptr<MeshAxis> new_val, bool fireResized = true);
 
    const shared_ptr<MeshAxis> getAxis0() const { return axis[0]; }
 
    void setAxis0(shared_ptr<MeshAxis> a0) { setAxis(0, a0); }
 
    const shared_ptr<MeshAxis> getAxis1() const { return axis[1]; }
 
    void setAxis1(shared_ptr<MeshAxis> a1) { setAxis(1, a1); }
 
 
 
    /*
     * Construct mesh with is based on given 1D meshes
     *
     * @param mesh0 mesh for the first coordinate, or constructor argument for the first coordinate mesh
     * @param mesh1 mesh for the second coordinate, or constructor argument for the second coordinate mesh
     * @param iterationOrder iteration order
     */
    /*template <typename Mesh0CtorArg, typename Mesh1CtorArg>
    RectangularMesh(Mesh0CtorArg&& mesh0, Mesh1CtorArg&& mesh1, IterationOrder iterationOrder = ORDER_10):
        axis0(std::forward<Mesh0CtorArg>(mesh0)), axis1(std::forward<Mesh1CtorArg>(mesh1)) elements(this) {
        axis0->owner = this; axis[1]->owner = this;
        setIterationOrder(iterationOrder); }*/
 
    const shared_ptr<MeshAxis>& tran() const { return axis[0]; }
 
    void setTran(shared_ptr<MeshAxis> a0) { setAxis0(a0); }
 
    const shared_ptr<MeshAxis>& vert() const { return axis[1]; }
 
    void setVert(shared_ptr<MeshAxis> a1) { setAxis1(a1); }
 
    const shared_ptr<MeshAxis>& ee_x() const { return axis[0]; }
 
    const shared_ptr<MeshAxis>& ee_y() const { return axis[1]; }
 
    const shared_ptr<MeshAxis>& rad_r() const { return axis[0]; }
 
    const shared_ptr<MeshAxis>& rad_z() const { return axis[1]; }
 
    const shared_ptr<MeshAxis>& getAxis(size_t n) const {
        if (n >= 2) throw Exception("bad axis number");
        return axis[n];
    }
 
    inline const shared_ptr<MeshAxis> majorAxis() const {
        return *major_axis;
    }
 
    inline std::size_t majorAxisIndex() const {
        return getIterationOrder() == IterationOrder::ORDER_01 ? 0 : 1;
    }
 
    inline const shared_ptr<MeshAxis> minorAxis() const {
        return *minor_axis;
    }
 
    inline std::size_t minorAxisIndex() const {
        return getIterationOrder() == IterationOrder::ORDER_01 ? 1 : 0;
    }
 
    bool operator==(const RectangularMesh2D& to_compare) const {
        if (empty()) return to_compare.empty();
        // here we know that axes are non-empty
        if (*axis[0] != *to_compare.axis[0] || *axis[1] != *to_compare.axis[1]) return false;
        // here we know that axes are equal
        return this->getIterationOrder() == to_compare.getIterationOrder() || (axis[0]->size() == 1) || (axis[1]->size() == 1);
    }
 
    bool operator!=(const RectangularMesh2D& to_compare) const {
        return !(*this == to_compare);
    }
 
    std::size_t size() const override { return axis[0]->size() * axis[1]->size(); }
 
    std::size_t getMaxSize() const { return std::max(axis[0]->size(), axis[1]->size()); }
 
    std::size_t getMinSize() const { return std::min(axis[0]->size(), axis[1]->size()); }
 
    void writeXML(XMLElement& object) const override;
 
    bool empty() const override { return axis[0]->empty() || axis[1]->empty(); }
 
    inline std::size_t index(std::size_t axis0_index, std::size_t axis1_index) const {
        return index_f(this, axis0_index, axis1_index);
    }
 
    inline std::size_t index(const Vec<2, std::size_t>& indexes) const {
        return index(indexes[0], indexes[1]);
    }
 
    inline std::size_t index0(std::size_t mesh_index) const {
        return index0_f(this, mesh_index);
    }
 
    inline std::size_t index1(std::size_t mesh_index) const {
        return index1_f(this, mesh_index);
    }
 
    inline Vec<2, std::size_t> indexes(std::size_t mesh_index) const {
        return Vec<2, std::size_t>(index0(mesh_index), index1(mesh_index));
    }
 
    inline std::size_t majorIndex(std::size_t mesh_index) const {
        return mesh_index / (*minor_axis)->size();
    }
 
    inline std::size_t minorIndex(std::size_t mesh_index) const {
        return mesh_index % (*minor_axis)->size();
    }
 
    inline Vec<2, double> at(std::size_t index0, std::size_t index1) const {
        return Vec<2, double>(axis[0]->at(index0), axis[1]->at(index1));
    }
 
    Vec<2, double> at(std::size_t index) const override {
        return Vec<2, double>(axis[0]->at(index0(index)), axis[1]->at(index1(index)));
    }
 
    inline Vec<2,double> operator[](std::size_t index) const {
        return Vec<2, double>(axis[0]->at(index0(index)), axis[1]->at(index1(index)));
    }
 
    inline Vec<2,double> operator()(std::size_t axis0_index, std::size_t axis1_index) const {
        return Vec<2, double>(axis[0]->at(axis0_index), axis[1]->at(axis1_index));
    }
 
    /*void clear() {
        axis0->clear();
        axis[1]->clear();
    }*/
 
    shared_ptr<RectangularMesh2D::ElementMesh> getElementMesh() const;
 
    template <typename RandomAccessContainer>
    auto interpolateLinear(const RandomAccessContainer& data, const Vec<2>& point, const InterpolationFlags& flags) const
        -> typename std::remove_reference<decltype(data[0])>::type
    {
        auto p = flags.wrap(point);
 
        size_t index0_lo, index0_hi;
        double left, right;
        bool invert_left, invert_right;
        prepareInterpolationForAxis(*axis[0], flags, p.c0, 0, index0_lo, index0_hi, left, right, invert_left, invert_right);
 
        size_t index1_lo, index1_hi;
        double bottom, top;
        bool invert_bottom, invert_top;
        prepareInterpolationForAxis(*axis[1], flags, p.c1, 1, index1_lo, index1_hi, bottom, top, invert_bottom, invert_top);
 
        typename std::remove_const<typename std::remove_reference<decltype(data[0])>::type>::type
            data_lb = data[index(index0_lo, index1_lo)],
            data_rb = data[index(index0_hi, index1_lo)],
            data_rt = data[index(index0_hi, index1_hi)],
            data_lt = data[index(index0_lo, index1_hi)];
        if (invert_left)   { data_lb = flags.reflect(0, data_lb); data_lt = flags.reflect(0, data_lt); }
        if (invert_right)  { data_rb = flags.reflect(0, data_rb); data_rt = flags.reflect(0, data_rt); }
        if (invert_top)    { data_lt = flags.reflect(1, data_lt); data_rt = flags.reflect(1, data_rt); }
        if (invert_bottom) { data_lb = flags.reflect(1, data_lb); data_rb = flags.reflect(1, data_rb); }
 
        return flags.postprocess(point, interpolation::bilinear(left, right, bottom, top, data_lb, data_rb, data_rt, data_lt, p.c0, p.c1));
    }
 
    template <typename RandomAccessContainer>
    auto interpolateNearestNeighbor(const RandomAccessContainer& data, const Vec<2>& point, const InterpolationFlags& flags) const
        -> typename std::remove_reference<decltype(data[0])>::type {
        auto p = flags.wrap(point);
        prepareNearestNeighborInterpolationForAxis(*axis[0], flags, p.c0, 0);
        prepareNearestNeighborInterpolationForAxis(*axis[1], flags, p.c1, 1);
        return flags.postprocess(point, data[this->index(axis[0]->findNearestIndex(p.c0), axis[1]->findNearestIndex(p.c1))]);
    }
 
    std::size_t getElementsCount0() const {
        const std::size_t s = axis[0]->size();
        return (s != 0)? s-1 : 0;
    }
 
    std::size_t getElementsCount1() const {
        const std::size_t s = axis[1]->size();
        return (s != 0)? s-1 : 0;
    }
 
    std::size_t getElementsCount() const {
        return (axis[0]->size() != 0 && axis[1]->size() != 0)?
            (axis[0]->size()-1) * (axis[1]->size()-1) : 0;
    }
 
    bool isLowIndexOfElement(std::size_t meshIndex) const {
        return index0(meshIndex) + 1 < axis[0]->size() && index1(meshIndex) + 1 < axis[1]->size();
    }
 
    std::size_t getElementIndexFromLowIndex(std::size_t mesh_index_of_el_bottom_left) const {
        return mesh_index_of_el_bottom_left - mesh_index_of_el_bottom_left / (*minor_axis)->size();
    }
 
    std::size_t getElementIndexFromLowIndexes(std::size_t axis0_index, std::size_t axis1_index) const {
        return getElementIndexFromLowIndex(index(axis0_index, axis1_index));
    }
 
    std::size_t getElementMeshLowIndex(std::size_t element_index) const {
        return element_index + (element_index / ((*minor_axis)->size()-1));
    }
 
    Vec<2, std::size_t> getElementMeshLowIndexes(std::size_t element_index) const {
        std::size_t bl_index = getElementMeshLowIndex(element_index);
        return Vec<2, std::size_t>(index0(bl_index), index1(bl_index));
    }
 
    double getElementArea(std::size_t index0, std::size_t index1) const {
        return (axis[0]->at(index0+1) - axis[0]->at(index0)) * (axis[1]->at(index1+1) - axis[1]->at(index1));
    }
 
    double getElementArea(std::size_t element_index) const {
        std::size_t bl_index = getElementMeshLowIndex(element_index);
        return getElementArea(index0(bl_index), index1(bl_index));
    }
 
    double getElementMidpoint0(std::size_t index0) const { return 0.5 * (axis[0]->at(index0) + axis[0]->at(index0+1)); }
 
    double getElementMidpoint1(std::size_t index1) const { return 0.5 * (axis[1]->at(index1) + axis[1]->at(index1+1)); }
 
    Vec<2, double> getElementMidpoint(std::size_t index0, std::size_t index1) const {
        return vec(getElementMidpoint0(index0), getElementMidpoint1(index1));
    }
 
    Vec<2, double> getElementMidpoint(std::size_t element_index) const {
        std::size_t bl_index = getElementMeshLowIndex(element_index);
        return getElementMidpoint(index0(bl_index), index1(bl_index));
    }
 
    Box2D getElementBox(std::size_t index0, std::size_t index1) const {
        return Box2D(axis[0]->at(index0), axis[1]->at(index1), axis[0]->at(index0+1), axis[1]->at(index1+1));
    }
 
    Box2D getElementBox(std::size_t element_index) const {
        std::size_t bl_index = getElementMeshLowIndex(element_index);
        return getElementBox(index0(bl_index), index1(bl_index));
    }
 
  protected:
 
    bool hasSameNodes(const MeshD<2> &to_compare) const override;
 
 
  private:
 
    // Common code for: left, right, bottom, top boundaries:
    struct BoundaryIteratorImpl: public BoundaryNodeSetImpl::IteratorImpl {
 
        const RectangularMesh2D &mesh;
 
        std::size_t line;
 
        std::size_t index;
 
        BoundaryIteratorImpl(const RectangularMesh2D& mesh, std::size_t line, std::size_t index): mesh(mesh), line(line), index(index) {}
 
        void increment() override { ++index; }
 
        bool equal(const typename BoundaryNodeSetImpl::IteratorImpl& other) const override {
            return index == static_cast<const BoundaryIteratorImpl&>(other).index;
        }
 
    };
 
    // iterator over vertical line (from bottom to top). for left and right boundaries
    struct VerticalIteratorImpl: public BoundaryIteratorImpl {
 
        VerticalIteratorImpl(const RectangularMesh2D& mesh, std::size_t line, std::size_t index): BoundaryIteratorImpl(mesh, line, index) {}
 
        std::size_t dereference() const override { return this->mesh.index(this->line, this->index); }
 
        std::unique_ptr<typename BoundaryNodeSetImpl::IteratorImpl> clone() const override {
            return std::unique_ptr<typename BoundaryNodeSetImpl::IteratorImpl>(new VerticalIteratorImpl(*this));
        }
    };
 
    // iterator over horizonstal line (from left to right), for bottom and top boundaries
    struct HorizontalIteratorImpl: public BoundaryIteratorImpl {
 
        HorizontalIteratorImpl(const RectangularMesh2D& mesh, std::size_t line, std::size_t index): BoundaryIteratorImpl(mesh, line, index) {}
 
        std::size_t dereference() const override { return this->mesh.index(this->index, this->line); }
 
        std::unique_ptr<typename BoundaryNodeSetImpl::IteratorImpl> clone() const override {
            return std::unique_ptr<typename BoundaryNodeSetImpl::IteratorImpl>(new HorizontalIteratorImpl(*this));
        }
    };
 
    struct VerticalBoundary: public BoundaryNodeSetWithMeshImpl<RectangularMesh2D> {
 
        typedef typename BoundaryNodeSetImpl::Iterator Iterator;
 
        std::size_t line;
 
        VerticalBoundary(const RectangularMesh2D& mesh, std::size_t line_axis0): BoundaryNodeSetWithMeshImpl<RectangularMesh2D>(mesh), line(line_axis0) {}
 
        //virtual LeftBoundary* clone() const { return new LeftBoundary(); }
 
        bool contains(std::size_t mesh_index) const override {
            return this->mesh.index0(mesh_index) == line;
        }
 
        Iterator begin() const override {
            return Iterator(new VerticalIteratorImpl(this->mesh, line, 0));
        }
 
        Iterator end() const override {
            return Iterator(new VerticalIteratorImpl(this->mesh, line, this->mesh.axis[1]->size()));
        }
 
        std::size_t size() const override {
            return this->mesh.axis[1]->size();
        }
    };
 
 
    struct VerticalBoundaryInRange: public BoundaryNodeSetWithMeshImpl<RectangularMesh2D> {
 
        typedef typename BoundaryNodeSetImpl::Iterator Iterator;
 
        std::size_t line, beginInLineIndex, endInLineIndex;
 
        VerticalBoundaryInRange(const RectangularMesh2D& mesh, std::size_t line_axis0, std::size_t beginInLineIndex, std::size_t endInLineIndex)
            : BoundaryNodeSetWithMeshImpl<RectangularMesh2D>(mesh), line(line_axis0), beginInLineIndex(beginInLineIndex), endInLineIndex(endInLineIndex) {}
 
        //virtual LeftBoundary* clone() const { return new LeftBoundary(); }
 
        bool contains(std::size_t mesh_index) const override {
            return this->mesh.index0(mesh_index) == line && in_range(this->mesh.index1(mesh_index), beginInLineIndex, endInLineIndex);
        }
 
        Iterator begin() const override {
            return Iterator(new VerticalIteratorImpl(this->mesh, line, beginInLineIndex));
        }
 
        Iterator end() const override {
            return Iterator(new VerticalIteratorImpl(this->mesh, line, endInLineIndex));
        }
 
        std::size_t size() const override {
            return endInLineIndex - beginInLineIndex;
        }
    };
 
    struct HorizontalBoundary: public BoundaryNodeSetWithMeshImpl<RectangularMesh2D> {
 
        typedef typename BoundaryNodeSetImpl::Iterator Iterator;
 
        std::size_t line;
 
        HorizontalBoundary(const RectangularMesh2D& mesh, std::size_t line_axis1): BoundaryNodeSetWithMeshImpl<RectangularMesh2D>(mesh), line(line_axis1) {}
 
        //virtual TopBoundary* clone() const { return new TopBoundary(); }
 
        bool contains(std::size_t mesh_index) const override {
            return this->mesh.index1(mesh_index) == line;
        }
 
        Iterator begin() const override {
            return Iterator(new HorizontalIteratorImpl(this->mesh, line, 0));
        }
 
        Iterator end() const override {
            return Iterator(new HorizontalIteratorImpl(this->mesh, line, this->mesh.axis[0]->size()));
        }
 
        std::size_t size() const override {
            return this->mesh.axis[0]->size();
        }
    };
 
    struct HorizontalBoundaryInRange: public BoundaryNodeSetWithMeshImpl<RectangularMesh2D> {
 
        typedef typename BoundaryNodeSetImpl::Iterator Iterator;
 
        std::size_t line, beginInLineIndex, endInLineIndex;
 
        HorizontalBoundaryInRange(const RectangularMesh2D& mesh, std::size_t line_axis1, std::size_t beginInLineIndex, std::size_t endInLineIndex)
            : BoundaryNodeSetWithMeshImpl<RectangularMesh2D>(mesh), line(line_axis1), beginInLineIndex(beginInLineIndex), endInLineIndex(endInLineIndex) {}
        //virtual TopBoundary* clone() const { return new TopBoundary(); }
 
        bool contains(std::size_t mesh_index) const override {
            return this->mesh.index1(mesh_index) == line && in_range(this->mesh.index0(mesh_index), beginInLineIndex, endInLineIndex);
        }
 
        Iterator begin() const override {
            return Iterator(new HorizontalIteratorImpl(this->mesh, line, beginInLineIndex));
        }
 
        Iterator end() const override {
            return Iterator(new HorizontalIteratorImpl(this->mesh, line, endInLineIndex));
        }
 
        std::size_t size() const override {
            return endInLineIndex - beginInLineIndex;
        }
    };
 
    //TODO
    /*struct HorizontalLineBoundary: public BoundaryLogicImpl<RectangularMesh2D> {
 
        double height;
 
        bool contains(const RectangularMesh &mesh, std::size_t mesh_index) const {
            return mesh.index1(mesh_index) == mesh.axis[1]->findNearestIndex(height);
        }
    };*/
 
  public:     // boundaries:
 
    BoundaryNodeSet createVerticalBoundaryAtLine(std::size_t line_nr_axis0) const override;
 
    BoundaryNodeSet createVerticalBoundaryAtLine(std::size_t line_nr_axis0, std::size_t indexBegin, std::size_t indexEnd) const override;
 
    BoundaryNodeSet createVerticalBoundaryNear(double axis0_coord) const override;
 
    BoundaryNodeSet createVerticalBoundaryNear(double axis0_coord, double from, double to) const override;
 
    BoundaryNodeSet createLeftBoundary() const override;
 
    BoundaryNodeSet createRightBoundary() const override;
 
    BoundaryNodeSet createLeftOfBoundary(const Box2D& box) const override;
 
    BoundaryNodeSet createRightOfBoundary(const Box2D& box) const override;
 
    BoundaryNodeSet createBottomOfBoundary(const Box2D& box) const override;
 
    BoundaryNodeSet createTopOfBoundary(const Box2D& box) const override;
 
    BoundaryNodeSet createHorizontalBoundaryAtLine(std::size_t line_nr_axis1) const override;
 
    BoundaryNodeSet createHorizontalBoundaryAtLine(std::size_t line_nr_axis1, std::size_t indexBegin, std::size_t indexEnd) const override;
 
    BoundaryNodeSet createHorizontalBoundaryNear(double axis1_coord) const override;
 
    BoundaryNodeSet createHorizontalBoundaryNear(double axis1_coord, double from, double to) const override;
 
    BoundaryNodeSet createTopBoundary() const override;
 
    BoundaryNodeSet createBottomBoundary() const override;
};
 
 
class RectangularMesh2D::ElementMesh: public RectangularMesh2D {
 
    const RectangularMesh2D* originalMesh;
 
  public:
    template <typename... Args>
    ElementMesh(const RectangularMesh2D* originalMesh, Args... args): RectangularMesh2D(args...), originalMesh(originalMesh) {}
 
    template <typename RandomAccessContainer>
    auto interpolateNearestNeighbor(const RandomAccessContainer& data, const Vec<2>& point, const InterpolationFlags& flags) const
        -> typename std::remove_reference<decltype(data[0])>::type {
        auto p = flags.wrap(point);
        prepareNearestNeighborInterpolationForAxis(*originalMesh->axis[0], flags, p.c0, 0);
        prepareNearestNeighborInterpolationForAxis(*originalMesh->axis[1], flags, p.c1, 1);
        size_t i0 = originalMesh->axis[0]->findUpIndex(p.c0), i1 = originalMesh->axis[1]->findUpIndex(p.c1);
        if (i0 == originalMesh->axis[0]->size()) --i0;
        if (i0 != 0) --i0;
        if (i1 == originalMesh->axis[1]->size()) --i1;
        if (i1 != 0) --i1;
        return flags.postprocess(point, data[this->index(i0, i1)]);
    }
 
protected:
 
    bool hasSameNodes(const MeshD<2> &to_compare) const override;
 
};
 
 
template <typename SrcT, typename DstT>
struct InterpolationAlgorithm<RectangularMesh2D, SrcT, DstT, INTERPOLATION_LINEAR> {
    static LazyData<DstT> interpolate(const shared_ptr<const RectangularMesh2D>& src_mesh, const DataVector<const SrcT>& src_vec,
                                      const shared_ptr<const MeshD<2>>& dst_mesh, const InterpolationFlags& flags) {
        if (src_mesh->axis[0]->size() == 0 || src_mesh->axis[1]->size() == 0)
            throw BadMesh("interpolate", "source mesh empty");
        return new LinearInterpolatedLazyDataImpl< DstT, RectangularMesh2D, SrcT >(src_mesh, src_vec, dst_mesh, flags);
    }
};
 
template <typename SrcT, typename DstT>
struct InterpolationAlgorithm<RectangularMesh2D, SrcT, DstT, INTERPOLATION_NEAREST> {
    static LazyData<DstT> interpolate(const shared_ptr<const RectangularMesh2D>& src_mesh, const DataVector<const SrcT>& src_vec,
                                      const shared_ptr<const MeshD<2>>& dst_mesh, const InterpolationFlags& flags) {
        if (src_mesh->axis[0]->size() == 0 || src_mesh->axis[1]->size() == 0)
            throw BadMesh("interpolate", "source mesh empty");
        return new NearestNeighborInterpolatedLazyDataImpl< DstT, RectangularMesh2D, SrcT >(src_mesh, src_vec, dst_mesh, flags);
    }
};
 
 
template <typename SrcT, typename DstT, InterpolationMethod method>
struct InterpolationAlgorithm<typename std::enable_if<method != INTERPOLATION_DEFAULT, RectangularMesh2D::ElementMesh>::type, SrcT, DstT, method> {
    static LazyData<DstT> interpolate(const shared_ptr<const RectangularMesh2D::ElementMesh>& src_mesh, const DataVector<const SrcT>& src_vec,
                                      const shared_ptr<const MeshD<2>>& dst_mesh, const InterpolationFlags& flags) {
        return InterpolationAlgorithm<RectangularMesh2D, SrcT, DstT, method>::interpolate(src_mesh, src_vec, dst_mesh, flags);
    }
};
 
template <typename SrcT, typename DstT>
struct InterpolationAlgorithm<RectangularMesh2D::ElementMesh, SrcT, DstT, INTERPOLATION_NEAREST> {
    static LazyData<DstT> interpolate(const shared_ptr<const RectangularMesh2D::ElementMesh>& src_mesh, const DataVector<const SrcT>& src_vec,
                                      const shared_ptr<const MeshD<2>>& dst_mesh, const InterpolationFlags& flags) {
        if (src_mesh->axis[0]->size() == 0 || src_mesh->axis[1]->size() == 0)
            throw BadMesh("interpolate", "source mesh empty");
        return new NearestNeighborInterpolatedLazyDataImpl<DstT, RectangularMesh2D::ElementMesh, SrcT>(src_mesh, src_vec, dst_mesh, flags);
    }
};
 
 
PLASK_API shared_ptr<RectangularMesh2D> make_rectangular_mesh(const RectangularMesh2D& to_copy);
inline shared_ptr<RectangularMesh2D> make_rectangular_mesh(shared_ptr<const RectangularMesh2D> to_copy) { return make_rectangular_mesh(*to_copy); }
 
} // namespace plask
 
#include "rectangular_spline.hpp"
 
#endif // PLASK__RECTANGULAR2D_H