rectangular_masked_common.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_MASKED_COMMON_H
#define PLASK__RECTANGULAR_MASKED_COMMON_H
 
#include <functional>
 
#include "rectangular.hpp"
#include "../utils/numbers_set.hpp"
 
#include <boost/thread.hpp>
 
 
namespace plask {
 
template <int DIM>
struct RectangularMaskedMeshBase: public RectangularMeshBase<DIM> {
 
    constexpr static double MIN_DISTANCE = 1e-9; // 1 femtometer
 
    RectangularMesh<DIM> fullMesh;
 
    using typename MeshD<DIM>::LocalCoords;
 
  protected:
 
    //typedef CompressedSetOfNumbers<std::uint32_t> Set;
    typedef CompressedSetOfNumbers<std::size_t> Set;
 
    Set nodeSet;
 
    Set elementSet;
 
    struct BoundaryIndexForAxis {
        std::size_t lo, up;
        void improveLo(std::size_t i) { if (i < lo) lo = i; }
        void improveUp(std::size_t i) { if (i > up) up = i; }
    };
    typedef BoundaryIndexForAxis BoundaryIndex[DIM];
    BoundaryIndex boundaryIndex;
 
    static void findIndexes(const MeshAxis& axis, double wrapped_point_coord, std::size_t& index_lo, std::size_t& index_hi) {
        index_hi = axis.findUpIndex(wrapped_point_coord);
        if (index_hi == axis.size()) --index_hi;    // p.c0 == axis0->at(axis0->size()-1)
        else if (index_hi == 0) index_hi = 1;
        index_lo = index_hi - 1;
    }
 
    static std::size_t nearest(double p, const MeshAxis& axis, std::size_t index_lo, std::size_t index_hi) {
        return p - axis.at(index_lo) <= axis.at(index_hi) - p ? index_lo : index_hi;
    }
 
    template <typename MaskedMeshType>
    struct ElementsBase {
 
        using Element = typename MaskedMeshType::Element;
 
        class const_iterator: public Set::ConstIteratorFacade<const_iterator, Element> {
 
            const MaskedMeshType* maskedMesh;
 
            Element dereference() const {
                return Element(*maskedMesh, this->getIndex(), this->getNumber());
            }
 
            friend class boost::iterator_core_access;
 
        public:
 
            template <typename... CtorArgs>
            explicit const_iterator(const MaskedMeshType& maskedMesh, CtorArgs&&... ctorArgs)
                : Set::ConstIteratorFacade<const_iterator, Element>(std::forward<CtorArgs>(ctorArgs)...), maskedMesh(&maskedMesh) {}
 
            const Set& set() const { return maskedMesh->elementSet; }
        };
 
        typedef const_iterator iterator;
 
        const MaskedMeshType* maskedMesh;
 
        explicit ElementsBase(const MaskedMeshType& maskedMesh): maskedMesh(&maskedMesh) {}
 
        std::size_t size() const { return maskedMesh->getElementsCount(); }
 
        const_iterator begin() const { return const_iterator(*maskedMesh, 0, maskedMesh->elementSet.segments.begin()); }
 
        const_iterator end() const { return const_iterator(*maskedMesh, size(), maskedMesh->elementSet.segments.end()); }
 
        Element operator[](std::size_t i) const { return Element(*maskedMesh, i); }
 
    };  // struct Elements
 
    template <typename MaskedMeshType>
    struct ElementMeshBase: MeshD<DIM> {
 
        using Element = typename MaskedMeshType::Element;
 
        class const_iterator: public Set::ConstIteratorFacade<const_iterator, LocalCoords> {
 
            const MaskedMeshType* originalMesh;
 
            LocalCoords dereference() const {
                return Element(*originalMesh, this->getIndex(), this->getNumber()).getMidpoint();
            }
 
            friend class boost::iterator_core_access;
 
        public:
 
            template <typename... CtorArgs>
            explicit const_iterator(const MaskedMeshType& originalMesh, CtorArgs&&... ctorArgs)
                : Set::ConstIteratorFacade<const_iterator, LocalCoords>(std::forward<CtorArgs>(ctorArgs)...), originalMesh(&originalMesh) {}
 
