rectangular_masked2d.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__RECTANGULAR_MASKED2D_H
#define PLASK__RECTANGULAR_MASKED2D_H
 
#include "rectangular_masked_common.hpp"
 
namespace plask {
 
struct PLASK_API RectangularMaskedMesh2D: public RectangularMaskedMeshBase<2> {
 
    typedef std::function<bool(const RectangularMesh2D::Element&)> Predicate;
 
    class PLASK_API Element {
 
        const RectangularMaskedMesh2D& maskedMesh;
 
        //std::uint32_t elementNumber;    ///< index of element in oryginal mesh
        std::size_t index0, index1; // probably this form allows to do most operation fastest in average, low indexes of element corner or just element indexes
 
        mutable std::size_t elementIndex;
 
        const RectangularMesh<2>& fullMesh() const { return maskedMesh.fullMesh; }
 
    public:
 
        enum: std::size_t { UNKNOWN_ELEMENT_INDEX = std::numeric_limits<std::size_t>::max() };
 
        Element(const RectangularMaskedMesh2D& maskedMesh, std::size_t elementIndex, std::size_t index0, std::size_t index1)
            : maskedMesh(maskedMesh), index0(index0), index1(index1), elementIndex(elementIndex)
        {
        }
 
        Element(const RectangularMaskedMesh2D& maskedMesh, std::size_t elementIndex, std::size_t elementIndexOfFullMesh)
            : maskedMesh(maskedMesh), elementIndex(elementIndex)
        {
            const std::size_t v = fullMesh().getElementMeshLowIndex(elementIndexOfFullMesh);
            index0 = fullMesh().index0(v);
            index1 = fullMesh().index1(v);
        }
 
        Element(const RectangularMaskedMesh2D& maskedMesh, std::size_t elementIndex)
            : Element(maskedMesh, elementIndex, maskedMesh.elementSet.at(elementIndex))
        {}
 
        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 fullMesh().axis[0]->at(index0); }
 
        inline double getLower1() const { return fullMesh().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 fullMesh().axis[0]->at(getUpperIndex0()); }
 
        inline double getUpper1() const { return fullMesh().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 maskedMesh.getElementMidpoint(index0, index1); }
 
        inline std::size_t getIndex() const {
            if (elementIndex == UNKNOWN_ELEMENT_INDEX)
                elementIndex = maskedMesh.getElementIndexFromLowIndexes(getLowerIndex0(), getLowerIndex1());
            return elementIndex;
        }
 
        inline Box2D toBox() const { return maskedMesh.getElementBox(index0, index1); }
 
        inline double getVolume() const { return getSize0() * getSize1(); }
 
        inline double getArea() const { return getVolume(); }
 
        inline std::size_t getLoLoIndex() const { return maskedMesh.index(getLowerIndex0(), getLowerIndex1()); }
 
        inline std::size_t getLoUpIndex() const { return maskedMesh.index(getLowerIndex0(), getUpperIndex1()); }
 
        inline std::size_t getUpLoIndex() const { return maskedMesh.index(getUpperIndex0(), getLowerIndex1()); }
 
        inline std::size_t getUpUpIndex() const { return maskedMesh.index(getUpperIndex0(), getUpperIndex1()); }
 
        inline Vec<2, double> getLoLo() const { return maskedMesh(getLowerIndex0(), getLowerIndex1()); }
 
        inline Vec<2, double> getLoUp() const { return maskedMesh(getLowerIndex0(), getUpperIndex1()); }
 
        inline Vec<2, double> getUpLo() const { return maskedMesh(getUpperIndex0(), getLowerIndex1()); }
 
        inline Vec<2, double> getUpUp() const { return maskedMesh(getUpperIndex0(), getUpperIndex1()); }
 
    };  // class Element
 
    struct PLASK_API Elements: ElementsBase<RectangularMaskedMesh2D> {
 
        explicit Elements(const RectangularMaskedMesh2D& mesh): ElementsBase(mesh) { mesh.ensureHasElements(); }
 
        Element operator()(std::size_t i0, std::size_t i1) const { return Element(*maskedMesh, Element::UNKNOWN_ELEMENT_INDEX, i0, i1); }
 
    };  // struct Elements
 
    struct PLASK_API ElementMesh: ElementMeshBase<RectangularMaskedMesh2D> {
 
