Skip to content

GitLab

Nektar
Nektar
Commit 372b8ffd authored by Dave Moxey's avatar Dave Moxey
Browse files
Options

Merge branch 'master' into fix/fieldconvert 

Conflicts:
	library/MultiRegions/DisContField2D.cpp
	library/MultiRegions/DisContField3D.cpp
parents 7d893399 2c7169fd
Hide whitespace changes
Inline Side-by-side
Showing with 27476 additions and 6470 deletions
library/Demos/LocalRegions/CMakeLists.txt
View file @ 372b8ffd
......@@ -14,7 +14,6 @@ ADD_NEKTAR_EXECUTABLE(LocProject1D demos LocProject1DSource)
TARGET_LINK_LIBRARIES(LocProject1D ${LinkLibraries})
SET_LAPACK_LINK_LIBRARIES(LocProject1D)
SET(LocProject2DSource LocProject2D.cpp)
ADD_NEKTAR_EXECUTABLE(LocProject2D demos LocProject2DSource)
TARGET_LINK_LIBRARIES(LocProject2D ${LinkLibraries})
......@@ -25,7 +24,6 @@ ADD_NEKTAR_EXECUTABLE(LocProject3D demos LocProject3DSource)
TARGET_LINK_LIBRARIES(LocProject3D ${LinkLibraries})
SET_LAPACK_LINK_LIBRARIES(LocProject3D)
SET(LocProject_Diff1DSource LocProject_Diff1D.cpp)
ADD_NEKTAR_EXECUTABLE(LocProject_Diff1D demos LocProject_Diff1DSource)
TARGET_LINK_LIBRARIES(LocProject_Diff1D ${LinkLibraries})
......@@ -41,46 +39,6 @@ ADD_NEKTAR_EXECUTABLE(LocProject_Diff3D demos LocProject_Diff3DSource)
TARGET_LINK_LIBRARIES(LocProject_Diff3D ${LinkLibraries})
SET_LAPACK_LINK_LIBRARIES(LocProject_Diff3D)
SET(LocTetExpDemoSource LocTetExpDemo.cpp)
ADD_NEKTAR_EXECUTABLE(LocTetExpDemo demos LocTetExpDemoSource)
TARGET_LINK_LIBRARIES(LocTetExpDemo ${LinkLibraries})
SET_LAPACK_LINK_LIBRARIES(LocTetExpDemo)
SET(LocPyrExpDemoSource LocPyrExpDemo.cpp)
ADD_NEKTAR_EXECUTABLE(LocPyrExpDemo demos LocPyrExpDemoSource)
TARGET_LINK_LIBRARIES(LocPyrExpDemo ${LinkLibraries})
SET_LAPACK_LINK_LIBRARIES(LocPyrExpDemo)
SET(LocPrismExpDemoSource LocPrismExpDemo.cpp)
ADD_NEKTAR_EXECUTABLE(LocPrismExpDemo demos LocPrismExpDemoSource)
TARGET_LINK_LIBRARIES(LocPrismExpDemo ${LinkLibraries})
SET_LAPACK_LINK_LIBRARIES(LocPrismExpDemo)
SET(LocHexExpDemoSource LocHexExpDemo.cpp)
ADD_NEKTAR_EXECUTABLE(LocHexExpDemo demos LocHexExpDemoSource)
TARGET_LINK_LIBRARIES(LocHexExpDemo ${LinkLibraries})
SET_LAPACK_LINK_LIBRARIES(LocHexExpDemo)
SET(LocHexExpDemoSource LocHexExpDemo_2.cpp)
ADD_NEKTAR_EXECUTABLE(LocHexExpDemo_2 demos LocHexExpDemoSource)
TARGET_LINK_LIBRARIES(LocHexExpDemo_2 ${LinkLibraries})
SET_LAPACK_LINK_LIBRARIES(LocHexExpDemo_2)
SET(LocPrismExpDemoSource LocPrismExpDemo_2.cpp)
ADD_NEKTAR_EXECUTABLE(LocPrismExpDemo_2 demos LocPrismExpDemoSource)
TARGET_LINK_LIBRARIES(LocPrismExpDemo_2 ${LinkLibraries})
SET_LAPACK_LINK_LIBRARIES(LocPrismExpDemo_2)
SET(LocPyrExpDemoSource LocPyrExpDemo_2.cpp)
ADD_NEKTAR_EXECUTABLE(LocPyrExpDemo_2 demos LocPyrExpDemoSource)
TARGET_LINK_LIBRARIES(LocPyrExpDemo_2 ${LinkLibraries})
SET_LAPACK_LINK_LIBRARIES(LocPyrExpDemo_2)
SET(LocTetExpDemoSource LocTetExpDemo_2.cpp)
ADD_NEKTAR_EXECUTABLE(LocTetExpDemo_2 demos LocTetExpDemoSource)
TARGET_LINK_LIBRARIES(LocTetExpDemo_2 ${LinkLibraries})
SET_LAPACK_LINK_LIBRARIES(LocTetExpDemo_2)
# Generate list of available subdirectories
FILE(GLOB dir_list "*")
FOREACH(dir ${dir_list})
......@@ -108,6 +66,7 @@ ADD_NEKTAR_TEST(LocProject3D_Tet_Ortho_Basis_P6_Q7)
ADD_NEKTAR_TEST(LocProject3D_Tet_Mod_Basis_P6_Q7)
ADD_NEKTAR_TEST(LocProject3D_Prism_Ortho_Basis_P6_Q7)
ADD_NEKTAR_TEST(LocProject3D_Prism_Mod_Basis_P6_Q7)
ADD_NEKTAR_TEST(LocProject3D_Pyr_Mod_Basis_P6_Q7)
ADD_NEKTAR_TEST(LocProject3D_Hex_Ortho_Basis_P6_Q7)
ADD_NEKTAR_TEST(LocProject3D_Hex_Mod_Basis_P6_Q7)
ADD_NEKTAR_TEST(LocProject3D_Hex_Lagrange_Basis_P6_Q7)
......@@ -123,6 +82,7 @@ ADD_NEKTAR_TEST(LocProject_Diff3D_Reg_Hex_Mod_Basis_P6_Q7)
ADD_NEKTAR_TEST(LocProject_Diff3D_Reg_Hex_Ortho_Basis_P6_Q7)
ADD_NEKTAR_TEST(LocProject_Diff3D_Reg_Prism_Mod_Basis_P6_Q7)
ADD_NEKTAR_TEST(LocProject_Diff3D_Reg_Prism_Ortho_Basis_P6_Q7)
ADD_NEKTAR_TEST(LocProject_Diff3D_Reg_Pyr_Mod_Basis_P6_Q7)
ADD_NEKTAR_TEST(LocProject_Diff3D_Reg_Tet_Mod_Basis_P6_Q7)
ADD_NEKTAR_TEST(LocProject_Diff3D_Reg_Tet_Ortho_Basis_P6_Q7)
......
library/Demos/LocalRegions/LocHexExpDemo.cpp deleted 100644 → 0
View file @ 7d893399
#include <StdRegions/StdExpansion3D.h>
#include <LocalRegions/HexExp.h>
#include <LibUtilities/Foundations/Foundations.hpp>
#include <LibUtilities/Foundations/Basis.h>
#include <SpatialDomains/MeshComponents.h>
#include <LibUtilities/BasicConst/NektarUnivTypeDefs.hpp>
#include <algorithm>
#include <iostream>
#include <cstdlib>
#include <cmath>
#include <iomanip>
#include <iosfwd>
using namespace std;
using namespace boost;
using namespace Nektar;
using namespace Nektar::LibUtilities;
using namespace Nektar::LocalRegions;
using namespace Nektar::StdRegions;
using namespace Nektar::SpatialDomains;
/// Defines the test solution which excites all modes in the Hex.
NekDouble Hex_sol(NekDouble x, NekDouble y, NekDouble z,
int order1, int order2, int order3,
LibUtilities::BasisType btype1,
LibUtilities::BasisType btype2,
LibUtilities::BasisType btype3);
/// Specialised pow function returning 1.0 for exponents < 0.
static double pow_loc(const double val, const int i)
{
return (i < 0)? 1.0: pow(val,i);
}
int main(int argc, char *argv[])
{
if( argc != 10 ) {
cerr << "Usage: HexDemo Type_x Type_y Type_z numModes_x numModes_y "
"numModes_z Qx Qy Qz" << endl;
cerr << "Where type is an interger value which dictates the basis as:"
<< endl;
cerr << "\t Ortho_A = 1\n";
cerr << "\t Modified_A = 4\n";
cerr << "\t Fourier = 7\n";
cerr << "\t Lagrange = 8\n";
cerr << "\t Legendre = 9\n";
cerr << "\t Chebyshev = 10\n";
cerr << "\n\n" << "Example: " << argv[0] << " 4 4 4 3 3 3 5 5 5"
<< endl;
cerr << endl;
exit(1);
}
// Set up the region shape and expansion types
int bType_x_val = atoi(argv[1]);
int bType_y_val = atoi(argv[2]);
int bType_z_val = atoi(argv[3]);
BasisType bType_x = static_cast<BasisType>( bType_x_val );
BasisType bType_y = static_cast<BasisType>( bType_y_val );
BasisType bType_z = static_cast<BasisType>( bType_z_val );
if( (bType_x_val == 13) || (bType_y_val == 13) || (bType_z_val == 13) )
{
bType_x = LibUtilities::eOrtho_A;
bType_y = LibUtilities::eOrtho_B;
bType_z = LibUtilities::eOrtho_C;
}
if( (bType_x == eOrtho_B) || (bType_x == eModified_B) ) {
NEKERROR(ErrorUtil::efatal,
"Basis 1 cannot be of type Ortho_B or Modified_B");
}
if( (bType_x == eOrtho_C) || (bType_x == eModified_C) ) {
NEKERROR(ErrorUtil::efatal,
"Basis 1 cannot be of type Ortho_C or Modified_C");
}
if( (bType_y == eOrtho_C) || (bType_y == eModified_C) ) {
NEKERROR(ErrorUtil::efatal,
"Basis 2 cannot be of type Ortho_C or Modified_C");
}
// Set up the number of quadrature points, order of bases, etc
int xModes = atoi(argv[4]);
int yModes = atoi(argv[5]);
int zModes = atoi(argv[6]);
int Qx = atoi(argv[7]);
int Qy = atoi(argv[8]);
int Qz = atoi(argv[9]);
int P = xModes - 1;
int Q = yModes - 1;
int R = zModes - 1;
const int three = 3;
Array<OneD, NekDouble> solution( Qx * Qy * Qz );
LibUtilities::PointsType Qtype_x, Qtype_y, Qtype_z;
Array<OneD,NekDouble> x = Array<OneD,NekDouble>( Qx * Qy * Qz );
Array<OneD,NekDouble> y = Array<OneD,NekDouble>( Qx * Qy * Qz );
Array<OneD,NekDouble> z = Array<OneD,NekDouble>( Qx * Qy * Qz );
if(bType_x != LibUtilities::eFourier)
{
Qtype_x = LibUtilities::eGaussLobattoLegendre;
}
else
{
Qtype_x = LibUtilities::eFourierEvenlySpaced;
}
if(bType_y != LibUtilities::eFourier)
{
Qtype_y = LibUtilities::eGaussLobattoLegendre;
}
else
{
Qtype_y = LibUtilities::eFourierEvenlySpaced;
}
if(bType_z != LibUtilities::eFourier)
{
Qtype_z = LibUtilities::eGaussLobattoLegendre;
}
else
{
Qtype_z = LibUtilities::eFourierEvenlySpaced;
}
//-----------------------------------------------
// Define a 3D expansion based on basis definition
StdRegions::StdExpansion3D *lhe = 0;
// ////////////////////////////////////////////////////////////////
// Set up Hexahedron vertex coordinates
// PointGeom (const int coordim, const int vid, double x,
// double y, double z)
const int nVerts = 8;
const double point[][3] = {
{0,0,0}, {1,0,0}, {1,1,0}, {0,1,0},
{0,0,1}, {1,0,1}, {1,1,1}, {0,1,1}
};
// Populate the list of verts
PointGeomSharedPtr verts[8];
for( int i = 0; i < nVerts; ++i ) {
verts[i] = MemoryManager<PointGeom>
::AllocateSharedPtr(three, i, point[i][0],
point[i][1], point[i][2]);
}
// ////////////////////////////////////////////////////////////////
// Set up Hexahedron Edges
// SegGeom (int id, const int coordim), EdgeComponent(id, coordim)
const int nEdges = 12;
const int vertexConnectivity[][2] = {
{0,1}, {1,2}, {2,3}, {0,3}, {0,4}, {1,5},
{2,6}, {3,7}, {4,5}, {5,6}, {6,7}, {4,7}
};
// Populate the list of edges
SegGeomSharedPtr edges[nEdges];
for( int i = 0; i < nEdges; ++i ) {
PointGeomSharedPtr vertsArray[2];
for( int j = 0; j < 2; ++j ) {
vertsArray[j] = verts[vertexConnectivity[i][j]];
}
edges[i] = MemoryManager<SegGeom>::
AllocateSharedPtr( i, three, vertsArray);
}
// ////////////////////////////////////////////////////////////////
// Set up Hexahedron faces
const int nFaces = 6;
const int edgeConnectivity[][4] = {
{0,1,2,3}, {0,5,8,4}, {1,6,9,5},
{2,7,10,6}, {3,7,11,4}, {8,9,10,11}
};
const bool isEdgeFlipped[][4] = {
{0,0,0,1}, {0,0,1,1}, {0,0,1,1},
{0,0,1,1}, {0,0,1,1}, {0,0,0,1}
};
// Populate the list of faces
QuadGeomSharedPtr faces[nFaces];
for( int i = 0; i < nFaces; ++i ) {
SegGeomSharedPtr edgeArray[4];
Orientation eorientArray[4];
for( int j = 0; j < 4; ++j ) {
edgeArray[j] = edges[edgeConnectivity[i][j]];
eorientArray[j] = isEdgeFlipped[i][j] ? eBackwards : eForwards;
}
faces[i] = MemoryManager<QuadGeom>::AllocateSharedPtr(i, edgeArray,
eorientArray);
}
const LibUtilities::PointsKey pkey1( Qx, Qtype_x );
const LibUtilities::PointsKey pkey2( Qy, Qtype_y );
const LibUtilities::PointsKey pkey3( Qz, Qtype_z );
const LibUtilities::BasisKey bkey1( bType_x, xModes, pkey1 );
const LibUtilities::BasisKey bkey2( bType_y, yModes, pkey2 );
const LibUtilities::BasisKey bkey3( bType_z, zModes, pkey3 );
Array<OneD, StdRegions::StdExpansion3DSharedPtr> xMap(3);
for(int i = 0; i < 3; ++i) {
xMap[i] = MemoryManager<StdRegions::StdHexExp>
::AllocateSharedPtr(bkey1, bkey2, bkey3);
}
SpatialDomains::HexGeomSharedPtr geom
= MemoryManager<SpatialDomains::HexGeom>::AllocateSharedPtr(faces);
geom->SetOwnData();
if( bType_x_val < 10 ) {
lhe = new LocalRegions::HexExp( bkey1, bkey2, bkey3, geom );
}
lhe->GetCoords(x,y,z);
//----------------------------------------------
// Define solution to be projected
for(int n = 0; n < Qx * Qy * Qz; ++n) {
solution[n] = Hex_sol( x[n], y[n], z[n], P, Q, R,
bType_x, bType_y, bType_z );
}
//----------------------------------------------
//---------------------------------------------
// Project onto Expansion
Array<OneD, NekDouble> phys (Qx * Qy * Qz);
Array<OneD, NekDouble> coeffs(lhe->GetNcoeffs());
lhe->FwdTrans( solution, coeffs );
//---------------------------------------------
//-------------------------------------------
// Backward Transform Solution to get projected values
lhe->BwdTrans( coeffs, phys );
//-------------------------------------------
//--------------------------------------------
// Calculate L_p error
cout << "L infinity error: " << lhe->Linf(phys, solution) << endl;
cout << "L 2 error: " << lhe->L2 (phys, solution) << endl;
//--------------------------------------------
//-------------------------------------------
// Evaulate solution at x = y = z = 0 and print error
Array<OneD, NekDouble> t = Array<OneD, NekDouble>(3);
t[0] = 0.5;
t[1] = 0.5;
t[2] = 0.5;
NekDouble numericSolution = lhe->PhysEvaluate(t, phys);
solution[0] = Hex_sol( t[0], t[1], t[2], P, Q, R,
bType_x, bType_y, bType_z );
cout << "Solution = " << solution[0] << endl;
cout << "Numeric = " << numericSolution << endl;
cout << "Interpolation difference from actual solution at x = ( "
<< t[0] << ", " << t[1] << ", " << t[2] << " ): "
<< numericSolution - solution[0] << endl;
//-------------------------------------------
return 0;
}
NekDouble Hex_sol(NekDouble x, NekDouble y, NekDouble z,
int order1, int order2, int order3,
LibUtilities::BasisType btype1,
LibUtilities::BasisType btype2,
LibUtilities::BasisType btype3 ) {
int i,j,k;
NekDouble sol = 0.0;
int Nx = (btype1 == LibUtilities::eFourier ? order1/2 : order1);
int Ny = (btype2 == LibUtilities::eFourier ? order2/2 : order2);
int Nz = (btype3 == LibUtilities::eFourier ? order3/2 : order3);
bool Fx = (btype1 == LibUtilities::eFourier);
bool Fy = (btype2 == LibUtilities::eFourier);
bool Fz = (btype3 == LibUtilities::eFourier);
NekDouble a;
for (i = 0; i < Nx; ++i)
{
for (j = 0; j < Ny; ++j)
{
for (k = 0; k < Nz; ++k)
{
a = (Fx ? sin(M_PI*i*x) + cos(M_PI*i*x) : pow_loc(x,i));
a *= (Fy ? sin(M_PI*j*y) + cos(M_PI*j*y) : pow_loc(y,j));
a *= (Fz ? sin(M_PI*k*z) + cos(M_PI*k*z) : pow_loc(z,k));
sol += a;
}
}
}
return sol;
}
library/Demos/LocalRegions/LocHexExpDemo_2.cpp deleted 100644 → 0
View file @ 7d893399
#include <cstdlib>
#include <LocalRegions/HexExp.h>
#include <SpatialDomains/MeshGraph3D.h>
#include <LibUtilities/BasicConst/NektarUnivTypeDefs.hpp>
using namespace Nektar;
/// modification to deal with exact solution. Return 1 if integer < 0
static double pow_loc(const double val, const int i)
{
return (i < 0)? 1.0: pow(val,i);
}
/// Solution to excite all modes in Hex expansion.
NekDouble Hex_sol( NekDouble x, NekDouble y, NekDouble z,
int order1, int order2, int order3,
LibUtilities::BasisType btype1,
LibUtilities::BasisType btype2,
LibUtilities::BasisType btype3);
/// Derivative of solution in Hex expansion.
NekDouble Hex_Diff_Sol( NekDouble x, NekDouble y, NekDouble z,
int P, int Q, int R,
LibUtilities::BasisType btype1,
LibUtilities::BasisType btype2,
LibUtilities::BasisType btype3, int dir);
int main(int argc, char *argv[])
{
if(argc != 8)
{
cerr << "usage: LocHexExpDemo MeshFile nummodes0 nummodes1 nummodes2 "
"nx ny nz" << endl;
exit(1);
}
int i;
string in(argv[1]);
int nummodes0 = atoi(argv[2]);
int nummodes1 = atoi(argv[3]);
int nummodes2 = atoi(argv[4]);
int nquad0 = atoi(argv[5]);
int nquad1 = atoi(argv[6]);
int nquad2 = atoi(argv[7]);
int P = nummodes0 - 1;
int Q = nummodes1 - 1;
int R = nummodes2 - 1;
int ntotquad = nquad0*nquad1*nquad2;
SpatialDomains::MeshGraph3D graph3D;
graph3D.ReadGeometry(in);
LibUtilities::BasisType bType = LibUtilities::eModified_A;
LibUtilities::PointsType qtype = LibUtilities::eGaussLobattoLegendre;
const LibUtilities::PointsKey pointsKey0( nquad0, qtype );
const LibUtilities::PointsKey pointsKey1( nquad1, qtype );
const LibUtilities::PointsKey pointsKey2( nquad2, qtype );
const LibUtilities::BasisKey basisKey0( bType, nummodes0, pointsKey0 );
const LibUtilities::BasisKey basisKey1( bType, nummodes1, pointsKey1 );
const LibUtilities::BasisKey basisKey2( bType, nummodes2, pointsKey2 );
SpatialDomains::HexGeomSharedPtr geom;
if(!(geom = boost::dynamic_pointer_cast<SpatialDomains::HexGeom>(
graph3D.GetCompositeItem(0,0))))
{
cerr << "Could not find HexGeom in input file" << endl;
exit(1);
}
LocalRegions::HexExpSharedPtr E = MemoryManager<LocalRegions::HexExp>
::AllocateSharedPtr(basisKey0,basisKey1,basisKey2,geom);
Array<OneD, NekDouble> sol( ntotquad );
Array<OneD, NekDouble> dx( ntotquad, 0.0 );
Array<OneD, NekDouble> dy( ntotquad, 0.0 );
Array<OneD, NekDouble> dz( ntotquad, 0.0 );
Array<OneD, NekDouble> diff_solution_x( ntotquad, 0.0 );
Array<OneD, NekDouble> diff_solution_y( ntotquad, 0.0 );
Array<OneD, NekDouble> diff_solution_z( ntotquad, 0.0 );
Array<OneD, NekDouble> derivatives( ntotquad, 0.0 );
Array<OneD, NekDouble> x( ntotquad );
Array<OneD, NekDouble> y( ntotquad );
Array<OneD, NekDouble> z( ntotquad );
E->GetCoords(x,y,z);
//----------------------------------------------
// Define solution to be projected
for(i = 0; i < ntotquad; i++)
{
sol[i] = Hex_sol(x[i],y[i],z[i],P,Q,R, bType, bType, bType);
}
//----------------------------------------------
//----------------------------------------------
// Define the derivative solution
for(int n = 0; n < ntotquad; ++n){
diff_solution_x[n] = Hex_Diff_Sol( x[n], y[n], z[n], P,Q,R,
bType, bType, bType, 1);
diff_solution_y[n] = Hex_Diff_Sol( x[n], y[n], z[n], P,Q,R,
bType, bType, bType, 2);
diff_solution_z[n] = Hex_Diff_Sol( x[n], y[n], z[n], P,Q,R,
bType, bType, bType, 3);
}
//---------------------------------------------
//---------------------------------------------
// Project onto Expansion
Array<OneD, NekDouble> coeffs(E->GetNcoeffs());
Array<OneD, NekDouble> phys (ntotquad);
E->FwdTrans( sol, coeffs );
//---------------------------------------------
//-------------------------------------------
// Backward Transform Solution to get projected values
E->BwdTrans( coeffs, phys );
//-------------------------------------------
//--------------------------------------------
// Calculate L_p error
cout << "L infinity error: " << E->Linf(phys, sol) << endl;
cout << "L 2 error: " << E->L2 (phys, sol) << endl;
//--------------------------------------------
//-------------------------------------------
// Evaulate solution at x = y = z = 0 and print error
Array<OneD, NekDouble> t = Array<OneD, NekDouble>(3);
t[0] = -0.39;
t[1] = -0.25;
t[2] = 0.5;
NekDouble exact_sol = 0.0;
exact_sol = Hex_sol( t[0], t[1], t[2], nummodes0, nummodes1, nummodes2,
bType, bType, bType);
NekDouble numericSolution = E->PhysEvaluate(t, phys);
cout << "Numeric solution = " << numericSolution << endl;
cout << "Exact solution = " << exact_sol << endl;
cout << "Difference at x = ( "
<< t[0] << ", " << t[1] << ", " << t[2] << " ): "
<< numericSolution - exact_sol << endl;
//-------------------------------------------
E->PhysDeriv( sol, dx, dy, dz);
double error_x = 0, error_y=0, error_z=0;
for( int n = 0; n < ntotquad; ++n ) {
error_x += fabs(diff_solution_x[n] - dx[n]);
error_y += fabs(diff_solution_y[n] - dy[n]);
error_z += fabs(diff_solution_z[n] - dz[n]);
}