            const Set& set() const { return originalMesh->elementSet; }
        };
 
        typedef const_iterator iterator;
 
        const_iterator begin() const { return const_iterator(*originalMesh, 0, originalMesh->elementSet.segments.begin()); }
 
        const_iterator end() const { return const_iterator(*originalMesh, size(), originalMesh->elementSet.segments.end()); }
 
        const MaskedMeshType* originalMesh;
        RectangularMesh<DIM> fullMesh;
 
        explicit ElementMeshBase(const MaskedMeshType* originalMesh):
            originalMesh(originalMesh), fullMesh(*originalMesh->fullMesh.getElementMesh()) {}
 
        explicit ElementMeshBase(const MaskedMeshType& originalMesh): ElementMeshBase(&originalMesh) {}
 
        std::size_t size() const override { return originalMesh->getElementsCount(); }
 
        LocalCoords at(std::size_t index) const override {
            return fullMesh.at(originalMesh->elementSet.at(index));
        }
 
        bool empty() const override { return originalMesh->elementSet.empty(); }
 
    };  // ElementMeshBase
 
    void resetBoundyIndex() {
        for (int d = 0; d < DIM; ++d) { // prepare for finding indexes by subclass constructor:
            boundaryIndex[d].lo = this->fullMesh.axis[d]->size()-1;
            boundaryIndex[d].up = 0;
        }
        boundaryIndexInitialized = false;
    }
 
    void reset() {
        nodeSet.clear();
        elementSet.clear();
        resetBoundyIndex();
    }
 
  public:
 
    enum:std::size_t { NOT_INCLUDED = Set::NOT_INCLUDED };
 
    RectangularMaskedMeshBase() = default;
 
    RectangularMaskedMeshBase(const RectangularMesh<DIM>& rectangularMesh, Set nodeSet, bool clone_axes = false)
        : fullMesh(rectangularMesh, clone_axes), nodeSet(std::move(nodeSet)), elementSetInitialized(false)
    {
        resetBoundyIndex();
    }
 
    RectangularMaskedMeshBase(const RectangularMesh<DIM>& rectangularMesh, bool clone_axes = false)
        : fullMesh(rectangularMesh, clone_axes)
    {
        resetBoundyIndex();
    }
 
    class const_iterator: public CompressedSetOfNumbers<std::size_t>::ConstIteratorFacade<const_iterator, LocalCoords> {
 
        friend class boost::iterator_core_access;
 
        const RectangularMaskedMeshBase* mesh;
 
        LocalCoords dereference() const {
            return mesh->fullMesh.at(this->getNumber());
        }
 
    public:
 
        template <typename... CtorArgs>
        explicit const_iterator(const RectangularMaskedMeshBase& mesh, CtorArgs&&... ctorArgs)
            : CompressedSetOfNumbers<std::size_t>::ConstIteratorFacade<const_iterator, LocalCoords>(std::forward<CtorArgs>(ctorArgs)...), mesh(&mesh) {}
 
        const CompressedSetOfNumbers<std::size_t>& set() const { return mesh->nodeSet; }
    };
 
    typedef const_iterator iterator;
 
    const_iterator begin() const { return const_iterator(*this, 0, nodeSet.segments.begin()); }
    const_iterator end() const { return const_iterator(*this, size(), nodeSet.segments.end()); }
 
    LocalCoords at(std::size_t index) const override {
        return fullMesh.at(nodeSet.at(index));
    }
 
    std::size_t size() const override { return nodeSet.size(); }
 
    bool empty() const override { return nodeSet.empty(); }
 
    bool full() const { return nodeSet.size() == fullMesh.size(); }
 
    void selectAll(const RectangularMesh<DIM>& rectangularMesh, bool clone_axes = false) {
        fullMesh.reset(rectangularMesh, clone_axes);
        selectAll();
    }
 
