TriGeom.cpp

Go to the documentation of this file.

////////////////////////////////////////////////////////////////////////////////
//
//  File: TriGeom.cpp
//
//  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:
//
//
////////////////////////////////////////////////////////////////////////////////
 
#include <LibUtilities/Foundations/Interp.h>
#include <LibUtilities/Foundations/ManagerAccess.h>
#include <SpatialDomains/TriGeom.h>
#include <StdRegions/StdNodalTriExp.h>
 
#include <SpatialDomains/GeomFactors.h>
#include <SpatialDomains/SegGeom.h>
#include <SpatialDomains/XmapFactory.hpp>
 
namespace Nektar::SpatialDomains
{
 
XmapFactory<StdRegions::StdTriExp, 2> &GetStdTriFactory()
{
    static XmapFactory<StdRegions::StdTriExp, 2> factory;
    return factory;
}
 
TriGeom::TriGeom()
{
    m_shapeType = LibUtilities::eTriangle;
}
 
TriGeom::TriGeom(const int id, std::array<SegGeom *, kNedges> edges,
                 Curve *curve)
    : Geometry2D(edges[0]->GetVertex(0)->GetCoordim(), curve)
{
    int j;
 
    m_shapeType = LibUtilities::eTriangle;
    m_globalID  = id;
 
    /// Copy the edge pointers
    m_edges = edges;
 
    for (j = 0; j < kNedges; ++j)
    {
        m_eorient[j] =
            SegGeom::GetEdgeOrientation(*edges[j], *edges[(j + 1) % kNedges]);
        m_verts[j] =
            edges[j]->GetVertex(m_eorient[j] == StdRegions::eForwards ? 0 : 1);
    }
 
    m_eorient[2] = m_eorient[2] == StdRegions::eBackwards
                       ? StdRegions::eForwards
                       : StdRegions::eBackwards;
 
    m_coordim = edges[0]->GetVertex(0)->GetCoordim();
    if (m_coordim <= 1)
    {
        NEKERROR(ErrorUtil::efatal, "Cannot call function with dim == 1");
    }
}
 
TriGeom::TriGeom(const TriGeom &in) : Geometry2D(in)
{
    // From Geometry
    m_shapeType = in.m_shapeType;
 
    // From TriFaceComponent
    m_globalID = in.m_globalID;
 
    // From TriGeom
    m_verts = in.m_verts;
    m_edges = in.m_edges;
    for (int i = 0; i < kNedges; i++)
    {
        m_eorient[i] = in.m_eorient[i];
    }
}
 
int TriGeom::v_AllLeftCheck(const Array<OneD, const NekDouble> &gloCoord)
{
    int nc = 1, d0 = m_manifold[0], d1 = m_manifold[1];
    if (0 == m_edgeNormal.size())
    {
        m_edgeNormal = Array<OneD, Array<OneD, NekDouble>>(m_verts.size());
        Array<OneD, Array<OneD, NekDouble>> x(2);
        x[0] = Array<OneD, NekDouble>(3);
        x[1] = Array<OneD, NekDouble>(3);
        m_verts[0]->GetCoords(x[0]);
        int i0 = 1, i1 = 0, direction = 1;
        for (size_t i = 0; i < m_verts.size(); ++i)
        {
            i0 ^= 1;
            i1 ^= 1;
            m_verts[(i + 1) % m_verts.size()]->GetCoords(x[i1]);
            if (m_edges[i]->GetXmap()->GetBasis(0)->GetNumModes() > 2)
            {
                continue;
            }
            m_edgeNormal[i]    = Array<OneD, NekDouble>(2);
            m_edgeNormal[i][0] = x[i0][d1] - x[i1][d1];
            m_edgeNormal[i][1] = x[i1][d0] - x[i0][d0];
        }
        if (m_coordim == 3)
        {
            for (size_t i = 0; i < m_verts.size(); ++i)
            {
                if (m_edgeNormal[i].size() == 2)
                {
                    m_verts[i]->GetCoords(x[0]);
                    m_verts[(i + 2) % m_verts.size()]->GetCoords(x[1]);
                    if (m_edgeNormal[i][0] * (x[1][d0] - x[0][d0]) <
                        m_edgeNormal[i][1] * (x[0][d1] - x[1][d1]))
                    {
                        direction = -1;
                    }
                    break;
                }
            }
        }
        if (direction == -1)
        {
            for (size_t i = 0; i < m_verts.size(); ++i)
            {
                if (m_edgeNormal[i].size() == 2)
                {
                    m_edgeNormal[i][0] = -m_edgeNormal[i][0];
                    m_edgeNormal[i][1] = -m_edgeNormal[i][1];
                }
            }
        }
    }
 
    Array<OneD, NekDouble> vertex(3);
    for (size_t i = 0; i < m_verts.size(); ++i)
    {
        int i1 = (i + 1) % m_verts.size();
        if (m_verts[i]->GetGlobalID() < m_verts[i1]->GetGlobalID())
        {
            m_verts[i]->GetCoords(vertex);
        }
        else
        {
            m_verts[i1]->GetCoords(vertex);
        }
        if (m_edgeNormal[i].size() == 0)
        {
            nc = 0; // not sure
            continue;
        }
        if (m_edgeNormal[i][0] * (gloCoord[d0] - vertex[d0]) <
            m_edgeNormal[i][1] * (vertex[d1] - gloCoord[d1]))
        {
            return -1; // outside
        }
    }
    // 3D manifold needs to check the distance
    if (m_coordim == 3)
    {
        nc = 0;
    }
    // nc: 1 (side element), 0 (maybe inside), -1 (outside)
    return nc;
}
 
NekDouble TriGeom::v_GetCoord(const int i,
                              const Array<OneD, const NekDouble> &Lcoord)
{
    if (m_state != ePtsFilled)
    {
        NEKERROR(ErrorUtil::ewarning, "Geometry is not in physical space");
    }
 
    Array<OneD, NekDouble> tmp(m_xmap->GetTotPoints());
    m_xmap->BwdTrans(m_coeffs[i], tmp);
 
    return m_xmap->PhysEvaluate(Lcoord, tmp);
}
 
StdRegions::Orientation TriGeom::GetFaceOrientation(const TriGeom &face1,
                                                    const TriGeom &face2,
                                                    bool doRot, int dir,
                                                    NekDouble angle,
                                                    NekDouble tol)
{
    return GetFaceOrientation(face1.m_verts, face2.m_verts, doRot, dir, angle,
                              tol);
}
 
StdRegions::Orientation TriGeom::GetFaceOrientation(
    std::array<PointGeom *, 3> face1, std::array<PointGeom *, 3> face2,
    bool doRot, int dir, NekDouble angle, NekDouble tol)
{
    int i, j, vmap[3] = {-1, -1, -1};
 
    if (doRot)
    {
        PointGeom rotPt;
 
        for (i = 0; i < 3; ++i)
        {
            rotPt.Rotate((*face1[i]), dir, angle);
            for (j = 0; j < 3; ++j)
            {
                if (rotPt.dist(*face2[j]) < tol)
                {
                    vmap[j] = i;
                    break;
                }
            }
        }
    }
    else
    {
 
        NekDouble x, y, z, x1, y1, z1, cx = 0.0, cy = 0.0, cz = 0.0;
 
        // For periodic faces, we calculate the vector between the centre
        // points of the two faces. (For connected faces this will be
        // zero). We can then use this to determine alignment later in the
        // algorithm.
        for (i = 0; i < 3; ++i)
        {
            cx += (*face2[i])(0) - (*face1[i])(0);
            cy += (*face2[i])(1) - (*face1[i])(1);
            cz += (*face2[i])(2) - (*face1[i])(2);
        }
        cx /= 3;
        cy /= 3;
        cz /= 3;
 
        // Now construct a mapping which takes us from the vertices of one
        // face to the other. That is, vertex j of face2 corresponds to
        // vertex vmap[j] of face1.
        for (i = 0; i < 3; ++i)
        {
            x = (*face1[i])(0);
            y = (*face1[i])(1);
            z = (*face1[i])(2);
            for (j = 0; j < 3; ++j)
            {
                x1 = (*face2[j])(0) - cx;
                y1 = (*face2[j])(1) - cy;
                z1 = (*face2[j])(2) - cz;
                if (sqrt((x1 - x) * (x1 - x) + (y1 - y) * (y1 - y) +
                         (z1 - z) * (z1 - z)) < 1e-8)
                {
                    vmap[j] = i;
                    break;
                }
            }
        }
    }
 
    if (vmap[1] == (vmap[0] + 1) % 3)
    {
        switch (vmap[0])
        {
            case 0:
                return StdRegions::eDir1FwdDir1_Dir2FwdDir2;
                break;
            case 1:
                return StdRegions::eDir1FwdDir2_Dir2BwdDir1;
                break;
            case 2:
                return StdRegions::eDir1BwdDir1_Dir2BwdDir2;
                break;
        }
    }
    else
    {
        switch (vmap[0])
        {
            case 0:
                return StdRegions::eDir1FwdDir2_Dir2FwdDir1;
                break;
            case 1:
                return StdRegions::eDir1BwdDir1_Dir2FwdDir2;
                break;
            case 2:
                return StdRegions::eDir1BwdDir2_Dir2BwdDir1;
                break;
        }
    }
    NEKERROR(ErrorUtil::efatal, "Unable to determine triangle orientation");
    return StdRegions::eNoOrientation;
}
 
GeomType TriGeom::v_CalcGeomType()
{
    if (!m_setupState)
    {
        TriGeom::v_Setup();
    }
    TriGeom::v_FillGeom();
 
    GeomType Gtype = eRegular;
 
    // check to see if expansions are linear
    if (m_xmap->GetBasisNumModes(0) != 2 || m_xmap->GetBasisNumModes(1) != 2)
    {
        Gtype = eDeformed;
    }
 
    m_manifold    = Array<OneD, int>(m_coordim);
    m_manifold[0] = 0;
    m_manifold[1] = 1;
    if (m_coordim == 3)
    {
        PointGeom e01, e21, norm;
        e01.Sub(*m_verts[0], *m_verts[1]);
        e21.Sub(*m_verts[2], *m_verts[1]);
        norm.Mult(e01, e21);
        int tmpi   = 0;
        double tmp = std::fabs(norm[0]);
        if (tmp < fabs(norm[1]))
        {
            tmp  = fabs(norm[1]);
            tmpi = 1;
        }
        if (tmp < fabs(norm[2]))
        {
            tmpi = 2;
        }
        m_manifold[0] = (tmpi + 1) % 3;
        m_manifold[1] = (tmpi + 2) % 3;
        m_manifold[2] = (tmpi + 3) % 3;
    }
    if (Gtype == eRegular)
    {
        Array<OneD, Array<OneD, NekDouble>> verts(m_verts.size());
        for (int i = 0; i < m_verts.size(); ++i)
        {
            verts[i] = Array<OneD, NekDouble>(3);
            m_verts[i]->GetCoords(verts[i]);
        }
        // a00 + a01 xi1 + a02 xi2
        // a10 + a11 xi1 + a12 xi2
        m_isoParameter = Array<OneD, Array<OneD, NekDouble>>(2);
        for (int i = 0; i < 2; ++i)
        {
            unsigned int d       = m_manifold[i];
            m_isoParameter[i]    = Array<OneD, NekDouble>(3, 0.);
            NekDouble A          = verts[0][d];
            NekDouble B          = verts[1][d];
            NekDouble C          = verts[2][d];
            m_isoParameter[i][0] = 0.5 * (B + C);  // 1
            m_isoParameter[i][1] = 0.5 * (-A + B); // xi1
            m_isoParameter[i][2] = 0.5 * (-A + C); // xi2
        }
    }
 
    if (Gtype == eRegular)
    {
        v_CalculateInverseIsoParam();
    }
 
    return Gtype;
}
 
GeomFactorsUniquePtr TriGeom::v_GenGeomFactors(
    LibUtilities::PointsKeyVector &keyTgt)
{
    GeomType Gtype = CalcGeomType();
 
    return ObjPoolManager<GeomFactors>::AllocateUniquePtr(
        Gtype, m_coordim, m_xmap, m_coeffs, keyTgt);
}
 
/**
 * Note verts and edges are listed according to anticlockwise
 * convention but points in _coeffs have to be in array format from
 * left to right.
 */
void TriGeom::v_FillGeom()
{
    // check to see if geometry structure is already filled
    if (m_state == ePtsFilled)
    {
        return;
    }
 
    int i, j, k;
    int nEdgeCoeffs = m_xmap->GetTraceNcoeffs(0);
 
    if (m_curve)
    {
        int pdim = LibUtilities::PointsManager()[LibUtilities::PointsKey(
                                                     2, m_curve->m_ptype)]
                       ->GetPointsDim();
 
        // Deal with 2D points type separately
        // (e.g. electrostatic or Fekete points) to 1D tensor
        // product.
        if (pdim == 2)
        {
            int N = m_curve->m_points.size();
            int nEdgePts =
                (-1 + (int)sqrt(static_cast<NekDouble>(8 * N + 1))) / 2;
 
            ASSERTL0(nEdgePts * (nEdgePts + 1) / 2 == N,
                     "NUMPOINTS should be a triangle number for"
                     " triangle curved face " +
                         std::to_string(m_globalID));
 
            // Sanity check 1: are curved vertices consistent with
            // triangle vertices?
            for (i = 0; i < 3; ++i)
            {
                NekDouble dist = m_verts[i]->dist(*(m_curve->m_points[i]));
                if (dist > NekConstants::kVertexTheSameDouble)
                {
                    std::stringstream ss;
                    ss << "Curved vertex " << i << " of triangle " << m_globalID
                       << " is separated from expansion vertex by"
                       << " more than " << NekConstants::kVertexTheSameDouble
                       << " (dist = " << dist << ")";
                    NEKERROR(ErrorUtil::ewarning, ss.str().c_str());
                }
            }
 
            // Sanity check 2: are curved edges from the face curvature
            // consistent with curved edges?
            for (i = 0; i < kNedges; ++i)
            {
                Curve *edgeCurve = m_edges[i]->GetCurve();
 
                ASSERTL0(edgeCurve->m_points.size() == nEdgePts,
                         "Number of edge points does not correspond "
                         "to number of face points in triangle " +
                             std::to_string(m_globalID));
 
                const int offset  = 3 + i * (nEdgePts - 2);
                NekDouble maxDist = 0.0;
 
                // Account for different ordering of nodal coordinates
                // vs. Cartesian ordering of element.
                StdRegions::Orientation orient = m_eorient[i];
 
                if (i == 2)
                {
                    orient = orient == StdRegions::eForwards
                                 ? StdRegions::eBackwards
                                 : StdRegions::eForwards;
                }
 
                if (orient == StdRegions::eForwards)
                {
                    for (j = 0; j < nEdgePts - 2; ++j)
                    {
                        NekDouble dist = m_curve->m_points[offset + j]->dist(
                            *(edgeCurve->m_points[j + 1]));
                        maxDist = dist > maxDist ? dist : maxDist;
                    }
                }
                else
                {
                    for (j = 0; j < nEdgePts - 2; ++j)
                    {
                        NekDouble dist = m_curve->m_points[offset + j]->dist(
                            *(edgeCurve->m_points[nEdgePts - 2 - j]));
                        maxDist = dist > maxDist ? dist : maxDist;
                    }
                }
 
                if (maxDist > NekConstants::kVertexTheSameDouble)
                {
                    std::stringstream ss;
                    ss << "Curved edge " << i << " of triangle " << m_globalID
                       << " has a point separated from edge interior"
                       << " points by more than "
                       << NekConstants::kVertexTheSameDouble
                       << " (maxdist = " << maxDist << ")";
                    NEKERROR(ErrorUtil::ewarning, ss.str().c_str());
                }
            }
 
            const LibUtilities::PointsKey P0(
                nEdgePts, LibUtilities::eGaussLobattoLegendre);
            const LibUtilities::PointsKey P1(
                nEdgePts, LibUtilities::eGaussRadauMAlpha1Beta0);
            const LibUtilities::BasisKey T0(LibUtilities::eOrtho_A, nEdgePts,
                                            P0);
            const LibUtilities::BasisKey T1(LibUtilities::eOrtho_B, nEdgePts,
                                            P1);
            Array<OneD, NekDouble> phys(
                std::max(nEdgePts * nEdgePts, m_xmap->GetTotPoints()));
            Array<OneD, NekDouble> tmp(nEdgePts * nEdgePts);
 
            for (i = 0; i < m_coordim; ++i)
            {
                // Create a StdNodalTriExp.
                StdRegions::StdNodalTriExpSharedPtr t =
                    MemoryManager<StdRegions::StdNodalTriExp>::
                        AllocateSharedPtr(T0, T1, m_curve->m_ptype);
 
                for (j = 0; j < N; ++j)
                {
                    phys[j] = (m_curve->m_points[j]->GetPtr())[i];
                }
 
                t->BwdTrans(phys, tmp);
 
                // Interpolate points to standard region.
                LibUtilities::Interp2D(
                    P0, P1, tmp, m_xmap->GetBasis(0)->GetPointsKey(),
                    m_xmap->GetBasis(1)->GetPointsKey(), phys);
 
                // Forwards transform to get coefficient space.
                m_xmap->FwdTrans(phys, m_coeffs[i]);
            }
        }
        else if (pdim == 1)
        {
            int npts     = m_curve->m_points.size();
            int nEdgePts = (int)sqrt(static_cast<NekDouble>(npts));
            Array<OneD, NekDouble> tmp(npts);
            Array<OneD, NekDouble> phys(m_xmap->GetTotPoints());
            LibUtilities::PointsKey curveKey(nEdgePts, m_curve->m_ptype);
 
            // Sanity checks:
            // - Curved faces should have square number of points;
            // - Each edge should have sqrt(npts) points.
            ASSERTL0(nEdgePts * nEdgePts == npts,
                     "NUMPOINTS should be a square number for"
                     " triangle " +
                         std::to_string(m_globalID));
 
            for (i = 0; i < kNedges; ++i)
            {
                ASSERTL0(m_edges[i]->GetXmap()->GetNcoeffs() == nEdgePts,
                         "Number of edge points does not correspond to "
                         "number of face points in triangle " +
                             std::to_string(m_globalID));
            }
 
            for (i = 0; i < m_coordim; ++i)
            {
                for (j = 0; j < npts; ++j)
                {
                    tmp[j] = (m_curve->m_points[j]->GetPtr())[i];
                }
 
                // Interpolate curve points to standard triangle
                // points.
                LibUtilities::Interp2D(curveKey, curveKey, tmp,
                                       m_xmap->GetBasis(0)->GetPointsKey(),
                                       m_xmap->GetBasis(1)->GetPointsKey(),
                                       phys);
 
                // Forwards transform to get coefficient space.
                m_xmap->FwdTrans(phys, m_coeffs[i]);
            }
        }
        else
        {
            NEKERROR(ErrorUtil::efatal, "Only 1D/2D points distributions "
                                        "supported.");
        }
    }
 
    Array<OneD, unsigned int> mapArray(nEdgeCoeffs);
    Array<OneD, int> signArray(nEdgeCoeffs);
 
    for (i = 0; i < kNedges; i++)
    {
        m_edges[i]->FillGeom();
        m_xmap->GetTraceToElementMap(i, mapArray, signArray, m_eorient[i]);
 
        nEdgeCoeffs = m_edges[i]->GetXmap()->GetNcoeffs();
 
        for (j = 0; j < m_coordim; j++)
        {
            for (k = 0; k < nEdgeCoeffs; k++)
            {
                m_coeffs[j][mapArray[k]] =
                    signArray[k] * m_edges[i]->GetCoeffs(j)[k];
            }
        }
    }
 
    m_state = ePtsFilled;
}
 
int TriGeom::v_GetDir(const int i, [[maybe_unused]] const int j) const
{
    return i == 0 ? 0 : 1;
}
 
void TriGeom::v_Reset(CurveMap &curvedEdges, CurveMap &curvedFaces)
{
    Geometry::v_Reset(curvedEdges, curvedFaces);
    CurveMap::iterator it = curvedFaces.find(m_globalID);
 
    if (it != curvedFaces.end())
    {
        m_curve = it->second.get();
    }
 
    for (int i = 0; i < 3; ++i)
    {
        m_edges[i]->Reset(curvedEdges, curvedFaces);
    }
 
    SetUpXmap();
    SetUpCoeffs(m_xmap->GetNcoeffs());
}
 
void TriGeom::v_Setup()
{
    if (!m_setupState)
    {
        for (int i = 0; i < 3; ++i)
        {
            m_edges[i]->Setup();
        }
        SetUpXmap();
        SetUpCoeffs(m_xmap->GetNcoeffs());
 
        // check to see if expansions are linear
        m_straightEdge = 1;
        if (m_xmap->GetBasisNumModes(0) != 2 ||
            m_xmap->GetBasisNumModes(1) != 2)
        {
            m_straightEdge = 0;
        }
 
        m_setupState = true;
    }
}
 
void TriGeom::SetUpXmap()
{
    int order0 = m_edges[0]->GetXmap()->GetBasis(0)->GetNumModes();
    int order1 = std::max(
        order0, std::max(m_edges[1]->GetXmap()->GetBasis(0)->GetNumModes(),
                         m_edges[2]->GetXmap()->GetBasis(0)->GetNumModes()));
 
    std::array<LibUtilities::BasisKey, 2> basis = {
        LibUtilities::BasisKey(
            LibUtilities::eModified_A, order0,
            LibUtilities::PointsKey(order0 + 1,
                                    LibUtilities::eGaussLobattoLegendre)),
        LibUtilities::BasisKey(
            LibUtilities::eModified_B, order1,
            LibUtilities::PointsKey(order1,
                                    LibUtilities::eGaussRadauMAlpha1Beta0))};
 
    m_xmap = GetStdTriFactory().CreateInstance(basis);
}
 
} // namespace Nektar::SpatialDomains