Commit d8fd5ab9 authored by Spencer Sherwin's avatar Spencer Sherwin

Merge branch 'master' into feature/MeshConvertPlyNormals

parents 7ad5f009 b6e1a185
......@@ -5,6 +5,9 @@
SET(MemoryManagerSources
MemoryManager.cpp)
SET(PartitionAnalyseSources
PartitionAnalyse.cpp)
SET(VecMatSources
VecMatExample.cpp)
......@@ -37,6 +40,9 @@ SET(TimeIntegrationDemoSources
ADD_NEKTAR_EXECUTABLE(Matrix demos VecMatSources)
SET_LAPACK_LINK_LIBRARIES(Matrix)
ADD_NEKTAR_EXECUTABLE(PartitionAnalyse demos PartitionAnalyseSources)
SET_LAPACK_LINK_LIBRARIES(PartitionAnalyse)
ADD_NEKTAR_EXECUTABLE(FoundationDemo demos FoundationSources )
SET_LAPACK_LINK_LIBRARIES(FoundationDemo)
......
///////////////////////////////////////////////////////////////////////////////
//
// File: PartitionAnalyse.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).
//
// 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: Small utility to export histogram of partition sizes.
//
///////////////////////////////////////////////////////////////////////////////
#include <LibUtilities/BasicUtils/SessionReader.h>
#include <LibUtilities/BasicUtils/MeshPartition.h>
#include <LibUtilities/Communication/CommSerial.h>
#include <iostream>
using namespace std;
using namespace Nektar;
using namespace Nektar::LibUtilities;
class FauxComm : public CommSerial
{
public:
FauxComm(int argc, char* argv[], int size)
: CommSerial(argc, argv)
{
m_size = size;
m_type = "Faux parallel";
}
FauxComm(int size) : CommSerial(0, NULL)
{
m_size = size;
m_type = "Faux parallel";
}
virtual ~FauxComm() {}
void v_SplitComm(int pRows, int pColumns)
{
m_commRow = boost::shared_ptr<FauxComm>(new FauxComm(pColumns));
m_commColumn = boost::shared_ptr<FauxComm>(new FauxComm(pRows));
}
};
int main(int argc, char *argv[])
{
if (argc < 3)
{
cerr << "Usage: PartitionAnalyse <nproc> <xml file1> [xml file 2..n]"
<< endl;
return 1;
}
int nParts = atoi(argv[1]);
vector<string> filenames(argv + 2, argv + argc);
CommSharedPtr vComm = boost::shared_ptr<FauxComm>(
new FauxComm(argc, argv, nParts));
char **new_argv = new char*[argc];
new_argv[0] = strdup("PartitionAnalyse");
new_argv[1] = strdup("--part-info");
for (int i = 0; i < argc-2; ++i)
{
new_argv[i+2] = strdup(filenames[i].c_str());
}
LibUtilities::SessionReaderSharedPtr vSession
= LibUtilities::SessionReader::CreateInstance(
argc, new_argv, filenames, vComm);
return 0;
}
......@@ -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)
......
This diff is collapsed.
#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]);
}
cout << "L 1 error of derivatives X = " << error_x << endl;
cout << "L 1 error of derivatives Y = " << error_y << endl;
cout << "L 1 error of derivatives Z = " << error_z << 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;
}
NekDouble Hex_Diff_Sol(NekDouble x, NekDouble y, NekDouble z,
int order1, int order2, int order3,
LibUtilities::BasisType btype1,
LibUtilities::BasisType btype2,
LibUtilities::BasisType btype3, int dir){
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)
{
switch (dir)
{
case 1:
a = i*(Fx ? M_PI*(cos(M_PI*i*x) - sin(M_PI*i*x))
: pow_loc(x,i-1));
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;
break;
case 2:
a = (Fx ? (sin(M_PI*i*x) + cos(M_PI*i*x)) : pow_loc(x,i));
a *= j*(Fy ? M_PI*(cos(M_PI*j*y) - sin(M_PI*j*y))
: pow_loc(y,j-1));
a *= (Fz ? (sin(M_PI*k*z) + cos(M_PI*k*z)) : pow_loc(z,k));
sol += a;
break;
case 3:
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 *= k*(Fz ? M_PI*(cos(M_PI*k*z) - sin(M_PI*k*z))
: pow_loc(z,k-1));
sol += a;
break;
}
}
}
}
return sol;
}
This diff is collapsed.
#include <cstdlib>
#include <LocalRegions/PrismExp.h>
#include <SpatialDomains/MeshGraph3D.h>
#include <LibUtilities/BasicConst/NektarUnivTypeDefs.hpp>
using namespace Nektar;
NekDouble Prism_sol(NekDouble x, NekDouble y, NekDouble z, int order1, int order2, int order3);
NekDouble Prism_sol2(NekDouble x, NekDouble y, NekDouble z, int P, int Q, int R,
LibUtilities::BasisType bType_x,
LibUtilities::BasisType bType_y,
LibUtilities::BasisType bType_z );
int main(int argc, char *argv[])
{
if(argc != 5)
{
cerr << "usage: LocPrismExpDemo MeshFile nummodes0 nummodes1 nummodes2" << 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 = nummodes0 + 1;
int nquad1 = nummodes1 + 1;
int nquad2 = nummodes2 + 1;
int ntotquad = nquad0*nquad1*nquad2;
SpatialDomains::MeshGraph3D graph3D;
graph3D.ReadGeometry(in);
LibUtilities::BasisType bTypeX = LibUtilities::eModified_A;
LibUtilities::BasisType bTypeY = LibUtilities::eModified_A;
LibUtilities::BasisType bTypeZ = LibUtilities::eModified_B;
LibUtilities::PointsType qtypeX = LibUtilities::eGaussLobattoLegendre;
LibUtilities::PointsType qtypeY = LibUtilities::eGaussLobattoLegendre;
LibUtilities::PointsType qtypeZ = LibUtilities::eGaussRadauMAlpha1Beta0;
const LibUtilities::PointsKey pointsKey0( nquad0, qtypeX );
const LibUtilities::PointsKey pointsKey1( nquad1, qtypeY );
const LibUtilities::PointsKey pointsKey2( nquad2, qtypeZ );
const LibUtilities::BasisKey basisKey0( bTypeX, nummodes0, pointsKey0 );
const LibUtilities::BasisKey basisKey1( bTypeY, nummodes1, pointsKey1 );
const LibUtilities::BasisKey basisKey2( bTypeZ, nummodes2, pointsKey2 );
SpatialDomains::PrismGeomSharedPtr geom;
if(!(geom = boost::dynamic_pointer_cast<SpatialDomains::PrismGeom>(graph3D.GetCompositeItem(0,0))))
{
cerr << "Could not find PrismGeom in input file" << endl;
exit(1);
}
LocalRegions::PrismExpSharedPtr E = MemoryManager<LocalRegions::PrismExp>::
AllocateSharedPtr(basisKey0,basisKey1,basisKey2,geom);
Array<OneD,NekDouble> x(ntotquad);
Array<OneD,NekDouble> y(ntotquad);
Array<OneD,NekDouble> z(ntotquad);
E->GetCoords(x,y,z);
Array<OneD,NekDouble> sol(ntotquad);
//----------------------------------------------
// Define solution to be projected
for(i = 0; i < ntotquad; i++)
{
sol[i] = Prism_sol(x[i],y[i],z[i],nummodes0,nummodes1,nummodes2);
}
//----------------------------------------------
//---------------------------------------------
// 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;
//--------------------------------------------
return 0;
}
NekDouble Prism_sol(NekDouble x, NekDouble y, NekDouble z, int nummodes1, int nummodes2, int nummodes3)