element_reference.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 <set>
 
#include "bitpit_CG.hpp"
#include "bitpit_containers.hpp"
#include "bitpit_operators.hpp"
 
#include "element_reference.hpp"
 
namespace bitpit {


ReferenceElementInfo::ReferenceElementInfo(int _dimension, ElementType _type, int _nVertices, int _nFaces, int _nEdges)
    : dimension(_dimension), type(_type),
      nVertices(_nVertices), nFaces(_nFaces), nEdges(_nEdges)
{
    faceTypeStorage.fill(ElementType::UNDEFINED);
    for (int i = 0; i < MAX_ELEM_FACES; ++i) {
        faceConnectStorage[i].fill(-1);
        faceEdgeStorage[i].fill(-1);
    }
 
    edgeTypeStorage.fill(ElementType::UNDEFINED);
    for (int i = 0; i < MAX_ELEM_EDGES; ++i) {
        edgeConnectStorage[i].fill(-1);
    }
}

bool ReferenceElementInfo::hasInfo(ElementType type)
{
 
    switch (type) {
 
    case (ElementType::VERTEX):
    case (ElementType::LINE):
    case (ElementType::TRIANGLE):
    case (ElementType::PIXEL):
    case (ElementType::QUAD):
    case (ElementType::TETRA):
    case (ElementType::VOXEL):
    case (ElementType::HEXAHEDRON):
    case (ElementType::PYRAMID):
    case (ElementType::WEDGE):
        return true;
 
    default:
        return false;
 
    }
}

const ReferenceElementInfo & ReferenceElementInfo::getInfo(ElementType type)
{
    switch (type) {
 
    case (ElementType::VERTEX):
        return ReferenceVertexInfo::info;
 
    case (ElementType::LINE):
        return ReferenceLineInfo::info;
 
    case (ElementType::TRIANGLE):
        return ReferenceTriangleInfo::info;
 
    case (ElementType::PIXEL):
        return ReferencePixelInfo::info;
 
    case (ElementType::QUAD):
        return ReferenceQuadInfo::info;
 
    case (ElementType::TETRA):
        return ReferenceTetraInfo::info;
 
    case (ElementType::VOXEL):
        return ReferenceVoxelInfo::info;
 
    case (ElementType::HEXAHEDRON):
        return ReferenceHexahedronInfo::info;
 
    case (ElementType::PYRAMID):
        return ReferencePyramidInfo::info;
 
    case (ElementType::WEDGE):
        return ReferenceWedgeInfo::info;
 
    default:
        BITPIT_UNREACHABLE("Unsupported element");
        throw std::runtime_error("Unsupported element");
 
    }
}

void ReferenceElementInfo::initializeFaceEdges(const std::vector<const ReferenceElementInfo *> &facesInfo,
                                               const std::vector<const ReferenceElementInfo *> &edgesInfo)
{
    for (int k = 0; k < nFaces; ++k) {
        const ReferenceElementInfo &faceInfo = *(facesInfo[k]);
 
        int nFaceEdges = faceInfo.nFaces;
        int faceEdgeCounter = 0;
        for (int i = 0; i < nFaceEdges; ++i) {
            const ReferenceElementInfo &faceEdgeInfo = *(edgesInfo[i]);
 
            // Connectivity of the edge associated to the face
            const int *localFaceEdgeConnect = faceInfo.faceConnectStorage[i].data();
 
            std::set<int> faceEdgeConnect;
            for (int n = 0; n < faceEdgeInfo.nVertices; ++n) {
                int localVertexId = localFaceEdgeConnect[n];
                int vertexId      = faceConnectStorage[k][localVertexId];
 
                faceEdgeConnect.insert(vertexId);
            }
 
            // Search the edge that has the same connectivity of the face edge
            for (int j = 0; j < nEdges; ++j) {
                const ReferenceElementInfo &guessEdgeInfo = *(edgesInfo[j]);
 
                // If face edge and the guess edge have a different type, the
                // two edge cannot be the same.
                if (guessEdgeInfo.type != faceEdgeInfo.type) {
                    continue;
                }
 
                // If the connecitivity of the face edge and the one of the
                // guess edge are the same, the two edges coincides.
                const std::set<int> guessEdgeConnect = std::set<int>(edgeConnectStorage[j].begin(), edgeConnectStorage[j].begin() + guessEdgeInfo.nVertices);
                if (faceEdgeConnect == guessEdgeConnect) {
                    faceEdgeStorage[k][faceEdgeCounter] = j;
                    ++faceEdgeCounter;
                }
            }
        }
 
        assert(faceEdgeCounter == nFaceEdges);
    }
}


Reference3DElementInfo::Reference3DElementInfo(ElementType type, int nVertices, int nFaces)
    : ReferenceElementInfo(3, type, nVertices, nFaces, nVertices + nFaces - 2)
{
}

double Reference3DElementInfo::evalSize(const std::array<double, 3> *vertexCoords) const
{
    double volume = evalVolume(vertexCoords);
 
    double faceMaxArea = 0;
    for (int face = 0; face < nFaces; ++face) {
        double faceArea = evalFaceArea(face, vertexCoords);
        faceMaxArea = std::max(faceArea, faceMaxArea);
    }
 
    double length = volume / faceMaxArea;
 
    return length;
}

double Reference3DElementInfo::evalSurfaceArea(const std::array<double, 3> *vertexCoords) const
{
    double area = 0;
    for (int i = 0; i < nFaces; ++i) {
        area += evalFaceArea(i, vertexCoords);
    }
 
    return area;
}

double Reference3DElementInfo::evalEdgePerimeter(const std::array<double, 3> *vertexCoords) const
{
    double perimeter = 0;
    for (int i = 0; i < nEdges; ++i) {
        perimeter += evalEdgeLength(i, vertexCoords);
    }
 
    return perimeter;
}

double Reference3DElementInfo::evalFaceArea(int face, const std::array<double, 3> *vertexCoords) const
{
    ElementType faceType = faceTypeStorage[face];
    const Reference2DElementInfo &faceInfo = static_cast<const Reference2DElementInfo &>(getInfo(faceType));
 
    std::array<std::array<double, 3>, MAX_ELEM_VERTICES> faceVertexCoords;
    for (int n = 0; n < faceInfo.nVertices; ++n) {
        faceVertexCoords[n] = vertexCoords[faceConnectStorage[face][n]];
    }
 
    return faceInfo.evalArea(faceVertexCoords.data());
}

double Reference3DElementInfo::evalEdgeLength(int edge, const std::array<double, 3> *vertexCoords) const
{
    const ReferenceLineInfo &edgeInfo = static_cast<const ReferenceLineInfo &>(getInfo(ElementType::LINE));
 
    std::array<std::array<double, 3>, MAX_ELEM_VERTICES> edgeVertexCoords;
    for (int n = 0; n < edgeInfo.nVertices; ++n) {
        edgeVertexCoords[n] = vertexCoords[edgeConnectStorage[edge][n]];
    }
 
    return edgeInfo.evalLength(edgeVertexCoords.data());
}

void Reference3DElementInfo::evalPointProjection(const std::array<double, 3> &point, const std::array<double, 3> *vertexCoords,
                                                 std::array<double, 3> *projection, double *distance) const
{
    std::array<std::array<double, 3>, MAX_ELEM_VERTICES> faceVertexCoords;
 
    *distance = std::numeric_limits<double>::max();
    for (int i = 0; i < nFaces; ++i) {
        ElementType faceType = faceTypeStorage[i];
        const Reference2DElementInfo &faceInfo = static_cast<const Reference2DElementInfo &>(getInfo(faceType));
        for (int n = 0; n < faceInfo.nVertices; ++n) {
            faceVertexCoords[n] = vertexCoords[faceConnectStorage[i][n]];
        }
 
        double faceDistance;
        std::array<double, 3> faceProjection;
        faceInfo.evalPointProjection(point, faceVertexCoords.data(), &faceProjection, &faceDistance);
 
        if (faceDistance < *distance) {
            *distance   = faceDistance;
            *projection = faceProjection;
        }
    }
}

double Reference3DElementInfo::evalPointDistance(const std::array<double, 3> &point, const std::array<double, 3> *vertexCoords) const
{
    std::array<std::array<double, 3>, MAX_ELEM_VERTICES> faceVertexCoords;
 
    double distance = std::numeric_limits<double>::max();
    for (int i = 0; i < nFaces; ++i) {
        ElementType faceType = faceTypeStorage[i];
        const Reference2DElementInfo &faceInfo = static_cast<const Reference2DElementInfo &>(getInfo(faceType));
        for (int n = 0; n < faceInfo.nVertices; ++n) {
            faceVertexCoords[n] = vertexCoords[faceConnectStorage[i][n]];
        }
 
        distance = std::min(faceInfo.evalPointDistance(point, faceVertexCoords.data()), distance);
    }
 
    return distance;
}

 
const ReferenceTetraInfo ReferenceTetraInfo::info;

ReferenceTetraInfo::ReferenceTetraInfo()
    : Reference3DElementInfo(ElementType::TETRA, 4, 4)
{
    const ReferenceLineInfo lineInfo;
    const ReferenceTriangleInfo triangleInfo;
 
    // Edge data
    std::vector<const ReferenceElementInfo *> edgesInfo(nEdges);
 
    for (int k = 0; k < nEdges; ++k) {
        edgesInfo[k] = &lineInfo;
 
        edgeTypeStorage[k] = LINE;
    }
 
    edgeConnectStorage[0][0] = 0;
    edgeConnectStorage[0][1] = 1;
 
    edgeConnectStorage[1][0] = 1;
    edgeConnectStorage[1][1] = 2;
 
    edgeConnectStorage[2][0] = 2;
    edgeConnectStorage[2][1] = 0;
 
    edgeConnectStorage[3][0] = 3;
    edgeConnectStorage[3][1] = 0;
 
    edgeConnectStorage[4][0] = 3;
    edgeConnectStorage[4][1] = 1;
 
    edgeConnectStorage[5][0] = 3;
    edgeConnectStorage[5][1] = 2;
 
    // Face data
    std::vector<const ReferenceElementInfo *> facesInfo(nFaces);
 
    for (int k = 0; k < nFaces; ++k) {
        facesInfo[k] = &triangleInfo;
 
        faceTypeStorage[k] = TRIANGLE;
    }
 
    faceConnectStorage[0][0] = 1;
    faceConnectStorage[0][1] = 0;
    faceConnectStorage[0][2] = 2;
 
    faceConnectStorage[1][0] = 0;
    faceConnectStorage[1][1] = 3;
    faceConnectStorage[1][2] = 2;
 
    faceConnectStorage[2][0] = 3;
    faceConnectStorage[2][1] = 1;
    faceConnectStorage[2][2] = 2;
 
    faceConnectStorage[3][0] = 0;
    faceConnectStorage[3][1] = 1;
    faceConnectStorage[3][2] = 3;
 
    initializeFaceEdges(facesInfo, edgesInfo);
}

double ReferenceTetraInfo::evalVolume(const std::array<double, 3> *vertexCoords) const
{
    const std::array<double, 3> &V_A = vertexCoords[0];
    const std::array<double, 3> &V_B = vertexCoords[1];
    const std::array<double, 3> &V_C = vertexCoords[2];
    const std::array<double, 3> &V_D = vertexCoords[3];
 
    double volume = std::abs(dotProduct(V_A - V_D, crossProduct(V_B - V_D, V_C - V_D))) / 6.;
 
    return volume;
}

double ReferenceTetraInfo::evalSize(const std::array<double, 3> *vertexCoords) const
{
    double volume = evalVolume(vertexCoords);
    double area   = evalSurfaceArea(vertexCoords);
 
    double inscribedRadius = 3 * volume / area;
 
    double length = 4 * inscribedRadius * sqrt(3. / 2.);
 
    return length;
}

 
const ReferenceVoxelInfo ReferenceVoxelInfo::info;

ReferenceVoxelInfo::ReferenceVoxelInfo()
    : Reference3DElementInfo(ElementType::VOXEL, 8, 6)
{
    const ReferenceLineInfo lineInfo;
    const ReferencePixelInfo pixelInfo;
 
    // Edge data
    std::vector<const ReferenceElementInfo *> edgesInfo(nEdges);
 
    for (int k = 0; k < nEdges; ++k) {
        edgesInfo[k] = &lineInfo;
 
        edgeTypeStorage[k] = LINE;
    }
 
    edgeConnectStorage[0][0] = 0;
    edgeConnectStorage[0][1] = 2;
 
    edgeConnectStorage[1][0] = 1;
    edgeConnectStorage[1][1] = 3;
 
    edgeConnectStorage[2][0] = 0;
    edgeConnectStorage[2][1] = 1;
 
    edgeConnectStorage[3][0] = 2;
    edgeConnectStorage[3][1] = 3;
 
    edgeConnectStorage[4][0] = 0;
    edgeConnectStorage[4][1] = 4;
 
    edgeConnectStorage[5][0] = 1;
    edgeConnectStorage[5][1] = 5;
 
    edgeConnectStorage[6][0] = 2;
    edgeConnectStorage[6][1] = 6;
 
    edgeConnectStorage[7][0] = 3;
    edgeConnectStorage[7][1] = 7;
 
    edgeConnectStorage[8][0] = 4;
    edgeConnectStorage[8][1] = 6;
 
    edgeConnectStorage[9][0] = 5;
    edgeConnectStorage[9][1] = 7;
 
    edgeConnectStorage[10][0] = 4;
    edgeConnectStorage[10][1] = 5;
 
    edgeConnectStorage[11][0] = 6;
    edgeConnectStorage[11][1] = 7;
 
    // Face data
    std::vector<const ReferenceElementInfo *> facesInfo(nFaces);
 
    for (int k = 0; k < nFaces; ++k) {
        facesInfo[k] = &pixelInfo;
 
        faceTypeStorage[k] = PIXEL;
    }
 
    faceConnectStorage[0][0] = 2;
    faceConnectStorage[0][1] = 0;
    faceConnectStorage[0][2] = 6;
    faceConnectStorage[0][3] = 4;
 
    faceConnectStorage[1][0] = 1;
    faceConnectStorage[1][1] = 3;
    faceConnectStorage[1][2] = 5;
    faceConnectStorage[1][3] = 7;
 
    faceConnectStorage[2][0] = 0;
    faceConnectStorage[2][1] = 1;
    faceConnectStorage[2][2] = 4;
    faceConnectStorage[2][3] = 5;
 
    faceConnectStorage[3][0] = 3;
    faceConnectStorage[3][1] = 2;
    faceConnectStorage[3][2] = 7;
    faceConnectStorage[3][3] = 6;
 
    faceConnectStorage[4][0] = 2;
    faceConnectStorage[4][1] = 3;
    faceConnectStorage[4][2] = 0;
    faceConnectStorage[4][3] = 1;
 
    faceConnectStorage[5][0] = 7;
    faceConnectStorage[5][1] = 6;
    faceConnectStorage[5][2] = 5;
    faceConnectStorage[5][3] = 4;
 
    initializeFaceEdges(facesInfo, edgesInfo);
}

double ReferenceVoxelInfo::evalVolume(const std::array<double, 3> *vertexCoords) const
{
    const std::array<double, 3> &V_A = vertexCoords[0];
    const std::array<double, 3> &V_B = vertexCoords[1];
    const std::array<double, 3> &V_C = vertexCoords[2];
    const std::array<double, 3> &V_D = vertexCoords[4];
 
    double volume = std::abs(dotProduct(V_D - V_A, crossProduct(V_B - V_A, V_C - V_A)));
 
    return volume;
}

double ReferenceVoxelInfo::evalSize(const std::array<double, 3> *vertexCoords) const
{
    const std::array<double, 3> &V_A = vertexCoords[0];
    const std::array<double, 3> &V_B = vertexCoords[1];
    const std::array<double, 3> &V_C = vertexCoords[2];
    const std::array<double, 3> &V_D = vertexCoords[4];
 
    double sideLength_x = norm2(V_B - V_A);
    double sideLength_y = norm2(V_C - V_A);
    double sideLength_z = norm2(V_D - V_A);
 
    double size = std::min({sideLength_x, sideLength_y, sideLength_z});
 
    return size;
}

 
const ReferenceHexahedronInfo ReferenceHexahedronInfo::info;

ReferenceHexahedronInfo::ReferenceHexahedronInfo()
    : Reference3DElementInfo(ElementType::HEXAHEDRON, 8, 6)
{
    const ReferenceLineInfo lineInfo;
    const ReferenceQuadInfo quadInfo;
 
    // Edge data
    std::vector<const ReferenceElementInfo *> edgesInfo(nEdges);
 
    for (int k = 0; k < nEdges; ++k) {
        edgesInfo[k] = &lineInfo;
 
        edgeTypeStorage[k] = LINE;
    }
 
    edgeConnectStorage[0][0] = 1;
    edgeConnectStorage[0][1] = 0;
 
    edgeConnectStorage[1][0] = 1;
    edgeConnectStorage[1][1] = 2;
 
    edgeConnectStorage[2][0] = 2;
    edgeConnectStorage[2][1] = 3;
 
    edgeConnectStorage[3][0] = 3;
    edgeConnectStorage[3][1] = 0;
 
    edgeConnectStorage[4][0] = 4;
    edgeConnectStorage[4][1] = 5;
 
    edgeConnectStorage[5][0] = 5;
    edgeConnectStorage[5][1] = 6;
 
    edgeConnectStorage[6][0] = 6;
    edgeConnectStorage[6][1] = 7;
 
    edgeConnectStorage[7][0] = 7;
    edgeConnectStorage[7][1] = 4;
 
    edgeConnectStorage[8][0] = 0;
    edgeConnectStorage[8][1] = 4;
 
    edgeConnectStorage[9][0] = 1;
    edgeConnectStorage[9][1] = 5;
 
    edgeConnectStorage[10][0] = 2;
    edgeConnectStorage[10][1] = 6;
 
    edgeConnectStorage[11][0] = 3;
    edgeConnectStorage[11][1] = 7;
 
    // Face data
    std::vector<const ReferenceElementInfo *> facesInfo(nFaces);
 
    for (int k = 0; k < nFaces; ++k) {
        facesInfo[k] = &quadInfo;
 
        faceTypeStorage[k] = QUAD;
    }
 
    faceConnectStorage[0][0] = 1;
    faceConnectStorage[0][1] = 0;
    faceConnectStorage[0][2] = 3;
    faceConnectStorage[0][3] = 2;
 
    faceConnectStorage[1][0] = 4;
    faceConnectStorage[1][1] = 5;
    faceConnectStorage[1][2] = 6;
    faceConnectStorage[1][3] = 7;
 
    faceConnectStorage[2][0] = 7;
    faceConnectStorage[2][1] = 3;
    faceConnectStorage[2][2] = 0;
    faceConnectStorage[2][3] = 4;
 
    faceConnectStorage[3][0] = 5;
    faceConnectStorage[3][1] = 1;
    faceConnectStorage[3][2] = 2;
    faceConnectStorage[3][3] = 6;
 
    faceConnectStorage[4][0] = 4;
    faceConnectStorage[4][1] = 0;
    faceConnectStorage[4][2] = 1;
    faceConnectStorage[4][3] = 5;
 
    faceConnectStorage[5][0] = 6;
    faceConnectStorage[5][1] = 2;
    faceConnectStorage[5][2] = 3;
    faceConnectStorage[5][3] = 7;
 
    initializeFaceEdges(facesInfo, edgesInfo);
}

double ReferenceHexahedronInfo::evalVolume(const std::array<double, 3> *vertexCoords) const
{
    std::array<std::array<double, 3>, MAX_ELEM_VERTICES> pyramidVertexCoords;
 
    const ReferenceQuadInfo    &quadInfo    = static_cast<const ReferenceQuadInfo &>(getInfo(ElementType::QUAD));
    const ReferencePyramidInfo &pyramidInfo = static_cast<const ReferencePyramidInfo &>(getInfo(ElementType::PYRAMID));
 
    // The centroid is the apex of all the pyramids
    pyramidVertexCoords[4] = vertexCoords[0];
    for (int i = 1; i < nVertices; ++i) {
        pyramidVertexCoords[4] += vertexCoords[i];
    }
    pyramidVertexCoords[4] = pyramidVertexCoords[4] / double(nVertices);
 
    // Sum the volume of the pyramids having a face as the base and the
    // centroid as the apex.
    double volume = 0.;
    for (int i = 0; i < nFaces; ++i) {
        for (int n = 0; n < quadInfo.nVertices; ++n) {
            pyramidVertexCoords[quadInfo.nVertices - n - 1] = vertexCoords[faceConnectStorage[i][n]];
        }
 
        volume += pyramidInfo.evalVolume(pyramidVertexCoords.data());
    }
 
    return volume;
}

 
const ReferencePyramidInfo ReferencePyramidInfo::info;

ReferencePyramidInfo::ReferencePyramidInfo()
    : Reference3DElementInfo(ElementType::PYRAMID, 5, 5)
{
    const ReferenceLineInfo lineInfo;
    const ReferenceTriangleInfo triangleInfo;
    const ReferenceQuadInfo quadInfo;
 
    // Edge data
    std::vector<const ReferenceElementInfo *> edgesInfo(nEdges);
 
    for (int k = 0; k < nEdges; ++k) {
        edgesInfo[k] = &lineInfo;
 
        edgeTypeStorage[k] = LINE;
    }
 
    edgeConnectStorage[0][0] = 0;
    edgeConnectStorage[0][1] = 1;
 
    edgeConnectStorage[1][0] = 1;
    edgeConnectStorage[1][1] = 2;
 
    edgeConnectStorage[2][0] = 2;
    edgeConnectStorage[2][1] = 3;
 
    edgeConnectStorage[3][0] = 3;
    edgeConnectStorage[3][1] = 0;
 
    edgeConnectStorage[4][0] = 4;
    edgeConnectStorage[4][1] = 0;
 
    edgeConnectStorage[5][0] = 4;
    edgeConnectStorage[5][1] = 1;
 
    edgeConnectStorage[6][0] = 4;
    edgeConnectStorage[6][1] = 2;
 
    edgeConnectStorage[7][0] = 4;
    edgeConnectStorage[7][1] = 3;
 
    // Face data
    std::vector<const ReferenceElementInfo *> facesInfo(nFaces);
 
    for (int k = 0; k < nFaces; ++k) {
        if (k == 0) {
            facesInfo[k] = &quadInfo;
 
            faceTypeStorage[k] = ElementType::QUAD;
        } else {
            facesInfo[k] = &triangleInfo;
 
            faceTypeStorage[k] = ElementType::TRIANGLE;
        }
    }
 
    faceConnectStorage[0][0] = 0;
    faceConnectStorage[0][1] = 3;
    faceConnectStorage[0][2] = 2;
    faceConnectStorage[0][3] = 1;
 
    faceConnectStorage[1][0] = 3;
    faceConnectStorage[1][1] = 0;
    faceConnectStorage[1][2] = 4;
 
    faceConnectStorage[2][0] = 0;
    faceConnectStorage[2][1] = 1;
    faceConnectStorage[2][2] = 4;
 
    faceConnectStorage[3][0] = 1;
    faceConnectStorage[3][1] = 2;
    faceConnectStorage[3][2] = 4;
 
    faceConnectStorage[4][0] = 2;
    faceConnectStorage[4][1] = 3;
    faceConnectStorage[4][2] = 4;
 
    initializeFaceEdges(facesInfo, edgesInfo);
}

double ReferencePyramidInfo::evalVolume(const std::array<double, 3> *vertexCoords) const
{
    const std::array<double, 3> RA = vertexCoords[0] - vertexCoords[4];
    const std::array<double, 3> DB = vertexCoords[3] - vertexCoords[1];
    const std::array<double, 3> AC = vertexCoords[2] - vertexCoords[0];
    const std::array<double, 3> AD = vertexCoords[1] - vertexCoords[0];
    const std::array<double, 3> AB = vertexCoords[3] - vertexCoords[0];
 
    double volume = dotProduct(RA, crossProduct(DB, AC)) / 6. + dotProduct(AC, crossProduct(AD, AB)) / 12.;
 
    return volume;
}

 
const ReferenceWedgeInfo ReferenceWedgeInfo::info;

ReferenceWedgeInfo::ReferenceWedgeInfo()
    : Reference3DElementInfo(ElementType::WEDGE, 6, 5)
{
    const ReferenceLineInfo lineInfo;
    const ReferenceTriangleInfo triangleInfo;
    const ReferenceQuadInfo quadInfo;
 
    // Edge data
    std::vector<const ReferenceElementInfo *> edgesInfo(nEdges);
 
    for (int k = 0; k < nEdges; ++k) {
        edgeTypeStorage[k]   = LINE;
        edgesInfo[k]   = &lineInfo;
    }
 
    edgeConnectStorage[0][0] = 1;
    edgeConnectStorage[0][1] = 0;
 
    edgeConnectStorage[1][0] = 1;
    edgeConnectStorage[1][1] = 2;
 
    edgeConnectStorage[2][0] = 2;
    edgeConnectStorage[2][1] = 0;
 
    edgeConnectStorage[3][0] = 3;
    edgeConnectStorage[3][1] = 4;
 
    edgeConnectStorage[4][0] = 4;
    edgeConnectStorage[4][1] = 5;
 
    edgeConnectStorage[5][0] = 5;
    edgeConnectStorage[5][1] = 3;
 
    edgeConnectStorage[6][0] = 3;
    edgeConnectStorage[6][1] = 0;
 
    edgeConnectStorage[7][0] = 4;
    edgeConnectStorage[7][1] = 1;
 
    edgeConnectStorage[8][0] = 5;
    edgeConnectStorage[8][1] = 2;
 
    // Face data
    std::vector<const ReferenceElementInfo *> facesInfo(nFaces);
 
    for (int k = 0; k < nFaces; ++k) {
        if (k == 0 || k == 1) {
            facesInfo[k] = &triangleInfo;
 
            faceTypeStorage[k] = TRIANGLE;
        } else {
            facesInfo[k] = &quadInfo;
 
            faceTypeStorage[k] = QUAD;
        }
    }
 
    faceConnectStorage[0][0] = 0;
    faceConnectStorage[0][1] = 1;
    faceConnectStorage[0][2] = 2;
 
    faceConnectStorage[1][0] = 4;
    faceConnectStorage[1][1] = 3;
    faceConnectStorage[1][2] = 5;
 
    faceConnectStorage[2][0] = 4;
    faceConnectStorage[2][1] = 1;
    faceConnectStorage[2][2] = 0;
    faceConnectStorage[2][3] = 3;
 
    faceConnectStorage[3][0] = 1;
    faceConnectStorage[3][1] = 4;
    faceConnectStorage[3][2] = 5;
    faceConnectStorage[3][3] = 2;
 
    faceConnectStorage[4][0] = 3;
    faceConnectStorage[4][1] = 0;
    faceConnectStorage[4][2] = 2;
    faceConnectStorage[4][3] = 5;
 
    initializeFaceEdges(facesInfo, edgesInfo);
}

double ReferenceWedgeInfo::evalVolume(const std::array<double, 3> *vertexCoords) const
{
    std::array<int, 2> triaFaces = {{0, 1}};
    std::array<int, 3> quadFaces = {{2, 3, 4}};
 
    std::array<std::array<double, 3>, MAX_ELEM_VERTICES> pyramidVertexCoords;
    std::array<std::array<double, 3>, MAX_ELEM_VERTICES> tetraVertexCoords;
 
    const ReferenceTriangleInfo &triaInfo    = static_cast<const ReferenceTriangleInfo &>(getInfo(ElementType::TRIANGLE));
    const ReferenceQuadInfo     &quadInfo    = static_cast<const ReferenceQuadInfo &>(getInfo(ElementType::QUAD));
    const ReferenceTetraInfo    &tetraInfo   = static_cast<const ReferenceTetraInfo &>(getInfo(ElementType::TETRA));
    const ReferencePyramidInfo  &pyramidInfo = static_cast<const ReferencePyramidInfo &>(getInfo(ElementType::PYRAMID));
 
    // The centroid is the apex of all the pyramids
    pyramidVertexCoords[4] = vertexCoords[0];
    for (int i = 1; i < nVertices; ++i) {
        pyramidVertexCoords[4] += vertexCoords[i];
    }
    pyramidVertexCoords[4] = pyramidVertexCoords[4] / double(nVertices);
 
    // The centroid is the apex of all the tetras
    tetraVertexCoords[3] = pyramidVertexCoords[4];
 
    // Sum the volume of the pyramids/tetra having a face as the base and the
    // centroid as the apex.
    double volume = 0.;
    for (int i : triaFaces) {
        for (int n = 0; n < triaInfo.nVertices; ++n) {
            tetraVertexCoords[triaInfo.nVertices - n - 1] = vertexCoords[faceConnectStorage[i][n]];
        }
 
        volume += tetraInfo.evalVolume(tetraVertexCoords.data());
    }
 
    for (int i : quadFaces) {
        for (int n = 0; n < quadInfo.nVertices; ++n) {
            pyramidVertexCoords[quadInfo.nVertices - n - 1] = vertexCoords[faceConnectStorage[i][n]];
        }
 
        volume += pyramidInfo.evalVolume(pyramidVertexCoords.data());
    }
 
    return volume;
}


Reference2DElementInfo::Reference2DElementInfo(ElementType type, int nVertices)
    : ReferenceElementInfo(2, type, nVertices, nVertices, nVertices)
{
}

double Reference2DElementInfo::evalSize(const std::array<double, 3> *vertexCoords) const
{
    double area = evalArea(vertexCoords);
 
    double faceMaxLength = 0;
    for (int face = 0; face < nFaces; ++face) {
        double faceLength = evalFaceLength(face, vertexCoords);
        faceMaxLength = std::max(faceLength, faceMaxLength);
    }
 
    double length = area / faceMaxLength;
 
    return length;
}

double Reference2DElementInfo::evalPerimeter(const std::array<double, 3> *vertexCoords) const
{
    double perimeter = 0;
    for (int i = 0; i < nFaces; ++i) {
        perimeter += evalFaceLength(i, vertexCoords);
    }
 
    return perimeter;
}

double Reference2DElementInfo::evalFaceLength(int face, const std::array<double, 3> *vertexCoords) const
{
    const ReferenceLineInfo &sideInfo = static_cast<const ReferenceLineInfo &>(getInfo(ElementType::LINE));
 
    std::array<std::array<double, 3>, MAX_ELEM_VERTICES> sideVertexCoords;
    for (int n = 0; n < sideInfo.nVertices; ++n) {
        sideVertexCoords[n] = vertexCoords[faceConnectStorage[face][n]];
    }
 
    return sideInfo.evalLength(sideVertexCoords.data());
}

void Reference2DElementInfo::evalPointProjection(const std::array<double, 3> &point, const std::array<double, 3> *vertexCoords,
                                                 std::array<double, 3> *projection, double *distance) const
{
    // Get vertices ordered counter-clockwise
    const std::array<double, 3> *ccwVertexCoords;
    std::array<std::array<double, 3>, MAX_ELEM_VERTICES> ccwVertexCoordsStorage;
    getCCWVertexCoords(vertexCoords, &ccwVertexCoords, ccwVertexCoordsStorage.data());
 
    // Evaluate projection
    int projectionFlag;
    *distance = CGElem::distancePointPolygon(point, nVertices, ccwVertexCoords, *projection, projectionFlag);
}

double Reference2DElementInfo::evalPointDistance(const std::array<double, 3> &point, const std::array<double, 3> *vertexCoords) const
{
    // Get vertices ordered counter-clockwise
    const std::array<double, 3> *ccwVertexCoords;
    std::array<std::array<double, 3>, MAX_ELEM_VERTICES> ccwVertexCoordsStorage;
    getCCWVertexCoords(vertexCoords, &ccwVertexCoords, ccwVertexCoordsStorage.data());
 
    // Evaluate distance
    return CGElem::distancePointPolygon(point, nVertices, ccwVertexCoords);
}

bool Reference2DElementInfo::areVerticesCCWOrdered() const
{
    return true;
}

int Reference2DElementInfo::getCCWOrderedVertex(int n) const
{
    return n;
}

bool Reference2DElementInfo::areFacesCCWOrdered() const
{
    return true;
}

int Reference2DElementInfo::getCCWOrderedFace(int n) const
{
    return n;
}

void Reference2DElementInfo::getCCWVertexCoords(const std::array<double, 3> *vertexCoords,
                                                const std::array<double, 3> **ccwVertexCoords,
                                                std::array<double, 3> *ccwVertexCoordsStorage) const
{
    // Early return if vertex are already ordered
    if (areVerticesCCWOrdered()) {
        *ccwVertexCoords = vertexCoords;
        return;
    }
 
    // Order vertices
    for (int i = 0; i < nVertices; ++i) {
        ccwVertexCoordsStorage[i] = vertexCoords[getCCWOrderedVertex(i)];
    }
    *ccwVertexCoords = ccwVertexCoordsStorage;
}

 
const ReferenceTriangleInfo ReferenceTriangleInfo::info;

ReferenceTriangleInfo::ReferenceTriangleInfo()
    : Reference2DElementInfo(ElementType::TRIANGLE, 3)
{
    const ReferenceVertexInfo vertexInfo;
    const ReferenceLineInfo lineInfo;
 
    // Edge data
    std::vector<const ReferenceElementInfo *> edgesInfo(nEdges);
 
    for (int k = 0; k < nEdges; ++k) {
        edgesInfo[k] = &vertexInfo;
 
        edgeTypeStorage[k]       = ElementType::VERTEX;
        edgeConnectStorage[k][0] = k;
    }
 
    // Face data
    std::vector<const ReferenceElementInfo *> facesInfo(nFaces);
 
    for (int k = 0; k < nFaces; ++k) {
        facesInfo[k] = &lineInfo;
 
        faceTypeStorage[k]       = LINE;
        faceConnectStorage[k][0] = k;
        faceConnectStorage[k][1] = (k + 1) % nVertices;
    }
 
    initializeFaceEdges(facesInfo, edgesInfo);
}

double ReferenceTriangleInfo::evalSize(const std::array<double, 3> *vertexCoords) const
{
    double area      = evalArea(vertexCoords);
    double perimeter = evalPerimeter(vertexCoords);
 
    double inscribedRadius = 2 * area / perimeter;
 
    double length = 3. * inscribedRadius;
 
    return length;
}

double ReferenceTriangleInfo::evalArea(const std::array<double, 3> *vertexCoords) const
{
    const std::array<double, 3> &V_A = vertexCoords[0];
    const std::array<double, 3> &V_B = vertexCoords[1];
    const std::array<double, 3> &V_C = vertexCoords[2];
 
    double area = 0.5 * norm2(crossProduct(V_B - V_A, V_C - V_A));
 
    return area;
}

std::array<double, 3> ReferenceTriangleInfo::evalNormal(const std::array<double, 3> *vertexCoords, const std::array<double, 3> &point) const
{
    BITPIT_UNUSED(point);
 
    const std::array<double, 3> &V_A = vertexCoords[0];
    const std::array<double, 3> &V_B = vertexCoords[1];
    const std::array<double, 3> &V_C = vertexCoords[2];
 
    std::array<double, 3> normal = crossProduct(V_B - V_A, V_C - V_A);
    normal = normal / norm2(normal);
 
    return normal;
}

void ReferenceTriangleInfo::evalPointProjection(const std::array<double, 3> &point, const std::array<double, 3> *vertexCoords,
                                                std::array<double, 3> *projection, double *distance) const
{
    *projection = CGElem::projectPointTriangle(point, vertexCoords[0], vertexCoords[1], vertexCoords[2]);
    *distance   = norm2(point - *projection);
}

double ReferenceTriangleInfo::evalPointDistance(const std::array<double, 3> &point, const std::array<double, 3> *vertexCoords) const
{
    return CGElem::distancePointTriangle(point, vertexCoords[0], vertexCoords[1], vertexCoords[2]);
}

 
const ReferencePixelInfo ReferencePixelInfo::info;

ReferencePixelInfo::ReferencePixelInfo()
    : Reference2DElementInfo(ElementType::PIXEL, 4)
{
    const ReferenceVertexInfo vertexInfo;
    const ReferenceLineInfo lineInfo;
 
    // Edge data
    std::vector<const ReferenceElementInfo *> edgesInfo(nEdges);
 
    for (int k = 0; k < nEdges; ++k) {
        edgesInfo[k] = &vertexInfo;
 
        edgeTypeStorage[k]       = ElementType::VERTEX;
        edgeConnectStorage[k][0] = k;
    }
 
    // Face data
    std::vector<const ReferenceElementInfo *> facesInfo(nFaces);
 
    for (int k = 0; k < nFaces; ++k) {
        facesInfo[k] = &lineInfo;
 
        faceTypeStorage[k] = LINE;
    }
 
    faceConnectStorage[0][0] = 2;
    faceConnectStorage[0][1] = 0;
 
    faceConnectStorage[1][0] = 1;
    faceConnectStorage[1][1] = 3;
 
    faceConnectStorage[2][0] = 0;
    faceConnectStorage[2][1] = 1;
 
    faceConnectStorage[3][0] = 3;
    faceConnectStorage[3][1] = 2;
 
    initializeFaceEdges(facesInfo, edgesInfo);
}

double ReferencePixelInfo::evalSize(const std::array<double, 3> *vertexCoords) const
{
    const std::array<double, 3> &V_A = vertexCoords[0];
    const std::array<double, 3> &V_B = vertexCoords[1];
    const std::array<double, 3> &V_C = vertexCoords[2];
 
    double sideLength_x = norm2(V_B - V_A);
    double sideLength_y = norm2(V_C - V_A);
 
    double size = std::min(sideLength_x, sideLength_y);
 
    return size;
}

double ReferencePixelInfo::evalArea(const std::array<double, 3> *vertexCoords) const
{
    const std::array<double, 3> &V_A = vertexCoords[0];
    const std::array<double, 3> &V_B = vertexCoords[1];
    const std::array<double, 3> &V_C = vertexCoords[2];
 
    double area = norm2(crossProduct(V_B - V_A, V_C - V_A));
 
    return area;
}

std::array<double, 3> ReferencePixelInfo::evalNormal(const std::array<double, 3> *vertexCoords, const std::array<double, 3> &point) const
{
    BITPIT_UNUSED(point);
 
    const std::array<double, 3> &V_A = vertexCoords[0];
    const std::array<double, 3> &V_B = vertexCoords[1];
    const std::array<double, 3> &V_C = vertexCoords[2];
 
    std::array<double, 3> normal = crossProduct(V_B - V_A, V_C - V_A);
    normal = normal / norm2(normal);
 
    return normal;
}

bool ReferencePixelInfo::areVerticesCCWOrdered() const
{
    return false;
}

int ReferencePixelInfo::getCCWOrderedVertex(int n) const
{
    static const std::array<int, 4> CCW_ORDERED_VERTICES = {{0, 1, 3, 2}};
 
    return CCW_ORDERED_VERTICES[n];
}

bool ReferencePixelInfo::areFacesCCWOrdered() const
{
    return false;
}

int ReferencePixelInfo::getCCWOrderedFace(int n) const
{
    static const std::array<int, 4> CCW_ORDERED_FACES = {{2, 1, 3, 0}};
 
    return CCW_ORDERED_FACES[n];
}

 
const ReferenceQuadInfo ReferenceQuadInfo::info;

ReferenceQuadInfo::ReferenceQuadInfo()
    : Reference2DElementInfo(ElementType::QUAD, 4)
{
    const ReferenceVertexInfo vertexInfo;
    const ReferenceLineInfo lineInfo;
 
    // Edge data
    std::vector<const ReferenceElementInfo *> edgesInfo(nEdges);
 
    for (int k = 0; k < nEdges; ++k) {
        edgesInfo[k] = &vertexInfo;
 
        edgeTypeStorage[k]       = ElementType::VERTEX;
        edgeConnectStorage[k][0] = k;
    }
 
    // Face data
    std::vector<const ReferenceElementInfo *> facesInfo(nFaces);
 
    for (int k = 0; k < nFaces; ++k) {
        facesInfo[k] = &lineInfo;
 
        faceTypeStorage[k]       = LINE;
        faceConnectStorage[k][0] = k;
        faceConnectStorage[k][1] = (k + 1) % nVertices;
    }
 
    initializeFaceEdges(facesInfo, edgesInfo);
}

double ReferenceQuadInfo::evalArea(const std::array<double, 3> *vertexCoords) const
{
    const std::array<double, 3> &V_A = vertexCoords[0];
    const std::array<double, 3> &V_B = vertexCoords[1];
    const std::array<double, 3> &V_C = vertexCoords[2];
    const std::array<double, 3> &V_D = vertexCoords[3];
 
    std::array<double, 3> mapping = crossProduct(V_B - V_A, V_D - V_A);
    mapping += 0.5 * crossProduct(V_B - V_A, V_C - V_D);
    mapping += 0.5 * crossProduct(V_C - V_B, V_D - V_A);
 
    double area = norm2(mapping);
 
    return area;
}

std::array<double, 3> ReferenceQuadInfo::evalNormal(const std::array<double, 3> *vertexCoords, const std::array<double, 3> &point) const
{
    const std::array<double, 3> AB = vertexCoords[3] - vertexCoords[0];
    const std::array<double, 3> AD = vertexCoords[1] - vertexCoords[0];
    const std::array<double, 3> BC = vertexCoords[2] - vertexCoords[3];
    const std::array<double, 3> DC = vertexCoords[2] - vertexCoords[1];
 
    const double csi = point[0];
    const double eta = point[1];
 
    std::array<double, 3> normal = crossProduct(AB, AD);
    normal += csi * crossProduct(AB, DC);
    normal += eta * crossProduct(BC, AD);
    normal = normal / (- norm2(normal));
 
    return normal;
}


Reference1DElementInfo::Reference1DElementInfo(ElementType type)
    : ReferenceElementInfo(1, type, 2, 2, 2)
{
}

 
const ReferenceLineInfo ReferenceLineInfo::info;

ReferenceLineInfo::ReferenceLineInfo()
    : Reference1DElementInfo(ElementType::LINE)
{
    const ReferenceVertexInfo vertexInfo;
 
    // Edge data
    std::vector<const ReferenceElementInfo *> edgesInfo(nEdges);
 
    for (int k = 0; k < nEdges; ++k) {
        edgesInfo[k] = &vertexInfo;
 
        edgeTypeStorage[k]       = ElementType::VERTEX;
        edgeConnectStorage[k][0] = k;
    }
 
    // Face data
    std::vector<const ReferenceElementInfo *> facesInfo(nFaces);
 
    for (int k = 0; k < nFaces; ++k) {
        facesInfo[k] = &vertexInfo;
 
        faceTypeStorage[k]       = ElementType::VERTEX;
        faceConnectStorage[k][0] = k;
    }
 
    initializeFaceEdges(facesInfo, edgesInfo);
}

double ReferenceLineInfo::evalSize(const std::array<double, 3> *vertexCoords) const
{
    return evalLength(vertexCoords);
}

double ReferenceLineInfo::evalLength(const std::array<double, 3> *vertexCoords) const
{
    const std::array<double, 3> &V_A = vertexCoords[0];
    const std::array<double, 3> &V_B = vertexCoords[1];
 
    double length = norm2(V_B - V_A);
 
    return length;
}

std::array<double, 3> ReferenceLineInfo::evalNormal(const std::array<double, 3> *vertexCoords,
                                                    const std::array<double, 3> &orientation,
                                                    const std::array<double, 3> &point) const
{
    BITPIT_UNUSED(point);
 
    const std::array<double, 3> tangent = vertexCoords[1] - vertexCoords[0];
 
    std::array<double, 3> normal = crossProduct(tangent, orientation);
    normal = normal / norm2(normal);
 
    return normal;
}

void ReferenceLineInfo::evalPointProjection(const std::array<double, 3> &point, const std::array<double, 3> *vertexCoords,
                                            std::array<double, 3> *projection, double *distance) const
{
    std::array<double, 2> lambda;
    *distance = CGElem::distancePointSegment(point, vertexCoords[0], vertexCoords[1], lambda);
 
    *projection = CGElem::reconstructPointFromBarycentricSegment(vertexCoords[0], vertexCoords[1], lambda);
}

double ReferenceLineInfo::evalPointDistance(const std::array<double, 3> &point, const std::array<double, 3> *vertexCoords) const
{
    return CGElem::distancePointSegment(point, vertexCoords[0], vertexCoords[1]);
}


Reference0DElementInfo::Reference0DElementInfo(ElementType type)
    : ReferenceElementInfo(0, type, 1, 1, 1)
{
}

 
const ReferenceVertexInfo ReferenceVertexInfo::info;

ReferenceVertexInfo::ReferenceVertexInfo()
    : Reference0DElementInfo(ElementType::VERTEX)
{
    // Edge data
    std::vector<const ReferenceElementInfo *> edgesInfo(nEdges);
 
    edgesInfo[0] = this;
 
    edgeTypeStorage[0] = ElementType::VERTEX;
 
    edgeConnectStorage[0][0] = 0;
 
    // Face data
    std::vector<const ReferenceElementInfo *> facesInfo(nFaces);
 
    facesInfo[0] = this;
 
    faceTypeStorage[0] = ElementType::VERTEX;
 
    faceConnectStorage[0][0] = 0;
 
    initializeFaceEdges(facesInfo, edgesInfo);
}

double ReferenceVertexInfo::evalSize(const std::array<double, 3> *vertexCoords) const
{
    BITPIT_UNUSED(vertexCoords);
 
    return 0.;
}

std::array<double, 3> ReferenceVertexInfo::evalNormal(const std::array<double, 3> *vertexCoords,
                                                      const std::array<double, 3> &orientation) const
{
    BITPIT_UNUSED(vertexCoords);
 
    std::array<double, 3> normal = orientation;
    normal = normal / norm2(normal);
 
    return normal;
}

void ReferenceVertexInfo::evalPointProjection(const std::array<double, 3> &point, const std::array<double, 3> *vertexCoords,
                                              std::array<double, 3> *projection, double *distance) const
{
    *distance   = evalPointDistance(point, vertexCoords);
    *projection = vertexCoords[0];
}

double ReferenceVertexInfo::evalPointDistance(const std::array<double, 3> &point, const std::array<double, 3> *vertexCoords) const
{
    return norm2(point - vertexCoords[0]);
}
 
}