    void selectAll() {
        this->nodeSet.assignRange(fullMesh.size());
        this->elementSet.assignRange(fullMesh.getElementsCount());
        elementSetInitialized = true;
        for (int d = 0; d < DIM; ++d) {
            boundaryIndex[d].lo = 0;
            boundaryIndex[d].up = fullMesh.axis[d]->size()-1;
        }
        boundaryIndexInitialized = true;
    }
 
    inline std::size_t index(const Vec<DIM, std::size_t>& indexes) const {
        return nodeSet.indexOf(fullMesh.index(indexes));
    }
 
    inline std::size_t index0(std::size_t mesh_index) const {
        return fullMesh.index0(nodeSet.at(mesh_index));
    }
 
    inline std::size_t index1(std::size_t mesh_index) const {
        return fullMesh.index1(nodeSet.at(mesh_index));
    }
 
    inline Vec<DIM, std::size_t> indexes(std::size_t mesh_index) const {
        return fullMesh.indexes(nodeSet.at(mesh_index));
    }
 
    inline std::size_t majorIndex(std::size_t mesh_index) const {
        return fullMesh.majorIndex(nodeSet.at(mesh_index));
    }
 
    inline std::size_t minorIndex(std::size_t mesh_index) const {
        return fullMesh.minorIndex(nodeSet.at(mesh_index));
    }
 
    std::size_t getElementsCount0() const {
        return fullMesh.getElementsCount0();
    }
 
    std::size_t getElementsCount1() const {
        return fullMesh.getElementsCount1();
    }
 
    std::size_t getElementsCount() const {
        return ensureHasElements().size();
    }
 
    std::size_t getElementIndexFromLowIndex(std::size_t mesh_index_of_el_bottom_left) const {
        return ensureHasElements().indexOf(fullMesh.getElementIndexFromLowIndex(nodeSet.at(mesh_index_of_el_bottom_left)));
    }
 
    std::size_t getElementMeshLowIndex(std::size_t element_index) const {
        return nodeSet.indexOf(fullMesh.getElementMeshLowIndex(ensureHasElements().at(element_index)));
    }
 
    Vec<DIM, std::size_t> getElementMeshLowIndexes(std::size_t element_index) const {
        return fullMesh.getElementMeshLowIndexes(ensureHasElements().at(element_index));
    }
 
    double getElementArea(std::size_t element_index) const {
        return fullMesh.getElementArea(ensureHasElements().at(element_index));
    }
 
    double getElementMidpoint0(std::size_t index0) const {
        return fullMesh.getElementMidpoint0(index0);
    }
 
    double getElementMidpoint1(std::size_t index1) const {
        return fullMesh.getElementMidpoint1(index1);
    }
 
    Vec<DIM, double> getElementMidpoint(std::size_t element_index) const {
        return fullMesh.getElementMidpoint(ensureHasElements().at(element_index));
    }
 
    typename Primitive<DIM>::Box getElementBox(std::size_t element_index) const {
        return fullMesh.getElementBox(ensureHasElements().at(element_index));
    }
 
  protected:    // constructing elementSet from nodes set (element is chosen when all its vertices are chosen) on-deamand
 
    DontCopyThisField<boost::mutex> writeMutex;
 
    bool elementSetInitialized = true;
 
    bool boundaryIndexInitialized;
 
  private:
 
    /*bool restVerticesIncluded(const RectangularMesh2D::Element& el) const {
        return //nodeSet.includes(el.getLoLoIndex()) &&
               nodeSet.includes(el.getUpLoIndex()) &&
               nodeSet.includes(el.getLoUpIndex()) &&
               nodeSet.includes(el.getUpUpIndex());
    }
    bool restVerticesIncluded(const RectangularMesh3D::Element& el) const {
        return //nodeSet.includes(el.getLoLoLoIndex()) &&
               nodeSet.includes(el.getUpLoLoIndex()) &&
               nodeSet.includes(el.getLoUpLoIndex()) &&
               nodeSet.includes(el.getLoLoUpIndex()) &&
               nodeSet.includes(el.getLoUpUpIndex()) &&
               nodeSet.includes(el.getUpLoUpIndex()) &&
               nodeSet.includes(el.getUpUpLoIndex()) &&
               nodeSet.includes(el.getUpUpUpIndex());
    }
    void calculateElements() {
        boost::lock_guard<boost::mutex> lock((boost::mutex&)writeElementSet);
        if (elementSetInitialized) return;  // another thread has initialized elementSet just when we waited for mutex
        // TODO faster implementation
        for (auto lowNode: nodeSet) {
            if (!fullMesh.isLowIndexOfElement(lowNode)) continue;
            auto elementIndex = fullMesh.getElementIndexFromLowIndex(lowNode);
            if (restVerticesIncluded(fullMesh.getElement(elementIndex))) elementSet.push_back(elementIndex);
        }
        elementSetInitialized = true;
    }*/ // ^ generic algorithm, slow but probably correct
 