        explicit ElementMesh(const RectangularMaskedMesh2D* originalMesh): ElementMeshBase<RectangularMaskedMesh2D>(originalMesh) {}
 
        inline std::size_t index(std::size_t axis0_index, std::size_t axis1_index) const {
            return originalMesh->elementSet.indexOf(fullMesh.index(axis0_index, axis1_index));
        }
 
 
        bool prepareInterpolation(const Vec<2>& point, Vec<2>& wrapped_point, std::size_t& index0_lo, std::size_t& index0_hi, std::size_t& index1_lo, std::size_t& index1_hi,
                                  const InterpolationFlags& flags) const {
            return originalMesh->prepareInterpolation(point, wrapped_point, index0_lo, index0_hi, index1_lo, index1_hi, flags);
        }
 
        // Convert to recctangular masked mesh
        RectangularMaskedMesh2D toMasked() const {
            return RectangularMaskedMesh2D(fullMesh, originalMesh->elementSet);
        }
 
        // Convert to recctangular masked mesh
        operator RectangularMaskedMesh2D() const {
            return RectangularMaskedMesh2D(fullMesh, originalMesh->elementSet);
        }
 
        template <typename RandomAccessContainer>
        auto interpolateLinear(const RandomAccessContainer& data, const Vec<2>& point, const InterpolationFlags& flags) const
            -> typename std::remove_reference<decltype(data[0])>::type {
            typedef typename std::remove_reference<decltype(data[0])>::type DataT;
            Vec<2> p;
            size_t index0, index0_hi, index1, index1_hi;
 
            if (!prepareInterpolation(point, p, index0, index0_hi, index1, index1_hi, flags))
                return NaN<decltype(data[0])>();
 
            Vec<2> pa = fullMesh.at(index0, index1);
 
            size_t step0 = (p.c0 < pa.c0)?
                (index0 == 0)? 0 : -1 :
                (index0_hi == fullMesh.axis[0]->size())? 0 : 1;
            size_t step1 = (p.c1 < pa.c1)?
                (index1 == 0)? 0 : -1 :
                (index1_hi == fullMesh.axis[1]->size())? 0 : 1;
 
            size_t index_aa = index(index0, index1), index_ab, index_ba, index_bb;
 
            typename std::remove_const<DataT>::type data_aa = data[index_aa], data_ab, data_ba, data_bb;
 
            if (step0 == 0 && step1 == 0) {
                index_ab = index_ba = index_bb = index_aa;
                data_ab = data_ba = data_bb = data_aa;
            } else {
                index_ab = index(index0, index1+step1);
                index_ba = index(index0+step0, index1);
                index_bb = index(index0+step0, index1+step1);
                data_ab = (index_ab != Element::UNKNOWN_ELEMENT_INDEX)? data[index_ab] : data_aa;
                data_ba = (index_ba != Element::UNKNOWN_ELEMENT_INDEX)? data[index_ba] : data_aa;
                data_bb = (index_bb != Element::UNKNOWN_ELEMENT_INDEX)? data[index_bb] : data_ab + data_ba - data_aa;
            }
 
            Vec<2> pb = fullMesh.at(index0+step0, index1+step1);
            if (step0 == 0) pb.c0 += 1.;
            if (step1 == 0) pb.c1 += 1.;
 
            return flags.postprocess(point,
                interpolation::bilinear(pa.c0, pb.c0, pa.c1, pb.c1,
                                        data_aa, data_ba, data_bb, data_ab, 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 {
            Vec<2> wrapped_point;
            std::size_t index0_lo, index0_hi, index1_lo, index1_hi;
 
            if (!originalMesh->prepareInterpolation(point, wrapped_point, index0_lo, index0_hi, index1_lo, index1_hi, flags))
                return NaN<decltype(data[0])>();
 
            return flags.postprocess(point, data[this->index(index0_lo, index1_lo)]);
        }
 
    };  // struct ElementMesh
 
    RectangularMaskedMesh2D() = default;
 
    void reset(const Predicate& predicate);
 
    RectangularMaskedMesh2D(const RectangularMesh<2>& fullMesh, const Predicate& predicate, bool clone_axes = false);
 
    void reset(const RectangularMesh<2>& fullMesh, const Predicate& predicate, bool clone_axes = false);
 
    RectangularMaskedMesh2D(const RectangularMesh<2>& fullMesh, const GeometryD<2>& geom, const std::function<bool(shared_ptr<const Material>)> materialPredicate, bool clone_axes = false)
        : RectangularMaskedMesh2D(fullMesh,
                                    [&](const RectangularMesh2D::Element& el) { return materialPredicate(geom.getMaterial(el.getMidpoint())); },
                                    clone_axes)
    {
    }
 
    void reset(const RectangularMesh<2>& rectangularMesh, const GeometryD<2>& geom, const std::function<bool(shared_ptr<const Material>)> materialPredicate, bool clone_axes = false) {
        reset(rectangularMesh,
              [&](const RectangularMesh2D::Element& el) { return materialPredicate(geom.getMaterial(el.getMidpoint())); },
              clone_axes);
    }
 
    RectangularMaskedMesh2D(const RectangularMesh<2>& rectangularMesh, const GeometryD<2>& geom, unsigned materialKinds, bool clone_axes = false)
        : RectangularMaskedMesh2D(rectangularMesh,
                                    [&](const RectangularMesh2D::Element& el) { return (geom.getMaterial(el.getMidpoint())->kind() & materialKinds) != 0; },
                                    clone_axes)
    {
    }
 
    void reset(const RectangularMesh<2>& rectangularMesh, const GeometryD<2>& geom, unsigned materialKinds, bool clone_axes = false) {
        reset(rectangularMesh,
             [&](const RectangularMesh2D::Element& el) { return (geom.getMaterial(el.getMidpoint())->kind() & materialKinds) != 0; },
             clone_axes);
    }
 
    RectangularMaskedMesh2D(const RectangularMesh<DIM>& rectangularMesh, Set nodeSet, bool clone_axes = false);
 
    Elements elements() const { return Elements(*this); }
    Elements getElements() const { return elements(); }
 
    Element element(std::size_t i0, std::size_t i1) const { ensureHasElements(); return Element(*this, Element::UNKNOWN_ELEMENT_INDEX, i0, i1); }
    Element getElement(std::size_t i0, std::size_t i1) const { return element(i0, i1); }
 
    Element element(std::size_t i) const { ensureHasElements(); return Element(*this, i); }
 
    Element getElement(std::size_t i) const { return element(i); }
 
    inline std::size_t index(std::size_t axis0_index, std::size_t axis1_index) const {
        return nodeSet.indexOf(fullMesh.index(axis0_index, axis1_index));
    }
 
    using RectangularMaskedMeshBase<2>::index;
    using RectangularMaskedMeshBase<2>::at;
 
    inline Vec<2, double> at(std::size_t index0, std::size_t index1) const {
        return fullMesh.at(index0, index1);
    }
 
    inline Vec<2,double> operator()(std::size_t axis0_index, std::size_t axis1_index) const {
        return fullMesh.operator()(axis0_index, axis1_index);
    }
 
    shared_ptr<RectangularMaskedMesh2D::ElementMesh> getElementMesh() const {
        return make_shared<RectangularMaskedMesh2D::ElementMesh>(this);
    }
 
  private:
 
    void initNodesAndElements(const RectangularMaskedMesh2D::Predicate &predicate);
 
    bool canBeIncluded(const Vec<2>& point) const {
        return
            fullMesh.axis[0]->at(0) - point[0] < MIN_DISTANCE && point[0] - fullMesh.axis[0]->at(fullMesh.axis[0]->size()-1) < MIN_DISTANCE &&
            fullMesh.axis[1]->at(0) - point[1] < MIN_DISTANCE && point[1] - fullMesh.axis[1]->at(fullMesh.axis[1]->size()-1) < MIN_DISTANCE;
    }
 
  public:
    bool prepareInterpolation(const Vec<2>& point, Vec<2>& wrapped_point,
                              std::size_t& index0_lo, std::size_t& index0_hi,
                              std::size_t& index1_lo, std::size_t& index1_hi,
                              const InterpolationFlags& flags) 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
    {
        Vec<2> wrapped_point;
        std::size_t index0_lo, index0_hi, index1_lo, index1_hi;
 
        if (!prepareInterpolation(point, wrapped_point, index0_lo, index0_hi, index1_lo, index1_hi, flags))
            return NaN<decltype(data[0])>();
 
        return flags.postprocess(point,
                                 interpolation::bilinear(
                                     fullMesh.axis[0]->at(index0_lo), fullMesh.axis[0]->at(index0_hi),
                                     fullMesh.axis[1]->at(index1_lo), fullMesh.axis[1]->at(index1_hi),
                                     data[index(index0_lo, index1_lo)],
                                     data[index(index0_hi, index1_lo)],
                                     data[index(index0_hi, index1_hi)],
                                     data[index(index0_lo, index1_hi)],
                                     wrapped_point.c0, wrapped_point.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
    {
        Vec<2> wrapped_point;
        std::size_t index0_lo, index0_hi, index1_lo, index1_hi;
 
        if (!prepareInterpolation(point, wrapped_point, index0_lo, index0_hi, index1_lo, index1_hi, flags))
            return NaN<decltype(data[0])>();
 
        return flags.postprocess(point,
                                 data[this->index(
                                     nearest(wrapped_point.c0, *fullMesh.axis[0], index0_lo, index0_hi),
                                     nearest(wrapped_point.c1, *fullMesh.axis[1], index1_lo, index1_hi)
                                 )]);
    }
 
    std::size_t getElementIndexFromLowIndexes(std::size_t axis0_index, std::size_t axis1_index) const {
        return ensureHasElements().indexOf(fullMesh.getElementIndexFromLowIndexes(axis0_index, axis1_index));
    }
 
    double getElementArea(std::size_t index0, std::size_t index1) const {
        return fullMesh.getElementArea(index0, index1);
    }
 
    Vec<2, double> getElementMidpoint(std::size_t index0, std::size_t index1) const {
        return fullMesh.getElementMidpoint(index0, index1);
    }
 
    Box2D getElementBox(std::size_t index0, std::size_t index1) const {
        return fullMesh.getElementBox(index0, index1);
    }
 
  protected:  // boundaries code:
 
    // Common code for: left, right, bottom, top boundaries:
    template <int CHANGE_DIR>
    struct BoundaryIteratorImpl: public plask::BoundaryNodeSetImpl::IteratorImpl {
 
        const RectangularMaskedMeshBase<2> &mesh;
 
        Vec<2, std::size_t> index;
 
        std::size_t endIndex;
 
        BoundaryIteratorImpl(const RectangularMaskedMeshBase<2>& mesh, Vec<2, std::size_t> index, std::size_t endIndex)
            : mesh(mesh), index(index), endIndex(endIndex)
        {
            // go to the first index existed in order to make dereference possible:
            while (this->index[CHANGE_DIR] < this->endIndex && mesh.index(this->index) == NOT_INCLUDED)
                ++this->index[CHANGE_DIR];
        }
 
        void increment() override {
            do {
                ++index[CHANGE_DIR];
            } while (index[CHANGE_DIR] < endIndex && mesh.index(index) == NOT_INCLUDED);
        }
 
        bool equal(const plask::BoundaryNodeSetImpl::IteratorImpl& other) const override {
            return index == static_cast<const BoundaryIteratorImpl&>(other).index;
        }
 
        std::size_t dereference() const override {
            return mesh.index(index);
        }
 
        std::unique_ptr<plask::BoundaryNodeSetImpl::IteratorImpl> clone() const override {
            return std::unique_ptr<plask::BoundaryNodeSetImpl::IteratorImpl>(new BoundaryIteratorImpl<CHANGE_DIR>(*this));
        }
 
    };
 
    template <int CHANGE_DIR>
    struct BoundaryNodeSetImpl: public BoundaryNodeSetWithMeshImpl<RectangularMaskedMeshBase<2>> {
 
        using typename BoundaryNodeSetWithMeshImpl<RectangularMaskedMeshBase<2>>::const_iterator;
 
        Vec<2, std::size_t> index;
 
        std::size_t endIndex;
 
        BoundaryNodeSetImpl(const RectangularMaskedMeshBase<2>& mesh, Vec<2, std::size_t> index, std::size_t endIndex)
            : BoundaryNodeSetWithMeshImpl<RectangularMaskedMeshBase<2>>(mesh), index(index), endIndex(endIndex) {}
 
        BoundaryNodeSetImpl(const RectangularMaskedMeshBase<2>& mesh, std::size_t index0, std::size_t index1, std::size_t endIndex)
            : BoundaryNodeSetWithMeshImpl<RectangularMaskedMeshBase<2>>(mesh), index(index0, index1), endIndex(endIndex) {}
 
        bool contains(std::size_t mesh_index) const override {
            if (mesh_index >= this->mesh.size()) return false;
            Vec<2, std::size_t> mesh_indexes = this->mesh.indexes(mesh_index);
            for (int i = 0; i < 2; ++i)
                if (i == CHANGE_DIR) {
                    if (mesh_indexes[i] < index[i] || mesh_indexes[i] >= endIndex) return false;
                } else
                    if (mesh_indexes[i] != index[i]) return false;
            return true;
        }
 
        const_iterator begin() const override {
            return Iterator(new BoundaryIteratorImpl<CHANGE_DIR>(this->mesh, index, endIndex));
        }
 
        const_iterator end() const override {
            Vec<2, std::size_t> index_end = index;
            index_end[CHANGE_DIR] = endIndex;
            return Iterator(new BoundaryIteratorImpl<CHANGE_DIR>(this->mesh, index_end, endIndex));
        }
    };
 
  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;
 
};
 
template <typename SrcT, typename DstT>
struct InterpolationAlgorithm<RectangularMaskedMesh2D, SrcT, DstT, INTERPOLATION_LINEAR> {
    static LazyData<DstT> interpolate(const shared_ptr<const RectangularMaskedMesh2D>& src_mesh, const DataVector<const SrcT>& src_vec,
                                      const shared_ptr<const MeshD<2>>& dst_mesh, const InterpolationFlags& flags) {
        if (src_mesh->empty()) throw BadMesh("interpolate", "source mesh empty");
        return new LinearInterpolatedLazyDataImpl<DstT, RectangularMaskedMesh2D, SrcT>(src_mesh, src_vec, dst_mesh, flags);
    }
};
 
template <typename SrcT, typename DstT>
struct InterpolationAlgorithm<RectangularMaskedMesh2D, SrcT, DstT, INTERPOLATION_NEAREST> {
    static LazyData<DstT> interpolate(const shared_ptr<const RectangularMaskedMesh2D>& src_mesh, const DataVector<const SrcT>& src_vec,
                                      const shared_ptr<const MeshD<2>>& dst_mesh, const InterpolationFlags& flags) {
        if (src_mesh->empty()) throw BadMesh("interpolate", "source mesh empty");
        return new NearestNeighborInterpolatedLazyDataImpl<DstT, RectangularMaskedMesh2D, SrcT>(src_mesh, src_vec, dst_mesh, flags);
    }
};
 
 
template <typename SrcT, typename DstT>
struct InterpolationAlgorithm<RectangularMaskedMesh2D::ElementMesh, SrcT, DstT, INTERPOLATION_LINEAR> {
    static LazyData<DstT> interpolate(const shared_ptr<const RectangularMaskedMesh2D::ElementMesh>& src_mesh, const DataVector<const SrcT>& src_vec,
                                      const shared_ptr<const MeshD<2>>& dst_mesh, const InterpolationFlags& flags) {
        if (src_mesh->empty()) throw BadMesh("interpolate", "source mesh empty");
        return new LinearInterpolatedLazyDataImpl<DstT, RectangularMaskedMesh2D::ElementMesh, SrcT>(src_mesh, src_vec, dst_mesh, flags);
    }
};
 
template <typename SrcT, typename DstT>
struct InterpolationAlgorithm<RectangularMaskedMesh2D::ElementMesh, SrcT, DstT, INTERPOLATION_NEAREST> {
    static LazyData<DstT> interpolate(const shared_ptr<const RectangularMaskedMesh2D::ElementMesh>& src_mesh, const DataVector<const SrcT>& src_vec,
                                      const shared_ptr<const MeshD<2>>& dst_mesh, const InterpolationFlags& flags) {
        if (src_mesh->empty()) throw BadMesh("interpolate", "source mesh empty");
        return new NearestNeighborInterpolatedLazyDataImpl<DstT, RectangularMaskedMesh2D::ElementMesh, SrcT>(src_mesh, src_vec, dst_mesh, flags);
    }
};
 
}   // namespace plask
 
#endif // PLASK__RECTANGULAR_MASKED2D_H