Commit 028490f7 authored by Michael Turner's avatar Michael Turner
Browse files

pick files from varopti_DNC

parent 051b5baa
......@@ -126,6 +126,7 @@ SET(FoundationHeaders
./Foundations/Interp.h
./Foundations/ManagerAccess.h
./Foundations/NodalPrismEvenlySpaced.h
./Foundations/NodalPrismElec.h
./Foundations/NodalTetElecData.h
./Foundations/NodalTetElec.h
./Foundations/NodalTetEvenlySpaced.h
......@@ -134,6 +135,11 @@ SET(FoundationHeaders
./Foundations/NodalTriEvenlySpaced.h
./Foundations/NodalTriFeketeData.h
./Foundations/NodalTriFekete.h
./Foundations/NodalTriSPIData.h
./Foundations/NodalTriSPI.h
./Foundations/NodalTetSPI.h
./Foundations/NodalPrismSPI.h
./Foundations/NodalTetSPIData.h
./Foundations/NodalUtil.h
./Foundations/PhysGalerkinProject.h
./Foundations/Points.h
......@@ -142,7 +148,7 @@ SET(FoundationHeaders
)
SET(FoundationSources
./Foundations/Basis.cpp
./Foundations/Basis.cpp
./Foundations/BLPoints.cpp
./Foundations/FourierPoints.cpp
./Foundations/FourierSingleModePoints.cpp
......@@ -151,14 +157,18 @@ SET(FoundationSources
./Foundations/InterpCoeff.cpp
./Foundations/Interp.cpp
./Foundations/ManagerAccess.cpp
./Foundations/NodalPrismEvenlySpaced.cpp
./Foundations/NodalPrismEvenlySpaced.cpp
./Foundations/NodalPrismElec.cpp
./Foundations/NodalTetElec.cpp
./Foundations/NodalTetEvenlySpaced.cpp
./Foundations/NodalTetEvenlySpaced.cpp
./Foundations/NodalTriElec.cpp
./Foundations/NodalTriEvenlySpaced.cpp
./Foundations/NodalTriEvenlySpaced.cpp
./Foundations/NodalTriFekete.cpp
./Foundations/NodalTriSPI.cpp
./Foundations/NodalTetSPI.cpp
./Foundations/NodalPrismSPI.cpp
./Foundations/NodalUtil.cpp
./Foundations/PhysGalerkinProject.cpp
./Foundations/PhysGalerkinProject.cpp
./Foundations/Points.cpp
./Foundations/PolyEPoints.cpp
)
......
......@@ -41,10 +41,10 @@
#include <string>
namespace Nektar
{
{
namespace LibUtilities
{
const char* const BasisTypeMap[] =
{
const char* const BasisTypeMap[] =
{
"NoBasisType",
"Ortho_A",
......@@ -66,7 +66,7 @@ namespace Nektar
const std::string kPointsTypeStr[] =
const std::string kPointsTypeStr[] =
{
"NoPointsType",
"GaussGaussLegendre",
......@@ -92,14 +92,17 @@ namespace Nektar
"BoundaryLayerPointsRev",
"NodalTriElec",
"NodalTriFekete",
"NodalTriSPI",
"NodalTriEvenlySpaced",
"NodalTetEvenlySpaced",
"NodalTetElec",
"NodalPrismEvenlySpaced"
"NodalTetSPI",
"NodalPrismEvenlySpaced",
"NodalPrismElec",
"NodalPrismSPI"
};
} // end of namespace
} // end of namespace
#endif //FOUNDATIONS_H
......@@ -59,6 +59,7 @@ namespace Nektar
class NodalTetEvenlySpaced;
class NodalTetElec;
class NodalPrismEvenlySpaced;
class NodalPrismElec;
class NodalTriElec;
class NodalTriEvenlySpaced;
class NodalTriFekete;
......
......@@ -28,8 +28,8 @@
// 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: C functions to provide access to managers.
//
// Description: C functions to provide access to managers.
//
///////////////////////////////////////////////////////////////////////////////
......@@ -42,16 +42,20 @@
#include <LibUtilities/Foundations/NodalTriElec.h>
#include <LibUtilities/Foundations/NodalTetElec.h>
#include <LibUtilities/Foundations/NodalTriFekete.h>
#include <LibUtilities/Foundations/NodalTriSPI.h>
#include <LibUtilities/Foundations/NodalTriEvenlySpaced.h>
#include <LibUtilities/Foundations/NodalTetEvenlySpaced.h>
#include <LibUtilities/Foundations/NodalTetSPI.h>
#include <LibUtilities/Foundations/NodalPrismSPI.h>
#include <LibUtilities/Foundations/NodalPrismEvenlySpaced.h>
#include <LibUtilities/Foundations/NodalPrismElec.h>
#include <LibUtilities/Foundations/Basis.h>
#include <LibUtilities/Foundations/Foundations.hpp>
#include <LibUtilities/Foundations/ManagerAccess.h>
namespace Nektar
{
namespace LibUtilities
namespace LibUtilities
{
// Register all points and basis creators.
namespace
......@@ -68,12 +72,12 @@ namespace Nektar
const bool gaussInited10 = PointsManager().RegisterCreator(PointsKey(0, eGaussRadauMAlpha0Beta2), GaussPoints::Create);
const bool gaussInited11 = PointsManager().RegisterCreator(PointsKey(0, eGaussRadauMAlpha1Beta0), GaussPoints::Create);
const bool gaussInited12 = PointsManager().RegisterCreator(PointsKey(0, eGaussRadauMAlpha2Beta0), GaussPoints::Create);
const bool gaussInited13 = PointsManager().RegisterCreator(PointsKey(0, eGaussKronrodLegendre), GaussPoints::Create);
const bool gaussInited14 = PointsManager().RegisterCreator(PointsKey(0, eGaussRadauKronrodMLegendre), GaussPoints::Create);
const bool gaussInited15 = PointsManager().RegisterCreator(PointsKey(0, eGaussRadauKronrodMAlpha1Beta0), GaussPoints::Create);
const bool gaussInited16 = PointsManager().RegisterCreator(PointsKey(0, eGaussLobattoKronrodLegendre), GaussPoints::Create);
const bool fourierInited0 = PointsManager().RegisterCreator(PointsKey(0, eFourierEvenlySpaced), FourierPoints::Create);
const bool fourierInitedSM0 = PointsManager().RegisterCreator(PointsKey(0, eFourierSingleModeSpaced), FourierSingleModePoints::Create);
const bool BLInited0 = PointsManager().RegisterCreator(PointsKey(0, eBoundaryLayerPoints), BLPoints::Create);
......@@ -81,10 +85,14 @@ namespace Nektar
const bool polyeInited10 = PointsManager().RegisterCreator(PointsKey(0, ePolyEvenlySpaced), PolyEPoints::Create);
const bool NodalTriInited0 = PointsManager().RegisterCreator(PointsKey(0, eNodalTriElec), NodalTriElec::Create);
const bool NodalTriInited1 = PointsManager().RegisterCreator(PointsKey(0, eNodalTriFekete), NodalTriFekete::Create);
const bool NodalTriInited2 = PointsManager().RegisterCreator(PointsKey(0, eNodalTriSPI), NodalTriSPI::Create);
const bool nodalTetElecInited = PointsManager().RegisterCreator(PointsKey(0, eNodalTetElec), NodalTetElec::Create);
const bool nodalTetElecInited1 = PointsManager().RegisterCreator(PointsKey(0, eNodalTetSPI), NodalTetSPI::Create);
const bool NodalTriEveInited = PointsManager().RegisterCreator(PointsKey(0, eNodalTriEvenlySpaced), NodalTriEvenlySpaced::Create);
const bool NodalTetEveInited = PointsManager().RegisterCreator(PointsKey(0, eNodalTetEvenlySpaced), NodalTetEvenlySpaced::Create);
const bool NodalPrismEveInited = PointsManager().RegisterCreator(PointsKey(0, eNodalPrismEvenlySpaced), NodalPrismEvenlySpaced::Create);
const bool NodalPrismElecInited = PointsManager().RegisterCreator(PointsKey(0, eNodalPrismElec), NodalPrismElec::Create);
const bool NodalPrismSPIInited = PointsManager().RegisterCreator(PointsKey(0, eNodalPrismSPI), NodalPrismSPI::Create);
const bool Basis_Inited = BasisManager().RegisterGlobalCreator(Basis::Create);
};
......@@ -101,5 +109,3 @@ namespace Nektar
} // end of namespace LibUtilities
} // end of namespace Nektar
///////////////////////////////////////////////////////////////////////////////
//
// File NodalPrismEleccpp
//
// 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).
//
// License for the specific language governing rights and limitations under
// 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: 3D Nodal Prism Evenly Spaced Point Definitions
//
///////////////////////////////////////////////////////////////////////////////
#include <LibUtilities/Foundations/NodalPrismElec.h>
#include <LibUtilities/Foundations/Points.h>
#include <LibUtilities/BasicUtils/ErrorUtil.hpp>
#include <LibUtilities/BasicUtils/SharedArray.hpp>
#include <vector>
namespace Nektar
{
namespace LibUtilities
{
namespace
{
bool isVertex(int t, int y, int npts)
{
return (t == 0 && y == 0 ) ||
(t == 1 && y == 0 ) ||
(t == 2 && y == 0 ) ||
(t == 0 && y == (npts-1)) ||
(t == 1 && y == (npts-1)) ||
(t == 2 && y == (npts-1));
}
bool isEdge_01(int t, int y, int npts)
{
return y == 0 && t > 2 && t <=npts;
}
bool isEdge_12(int t, int y, int npts)
{
return t == 1;
}
bool isEdge_23(int t, int y, int npts)
{
return y == (npts-1) && t > 2 && t <=npts;
}
bool isEdge_30(int t, int y, int npts)
{
return t == 0;
}
bool isEdge_04(int t, int y, int npts)
{
return y == 0 && t >= 3+2*(npts-2) && t < 3+3*(npts-2);
}
bool isEdge_14(int t, int y, int npts)
{
return y == 0 && t >= 3+(npts-2) && t < 3+2*(npts-2);
}
bool isEdge_25(int t, int y, int npts)
{
return y == npts-1 && t >= 3+(npts-2) && t < 3+2*(npts-2);
}
bool isEdge_35(int t, int y, int npts)
{
return y == npts-1 && t >= 3+2*(npts-2) && t < 3+3*(npts-2);
}
bool isEdge_45(int t, int y, int npts)
{
return t == 2;
}
bool isEdge(int t, int y, int npts){
return isEdge_01(t,y,npts) || isEdge_12(t,y,npts) ||
isEdge_23(t,y,npts) || isEdge_30(t,y,npts) ||
isEdge_04(t,y,npts) || isEdge_14(t,y,npts) ||
isEdge_25(t,y,npts) || isEdge_35(t,y,npts) ||
isEdge_45(t,y,npts);
}
bool isFace_0123(int t, int y, int npts)
{
return t < 3 + (npts-2);
}
bool isFace_014(int t, int y, int npts)
{
return y == 0;
}
bool isFace_1254(int t, int y, int npts)
{
return t < 3 + 2*(npts-2) && t >= 3 + (npts-2);
}
bool isFace_325(int t, int y, int npts)
{
return y == (npts-1);
}
bool isFace_0354(int t, int y, int npts)
{
return t < 3 + 3*(npts-2) && t >= 3 + 2*(npts-2);
}
bool isFace(int t, int y, int npts){
return isFace_0123(t,y,npts) || isFace_014(t,y,npts) ||
isFace_1254(t,y,npts) || isFace_325(t,y,npts) ||
isFace_0354(t,y,npts);
}
}
// Calculate evenly spaced number of points
void NodalPrismElec::CalculatePoints()
{
// Allocate the storage for points
PointsBaseType::CalculatePoints();
// Populate m_points
unsigned int npts = GetNumPoints();
LibUtilities::PointsKey pkey1(npts, LibUtilities::eNodalTriElec);
Array<OneD, NekDouble> u1, v1;
LibUtilities::PointsManager()[pkey1]->GetPoints(u1, v1);
LibUtilities::PointsKey pkey2(npts, LibUtilities::eGaussLobattoLegendre);
Array<OneD, NekDouble> u;
LibUtilities::PointsManager()[pkey2]->GetPoints(u);
for(int y = 0, index = 0; y < npts; y++)
{
for(int t = 0; t < u1.num_elements(); t++, index++)
{
m_points[0][index] = u1[t];
m_points[1][index] = u[y];
m_points[2][index] = v1[t];
}
}
NodalPointReorder3d();
m_util = MemoryManager<NodalUtilPrism>::AllocateSharedPtr(
npts - 1, m_points[0], m_points[1], m_points[2]);
}
void NodalPrismElec::NodalPointReorder3d()
{
unsigned int npts = GetNumPoints();
using std::vector;
vector<int> vertex;
vector<int> iEdge_01; // interior edge 0
vector<int> iEdge_12; // interior edge 1
vector<int> iEdge_23; // interior edge 2
vector<int> iEdge_30; // interior edge 3
vector<int> iEdge_04; // interior edge 4
vector<int> iEdge_14; // interior edge 5
vector<int> iEdge_25; // interior edge 6
vector<int> iEdge_35; // interior edge 7
vector<int> iEdge_45; // interior edge 8
vector<int> iFace_0123; // interior face 0
vector<int> iFace_014; // interior face 1
vector<int> iFace_1254; // interior face 2
vector<int> iFace_325; // interior face 3
vector<int> iFace_0354; // interior face 4
vector<int> interiorVolumePoints; // interior volume points
vector<int> map;
// Build the lattice prism left to right - bottom to top
for(int y = 0, index = 0; y < npts; y++)
{
for(int t = 0; t < npts*(npts+1)/2; t++, index++)
{
if (isVertex(t,y,npts))
{
vertex.push_back(index);
}
else if (isEdge(t,y,npts))
{
if (isEdge_01(t,y,npts))
{
iEdge_01.push_back(index);
}
else if (isEdge_12(t,y,npts))
{
iEdge_12.push_back(index);
}
else if (isEdge_23(t,y,npts))
{
iEdge_23.push_back(index);
}
else if (isEdge_30(t,y,npts))
{
iEdge_30.push_back(index);
}
else if (isEdge_04(t,y,npts))
{
iEdge_04.push_back(index);
}
else if (isEdge_14(t,y,npts))
{
iEdge_14.push_back(index);
}
else if (isEdge_25(t,y,npts))
{
iEdge_25.push_back(index);
}
else if (isEdge_35(t,y,npts))
{
iEdge_35.push_back(index);
}
else if (isEdge_45(t,y,npts))
{
iEdge_45.push_back(index);
}
}
else if (isFace(t,y,npts))
{
if (isFace_0123(t,y,npts))
{
iFace_0123.push_back(index);
}
else if (isFace_014(t,y,npts))
{
iFace_014.push_back(index);
}
else if (isFace_1254(t,y,npts))
{
iFace_1254.push_back(index);
}
else if (isFace_325(t,y,npts))
{
iFace_325.push_back(index);
}
else if (isFace_0354(t,y,npts))
{
iFace_0354.push_back(index);
}
}
else
{
interiorVolumePoints.push_back(index);
}
}
}
//sort vertices
swap(vertex[2],vertex[4]);
//sort edges
reverse(iEdge_23.begin(),iEdge_23.end());
reverse(iEdge_30.begin(),iEdge_30.end());
reverse(iEdge_04.begin(),iEdge_04.end());
reverse(iEdge_35.begin(),iEdge_35.end());
//faces
for(int i = 0; i < npts-2; i++)
{
for(int j = i+1; j < npts-2; j++)
{
swap(iFace_1254[i*(npts-2)+j],iFace_1254[j*(npts-2)+i]);
}
}
for(int i = 0; i < npts-2; i++)
{
reverse(iFace_0354.begin()+i*(npts-2), iFace_0354.begin()+i*(npts-2)+npts-2);
}
for(int i = 0; i < npts-2; i++)
{
for(int j = i+1; j < npts-2; j++)
{
swap(iFace_0354[i*(npts-2)+j],iFace_0354[j*(npts-2)+i]);
}
}
for (unsigned int n=0; n<vertex.size(); ++n)
{
map.push_back(vertex[n]);
}
for (unsigned int n=0; n<iEdge_01.size(); ++n)
{
map.push_back(iEdge_01[n]);
}
for (unsigned int n=0; n<iEdge_12.size(); ++n)
{
map.push_back(iEdge_12[n]);
}
for (unsigned int n=0; n<iEdge_23.size(); ++n)
{
map.push_back(iEdge_23[n]);
}
for (unsigned int n=0; n<iEdge_30.size(); ++n)
{
map.push_back(iEdge_30[n]);
}
for (unsigned int n=0; n<iEdge_04.size(); ++n)
{
map.push_back(iEdge_04[n]);
}
for (unsigned int n=0; n<iEdge_14.size(); ++n)
{
map.push_back(iEdge_14[n]);
}
for (unsigned int n=0; n<iEdge_25.size(); ++n)
{
map.push_back(iEdge_25[n]);
}
for (unsigned int n=0; n<iEdge_35.size(); ++n)
{
map.push_back(iEdge_35[n]);
}
for (unsigned int n=0; n<iEdge_45.size(); ++n)
{
map.push_back(iEdge_45[n]);
}
for (unsigned int n=0; n<iFace_0123.size(); ++n)
{
map.push_back(iFace_0123[n]);
}
for (unsigned int n=0; n<iFace_014.size(); ++n)
{
map.push_back(iFace_014[n]);
}
for(unsigned int n=0; n<iFace_1254.size(); ++n)
{
map.push_back(iFace_1254[n]);
}
for(unsigned int n=0; n<iFace_325.size(); ++n)
{
map.push_back(iFace_325[n]);
}
for(unsigned int n=0; n<iFace_0354.size(); ++n)
{
map.push_back(iFace_0354[n]);
}
for(unsigned int n=0; n<interiorVolumePoints.size(); ++n)
{
map.push_back(interiorVolumePoints[n]);
}
Array<OneD, NekDouble> points[3];
points[0] = Array<OneD, NekDouble>(GetTotNumPoints());
points[1] = Array<OneD, NekDouble>(GetTotNumPoints());
points[2] = Array<OneD, NekDouble>(GetTotNumPoints());
for(unsigned int index=0; index<map.size(); ++index)
{
points[0][index] = m_points[0][index];
points[1][index] = m_points[1][index];
points[2][index] = m_points[2][index];
}
for(unsigned int index=0; index<map.size(); ++index)
{
m_points[0][index] = points[0][map[index]];
m_points[1][index] = points[1][map[index]];
m_points[2][index] = points[2][map[index]];
}
}
void NodalPrismElec::CalculateWeights()
{
// Allocate the storage for points
PointsBaseType::CalculateWeights();
typedef DataType T;
// Solve the Vandermonde system of integrals for the weight vector
NekVector<T> w = m_util->GetWeights();
m_weights = Array<OneD,T>( w.GetRows(), w.GetPtr() );
}
// ////////////////////////////////////////
// CalculateInterpMatrix()
void NodalPrismElec::CalculateInterpMatrix(const Array<OneD, const NekDouble>& xia,
const Array<OneD, const NekDouble>& yia,
const Array<OneD, const NekDouble>& zia,