    template <int d = DIM>
    typename std::enable_if<d == 2>::type calculateElements() {
        boost::lock_guard<boost::mutex> lock((boost::mutex&)writeMutex);
        if (elementSetInitialized) return;  // another thread has initialized elementSet just when we waited for mutex
 
        if (fullMesh.axis[0]->size() <= 1 || fullMesh.axis[1]->size() <= 1) {
            elementSetInitialized = true;
            return;
        }
 
        elementSet = nodeSet.transformed([] (std::size_t&, std::size_t& e) { --e; });   // same as nodeSet.intersected(nodeSet.shiftedLeft(1))
        auto minor_axis_size = fullMesh.minorAxis()->size();
        elementSet = elementSet.intersection(elementSet.shiftedLeft(minor_axis_size));
        // now elementSet includes all low indexes which have other corners of elements (plus some indexes in the last column)
        // we have to transform low indexes to indexes of elements:
        elementSet = elementSet.transformed([&, minor_axis_size] (std::size_t& b, std::size_t& e) {  // here: 0 <= b < e
            if (e % minor_axis_size == 0) --e;  // end of segment cannot lie at the last column, as getElementIndexFromLowIndex confuses last column with the first element in the next row
            // no need to fix b (if b % minor_axis_size == 0) as rounding in getElementIndexFromLowIndex will give good value (as getElementIndexFromLowIndex(b)==getElementIndexFromLowIndex(b+1))
            b = fullMesh.getElementIndexFromLowIndex(b);
            e = fullMesh.getElementIndexFromLowIndex(e);
        });
 
        elementSetInitialized = true;
    }
 
    template <int d = DIM>
    typename std::enable_if<d == 3>::type calculateElements() {
        boost::lock_guard<boost::mutex> lock((boost::mutex&)writeMutex);
        if (elementSetInitialized) return;  // another thread has initialized elementSet just when we waited for mutex
 
        if (fullMesh.axis[0]->size() <= 1 || fullMesh.axis[1]->size() <= 1 || fullMesh.axis[2]->size() <= 1) {
            elementSetInitialized = true;
            return;
        }
 
        elementSet = nodeSet.transformed([] (std::size_t&, std::size_t& e) { --e; });   // same as nodeSet.intersected(nodeSet.shiftedLeft(1))
        auto minor_axis_size = fullMesh.minorAxis()->size();
        elementSet = elementSet.intersection(elementSet.shiftedLeft(minor_axis_size));
        auto medium_axis_size = fullMesh.mediumAxis()->size();
        elementSet = elementSet.intersection(elementSet.shiftedLeft(minor_axis_size*medium_axis_size));
        // now elementSet includes all low indexes which have other corners of elements (plus some indexes in the last column)
        // we have to transform low indexes to indexes of elements:
        elementSet = elementSet.transformed([&, minor_axis_size, medium_axis_size] (std::size_t& b, std::size_t& e) {  // here: 0 <= b < e
            const std::size_t b_div = b / minor_axis_size;
            if (b_div % medium_axis_size == (medium_axis_size-1))
                b = (b_div+1) * minor_axis_size;    // first index in the next plane
            b = fullMesh.getElementIndexFromLowIndex(b);
 
            const std::size_t e_div = (e-1) / minor_axis_size;    // e-1 is the last index of the range
            if (e_div % medium_axis_size == (medium_axis_size-1))
                e = e_div * minor_axis_size - 1;    // -1 to not be divisible by minor_axis_size
            else if (e % minor_axis_size == 0) --e;
            e = fullMesh.getElementIndexFromLowIndex(e);
        });
 
        elementSetInitialized = true;
    }
 
