element.cpp

/*---------------------------------------------------------------------------*\
 *
 *  bitpit
 *
 *  Copyright (C) 2015-2021 OPTIMAD engineering Srl
 *
 *  -------------------------------------------------------------------------
 *  License
 *  This file is part of bitpit.
 *
 *  bitpit is free software: you can redistribute it and/or modify it
 *  under the terms of the GNU Lesser General Public License v3 (LGPL)
 *  as published by the Free Software Foundation.
 *
 *  bitpit 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 Lesser General Public
 *  License for more details.
 *
 *  You should have received a copy of the GNU Lesser General Public License
 *  along with bitpit. If not, see <http://www.gnu.org/licenses/>.
 *
\*---------------------------------------------------------------------------*/
 
#include <cassert>
#include <limits>
#include <set>
 
#include "bitpit_common.hpp"
#include "bitpit_CG.hpp"
#include "bitpit_operators.hpp"
 
#include "element.hpp"
 
namespace bitpit {
 
IBinaryStream& operator>>(IBinaryStream &buffer, Element &element)
{
    // Initialize the element
    ElementType type;
    buffer >> type;
 
    long id;
    buffer >> id;
 
    int connectSize;
    if (ReferenceElementInfo::hasInfo(type)) {
        element._initialize(id, type);
        connectSize = element.getConnectSize();
    } else {
        buffer >> connectSize;
        element._initialize(id, type, connectSize);
    }
 
    // Set connectivity
    buffer.read(reinterpret_cast<char *>(element.m_connect.get()), connectSize * sizeof(long));
 
    // Set PID
    int pid;
    buffer >> pid;
 
    element.setPID(pid);
 
    return buffer;
}
 
OBinaryStream& operator<<(OBinaryStream  &buffer, const Element &element)
{
    buffer << element.getType();
    buffer << element.getId();
 
    int connectSize = element.getConnectSize();
    if (!ReferenceElementInfo::hasInfo(element.m_type)) {
        buffer << connectSize;
    }
 
    buffer.write(reinterpret_cast<const char *>(element.m_connect.get()), connectSize * sizeof(long));
 
    buffer << element.getPID();
 
    return buffer;
}
 
Element::Tesselation::Tesselation()
    : m_nTiles(0)
{
}
 
int Element::Tesselation::importVertexCoordinates(const std::array<double, 3> &coordinates)
{
    return importVertexCoordinates(std::array<double, 3>(coordinates));
}
 
int Element::Tesselation::importVertexCoordinates(std::array<double, 3> &&coordinates)
{
    m_coordinates.push_back(std::move(coordinates));
 
    return (m_coordinates.size() - 1);
}
 
std::vector<int> Element::Tesselation::importVertexCoordinates(const std::array<double, 3> *coordinates, int nVertices)
{
    int nStoredVertices = m_coordinates.size();
    m_coordinates.reserve(nStoredVertices + nVertices);
 
    std::vector<int> ids(nVertices);
    for (int k = 0; k < nVertices; ++k) {
        m_coordinates.push_back(coordinates[k]);
        ids[k] = nStoredVertices++;
    }
 
    return ids;
}
 
void Element::Tesselation::importPolygon(const std::vector<int> &vertexIds)
{
    int nVertices = vertexIds.size();
    if (nVertices == 3 || nVertices == 4) {
        m_nTiles++;
        m_types.push_back(nVertices == 3 ? ElementType::TRIANGLE : ElementType::QUAD);
        m_connects.push_back(vertexIds);
 
        return;
    }
 
    // Add the centroid
    std::array<double, 3> centroid = {{0., 0., 0.}};
    for (int k = 0; k < nVertices; ++k) {
        centroid += m_coordinates[vertexIds[k]];
    }
    centroid = centroid / double(nVertices);
 
    int centroidId = importVertexCoordinates(std::move(centroid));
 
    // Decompose the polygon in triangles
    //
    // Each triangle is composed by the two vertices of a side and by the
    // centroid.
    ElementType tileType = ElementType::TRIANGLE;
    int nTileVertices = ReferenceElementInfo::getInfo(tileType).nVertices;
    int nSideVertices = ReferenceElementInfo::getInfo(ElementType::LINE).nVertices;
 
    int nSides = nVertices;
    m_types.resize(m_nTiles + nSides, tileType);
    m_connects.resize(m_nTiles + nSides, std::vector<int>(nTileVertices));
    for (int i = 0; i < nSides; ++i) {
        m_nTiles++;
        for (int k = 0; k < nSideVertices; ++k) {
            m_connects[m_nTiles - 1][k] = vertexIds[(i + k) % nVertices];
        }
        m_connects[m_nTiles - 1][nSideVertices] = centroidId;
    }
}
 
void Element::Tesselation::importPolyhedron(const std::vector<int> &vertexIds, const std::vector<std::vector<int>> &faceVertexIds)
{
    int nFaces    = faceVertexIds.size();
    int nVertices = vertexIds.size();
 
    // Generate the tesselation of the surface
    int nInitialTiles = m_nTiles;
    for (int i = 0; i < nFaces; ++i) {
        importPolygon(faceVertexIds[i]);
    }
    int nFinalTiles = m_nTiles;
 
    // Add the centroid of the element to the tesselation
    std::array<double, 3> centroid = {{0., 0., 0.}};
    for (int k = 0; k < nVertices; ++k) {
        centroid += m_coordinates[vertexIds[k]];
    }
    centroid = centroid / double(nVertices);
 
    int centroidTesselationId = importVertexCoordinates(std::move(centroid));
 
    // Decompose the polyhedron in prisms and pyramids.
    //
    // Each prism/pyramid has a tile of the surface tesselation as the
    // base and the centroid of the element as the apex. Since we have
    // already generated the surface tesselation, we can "convert" the
    // two-dimensional tiles of that tesselation in three-dimensional
    // tiles to obtain the volume tesselation.
    for (int i = nInitialTiles; i < nFinalTiles; ++i) {
        // Change the tile type
        ElementType &tileType = m_types[i];
        if (tileType == ElementType::TRIANGLE) {
            tileType = ElementType::TETRA;
        } else if (tileType == ElementType::QUAD) {
            tileType = ElementType::PYRAMID;
        } else {
            BITPIT_UNREACHABLE("Unsupported tile");
            throw std::runtime_error ("Unsupported tile");
        }
 
        // Fix the order of the connectivity the match the order of the
        // three-dimensional element
        if (tileType == ElementType::TETRA) {
            std::swap(m_connects[i][0], m_connects[i][2]);
        } else if (tileType == ElementType::PYRAMID) {
            std::swap(m_connects[i][1], m_connects[i][3]);
        }
 
        // Add the centroid to the tile connectivity
        m_connects[i].push_back(centroidTesselationId);
    }
}
 
int Element::Tesselation::getTileCount() const
{
    return m_nTiles;
}
 
ElementType Element::Tesselation::getTileType(int tile) const
{
    return m_types[tile];
}
 
std::vector<std::array<double, 3>> Element::Tesselation::getTileVertexCoordinates(int tile) const
{
    const ElementType tileType = getTileType(tile);
    const int nTileVertices = ReferenceElementInfo::getInfo(tileType).nVertices;
    const std::vector<int> &tileConnect = m_connects[tile];
 
    std::vector<std::array<double, 3>> coordinates(nTileVertices);
    for (int i = 0; i < nTileVertices; ++i) {
        coordinates[i] = m_coordinates[tileConnect[i]];
    }
 
    return coordinates;
}
 
const long Element::NULL_ID = std::numeric_limits<long>::min();
 
Element::Element()
{
    _initialize(NULL_ID, ElementType::UNDEFINED);
}
 
Element::Element(long id, ElementType type, int connectSize)
{
    _initialize(id, type, connectSize);
}
 
Element::Element(long id, ElementType type, std::unique_ptr<long[]> &&connectStorage)
{
    _initialize(id, type, std::move(connectStorage));
}
 
Element::Element(const Element &other)
{
    int connectSize;
    if (other.m_connect) {
        connectSize = other.getConnectSize();
    } else {
        connectSize = 0;
    }
 
    _initialize(other.m_id, other.m_type, connectSize);
 
    m_pid = other.m_pid;
 
    if (other.m_connect) {
        std::copy(other.m_connect.get(), other.m_connect.get() + connectSize, m_connect.get());
    }
}
 
Element & Element::operator=(const Element &other)
{
    Element tmp(other);
    swap(tmp);
 
    return *this;
}
 
void Element::swap(Element &other) noexcept
{
    std::swap(other.m_id, m_id);
    std::swap(other.m_type, m_type);
    std::swap(other.m_pid, m_pid);
    std::swap(other.m_connect, m_connect);
}
 
void Element::initialize(long id, ElementType type, int connectSize)
{
    _initialize(id, type, connectSize);
}
 
void Element::initialize(long id, ElementType type, std::unique_ptr<long[]> &&connectStorage)
{
    _initialize(id, type, std::move(connectStorage));
}
 
void Element::_initialize(long id, ElementType type, int connectSize)
{
    // Get previous connect size
    int previousConnectSize = 0;
    if (m_connect) {
        if (hasInfo()) {
            previousConnectSize = getInfo().nVertices;
        }
    }
 
    // Initialize connectivity storage
    if (ReferenceElementInfo::hasInfo(type)) {
        connectSize = ReferenceElementInfo::getInfo(type).nVertices;
    }
 
    std::unique_ptr<long[]> connectStorage;
    if (connectSize != previousConnectSize) {
        connectStorage = std::unique_ptr<long[]>(new long[connectSize]);
    } else {
        connectStorage = std::move(m_connect);
    }
 
    // Initialize element
    _initialize(id, type, std::move(connectStorage));
}
 
void Element::_initialize(long id, ElementType type, std::unique_ptr<long[]> &&connectStorage)
{
    // Set the id
    setId(id);
 
    // Set type
    setType(type);
 
    // Initialize PID
    setPID(0);
 
    // Initialize connectivity
    setConnect(std::move(connectStorage));
}
 
void Element::setId(long id)
{
    m_id = id;
}
 
long Element::getId() const
{
    return m_id;
}
 
bool Element::hasInfo() const
{
    return ReferenceElementInfo::hasInfo(m_type);
}
 
const ReferenceElementInfo & Element::getInfo() const
{
    return ReferenceElementInfo::getInfo(m_type);
}
 
void Element::setType(ElementType type)
{
    m_type = type;
}
 
ElementType Element::getType() const
{
    return m_type;
}
 
void Element::setPID(int pid)
{
    m_pid = pid;
}
 
int Element::getPID() const
{
    return m_pid;
}
 
void Element::setConnect(std::unique_ptr<long[]> &&connect)
{
    m_connect = std::move(connect);
}
 
void Element::unsetConnect()
{
    m_connect.reset(nullptr);
}
 
const long * Element::getConnect() const
{
    return m_connect.get();
}
 
long * Element::getConnect()
{
    return m_connect.get();
}
 
int Element::getConnectSize() const
{
    switch (m_type) {
 
    case (ElementType::POLYGON):
        return 1 + getVertexCount();
 
    case (ElementType::POLYHEDRON):
        return getFaceStreamSize();
 
    default:
        assert(m_type != ElementType::UNDEFINED);
 
        return getVertexCount();
 
    }
}
 
bool Element::hasSameConnect(const Element &other) const
{
    int cellConnectSize = getConnectSize();
    if (other.getConnectSize() != cellConnectSize) {
        return false;
    }
 
    const long *cellConnect  = getConnect();
    const long *otherConnect = other.getConnect();
    for (int k = 0; k < cellConnectSize; ++k) {
        if (cellConnect[k] != otherConnect[k]) {
            return false;
        }
    }
 
    return true;
}
 
int Element::getFaceCount() const
{
    switch (m_type) {
 
    case (ElementType::POLYGON):
        return countPolygonFaces(getConnect());
 
    case (ElementType::POLYHEDRON):
        return countPolyhedronFaces(getConnect());
 
    default:
        assert(m_type != ElementType::UNDEFINED);
 
        return getInfo().nFaces;
 
    }
}
 
ElementType Element::getFaceType(int face) const
{
    switch (m_type) {
 
    case (ElementType::POLYGON):
    {
        return ElementType::LINE;
    }
 
    case (ElementType::POLYHEDRON):
    {
        int nFaceVertices = getFaceVertexCount(face);
        switch (nFaceVertices) {
 
        case 3:
            return ElementType::TRIANGLE;
 
        case 4:
            return ElementType::QUAD;
 
        default:
            return ElementType::POLYGON;
 
        }
    }
 
    default:
    {
        assert(m_type != ElementType::UNDEFINED);
 
        return getInfo().faceTypeStorage[face];
    }
 
    }
}
 
int Element::getFaceVertexCount(int face) const
{
    switch (m_type) {
 
    case (ElementType::POLYHEDRON):
    {
        const long *connectivity = getConnect();
        int facePos = getFaceStreamPosition(connectivity, face);
 
        return connectivity[facePos];
    }
 
    default:
    {
        assert(m_type != ElementType::UNDEFINED);
 
        ElementType faceType = getFaceType(face);
 
        return ReferenceElementInfo::getInfo(faceType).nVertices;
    }
 
    }
}
 
ConstProxyVector<int> Element::getFaceLocalConnect(int face) const
{
    switch (m_type) {
 
    case (ElementType::POLYGON):
    {
        int nVertices = getVertexCount();
 
        int faceConnectSize = getFaceVertexCount(face);
 
        ConstProxyVector<int> localFaceConnect(ConstProxyVector<int>::INTERNAL_STORAGE, faceConnectSize);
        ConstProxyVector<int>::storage_pointer localFaceConnectStorage = localFaceConnect.storedData();
        for (int i = 0; i < faceConnectSize; ++i) {
            localFaceConnectStorage[i] = (face + i) % nVertices;
        }
 
        return localFaceConnect;
    }
 
    case (ElementType::POLYHEDRON):
    {
        // Get face information
        ElementType faceType = getFaceType(face);
        bool faceHasReferenceInfo = ReferenceElementInfo::hasInfo(faceType);
 
        // Get face vertices
        ConstProxyVector<long> faceVertexIds = getFaceVertexIds(face);
        int nFaceVertices = faceVertexIds.size();
 
        // Get element vertices
        ConstProxyVector<long> vertexIds = getVertexIds();
 
        // Build list of local face verties
        int faceConnectSize  = nFaceVertices;
        int localVertexOffset = 0;
        if (!faceHasReferenceInfo) {
            ++faceConnectSize;
            ++localVertexOffset;
        }
 
        ConstProxyVector<int> localFaceConnect(ConstProxyVector<int>::INTERNAL_STORAGE, faceConnectSize);
        ConstProxyVector<int>::storage_pointer localFaceConnectStorage = localFaceConnect.storedData();
        if (!faceHasReferenceInfo) {
            localFaceConnectStorage[0] = nFaceVertices;
        }
 
        for (int k = 0; k < nFaceVertices; ++k) {
            int vertexId = faceVertexIds[k];
            auto localVertexIdItr = std::find(vertexIds.begin(), vertexIds.end(), vertexId);
            assert(localVertexIdItr != vertexIds.end() && *localVertexIdItr == vertexId);
            localFaceConnectStorage[localVertexOffset + k] = static_cast<int>(std::distance(vertexIds.begin(), localVertexIdItr));
        }
 
        return localFaceConnect;
    }
 
    default:
    {
        assert(m_type != ElementType::UNDEFINED);
 
        const int localConnectSize = getFaceVertexCount(face);
        const int *localFaceConnect = getInfo().faceConnectStorage[face].data();
 
        return ConstProxyVector<int>(localFaceConnect, localConnectSize);
    }
 
    }
}
 
ConstProxyVector<long> Element::getFaceConnect(int face) const
{
    const long *connectivity = getConnect();
 
    switch (m_type) {
 
    case (ElementType::POLYGON):
    {
        int connectSize = getConnectSize();
 
        int facePos         = 1 + face;
        int faceConnectSize = getFaceVertexCount(face);
        if (facePos + faceConnectSize <= connectSize) {
            return ConstProxyVector<long>(connectivity + facePos, faceConnectSize);
        } else {
            ConstProxyVector<long> faceConnect(ConstProxyVector<long>::INTERNAL_STORAGE, faceConnectSize);
            ConstProxyVector<long>::storage_pointer faceConnectStorage = faceConnect.storedData();
            for (int i = 0; i < faceConnectSize; ++i) {
                int position = facePos + i;
                if (position >= connectSize) {
                    position = position % connectSize + 1;
                }
 
                faceConnectStorage[i] = connectivity[position];
            }
 
            return faceConnect;
        }
    }
 
    case (ElementType::POLYHEDRON):
    {
        ElementType faceType = getFaceType(face);
 
        int facePos          = getFaceStreamPosition(connectivity, face);
        int faceConnectSize  = connectivity[facePos];
        int faceConnectBegin = facePos + 1;
        if (!ReferenceElementInfo::hasInfo(faceType)) {
            faceConnectSize++;
            faceConnectBegin--;
        }
 
        return ConstProxyVector<long>(connectivity + faceConnectBegin, faceConnectSize);
    }
 
    default:
    {
        assert(m_type != ElementType::UNDEFINED);
 
        // If we are here, the element has a reference element, therefore we
        // can retrieve the local face connectivity directly form the info
        // associated the to referenc element.
        int faceConnectSize = getFaceVertexCount(face);
        const int *localFaceConnect = getInfo().faceConnectStorage[face].data();
 
        ConstProxyVector<long> faceConnect(ConstProxyVector<long>::INTERNAL_STORAGE, faceConnectSize);
        ConstProxyVector<long>::storage_pointer faceConnectStorage = faceConnect.storedData();
        for (int k = 0; k < faceConnectSize; ++k) {
            int localVertexId = localFaceConnect[k];
            long vertexId = connectivity[localVertexId];
            faceConnectStorage[k] = vertexId;
        }
 
        return faceConnect;
    }
 
    }
}
 
int Element::getEdgeCount() const
{
    switch (m_type) {
 
    case (ElementType::POLYGON):
    {
        return getVertexCount();
    }
 
    case (ElementType::POLYHEDRON):
    {
        int nVertices = getVertexCount();
        int nFaces    = getFaceCount();
        int nEdges    = nVertices + nFaces - 2;
 
        return nEdges;
    }
 
    default:
    {
        assert(m_type != ElementType::UNDEFINED);
 
        return getInfo().nEdges;
    }
 
    }
}
 
ElementType Element::getEdgeType(int edge) const
{
    BITPIT_UNUSED(edge);
 
    int dimension = getDimension();
    switch (dimension) {
 
    case 0:
        return ElementType::UNDEFINED;
 
    case 1:
    case 2:
        return ElementType::VERTEX;
 
    default:
        return ElementType::LINE;
 
    }
}
 
int Element::getEdgeVertexCount(int edge) const
{
    ElementType edgeType = getEdgeType(edge);
 
    return ReferenceElementInfo::getInfo(edgeType).nVertices;
}
 
ConstProxyVector<int> Element::getEdgeLocalConnect(int edge) const
{
    switch (m_type) {
 
    case (ElementType::POLYGON):
    {
        int nEdgeVertices = ReferenceElementInfo::getInfo(ElementType::VERTEX).nVertices;
 
        ConstProxyVector<int> localEdgeConnect(ConstProxyVector<int>::INTERNAL_STORAGE, nEdgeVertices);
        ConstProxyVector<int>::storage_pointer localEdgeConnectStorage = localEdgeConnect.storedData();
        for (int k = 0; k < nEdgeVertices; ++k) {
            localEdgeConnectStorage[k] = k;
        }
 
        return localEdgeConnect;
    }
 
    case (ElementType::POLYHEDRON):
    {
        // Get edge vertices
        ConstProxyVector<long> edgeVertexIds = getEdgeVertexIds(edge);
        int nEdgeVertices = edgeVertexIds.size();
 
        // Get element vertices
        ConstProxyVector<long> vertexIds = getVertexIds();
 
        // Build local edge connectivity
        ConstProxyVector<int> localEdgeConnect(ConstProxyVector<int>::INTERNAL_STORAGE, nEdgeVertices);
        ConstProxyVector<int>::storage_pointer localEdgeConnectStorage = localEdgeConnect.storedData();
        for (int k = 0; k < nEdgeVertices; ++k) {
            int vertexId = edgeVertexIds[k];
            auto localVertexIdItr = std::find(vertexIds.begin(), vertexIds.end(), vertexId);
            assert(localVertexIdItr != vertexIds.end() && *localVertexIdItr == vertexId);
            localEdgeConnectStorage[k] = static_cast<int>(std::distance(vertexIds.begin(), localVertexIdItr));
        }
 
        return localEdgeConnect;
    }
 
    default:
    {
        assert(m_type != ElementType::UNDEFINED);
 
        const int localConnectSize = getEdgeVertexCount(edge);
        const int *localEdgeConnect = getInfo().edgeConnectStorage[edge].data();
 
        return ConstProxyVector<int>(localEdgeConnect, localConnectSize);
    }
 
    }
}
 
ConstProxyVector<long> Element::getEdgeConnect(int edge) const
{
    const long *connectivity = getConnect();
 
    switch (m_type) {
 
    case (ElementType::POLYGON):
    {
        return ConstProxyVector<long>(connectivity + 1 + edge, 1);
    }
 
    case (ElementType::POLYHEDRON):
    {
        std::vector<ConstProxyVector<long>> edgeConnectStorage = evalEdgeConnects(edge + 1);
 
        return edgeConnectStorage[edge];
    }
 
    default:
    {
        assert(m_type != ElementType::UNDEFINED);
 
        ConstProxyVector<int> localEdgeConnect = getEdgeLocalConnect(edge);
        int nEdgeVertices = localEdgeConnect.size();
 
        ConstProxyVector<long> edgeConnect(ConstProxyVector<long>::INTERNAL_STORAGE, nEdgeVertices);
        ConstProxyVector<long>::storage_pointer edgeConnectStorage = edgeConnect.storedData();
        for (int k = 0; k < nEdgeVertices; ++k) {
            int localVertexId = localEdgeConnect[k];
            edgeConnectStorage[k] = connectivity[localVertexId];
        }
 
        return edgeConnect;
    }
 
    }
}
 
int Element::getDimension(ElementType type)
{
    switch (type) {
 
    case ElementType::POLYGON:
        return 2;
 
    case ElementType::POLYHEDRON:
        return 3;
 
    default:
        assert(type != ElementType::UNDEFINED);
 
        return ReferenceElementInfo::getInfo(type).dimension;
 
    }
}
 
int Element::getDimension() const
{
    return getDimension(m_type);
}
 
bool Element::isThreeDimensional(ElementType type)
{
    return (getDimension(type) == 3);
}
 
bool Element::isThreeDimensional() const
{
    return (getDimension() == 3);
}
 
int Element::getVertexCount() const
{
    switch (m_type) {
 
    case (ElementType::POLYGON):
    {
        return countPolygonVertices(getConnect());
    }
 
    case (ElementType::POLYHEDRON):
    {
        return getVertexIds().size();
    }
 
    default:
    {
        assert(m_type != ElementType::UNDEFINED);
 
        return getInfo().nVertices;
    }
 
    }
}
 
ConstProxyVector<long> Element::getVertexIds() const
{
    return getVertexIds(m_type, getConnect());
}
 
ConstProxyVector<long> Element::getVertexIds(ElementType type, const long *connectivity)
{
    switch (type) {
 
    case (ElementType::POLYGON):
    {
        return ConstProxyVector<long>(connectivity + 1, countPolygonVertices(connectivity));
    }
 
    case (ElementType::POLYHEDRON):
    {
        int nFaces = countPolyhedronFaces(connectivity);
 
        // Identify unique vertices
        //
        // We need to keep track of the order in which unique vertices are
        // identified (i.e., the order in which the vertices first appear in
        // the face stream), because the list of vertex ids should be filled
        // using the same order. This guarantees that vertices will be sorted
        // in the same order regardless of the id of the vertices. In this way
        // given an element, the list of its vertices generated by different
        // processes will iterate the vertices in the same order.
        std::unordered_map<long, std::size_t> uniqueVertexIds;
        for (int i = 0; i < nFaces; ++i) {
            int facePos = getFaceStreamPosition(connectivity, i);
 
            int beginVertexPos = facePos + 1;
            int endVertexPos   = facePos + 1 + connectivity[facePos];
            for (int vertexPos = beginVertexPos; vertexPos < endVertexPos; ++vertexPos) {
                long vertexId = connectivity[vertexPos];
                if (uniqueVertexIds.count(vertexId) != 0) {
                    continue;
                }
 
                std::size_t vertexSortIndex = uniqueVertexIds.size();
                uniqueVertexIds.insert({vertexId, vertexSortIndex});
            }
        }
 
        // Fill list of element's vertices
        //
        // The list of unique vertices should contain the vertices sorted in
        // the same order they first appear in the face stream.
        std::size_t nVertices = uniqueVertexIds.size();
 
        ConstProxyVector<long> vertexIds(ConstProxyVector<long>::INTERNAL_STORAGE, nVertices);
        ConstProxyVector<long>::storage_pointer vertexIdsStorage = vertexIds.storedData();
        for (const auto &vertexEntry : uniqueVertexIds) {
            long vertexId = vertexEntry.first;
            std::size_t vertexSortIndex = vertexEntry.second;
            vertexIdsStorage[vertexSortIndex] = vertexId;
        }
 
        return vertexIds;
    }
 
    default:
    {
        assert(type != ElementType::UNDEFINED);
 
        return ConstProxyVector<long>(connectivity, ReferenceElementInfo::getInfo(type).nVertices);
    }
 
    }
}
 
long Element::getVertexId(int vertex) const
{
    switch (m_type) {
 
    case (ElementType::POLYGON):
    case (ElementType::POLYHEDRON):
    {
        ConstProxyVector<long> vertexIds = getVertexIds();
 
        return vertexIds[vertex];
    }
 
    default:
    {
        assert(m_type != ElementType::UNDEFINED);
 
        const long *connectivity = getConnect();
 
        return connectivity[vertex];
    }
 
    }
}
 
int Element::findVertex(long vertexId) const
{
    ConstProxyVector<long> cellVertexIds = getVertexIds();
 
    auto localVertexItr = std::find(cellVertexIds.begin(), cellVertexIds.end(), vertexId);
    if (localVertexItr == cellVertexIds.end()) {
        return -1;
    }
 
    return static_cast<int>(std::distance(cellVertexIds.begin(), localVertexItr));
}
 
ConstProxyVector<long> Element::getFaceVertexIds(int face) const
{
    ConstProxyVector<long> vertexIds = getFaceConnect(face);
    if  (m_type == ElementType::POLYHEDRON) {
        ElementType faceType = getFaceType(face);
        if (faceType == ElementType::POLYGON) {
            assert(!vertexIds.storedData());
            vertexIds.set(vertexIds.data() + 1, vertexIds.size() - 1);
        }
    }
 
    return vertexIds;
}
 
long Element::getFaceVertexId(int face, int vertex) const
{
    switch (m_type) {
 
    case (ElementType::POLYGON):
    case (ElementType::POLYHEDRON):
    {
        ConstProxyVector<long> faceVertexIds = getFaceVertexIds(face);
 
        return faceVertexIds[vertex];
    }
 
    default:
    {
        ConstProxyVector<long> cellVertexIds = getVertexIds();
        ConstProxyVector<int> faceLocalVertexIds = getFaceLocalVertexIds(face);
 
        return cellVertexIds[faceLocalVertexIds[vertex]];
    }
 
    }
}
 
ConstProxyVector<int> Element::getFaceLocalVertexIds(int face) const
{
    switch (m_type) {
 
    case (ElementType::POLYHEDRON):
    {
        ElementType faceType = getFaceType(face);
        if (faceType != ElementType::POLYGON) {
            return getFaceLocalConnect(face);
        }
 
        ConstProxyVector<int> faceLocalConnect = getFaceLocalConnect(face);
        std::size_t faceLocalConnectSize = faceLocalConnect.size();
 
        std::size_t nFaceVertices = faceLocalConnectSize - 1;
        ConstProxyVector<int> faceLocalVertexIds(ConstProxyVector<int>::INTERNAL_STORAGE, nFaceVertices);
        ConstProxyVector<int>::storage_pointer faceLocalVertexIdsStorage = faceLocalVertexIds.storedData();
        for (std::size_t i = 0; i < nFaceVertices; ++i) {
            faceLocalVertexIdsStorage[i] = faceLocalConnect[i + 1];
        }
 
        return faceLocalVertexIds;
    }
 
    default:
    {
        return getFaceLocalConnect(face);
    }
 
    }
}
 
ConstProxyVector<long> Element::getEdgeVertexIds(int edge) const
{
    return getEdgeConnect(edge);
}
 
long Element::getEdgeVertexId(int edge, int vertex) const
{
    ConstProxyVector<long> edgeVertexIds = getEdgeVertexIds(edge);
 
    return edgeVertexIds[vertex];
}
 
ConstProxyVector<int> Element::getEdgeLocalVertexIds(int edge) const
{
    return getEdgeLocalConnect(edge);
}
 
int Element::renumberVertices(const std::unordered_map<long, long> &map)
{
    int nRenumberedVertices = 0;
    switch (m_type) {
 
    case ElementType::POLYGON:
    {
        int nVertices = getVertexCount();
        long *connectivity = getConnect();
        for (int k = 1; k < nVertices + 1; ++k) {
            auto mapItr = map.find(connectivity[k]);
            if (mapItr != map.end()) {
                connectivity[k] = mapItr->second;
                ++nRenumberedVertices;
            }
        }
 
        break;
    }
 
    case ElementType::POLYHEDRON:
    {
        int nFaces = getFaceCount();
        long *connectivity = getConnect();
 
        for (int i = 0; i < nFaces; ++i) {
            int facePos = getFaceStreamPosition(connectivity, i);
 
            int beginVertexPos = facePos + 1;
            int endVertexPos   = facePos + 1 + connectivity[facePos];
            for (int k = beginVertexPos; k < endVertexPos; ++k) {
                auto mapItr = map.find(connectivity[k]);
                if (mapItr != map.end()) {
                    connectivity[k] = mapItr->second;
                    ++nRenumberedVertices;
                }
            }
        }
 
        break;
    }
 
    default:
    {
        assert(m_type != ElementType::UNDEFINED);
 
        int nVertices = getVertexCount();
        long *connectivity = getConnect();
        for (int k = 0; k < nVertices; ++k) {
            auto mapItr = map.find(connectivity[k]);
            if (mapItr != map.end()) {
                connectivity[k] = mapItr->second;
                ++nRenumberedVertices;
            }
        }
 
        break;
    }
 
    }
 
    return nRenumberedVertices;
}
 
std::array<double, 3> Element::evalCentroid(const std::array<double, 3> *coordinates) const
{
    int nVertices = getVertexCount();
    if (nVertices == 0) {
        return {{0., 0., 0.}};
    }
 
    std::array<double, 3> centroid = coordinates[0];
    for (int i = 1; i < nVertices; ++i) {
        const std::array<double, 3> &vertexCoordinates = coordinates[i];
        for (int k = 0; k < 3; ++k) {
            centroid[k] += vertexCoordinates[k];
        }
    }
 
    for (int k = 0; k < 3; ++k) {
        centroid[k] /= nVertices;
    }
 
    return centroid;
}
 
double Element::evalSize(const std::array<double, 3> *coordinates) const
{
    switch (m_type) {
 
    case ElementType::POLYHEDRON:
    {
        // The characteristics size of a polyhedron is evaluated as the
        // volume divided by the area of the largest face.
        int nFaces = getFaceCount();
        BITPIT_CREATE_WORKSPACE(faceVertexCoords, std::array<double BITPIT_COMMA 3>, getVertexCount(), ReferenceElementInfo::MAX_ELEM_VERTICES);
 
        double volume = evalVolume(coordinates);
 
        double faceMaxArea = 0;
        for (int face = 0; face < nFaces; ++face) {
            ElementType faceType = getFaceType(face);
            bool faceHasReferenceInfo = ReferenceElementInfo::hasInfo(faceType);
 
            ConstProxyVector<int> faceLocalVertexIds = getFaceLocalVertexIds(face);
            for (int n = 0; n < getFaceVertexCount(face); ++n) {
                faceVertexCoords[n] = coordinates[faceLocalVertexIds[n]];
            }
 
            double faceArea;
            if (faceHasReferenceInfo) {
                const Reference2DElementInfo &faceInfo = static_cast<const Reference2DElementInfo &>(ReferenceElementInfo::getInfo(faceType));
                faceArea = faceInfo.evalArea(faceVertexCoords);
            } else {
                ConstProxyVector<long> faceConnect = getFaceConnect(face);
                std::size_t faceConnectSize = faceConnect.size();
 
                Element faceElement(0, faceType, faceConnect.size());
                long *faceElementConnect = faceElement.getConnect();
                for (std::size_t i = 0; i < faceConnectSize; ++i) {
                    faceElementConnect[i] = faceConnect[i];
                }
 
                faceArea = faceElement.evalArea(faceVertexCoords);
            }
            faceMaxArea = std::max(faceArea, faceMaxArea);
        }
 
        double size = volume / faceMaxArea;
 
        return size;
    }
 
    case ElementType::POLYGON:
    {
        // The characteristics size of a polygon is evaluated as the
        // area divided by the length of the longest face.
        int nFaces = getFaceCount();
        BITPIT_CREATE_WORKSPACE(faceVertexCoords, std::array<double BITPIT_COMMA 3>, getVertexCount(), ReferenceElementInfo::MAX_ELEM_VERTICES);
 
        double area = evalArea(coordinates);
 
        double faceMaxLength = 0;
        for (int face = 0; face < nFaces; ++face) {
            ElementType faceType = getFaceType(face);
            assert(ReferenceElementInfo::hasInfo(faceType));
 
            ConstProxyVector<int> faceLocalVertexIds = getFaceLocalVertexIds(face);
            for (int n = 0; n < getFaceVertexCount(face); ++n) {
                faceVertexCoords[n] = coordinates[faceLocalVertexIds[n]];
            }
 
            const Reference1DElementInfo &faceInfo = static_cast<const Reference1DElementInfo &>(ReferenceElementInfo::getInfo(faceType));
            faceMaxLength = std::max(faceInfo.evalLength(faceVertexCoords), faceMaxLength);
        }
 
        double size = area / faceMaxLength;
 
        return size;
    }
 
    case ElementType::UNDEFINED:
    {
        return 0.;
    }
 
    default:
    {
        const ReferenceElementInfo &referenceInfo = static_cast<const ReferenceElementInfo &>(getInfo());
 
        return referenceInfo.evalSize(coordinates);
    }
 
    }
}
 
 
double Element::evalVolume(const std::array<double, 3> *coordinates) const
{
    switch (m_type) {
 
    case ElementType::POLYHEDRON:
    {
        Tesselation tesselation = generateTesselation(coordinates);
        int nTiles = tesselation.getTileCount();
 
        double volume = 0.;
        for (int i = 0; i < nTiles; ++i) {
            ElementType tileType = tesselation.getTileType(i);
            const std::vector<std::array<double, 3>> tileCoordinates = tesselation.getTileVertexCoordinates(i);
            const Reference3DElementInfo &referenceInfo = static_cast<const Reference3DElementInfo &>(ReferenceElementInfo::getInfo(tileType));
 
            volume += referenceInfo.evalVolume(tileCoordinates.data());
        }
 
        return volume;
    }
 
    default:
    {
        assert(getDimension() == 3);
 
        const Reference3DElementInfo &referenceInfo = static_cast<const Reference3DElementInfo &>(getInfo());
 
        return referenceInfo.evalVolume(coordinates);
    }
 
    }
}
 
double Element::evalArea(const std::array<double, 3> *coordinates) const
{
    switch (m_type) {
 
    case ElementType::POLYGON:
    {
        Tesselation tesselation = generateTesselation(coordinates);
        int nTiles = tesselation.getTileCount();
 
        double area = 0.;
        for (int i = 0; i < nTiles; ++i) {
            ElementType tileType = tesselation.getTileType(i);
            const std::vector<std::array<double, 3>> tileCoordinates = tesselation.getTileVertexCoordinates(i);
            const Reference2DElementInfo &referenceInfo = static_cast<const Reference2DElementInfo &>(ReferenceElementInfo::getInfo(tileType));
 
            area += referenceInfo.evalArea(tileCoordinates.data());
        }
 
        return area;
    }
 
    default:
    {
        assert(getDimension() == 2);
 
        const Reference2DElementInfo &referenceInfo = static_cast<const Reference2DElementInfo &>(getInfo());
 
        return referenceInfo.evalArea(coordinates);
    }
 
    }
}
 
double Element::evalLength(const std::array<double, 3> *coordinates) const
{
    switch (m_type) {
 
    case ElementType::POLYGON:
    case ElementType::POLYHEDRON:
    case ElementType::UNDEFINED:
    {
        return 0.;
    }
 
    default:
    {
        assert(getDimension() == 1);
 
        const Reference1DElementInfo &referenceInfo = static_cast<const Reference1DElementInfo &>(getInfo());
 
        return referenceInfo.evalLength(coordinates);
    }
 
    }
}
 
std::array<double, 3> Element::evalNormal(const std::array<double, 3> *coordinates,
                                          const std::array<double, 3> &orientation,
                                          const std::array<double, 3> &point) const
{
    switch (m_type) {
 
    case ElementType::POLYGON:
    {
        // The normal of a polygonal element is evaluated as the weighted average of the
        // normals of its tiles. The resulting vector should be normalized in order to
        // obtain a versor (because of triangular inequality it is not possible
        // to evaluate the weights in order to automatically obtain a versor).
        int dimension = getDimension();
 
        Tesselation tesselation = generateTesselation(coordinates);
        int nTiles = tesselation.getTileCount();
 
        std::array<double, 3> normal = {{0., 0., 0.}};
        for (int i = 0; i < nTiles; ++i) {
            ElementType tileType = tesselation.getTileType(i);
            const std::vector<std::array<double, 3>> tileCoordinates = tesselation.getTileVertexCoordinates(i);
            const Reference2DElementInfo &referenceInfo = static_cast<const Reference2DElementInfo &>(ReferenceElementInfo::getInfo(tileType));
 
            double tileArea = referenceInfo.evalArea(tileCoordinates.data());
            std::array<double, 3> tileNormal = {{0., 0., 0.}};
            if (dimension == 2) {
                tileNormal = referenceInfo.evalNormal(tileCoordinates.data(), point);
            } else if (dimension == 1) {
                const Reference1DElementInfo &referenceInfo1D = static_cast<const Reference1DElementInfo &>(ReferenceElementInfo::getInfo(tileType));
                tileNormal = referenceInfo1D.evalNormal(tileCoordinates.data(), orientation, point);
            } else if (dimension == 0) {
                tileNormal = orientation;
            }
 
            normal += tileArea * tileNormal;
        }
        normal /= norm2(normal);
 
        return normal;
    }
 
    default:
    {
        assert(getDimension() != 3);
 
        int dimension = getDimension();
        if (dimension == 2) {
            const Reference2DElementInfo &referenceInfo = static_cast<const Reference2DElementInfo &>(getInfo());
 
            return referenceInfo.evalNormal(coordinates, point);
        } else if (dimension == 1) {
            const Reference1DElementInfo &referenceInfo = static_cast<const Reference1DElementInfo &>(getInfo());
 
            return referenceInfo.evalNormal(coordinates, orientation, point);
        } else {
            return orientation;
        }
    }
 
    }
}
 
double Element::evalPointDistance(const std::array<double, 3> &point, const std::array<double, 3> *coordinates) const
{
    double distance;
    std::array<double, 3> projection;
    evalPointProjection(point, coordinates, &projection, &distance);
 
    return distance;
}
 
void Element::evalPointProjection(const std::array<double, 3> &point, const std::array<double, 3> *coordinates,
                                  std::array<double, 3> *projection, double *distance) const
{
    switch (m_type) {
 
    case ElementType::POLYGON:
    {
        int projectionFlag;
        *distance = CGElem::distancePointPolygon(point, getVertexCount(), coordinates, *projection, projectionFlag);
 
        break;
    }
 
    case ElementType::POLYHEDRON:
    {
        *distance = std::numeric_limits<double>::max();
 
        int nFaces = getFaceCount();
        std::vector<std::array<double, 3>> faceCoordinates;
        for (int i = 0; i < nFaces; ++i) {
            ElementType faceType = getFaceType(i);
            ConstProxyVector<int> faceVertexIds = getFaceLocalVertexIds(i);
            int nFaceVertices = faceVertexIds.size();
            faceCoordinates.resize(nFaceVertices);
            for (int k = 0; k < nFaceVertices; ++k) {
                faceCoordinates[k] = coordinates[faceVertexIds[k]];
            }
 
            double faceDistance;
            std::array<double, 3> faceProjection;
            bool faceHasReferenceInfo = ReferenceElementInfo::hasInfo(faceType);
            if (faceHasReferenceInfo) {
                ReferenceElementInfo::getInfo(faceType).evalPointProjection(point, faceCoordinates.data(), &faceProjection, &faceDistance);
            } else {
                int faceProjectionFlag;
                faceDistance = CGElem::distancePointPolygon(point, nFaceVertices, faceCoordinates.data(), faceProjection, faceProjectionFlag);
            }
 
            if (faceDistance < *distance) {
                *distance   = faceDistance;
                *projection = faceProjection;
            }
        }
 
        break;
    }
 
    default:
    {
        assert(ReferenceElementInfo::hasInfo(m_type));
 
        getInfo().evalPointProjection(point, coordinates, projection, distance);
 
        break;
    }
 
    }
}
 
Element::Tesselation Element::generateTesselation(const std::array<double, 3> *coordinates) const
{
    Tesselation tesselation;
 
    // Add the coordinates of the vertices to the tesselation
    int nVertices = getVertexCount();
    std::vector<int> vertexTesselationIds = tesselation.importVertexCoordinates(coordinates, nVertices);
 
    // Generate the tesselation
    ElementType type = getType();
    switch(type) {
 
    case ElementType::POLYGON:
    {
        tesselation.importPolygon(vertexTesselationIds);
 
        break;
    }
 
    case ElementType::POLYHEDRON:
    {
        int nFaces = getFaceCount();
        std::vector<std::vector<int>> faceTesselationIds(nFaces);
        for (int i = 0; i < nFaces; ++i) {
            ConstProxyVector<int> localVertexIds = getFaceLocalVertexIds(i);
            int nFaceVertices = localVertexIds.size();
            faceTesselationIds[i].resize(nFaceVertices);
            for (int k = 0; k < nFaceVertices; ++k) {
                faceTesselationIds[i][k] = vertexTesselationIds[localVertexIds[k]];
            }
        }
 
        tesselation.importPolyhedron(vertexTesselationIds, faceTesselationIds);
 
        break;
    }
 
    default:
    {
        assert(ReferenceElementInfo::hasInfo(type));
 
        tesselation.m_nTiles = 1;
        tesselation.m_types.push_back(type);
        tesselation.m_connects.push_back(vertexTesselationIds);
 
        break;
    }
 
    }
 
    return tesselation;
}
 
int Element::getFaceStreamSize() const
{
    int nFaces = getFaceCount();
    int size   = 1 + (getFaceStreamPosition(nFaces) - 1);
 
    return size;
}
 
std::vector<long> Element::getFaceStream() const
{
    int nFaces = getFaceCount();
    int faceStreamSize = getFaceStreamSize();
    std::vector<long> faceStream(faceStreamSize);
 
    int pos = 0;
    faceStream[pos] = nFaces;
    for (int i = 0; i < nFaces; ++i) {
        ConstProxyVector<long> faceVertexIds = getFaceVertexIds(i);
        int nFaceVertices = faceVertexIds.size();
 
        ++pos;
        faceStream[pos] = getFaceVertexCount(i);
        for (int k = 0; k < nFaceVertices; ++k) {
            ++pos;
            faceStream[pos] = faceVertexIds[k];
        }
    }
 
    return faceStream;
}
 
void Element::renumberFaceStream(const PiercedStorage<long, long> &map, std::vector<long> *faceStream)
{
    int pos = 0;
    int nFaces = (*faceStream)[pos];
    for (int i = 0; i < nFaces; ++i) {
        ++pos;
        int nFaceVertices = (*faceStream)[pos];
        for (int k = 0; k < nFaceVertices; ++k) {
            ++pos;
            auto mapItr = map.find((*faceStream)[pos]);
            if (mapItr != map.end()) {
                (*faceStream)[pos] = *mapItr;
            }
        }
    }
}
 
int Element::getFaceStreamPosition(int face) const
{
    return getFaceStreamPosition(getConnect(), face);
}
 
int Element::getFaceStreamPosition(const long *connectivity,  int face)
{
    int position = 1;
    for (int i = 0; i < face; ++i) {
        position += 1 + connectivity[position];
    }
 
    return position;
}
 
int Element::countPolygonVertices(const long *connectivity)
{
    return connectivity[0];
}
 
int Element::countPolygonFaces(const long *connectivity)
{
    return countPolygonVertices(connectivity);
}
 
int Element::countPolyhedronFaces(const long *connectivity)
{
    return connectivity[0];
}
 
std::vector<ConstProxyVector<long>> Element::evalEdgeConnects(int nRequestedEdges) const
{
    if (nRequestedEdges == -1) {
        nRequestedEdges = getEdgeCount();
    }
    assert(nRequestedEdges <= getEdgeCount());
 
    std::set<std::pair<long, long>> edgeSet;
    std::vector<ConstProxyVector<long>> edgeConnectStorage(nRequestedEdges, ConstProxyVector<long>(ConstProxyVector<long>::INTERNAL_STORAGE, 2));
 
    int nFaces = getFaceCount();
    for (int i = 0; i < nFaces; ++i) {
        ConstProxyVector<long> faceVertexIds = getFaceVertexIds(i);
        int nFaceVertices = faceVertexIds.size();
        for (int k = 0; k < nFaceVertices; ++k) {
            long vertex_A = faceVertexIds[k];
            long vertex_B = faceVertexIds[(k + 1) % nFaceVertices];
            if (vertex_A > vertex_B) {
                std::swap(vertex_A, vertex_B);
            }
 
            std::pair<long, long> edgePair = std::pair<long, long>(vertex_A, vertex_B);
            std::pair<std::set<std::pair<long, long>>::iterator, bool> insertResult = edgeSet.insert(edgePair);
            if (insertResult.second) {
                ConstProxyVector<long>::storage_pointer connectStorage = edgeConnectStorage[edgeSet.size() - 1].storedData();
                connectStorage[0] = edgePair.first;
                connectStorage[1] = edgePair.second;
 
                if (edgeSet.size() == (std::size_t) nRequestedEdges) {
                    return edgeConnectStorage;
                }
            }
        }
    }
 
    assert((int) edgeSet.size() == nRequestedEdges);
 
    return edgeConnectStorage;
}
 
unsigned int Element::getBinarySize() const
{
    unsigned int binarySize = sizeof(m_type) + sizeof(m_id) + getConnectSize() * sizeof(long) + sizeof(m_pid);
    if (!bitpit::ReferenceElementInfo::hasInfo(m_type)) {
        binarySize += sizeof(int);
    }
 
    return binarySize;
}
 
// Explicit instantiation of the Element containers
template class PiercedVector<Element>;
 
}