PointGeom.cpp

Go to the documentation of this file.

////////////////////////////////////////////////////////////////////////////////
//
//  File: PointGeom.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: Point geometry information
//
////////////////////////////////////////////////////////////////////////////////
 
#include <fstream>
 
#include <SpatialDomains/PointGeom.h>
#include <SpatialDomains/SegGeom.h>
#include <StdRegions/StdPointExp.h>
#include <StdRegions/StdRegions.hpp>
 
namespace Nektar::SpatialDomains
{
 
PointGeom::PointGeom() : NekPoint<NekDouble>(0.0, 0.0, 0.0)
{
    m_shapeType = LibUtilities::ePoint;
    m_coordim   = 0;
    m_globalID  = 0;
}
 
PointGeom::PointGeom(const int coordim, const int vid, NekDouble x, NekDouble y,
                     NekDouble z)
    : NekPoint<NekDouble>(x, y, z)
{
    m_shapeType = LibUtilities::ePoint;
    m_coordim   = coordim;
    m_globalID  = vid;
}
 
// copy constructor
PointGeom::PointGeom(const PointGeom &T) : Geometry0D(T), NekPoint<NekDouble>(T)
{
    m_shapeType = T.m_shapeType;
    m_globalID  = T.m_globalID;
    m_coordim   = T.m_coordim;
}
 
void PointGeom::GetCoords(NekDouble &x, NekDouble &y, NekDouble &z)
{
    switch (m_coordim)
    {
        case 3:
            z = (*this)(2);
            /* Falls through. */
        case 2:
            y = (*this)(1);
            /* Falls through. */
        case 1:
            x = (*this)(0);
            break;
    }
}
 
void PointGeom::GetCoords(Array<OneD, NekDouble> &coords)
{
    switch (m_coordim)
    {
        case 3:
            coords[2] = (*this)(2);
            /* Falls through. */
        case 2:
            coords[1] = (*this)(1);
            /* Falls through. */
        case 1:
            coords[0] = (*this)(0);
            break;
    }
}
 
void PointGeom::UpdatePosition(NekDouble x, NekDouble y, NekDouble z)
{
    (*this)(0) = x;
    (*this)(1) = y;
    (*this)(2) = z;
}
 
// _this = a + b
void PointGeom::Add(PointGeom &a, PointGeom &b)
{
    (*this)(0) = a[0] + b[0];
    (*this)(1) = a[1] + b[1];
    (*this)(2) = a[2] + b[2];
    m_coordim  = std::max(a.GetCoordim(), b.GetCoordim());
}
 
// _this = a + b
void PointGeom::Sub(PointGeom &a, PointGeom &b)
{
    (*this)(0) = a[0] - b[0];
    (*this)(1) = a[1] - b[1];
    (*this)(2) = a[2] - b[2];
    m_coordim  = std::max(a.GetCoordim(), b.GetCoordim());
}
 
/// \brief _this = a x b
void PointGeom::Mult(PointGeom &a, PointGeom &b)
{
    (*this)(0) = a[1] * b[2] - a[2] * b[1];
    (*this)(1) = a[2] * b[0] - a[0] * b[2];
    (*this)(2) = a[0] * b[1] - a[1] * b[0];
    m_coordim  = 3;
}
 
/// \brief _this  = rotation of a by angle 'angle' around axis dir
void PointGeom::Rotate(PointGeom &a, int dir, NekDouble angle)
{
    switch (dir)
    {
        case 0:
        {
            NekDouble yrot = cos(angle) * a.y() - sin(angle) * a.z();
            NekDouble zrot = sin(angle) * a.y() + cos(angle) * a.z();
 
            (*this)(0) = a.x();
            (*this)(1) = yrot;
            (*this)(2) = zrot;
        }
        break;
        case 1:
        {
            NekDouble zrot = cos(angle) * a.z() - sin(angle) * a.x();
            NekDouble xrot = sin(angle) * a.z() + cos(angle) * a.x();
 
            (*this)(0) = xrot;
            (*this)(1) = a.y();
            (*this)(2) = zrot;
        }
        break;
        case 2:
        {
            NekDouble xrot = cos(angle) * a.x() - sin(angle) * a.y();
            NekDouble yrot = sin(angle) * a.x() + cos(angle) * a.y();
 
            (*this)(0) = xrot;
            (*this)(1) = yrot;
            (*this)(2) = a.z();
        }
        break;
    }
}
 
///  \brief return distance between this and input a
NekDouble PointGeom::dist(PointGeom &a)
{
    return sqrt((x() - a.x()) * (x() - a.x()) + (y() - a.y()) * (y() - a.y()) +
                (z() - a.z()) * (z() - a.z()));
}
 
/// \brief retun the dot product between this and input a
NekDouble PointGeom::dot(PointGeom &a)
{
    return (x() * a.x() + y() * a.y() + z() * a.z());
}
 
/// Determine equivalence by the ids.  No matter what the position,
/// if the ids are the same, then they are equivalent, and vice versa.
bool operator==(const PointGeom &x, const PointGeom &y)
{
    return (x.m_globalID == y.m_globalID);
}
 
bool operator==(const PointGeom &x, const PointGeom *y)
{
    return (x.m_globalID == y->m_globalID);
}
 
bool operator==(const PointGeom *x, const PointGeom &y)
{
    return (x->m_globalID == y.m_globalID);
}
 
bool operator!=(const PointGeom &x, const PointGeom &y)
{
    return (x.m_globalID != y.m_globalID);
}
 
bool operator!=(const PointGeom &x, const PointGeom *y)
{
    return (x.m_globalID != y->m_globalID);
}
 
bool operator!=(const PointGeom *x, const PointGeom &y)
{
    return (x->m_globalID != y.m_globalID);
}
 
void PointGeom::v_GenGeomFactors()
{
    if (!m_setupState)
    {
        m_xmap = MemoryManager<StdRegions::StdPointExp>::AllocateSharedPtr();
        SetUpCoeffs(m_xmap->GetNcoeffs());
        m_setupState = true;
    }
 
    if (m_geomFactorsState != ePtsFilled)
    {
        SpatialDomains::GeomType gType = eRegular;
        m_geomFactors = MemoryManager<GeomFactors>::AllocateSharedPtr(
            gType, m_coordim, m_xmap, m_coeffs);
        m_geomFactorsState = ePtsFilled;
    }
}
 
} // namespace Nektar::SpatialDomains