Skip to content

GitLab

Commit 29f1a54d authored by Gianmarco Mengaldo's avatar Gianmarco Mengaldo
Browse files
Options

Merge branch 'feature/cfs_cleaning' of gitlab.nektar.info:nektar into feature/cfs_cleaning

parents 7e534371 4abf4651
Expand all Hide whitespace changes
Inline Side-by-side
Showing with 4088 additions and 669 deletions
dist-exclude
View file @ 29f1a54d
......@@ -13,6 +13,7 @@ solvers/CardiacEPSolver
solvers/FluxReconstruction
solvers/VortexWaveInteraction
solvers/ImageWarpingSolver
solvers/PulseWaveSolver
docs/*.doc
docs/arch
docs/emacs
......
library/Demos/LocalRegions/LocHexExpDemo.cpp
View file @ 29f1a54d
......@@ -258,8 +258,8 @@ int main(int argc, char *argv[])
//-------------------------------------------
// 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[0] = 0.5;
t[1] = 0.5;
t[2] = 0.5;
NekDouble numericSolution = lhe->PhysEvaluate(t);
......
library/Demos/LocalRegions/LocPrismExpDemo.cpp
View file @ 29f1a54d
......@@ -247,11 +247,11 @@ int main(int argc, char *argv[])
//--------------------------------------------
//-------------------------------------------
// Evaulate solution at x = y = z = 0 and print error
// Evaulate solution at mid point and print error
Array<OneD, NekDouble> t = Array<OneD, NekDouble>(3);
t[0] = -0.39;
t[1] = -0.25;
t[2] = -0.50;
t[0] = 0.5;
t[1] = 0.5;
t[2] = 0.2;
if( regionShape == StdRegions::ePrism ) {
solution[0] = Prism_sol( t[0], t[1], t[2], P, Q, R, bType_x, bType_y, bType_z );
......
library/Demos/LocalRegions/LocProject2D.cpp
View file @ 29f1a54d
......@@ -19,9 +19,8 @@ NekDouble Quad_sol(NekDouble x, NekDouble y, int order1, int order2,
int main(int argc, char *argv[])
{
int i;
int order1,order2, nq1,nq2,NodalTri = 0;
LibUtilities::PointsType Qtype1,Qtype2;
LibUtilities::BasisType btype1,btype2;
......@@ -31,23 +30,23 @@ int main(int argc, char *argv[])
Array<OneD, NekDouble> sol;
Array<OneD, NekDouble> coords(8);
StdRegions::Orientation edgeDir = StdRegions::eForwards;
if((argc != 16)&&(argc != 14))
{
// arg[0] arg[1] arg[2] arg[3] arg[4] arg[5] arg[6] arg[7] arg[8] arg[9] arg[10] arg[11] arg[12] arg[13]
fprintf(stderr,"Usage: Project2D RegionShape Type1 Type2 order1 order2 nq1 nq2 x1, y1, x2, y2, x3, y3 \n");
fprintf(stderr,"Example : ./LocProject2D-g 2 4 5 3 3 4 4 .1 -.5 .6 .1 .3 .2 \n" );
fprintf(stderr,"Where RegionShape is an integer value which "
"dictates the region shape:\n");
fprintf(stderr,"\t Triangle = 2\n");
fprintf(stderr,"\t Quadrilateral = 3\n");
fprintf(stderr,"Where type is an integer value which "
"dictates the basis as:\n");
fprintf(stderr,"\t Ortho_A = 1\n");
fprintf(stderr,"\t Ortho_B = 2\n");
fprintf(stderr,"\t Modified_A = 4\n");
......@@ -58,24 +57,24 @@ int main(int argc, char *argv[])
fprintf(stderr,"\t Chebyshev = 10\n");
fprintf(stderr,"\t Nodal tri (Electro) = 11\n");
fprintf(stderr,"\t Nodal tri (Fekete) = 12\n");
fprintf(stderr,"Note type = 3,6 are for three-dimensional basis\n");
fprintf(stderr,"The last series of values are the coordinates\n");
exit(1);
}
regionshape = (StdRegions::ExpansionType) atoi(argv[1]);
// Check to see if 2D region
if((regionshape != StdRegions::eTriangle)&&(regionshape != StdRegions::eQuadrilateral))
{
NEKERROR(ErrorUtil::efatal,"This shape is not a 2D region");
}
int btype1_val = atoi(argv[2]);
int btype2_val = atoi(argv[3]);
if(( btype1_val <= 10)&&( btype2_val <= 10))
{
btype1 = (LibUtilities::BasisType) btype1_val;
......@@ -85,7 +84,7 @@ int main(int argc, char *argv[])
{
btype1 = LibUtilities::eOrtho_A;
btype2 = LibUtilities::eOrtho_B;
if(btype1_val == 11)
{
NodalType = LibUtilities::eNodalTriElec;
......@@ -94,49 +93,47 @@ int main(int argc, char *argv[])
{
NodalType = LibUtilities::eNodalTriFekete;
}
}
// Check to see that correct Expansions are used
switch(regionshape)
{
case StdRegions::eTriangle:
if((btype1 == LibUtilities::eOrtho_B)||(btype1 == LibUtilities::eModified_B))
{
NEKERROR(ErrorUtil::efatal,
"Basis 1 cannot be of type Ortho_B or Modified_B");
}
break;
case StdRegions::eQuadrilateral:
if((btype1 == LibUtilities::eOrtho_B)||(btype1 == LibUtilities::eOrtho_C)||
(btype1 == LibUtilities::eModified_B)||(btype1 == LibUtilities::eModified_C))
{
NEKERROR(ErrorUtil::efatal,
"Basis 1 is for 2 or 3D expansions");
}
if((btype2 == LibUtilities::eOrtho_B)||(btype2 == LibUtilities::eOrtho_C)||
(btype2 == LibUtilities::eModified_B)||(btype2 == LibUtilities::eModified_C))
{
NEKERROR(ErrorUtil::efatal,
"Basis 2 is for 2 or 3D expansions");
}
break;
default:
ASSERTL0(false, "Not a 2D expansion.");
break;
case StdRegions::eTriangle:
if((btype1 == LibUtilities::eOrtho_B)||(btype1 == LibUtilities::eModified_B))
{
NEKERROR(ErrorUtil::efatal,
"Basis 1 cannot be of type Ortho_B or Modified_B");
}
break;
case StdRegions::eQuadrilateral:
if((btype1 == LibUtilities::eOrtho_B)||(btype1 == LibUtilities::eOrtho_C)||
(btype1 == LibUtilities::eModified_B)||(btype1 == LibUtilities::eModified_C))
{
NEKERROR(ErrorUtil::efatal,
"Basis 1 is for 2 or 3D expansions");
}
if((btype2 == LibUtilities::eOrtho_B)||(btype2 == LibUtilities::eOrtho_C)||
(btype2 == LibUtilities::eModified_B)||(btype2 == LibUtilities::eModified_C))
{
NEKERROR(ErrorUtil::efatal,
"Basis 2 is for 2 or 3D expansions");
}
break;
default:
NEKERROR(ErrorUtil::efatal, "Not a valid 2D expansion.");
}
order1 = atoi(argv[4]);
order2 = atoi(argv[5]);
nq1 = atoi(argv[6]);
nq2 = atoi(argv[7]);
sol = Array<OneD, NekDouble>(nq1*nq2);
if(btype1 != LibUtilities::eFourier)
{
Qtype1 = LibUtilities::eGaussLobattoLegendre;
......@@ -145,7 +142,7 @@ int main(int argc, char *argv[])
{
Qtype1 = LibUtilities::eFourierEvenlySpaced;
}
if(btype2 != LibUtilities::eFourier)
{
if (regionshape == StdRegions::eTriangle) {
......@@ -154,28 +151,28 @@ int main(int argc, char *argv[])
else
{
Qtype2 = LibUtilities::eGaussLobattoLegendre;
}
}
}
else
{
Qtype2 = LibUtilities::eFourierEvenlySpaced;
}
//-----------------------------------------------
// Define a 2D expansion based on basis definition
switch(regionshape)
{
case StdRegions::eTriangle:
case StdRegions::eTriangle:
{
coords[0] = atof(argv[8]);
coords[1] = atof(argv[9]);
coords[2] = atof(argv[10]);
coords[3] = atof(argv[11]);
coords[4] = atof(argv[12]);
coords[5] = atof(argv[13]);
// Set up coordinates
SpatialDomains::VertexComponentSharedPtr verts[3];
const int zero = 0;
......@@ -185,21 +182,21 @@ int main(int argc, char *argv[])
verts[0] = MemoryManager<SpatialDomains::VertexComponent>::AllocateSharedPtr(two,zero,coords[0],coords[1],dZero);
verts[1] = MemoryManager<SpatialDomains::VertexComponent>::AllocateSharedPtr(two,one,coords[2],coords[3],dZero);
verts[2] = MemoryManager<SpatialDomains::VertexComponent>::AllocateSharedPtr(two,two,coords[4],coords[5],dZero);
// Set up Edges
SpatialDomains::SegGeomSharedPtr edges[3];
edges[0] = MemoryManager<SpatialDomains::SegGeom>::AllocateSharedPtr(zero,verts[0],verts[1]);
edges[1] = MemoryManager<SpatialDomains::SegGeom>::AllocateSharedPtr(one,verts[1],verts[2]);
edges[2] = MemoryManager<SpatialDomains::SegGeom>::AllocateSharedPtr(two,verts[2],verts[0]);
StdRegions::Orientation eorient[3];
eorient[0] = edgeDir;
eorient[1] = edgeDir;
eorient[2] = edgeDir;
SpatialDomains::TriGeomSharedPtr geom = MemoryManager<SpatialDomains::TriGeom>::AllocateSharedPtr(zero,verts,edges,eorient);
geom->SetOwnData();
const LibUtilities::PointsKey Pkey1(nq1,Qtype1);
const LibUtilities::PointsKey Pkey2(nq2,Qtype2);
const LibUtilities::BasisKey Bkey1(btype1,order1,Pkey1);
......@@ -213,11 +210,11 @@ int main(int argc, char *argv[])
{
E = new LocalRegions::TriExp(Bkey1,Bkey2,geom);
}
Array<OneD,NekDouble> x = Array<OneD,NekDouble>(nq1*nq2);
Array<OneD,NekDouble> y = Array<OneD,NekDouble>(nq1*nq2);
E->GetCoords(x,y);
//----------------------------------------------
// Define solution to be projected
for(i = 0; i < nq1*nq2; ++i)
......@@ -225,10 +222,10 @@ int main(int argc, char *argv[])
sol[i] = Tri_sol(x[i],y[i],order1,order2);
}
//----------------------------------------------
}
break;
case StdRegions::eQuadrilateral:
case StdRegions::eQuadrilateral:
{
// Gather coordinates
coords[0] = atof(argv[8]);
......@@ -239,7 +236,7 @@ int main(int argc, char *argv[])
coords[5] = atof(argv[13]);
coords[6] = atof(argv[14]);
coords[7] = atof(argv[15]);
// Set up coordinates
const int zero=0;
const int one=1;
......@@ -251,79 +248,77 @@ int main(int argc, char *argv[])
verts[1] = MemoryManager<SpatialDomains::VertexComponent>::AllocateSharedPtr(two,one,coords[2],coords[3],dZero);
verts[2] = MemoryManager<SpatialDomains::VertexComponent>::AllocateSharedPtr(two,two,coords[4],coords[5],dZero);
verts[3] = MemoryManager<SpatialDomains::VertexComponent>::AllocateSharedPtr(two,three,coords[6],coords[7],dZero);
// Set up Edges
SpatialDomains::SegGeomSharedPtr edges[4];
edges[0] = MemoryManager<SpatialDomains::SegGeom>::AllocateSharedPtr(zero,verts[0],verts[1]);
edges[1] = MemoryManager<SpatialDomains::SegGeom>::AllocateSharedPtr(one,verts[1],verts[2]);
edges[2] = MemoryManager<SpatialDomains::SegGeom>::AllocateSharedPtr(two,verts[2],verts[3]);
edges[3] = MemoryManager<SpatialDomains::SegGeom>::AllocateSharedPtr(three,verts[3],verts[0]);
StdRegions::Orientation eorient[4];
eorient[0] = edgeDir;
eorient[1] = edgeDir;
eorient[2] = edgeDir;
eorient[3] = edgeDir;
SpatialDomains::QuadGeomSharedPtr geom = MemoryManager<SpatialDomains::QuadGeom>::AllocateSharedPtr(zero,verts,edges,eorient);
geom->SetOwnData();
const LibUtilities::PointsKey Pkey1(nq1,Qtype1);
const LibUtilities::PointsKey Pkey2(nq2,Qtype2);
const LibUtilities::BasisKey Bkey1(btype1,order1,Pkey1);
const LibUtilities::BasisKey Bkey2(btype2,order2,Pkey2);
E = new LocalRegions::QuadExp(Bkey1,Bkey2,geom);
//----------------------------------------------
// Define solution to be projected
Array<OneD, NekDouble> x = Array<OneD, NekDouble>(nq1*nq2);
Array<OneD, NekDouble> y = Array<OneD, NekDouble>(nq1*nq2);
E->GetCoords(x,y);
for(i = 0; i < nq1*nq2; ++i)
{
sol[i] = Quad_sol(x[i],y[i],order1,order2,btype1,btype2);
}
//---------------------------------------------
}
break;
default:
ASSERTL0(false, "Not a 2D expansion.");
break;
default:
NEKERROR(ErrorUtil::efatal, "Not a valid 2D expansion.");
}
//---------------------------------------------
// Project onto Expansion
E->FwdTrans(sol,E->UpdateCoeffs());
//---------------------------------------------
//-------------------------------------------
// Backward Transform Solution to get projected values
E->BwdTrans(E->GetCoeffs(),E->UpdatePhys());
//-------------------------------------------
//--------------------------------------------
// Write solution
ofstream outfile("ProjectFile2D.dat");
E->WriteToFile(outfile,eTecplot);
outfile.close();
//-------------------------------------------
//--------------------------------------------
// Calculate L_inf error
cout << "L infinity error: " << E->Linf(sol) << endl;
cout << "L 2 error: " << E->L2 (sol) << endl;
//--------------------------------------------
//-------------------------------------------
// Evaulate solution at x = y =0 and print error
Array<OneD, NekDouble> x(2);
x[0] = (coords[0] + coords[2])*0.5;
x[1] = (coords[1] + coords[5])*0.5;
if(regionshape == StdRegions::eTriangle)
{
sol[0] = Tri_sol(x[0],x[1],order1,order2);
......@@ -332,11 +327,10 @@ int main(int argc, char *argv[])
{
sol[0] = Quad_sol(x[0],x[1],order1,order2,btype1,btype2);
}
Array<OneD,NekDouble> lcoord(2,0.0);
NekDouble nsol = E->PhysEvaluate(lcoord);
NekDouble nsol = E->PhysEvaluate(x);
cout << "error at x = (" <<x[0] <<","<<x[1] <<"): " << nsol - sol[0] << endl;
return 0;
}
......
library/Demos/MultiRegions/Helmholtz1D.cpp
View file @ 29f1a54d
......@@ -100,7 +100,9 @@ int main(int argc, char *argv[])
//----------------------------------------------
//----------------------------------------------
// Helmholtz solution taking physical forcing
//Helmholtz solution taking physical forcing after setting
//initial condition to zero
Vmath::Zero(Exp->GetNcoeffs(),Exp->UpdateCoeffs(),1);
Exp->HelmSolve(Fce->GetPhys(), Exp->UpdateCoeffs(), NullFlagList, factors);
//----------------------------------------------
......
library/Demos/MultiRegions/Helmholtz2D.cpp
View file @ 29f1a54d
......@@ -128,8 +128,10 @@ int main(int argc, char *argv[])
//----------------------------------------------
Timing("Define forcing ..");
//----------------------------------------------
// Helmholtz solution taking physical forcing
//----------------------------------------------
//Helmholtz solution taking physical forcing after setting
//initial condition to zero
Vmath::Zero(Exp->GetNcoeffs(),Exp->UpdateCoeffs(),1);
Exp->HelmSolve(Fce->GetPhys(), Exp->UpdateCoeffs(), flags, factors, varcoeffs);
//----------------------------------------------
Timing("Helmholtz Solve ..");
......
library/Demos/MultiRegions/Helmholtz3D.cpp
View file @ 29f1a54d
......@@ -140,8 +140,10 @@ int main(int argc, char *argv[])
Fce->SetPhys(fce);
//----------------------------------------------
//----------------------------------------------
// Helmholtz solution taking physical forcing
//----------------------------------------------
//Helmholtz solution taking physical forcing after setting
//initial condition to zero
Vmath::Zero(Exp->GetNcoeffs(),Exp->UpdateCoeffs(),1);
Exp->HelmSolve(Fce->GetPhys(), Exp->UpdateCoeffs(), flags, factors);
//----------------------------------------------
......
library/LibUtilities/BasicUtils/DBUtils.hpp 0 → 100644
View file @ 29f1a54d
///////////////////////////////////////////////////////////////////////////////
//
// File DBUtils.hpp
//
// 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: output function for use in debugging
//
///////////////////////////////////////////////////////////////////////////////
#ifndef NEKTAR_LIB_LIBUTILITIES_DBUTILS_HPP
#define NEKTAR_LIB_LIBUTILITIES_DBUTILS_HPP
#include <LibUtilities/BasicUtils/VmathArray.hpp>
namespace DBUtils
{
using namespace Nektar;
const int StopDefault = -99;
template<class T> void Output1DArray(const Array<OneD, const T> &in, const int start = 0,
const int stop = StopDefault)
{
int i;
ASSERTL1(start < in.num_elements(), "Start value is outside array range ");
if(stop == StopDefault)
{
for(i = start; i < in.num_elements(); ++i)
{
cout << in[i] << endl;
}
}
else
{
ASSERTL1(stop <= in.num_elements(), "Stop value is outside array range ");
for(i = start; i < in.num_elements(); ++i)
{
cout << in[i] << endl;
}
}
}
template<class T> void Output1DArray(const Array<OneD, const T> &in, std::string outfile,
const int start = 0,
const int stop = StopDefault)
{
int i;
ASSERTL1(start < in.num_elements(), "Start value is outside array range ");
ofstream ofile(outfile.c_str());
if(stop == StopDefault)
{
for(i = start; i < in.num_elements(); ++i)
{
ofile << in[i] << endl;
}
}
else
{
ASSERTL1(stop <= in.num_elements(), "Stop value is outside array range ");
for(i = start; i < stop; ++i)
{
ofile << in[i] << endl;
}
}
}
template<class T> void Output1DArray(const Array<OneD, const T> &in, ofstream &ofile,
const int start = 0,
const int stop = StopDefault)
{
int i;
ASSERTL1(start < in.num_elements(), "Start value is outside array range ");
if(stop == StopDefault)
{
for(i = start; i < in.num_elements(); ++i)
{
ofile << in[i] << endl;
}
}
else
{
ASSERTL1(stop <= in.num_elements(), "Stop value is outside array range ");
for(i = start; i < stop; ++i)
{