ShapeType.hpp

Go to the documentation of this file.

////////////////////////////////////////////////////////////////////////////////
//
//  File: ShapeType.hpp
//
//  For more information, please see: http://www.nektar.info/
//
//  The MIT License
//
//  Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
//  Department of Aeronautics, Imperial College London (UK), and Scientific
//  Computing and Imaging Institute, University of Utah (USA).
//
//  Permission is hereby granted, free of charge, to any person obtaining a
//  copy of this software and associated documentation files (the "Software"),
//  to deal in the Software without restriction, including without limitation
//  the rights to use, copy, modify, merge, publish, distribute, sublicense,
//  and/or sell copies of the Software, and to permit persons to whom the
//  Software is furnished to do so, subject to the following conditions:
//
//  The above copyright notice and this permission notice shall be included
//  in all copies or substantial portions of the Software.
//
//  THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
//  OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
//  FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
//  THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
//  LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
//  FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
//  DEALINGS IN THE SOFTWARE.
//
//  Description:
//
////////////////////////////////////////////////////////////////////////////////
 
#ifndef NEKTAR_LIB_UTILITIES_BASIC_UTILS_SHAPE_TYPE_H
#define NEKTAR_LIB_UTILITIES_BASIC_UTILS_SHAPE_TYPE_H
 
#include <algorithm>
#include <vector>
 
#include <LibUtilities/BasicUtils/ErrorUtil.hpp>
 
#ifdef min
#undef min
#endif
 
namespace Nektar::LibUtilities
{
 
// Types of geometry types.
enum ShapeType
{
    eNoShapeType,
    ePoint,
    eSegment,
    eTriangle,
    eQuadrilateral,
    eTetrahedron,
    ePyramid,
    ePrism,
    eHexahedron,
    eNodalTri,
    eNodalTet,
    eNodalPrism,
    SIZE_ShapeType,
 