    template <int d = DIM>
    typename std::enable_if<d == 2>::type calculateBoundaryIndex() {
        boost::lock_guard<boost::mutex> lock((boost::mutex&)writeMutex);
        if (boundaryIndexInitialized) return;  // another thread has initialized boundaryIndex just when we waited for mutex
 
        const auto minor = fullMesh.minorAxisIndex();
        const auto major = fullMesh.majorAxisIndex();
        nodeSet.forEachSegment([this, minor, major] (std::size_t b, std::size_t e) {
            const auto indexes_f = fullMesh.indexes(b);
            const auto indexes_l = fullMesh.indexes(e-1);
            if (indexes_f[major] != indexes_l[major]) {
                boundaryIndex[minor].lo = 0;
                boundaryIndex[minor].up = fullMesh.minorAxis()->size()-1;
            } else {
                boundaryIndex[minor].improveLo(indexes_f[minor]);
                boundaryIndex[minor].improveUp(indexes_l[minor]);
            }
            boundaryIndex[major].improveLo(indexes_f[major]);
            boundaryIndex[major].improveUp(indexes_l[major]);
        });
 
        boundaryIndexInitialized = true;
    }
 
    template <int d = DIM>
    typename std::enable_if<d == 3>::type calculateBoundaryIndex() {
        boost::lock_guard<boost::mutex> lock((boost::mutex&)writeMutex);
        if (boundaryIndexInitialized) return;  // another thread has initialized boundaryIndex just when we waited for mutex
 
        const auto minor = fullMesh.minorAxisIndex();
        const auto medium = fullMesh.mediumAxisIndex();
        const auto major = fullMesh.majorAxisIndex();
        nodeSet.forEachSegment([this, minor, medium, major] (std::size_t b, std::size_t e) {
            const auto indexes_f = fullMesh.indexes(b);
            const auto indexes_l = fullMesh.indexes(e-1);
            if (indexes_f[major] != indexes_l[major]) {
                boundaryIndex[minor].lo = 0;
                boundaryIndex[minor].up = fullMesh.minorAxis()->size()-1;
                boundaryIndex[medium].lo = 0;
                boundaryIndex[medium].up = fullMesh.mediumAxis()->size()-1;
            } else {
                if (indexes_f[medium] != indexes_l[medium]) {
                    boundaryIndex[minor].lo = 0;
                    boundaryIndex[minor].up = fullMesh.minorAxis()->size()-1;
                } else {   // here:   indexes_f[minor] <= indexes_l[minor]
                    boundaryIndex[minor].improveLo(indexes_f[minor]);
                    boundaryIndex[minor].improveUp(indexes_l[minor]);
                }
                // indexes_f[major] == indexes_l[major]   =>   indexes_f[medium] <= indexes_l[medium]
                boundaryIndex[medium].improveLo(indexes_f[medium]);
                boundaryIndex[medium].improveUp(indexes_l[medium]);
            }
            boundaryIndex[major].improveLo(indexes_f[major]);
            boundaryIndex[major].improveUp(indexes_l[major]);
        });
 
        boundaryIndexInitialized = true;
    }
 
  protected:
    const Set& ensureHasElements() const {
        if (!elementSetInitialized) const_cast<RectangularMaskedMeshBase<DIM>*>(this)->calculateElements();
        return elementSet;
    }
 
    const BoundaryIndex& ensureHasBoundaryIndex() const {
        if (!boundaryIndexInitialized) const_cast<RectangularMaskedMeshBase<DIM>*>(this)->calculateBoundaryIndex();
        return boundaryIndex;
    }
};
 
 
}   // namespace plask
 
#endif // PLASK__RECTANGULAR_MASKED_COMMON_H