    // These are the short names used for MatrixFree operators
    Point      = ePoint,
    Seg        = eSegment,
    Tri        = eTriangle,
    Quad       = eQuadrilateral,
    Tet        = eTetrahedron,
    Pyr        = ePyramid,
    Prism      = ePrism,
    Hex        = eHexahedron,
    NodalTri   = eNodalTri,
    NodalTet   = eNodalTet,
    NodalPrism = eNodalPrism
};
 
const char *const ShapeTypeMap[SIZE_ShapeType] = {
    "NoGeomShapeType", "Point",       "Segment",  "Triangle",
    "Quadrilateral",   "Tetrahedron", "Pyramid",  "Prism",
    "Hexahedron",      "NodalTri",    "NodalTet", "NodalPrism"};
 
// Hold the dimension of each of the types of shapes.
constexpr unsigned int ShapeTypeDimMap[SIZE_ShapeType] = {
    0, // Unknown
    0, // ePoint
    1, // eSegment
    2, // eTriangle
    2, // eQuadrilateral
    3, // eTetrahedron
    3, // ePyramid
    3, // ePrism
    3, // eHexahedron
    2, // eNodalTtri
    3, // eNodalTet
    3, // eNodalPrism
};
 
namespace StdSegData
{
inline constexpr int getNumberOfCoefficients(int Na)
{
    return Na;
}
 
inline constexpr int getNumberOfBndCoefficients([[maybe_unused]] int Na)
{
    return 2;
}
} // namespace StdSegData
 
// Dimensions of coefficients for each space
namespace StdTriData
{
inline constexpr int getNumberOfCoefficients(int Na, int Nb)
{
    // Note these assertions have been set to > 0 because
    // it can also be used to evaluate face expansion
    // order
    ASSERTL2(Na > 0, "Order in 'a' direction must be > 0.");
    ASSERTL2(Nb > 0, "Order in 'b' direction must be > 0.");
    ASSERTL1(Na <= Nb, "order in 'a' direction is higher "
                       "than order in 'b' direction");
    return Na * (Na + 1) / 2 + Na * (Nb - Na);
}
 
inline constexpr int getNumberOfBndCoefficients(int Na, int Nb)
{
    ASSERTL2(Na > 1, "Order in 'a' direction must be > 1.");
    ASSERTL2(Nb > 1, "Order in 'b' direction must be > 1.");
    ASSERTL1(Na <= Nb, "order in 'a' direction is higher "
                       "than order in 'b' direction");
    return (Na - 1) + 2 * (Nb - 1);
}
} // namespace StdTriData
 
// Dimensions of coefficients for each space
namespace StdNodalTriData
{
inline constexpr int getNumberOfCoefficients(int Na, int Nb)
{
    ASSERTL1(Na == Nb, "order in 'a' direction needs to be the same as "
                       "b' direction for nodaltri");
    return Na * (Nb + 1) / 2;
}
 
inline constexpr int getNumberOfBndCoefficients(int Na, int Nb)
{
    ASSERTL1(Na == Nb, "order in 'a' direction needs to be the same as "
                       "b' direction for nodaltri");
    return (Na - 1) + 2 * (Nb - 1);
}
} // namespace StdNodalTriData
 
namespace StdQuadData
{
inline constexpr int getNumberOfCoefficients(int Na, int Nb)
{
    // Note these assertions have been set to > 0 because
    // it can also be used to evaluate face expansion
    // order
    ASSERTL2(Na > 0, "Order in 'a' direction must be > 0.");
    ASSERTL2(Nb > 0, "Order in 'b' direction must be > 0.");
    return Na * Nb;
}
 
inline constexpr int getNumberOfBndCoefficients(int Na, int Nb)
{
    ASSERTL2(Na > 1, "Order in 'a' direction must be > 1.");
    ASSERTL2(Nb > 1, "Order in 'b' direction must be > 1.");
    return 2 * (Na - 1) + 2 * (Nb - 1);
}
} // namespace StdQuadData
 
namespace StdHexData
{
inline constexpr int getNumberOfCoefficients(int Na, int Nb, int Nc)
{
    ASSERTL2(Na > 1, "Order in 'a' direction must be > 1.");
    ASSERTL2(Nb > 1, "Order in 'b' direction must be > 1.");
    ASSERTL2(Nc > 1, "Order in 'c' direction must be > 1.");
    return Na * Nb * Nc;
}
 
inline constexpr int getNumberOfBndCoefficients(int Na, int Nb, int Nc)
{
    ASSERTL2(Na > 1, "Order in 'a' direction must be > 1.");
    ASSERTL2(Nb > 1, "Order in 'b' direction must be > 1.");
    ASSERTL2(Nc > 1, "Order in 'c' direction must be > 1.");
    return 2 * Na * Nb + 2 * Na * Nc + 2 * Nb * Nc - 4 * (Na + Nb + Nc) + 8;
}
} // namespace StdHexData
 
namespace StdTetData
{
/**
 * Adds up the number of cells in a truncated Nc by Nc by Nc
 * pyramid, where the longest Na rows and longest Nb columns are
 * kept. Example: (Na, Nb, Nc) = (3, 4, 5); The number of
 * coefficients is the sum of the elements of the following
 * matrix:
 *
 * |5  4  3  2  0|
 * |4  3  2  0   |
 * |3  2  0      |
 * |0  0         |
 * |0            |
 *
 * Sum = 28 = number of tet coefficients.
 */
inline constexpr int getNumberOfCoefficients(int Na, int Nb, int Nc)
{
    ASSERTL2(Na > 1, "Order in 'a' direction must be > 1.");
    ASSERTL2(Nb > 1, "Order in 'b' direction must be > 1.");
    ASSERTL2(Nc > 1, "Order in 'c' direction must be > 1.");
    ASSERTL1(Na <= Nc, "order in 'a' direction is higher "
                       "than order in 'c' direction");
    ASSERTL1(Nb <= Nc, "order in 'b' direction is higher "
                       "than order in 'c' direction");
    int nCoef = Na * Nb * Nc - Na * Nc * (Na - 1) / 2 - Na * Nb * (Nb - 1) / 2 +
                Na * (Na - 1) * (Na - 2) / 6;
    return nCoef;
}
inline constexpr int getNumberOfBndCoefficients(int Na, int Nb, int Nc)
{
    ASSERTL2(Na > 1, "Order in 'a' direction must be > 1.");
    ASSERTL2(Nb > 1, "Order in 'b' direction must be > 1.");
    ASSERTL2(Nc > 1, "Order in 'c' direction must be > 1.");
    ASSERTL1(Na <= Nc, "order in 'a' direction is higher "
                       "than order in 'c' direction");
    ASSERTL1(Nb <= Nc, "order in 'b' direction is higher "
                       "than order in 'c' direction");
 
    int nCoef = Na * (Na + 1) / 2 + (Nb - Na) * Na         // base
                + Na * (Na + 1) / 2 + (Nc - Na) * Na       // front
                + 2 * (Nb * (Nb + 1) / 2 + (Nc - Nb) * Nb) // 2 other sides
                - Na - 2 * Nb - 3 * Nc                     // less edges
                + 4;                                       // plus vertices
 
    return nCoef;
}
} // namespace StdTetData
 
namespace StdNodalTetData
{
/**
 * Adds up the number of cells in a truncated Nc by Nc by Nc
 * pyramid, where the longest Na rows and longest Nb columns are
 * kept. Example: (Na, Nb, Nc) = (3, 4, 5); The number of
 * coefficients is the sum of the elements of the following
 * matrix:
 *
 * |5  4  3  2  0|
 * |4  3  2  0   |
 * |3  2  0      |
 * |0  0         |
 * |0            |
 *
 * Sum = 28 = number of tet coefficients.
 */
inline constexpr int getNumberOfCoefficients(int Na, int Nb, int Nc)
{
    ASSERTL1(Na == Nb, "order in 'a' direction needs to be the same as "
                       "b' direction for nodaltet");
    ASSERTL1(Nb == Nc, "order in 'a' direction needs to be the same as "
                       "c' direction for nodaltet");
    return Na * (Nb + 1) * (Nc + 2) / 6;
}
 
inline constexpr int getNumberOfBndCoefficients(int Na, int Nb, int Nc)
{
    ASSERTL1(Na == Nb, "order in 'a' direction needs to be the same as "
                       "b' direction for nodaltet");
    ASSERTL1(Nb == Nc, "order in 'a' direction needs to be the same as "
                       "c' direction for nodaltet");
    int nCoef = Na * (Na + 1) / 2 + (Nb - Na) * Na         // base
                + Na * (Na + 1) / 2 + (Nc - Na) * Na       // front
                + 2 * (Nb * (Nb + 1) / 2 + (Nc - Nb) * Nb) // 2 other sides
                - Na - 2 * Nb - 3 * Nc                     // less edges
                + 4;                                       // plus vertices
 
    return nCoef;
}
} // namespace StdNodalTetData
 
namespace StdPyrData
{
inline constexpr int getNumberOfCoefficients(int Na, int Nb, int Nc)
{
    ASSERTL1(Na > 1, "Order in 'a' direction must be > 1.");
    ASSERTL1(Nb > 1, "Order in 'b' direction must be > 1.");
    ASSERTL1(Nc > 1, "Order in 'c' direction must be > 1.");
    ASSERTL1(Na <= Nc, "Order in 'a' direction is higher "
                       "than order in 'c' direction.");
    ASSERTL1(Nb <= Nc, "Order in 'b' direction is higher "
                       "than order in 'c' direction.");
 
    // Count number of coefficients explicitly.
    int nCoeff = 0;
 
    // Count number of interior tet modes
    for (int a = 0; a < Na; ++a)
    {
        for (int b = 0; b < Nb; ++b)
        {
            for (int c = 0; c < Nc - std::max(a, b); ++c)
            {
                ++nCoeff;
            }
        }
    }
    return nCoeff;
}
 
inline constexpr int getNumberOfBndCoefficients(int Na, int Nb, int Nc)
{
    ASSERTL1(Na > 1, "Order in 'a' direction must be > 1.");
    ASSERTL1(Nb > 1, "Order in 'b' direction must be > 1.");
    ASSERTL1(Nc > 1, "Order in 'c' direction must be > 1.");
    ASSERTL1(Na <= Nc, "Order in 'a' direction is higher "
                       "than order in 'c' direction.");
    ASSERTL1(Nb <= Nc, "Order in 'b' direction is higher "
                       "than order in 'c' direction.");
 
    return Na * Nb                                    // base
           + 2 * (Na * (Na + 1) / 2 + (Nc - Na) * Na) // front and back
           + 2 * (Nb * (Nb + 1) / 2 + (Nc - Nb) * Nb) // sides
           - 2 * Na - 2 * Nb - 4 * Nc                 // less edges
           + 5;                                       // plus vertices
}
} // namespace StdPyrData
 
namespace StdPrismData
{
inline constexpr int getNumberOfCoefficients(int Na, int Nb, int Nc)
{
    ASSERTL1(Na > 1, "Order in 'a' direction must be > 1.");
    ASSERTL1(Nb > 1, "Order in 'b' direction must be > 1.");
    ASSERTL1(Nc > 1, "Order in 'c' direction must be > 1.");
    ASSERTL1(Na <= Nc, "Order in 'a' direction is higher "
                       "than order in 'c' direction.");
 
    return Nb * StdTriData::getNumberOfCoefficients(Na, Nc);
}
 
inline constexpr int getNumberOfBndCoefficients(int Na, int Nb, int Nc)
{
    ASSERTL1(Na > 1, "Order in 'a' direction must be > 1.");
    ASSERTL1(Nb > 1, "Order in 'b' direction must be > 1.");
    ASSERTL1(Nc > 1, "Order in 'c' direction must be > 1.");
    ASSERTL1(Na <= Nc, "Order in 'a' direction is higher "
                       "than order in 'c' direction.");
 
    return Na * Nb + 2 * Nb * Nc                      // rect faces
           + 2 * (Na * (Na + 1) / 2 + (Nc - Na) * Na) // tri faces
           - 2 * Na - 3 * Nb - 4 * Nc                 // less edges
           + 6;                                       // plus vertices
}
} // namespace StdPrismData
 
namespace StdNodalPrismData
{
inline constexpr int getNumberOfCoefficients(int Na, int Nb, int Nc)
{
    ASSERTL1(Na > 1, "Order in 'a' direction must be > 1.");
    ASSERTL1(Nb > 1, "Order in 'b' direction must be > 1.");
    ASSERTL1(Nc > 1, "Order in 'c' direction must be > 1.");
    ASSERTL1(Na <= Nc, "Order in 'a' direction is higher "
                       "than order in 'c' direction.");
 
    return Nb * StdTriData::getNumberOfCoefficients(Na, Nc);
}
 
inline constexpr int getNumberOfBndCoefficients(int Na, int Nb, int Nc)
{
    ASSERTL1(Na > 1, "Order in 'a' direction must be > 1.");
    ASSERTL1(Nb > 1, "Order in 'b' direction must be > 1.");
    ASSERTL1(Nc > 1, "Order in 'c' direction must be > 1.");
    ASSERTL1(Na <= Nc, "Order in 'a' direction is higher "
                       "than order in 'c' direction.");
 
    return Na * Nb + 2 * Nb * Nc                      // rect faces
           + 2 * (Na * (Na + 1) / 2 + (Nc - Na) * Na) // tri faces
           - 2 * Na - 3 * Nb - 4 * Nc                 // less edges
           + 6;                                       // plus vertices
}
} // namespace StdNodalPrismData
 
inline constexpr int GetNumberOfCoefficients(ShapeType shape,
                                             std::vector<unsigned int> &modes,
                                             int offset = 0)
{
    int returnval = 0;
    switch (shape)
    {
        case eSegment:
            returnval = modes[offset];
            break;
        case eTriangle:
            returnval = StdTriData::getNumberOfCoefficients(modes[offset],
                                                            modes[offset + 1]);
            break;
        case eNodalTri:
            returnval = StdNodalTriData::getNumberOfCoefficients(
                modes[offset], modes[offset + 1]);
            break;
        case eQuadrilateral:
            returnval = modes[offset] * modes[offset + 1];
            break;
        case eTetrahedron:
            returnval = StdTetData::getNumberOfCoefficients(
                modes[offset], modes[offset + 1], modes[offset + 2]);
            break;
        case ePyramid:
            returnval = StdPyrData::getNumberOfCoefficients(
                modes[offset], modes[offset + 1], modes[offset + 2]);
            break;
        case ePrism:
            returnval = StdPrismData::getNumberOfCoefficients(
                modes[offset], modes[offset + 1], modes[offset + 2]);
            break;
        case eHexahedron:
            returnval = modes[offset] * modes[offset + 1] * modes[offset + 2];
            break;
        default:
            NEKERROR(ErrorUtil::efatal, "Unknown Shape Type");
            break;
    }
 
    return returnval;
}
 
inline constexpr int GetNumberOfCoefficients(ShapeType shape, int na,
                                             int nb = 0, int nc = 0)
{
    int returnval = 0;
    switch (shape)
    {
        case eSegment:
            returnval = na;
            break;
        case eTriangle:
            returnval = StdTriData::getNumberOfCoefficients(na, nb);
            break;
        case eNodalTri:
            returnval = StdNodalTriData::getNumberOfCoefficients(na, nb);
            break;
        case eQuadrilateral:
            returnval = na * nb;
            break;
        case eTetrahedron:
            returnval = StdTetData::getNumberOfCoefficients(na, nb, nc);
            break;
        case eNodalTet:
            returnval = StdNodalTetData::getNumberOfCoefficients(na, nb, nc);
            break;
        case ePyramid:
            returnval = StdPyrData::getNumberOfCoefficients(na, nb, nc);
            break;
        case ePrism:
            returnval = StdPrismData::getNumberOfCoefficients(na, nb, nc);
            break;
        case eNodalPrism:
            returnval = StdNodalPrismData::getNumberOfCoefficients(na, nb, nc);
            break;
        case eHexahedron:
            returnval = na * nb * nc;
            break;
        default:
            NEKERROR(ErrorUtil::efatal, "Unknown Shape Type");
            break;
    }
 
    return returnval;
}
} // namespace Nektar::LibUtilities
 
#endif