Commit db95f709 authored by Dave Moxey's avatar Dave Moxey
Browse files

Lots of temporary code

parent e78bc9c2
......@@ -10,6 +10,64 @@
using namespace Nektar;
void printSolution(StdRegions::StdExpansion *F,
std::string name,
Array<OneD, NekDouble> &x,
Array<OneD, NekDouble> &y,
Array<OneD, NekDouble> &phys)
{
int nPts = F->GetTotPoints();
ofstream outf(name);
int numBlocks = 1;
for (int j = 0; j < 2; ++j)
{
numBlocks *= F->GetNumPoints(j)-1;
}
outf << "VARIABLES = x, y, m" << endl;
outf << "Zone, N=" << nPts << ", E="
<< numBlocks << ", F=FEBlock, ET=QUADRILATERAL" << endl;
const int np0 = F->GetNumPoints(0);
const int np1 = F->GetNumPoints(1);
for (int j = 0; j < nPts; ++j)
{
outf << x[j] << " ";
if (j % 100 == 0 && j > 0)
outf << endl;
}
for (int j = 0; j < nPts; ++j)
{
outf << y[j] << " ";
if (j % 100 == 0 && j > 0)
outf << endl;
}
for (int j = 0; j < nPts; ++j)
{
outf << phys[j] << " ";
if (j % 100 == 0 && j > 0)
outf << endl;
}
outf << endl;
for(int j = 1; j < np1; ++j)
{
for(int k = 1; k < np0; ++k)
{
outf << (j-1)*np0 + k << " ";
outf << (j-1)*np0 + k+1 << " ";
outf << j *np0 + k+1 << " ";
outf << j *np0 + k << endl;
}
}
outf.close();
}
NekDouble Tri_sol(NekDouble x, NekDouble y, int order1, int order2);
NekDouble Quad_sol(NekDouble x, NekDouble y, int order1, int order2,
LibUtilities::BasisType btype1, LibUtilities::BasisType btype2);
......@@ -225,6 +283,13 @@ int main(int argc, char *argv[])
{
sol[i] = Quad_sol(x[i],y[i],order1,order2,btype1,btype2);
}
int nCoeffs = E->GetNcoeffs();
Array<OneD, NekDouble> tmp(nCoeffs, 0.0);
tmp[24] = 1.0;
E->BwdTrans(tmp, sol);
printSolution(E, "blah.dat", x, y, sol);
//---------------------------------------------
}
break;
......
......@@ -17,6 +17,79 @@ using namespace Nektar;
using namespace Nektar::LibUtilities;
using namespace Nektar::StdRegions;
void printSolution(StdRegions::StdExpansion *F,
std::string name,
Array<OneD, NekDouble> &x,
Array<OneD, NekDouble> &y,
Array<OneD, NekDouble> &z,
Array<OneD, NekDouble> &phys)
{
int nPts = F->GetTotPoints();
ofstream outf(name);
int numBlocks = 1;
for (int j = 0; j < 3; ++j)
{
numBlocks *= F->GetNumPoints(j)-1;
}
outf << "VARIABLES = x, y, z, m" << endl;
outf << "Zone, N=" << nPts << ", E="
<< numBlocks << ", F=FEBlock, ET=BRICK" << endl;
const int np0 = F->GetNumPoints(0);
const int np1 = F->GetNumPoints(1);
const int np2 = F->GetNumPoints(2);
const int np01 = np0*np1;
for (int j = 0; j < nPts; ++j)
{
outf << x[j] << " ";
if (j % 100 == 0 && j > 0)
outf << endl;
}
for (int j = 0; j < nPts; ++j)
{
outf << y[j] << " ";
if (j % 100 == 0 && j > 0)
outf << endl;
}
for (int j = 0; j < nPts; ++j)
{
outf << z[j] << " ";
if (j % 100 == 0 && j > 0)
outf << endl;
}
for (int j = 0; j < nPts; ++j)
{
outf << phys[j] << " ";
if (j % 100 == 0 && j > 0)
outf << endl;
}
for(int j = 1; j < np2; ++j)
{
for(int k = 1; k < np1; ++k)
{
for(int l = 1; l < np0; ++l)
{
outf << (j-1)*np01 + (k-1)*np0 + l << " ";
outf << (j-1)*np01 + (k-1)*np0 + l+1 << " ";
outf << (j-1)*np01 + k *np0 + l+1 << " ";
outf << (j-1)*np01 + k *np0 + l << " ";
outf << j *np01 + (k-1)*np0 + l << " ";
outf << j *np01 + (k-1)*np0 + l+1 << " ";
outf << j *np01 + k *np0 + l+1 << " ";
outf << j *np01 + k *np0 + l << endl;
}
}
}
outf.close();
}
/// Defines a solution which excites all modes in a StdTet expansion.
NekDouble Tet_sol(NekDouble x, NekDouble y, NekDouble z,
int order1, int order2, int order3);
......@@ -332,7 +405,7 @@ int main(int argc, char *argv[]){
StdRegions::StdPyrExp *F = new StdRegions::StdPyrExp(Bkey1,Bkey2,Bkey3);
E = F;
pyrIdx = F->GetMap();
vector<int> &rmap = F->GetRMap();
E->GetCoords(x,y,z);
//----------------------------------------------
......@@ -346,85 +419,98 @@ int main(int argc, char *argv[]){
#if 0
int nCoeffs = F->GetNcoeffs();
int nPts = F->GetTotPoints();
Array<OneD, NekDouble> blah(nCoeffs), blah2(nPts);
for (i = 0; i < nCoeffs; ++i)
{
Vmath::Zero(nCoeffs, blah, 1);
blah[i] = 1.0;
F->BwdTrans(blah, blah2);
boost::format pad("mode-%03d.dat");
pad % i;
ofstream outf(pad.str());
int numBlocks = 1;
for (int j = 0; j < 3; ++j)
{
numBlocks *= F->GetNumPoints(j)-1;
}
outf << "VARIABLES = x, y, z, m" << endl;
outf << "Zone, N=" << nPts << ", E="
<< numBlocks << ", F=FEBlock, ET=BRICK" << endl;
const int np0 = F->GetNumPoints(0);
const int np1 = F->GetNumPoints(1);
const int np2 = F->GetNumPoints(2);
const int np01 = np0*np1;
const Array<OneD, const NekDouble> &bx = F->GetBase()[0]->GetBdata();
const Array<OneD, const NekDouble> &by = F->GetBase()[1]->GetBdata();
const Array<OneD, const NekDouble> &bz = F->GetBase()[2]->GetBdata();
for (int j = 0; j < nPts; ++j)
for (int cnt = 0; cnt < nCoeffs; ++cnt)
{
StdPyrExp::triple &idx = pyrIdx[cnt];
const int p = boost::get<0>(idx);
const int q = boost::get<1>(idx);
const int r = boost::get<2>(idx);
if (r == 0 && p == 4 && q == 4)
{
outf << x[j] << " ";
if (j % 100 == 0 && j > 0)
outf << endl;
Array<OneD, NekDouble> asd(nPts);
for (int k = 0; k < nq3; ++k)
{
cout << z[k*nq1*nq2] << " " << bz[k + nq3*F->GetTetMode(3,2,1)] << endl;
/*
for (int j = 0; j < nq2; ++j)
{
for (int i = 0; i < nq1; ++i)
{
int s = i + nq1*(j + nq2*k);
asd[s] =
bx[i + nq1*p]*
by[j + nq2*q]*
bz[k + nq3*rmap[cnt]]*
bz[k + nq3*F->GetTetMode(q-1,0,0)];
}
}
*/
}
printSolution(F, "quadmode.dat", x, y, z, asd);
exit(0);
}
}
#endif
#if 0
int nCoeffs = F->GetNcoeffs();
int nPts = F->GetTotPoints();
Array<OneD, NekDouble> blah(nCoeffs), blah2(nPts);
for (i = 0; i < nCoeffs; ++i)
{
int p = boost::get<0>(pyrIdx[i]);
int q = boost::get<1>(pyrIdx[i]);
int r = boost::get<2>(pyrIdx[i]);
for (int j = 0; j < nPts; ++j)
if (r == 0 && p >= 2 && q >= 2)
{
outf << y[j] << " ";
if (j % 100 == 0 && j > 0)
outf << endl;
}
Vmath::Zero(nCoeffs, blah, 1);
blah[i] = 1.0;
F->BwdTrans(blah, blah2);
for (int j = 0; j < nPts; ++j)
{
outf << z[j] << " ";
if (j % 100 == 0 && j > 0)
outf << endl;
}
boost::format pad("mode-%03d.dat");
pad % i;
for (int j = 0; j < nPts; ++j)
{
outf << blah2[j] << " ";
if (j % 100 == 0 && j > 0)
outf << endl;
printSolution(F, pad.str(), x, y, z, blah2);
}
}
exit(0);
#endif
for(int j = 1; j < np2; ++j)
#if 0
int nCoeffs = F->GetNcoeffs();
int nPts = F->GetTotPoints();
Array<OneD, NekDouble> blah(nCoeffs), blah2(nPts);
for (i = 0; i < nCoeffs; ++i)
{
Vmath::Zero(nCoeffs, blah, 1);
blah[i] = 1.0;
F->BwdTrans(blah, blah2);
F->FwdTrans(blah2, blah);
int zeroCount = 0;
int oneCount = 0;
for (int j = 0; j < nCoeffs; ++j)
{
for(int k = 1; k < np1; ++k)
if (fabs(blah[j]) < 1e-8)
{
for(int l = 1; l < np0; ++l)
{
outf << (j-1)*np01 + (k-1)*np0 + l << " ";
outf << (j-1)*np01 + (k-1)*np0 + l+1 << " ";
outf << (j-1)*np01 + k *np0 + l+1 << " ";
outf << (j-1)*np01 + k *np0 + l << " ";
outf << j *np01 + (k-1)*np0 + l << " ";
outf << j *np01 + (k-1)*np0 + l+1 << " ";
outf << j *np01 + k *np0 + l+1 << " ";
outf << j *np01 + k *np0 + l << endl;
}
zeroCount++;
}
if (fabs(blah[j]-1.0) < 1e-8)
{
oneCount++;
}
}
outf.close();
//F->FwdTrans(sol, blah);
//cout << Vmath::Vmin(nCoeffs, blah, 1) << " " << Vmath::Vmax(nCoeffs, blah, 1) << endl;
cout << i << " " << zeroCount << " " << oneCount << endl;
}
exit(0);
#endif
//----------------------------------------------
}
break;
......@@ -500,6 +586,12 @@ int main(int argc, char *argv[]){
E->BwdTrans(coeffs,phys);
//-------------------------------------------
int nPts = E->GetTotPoints();
Array<OneD, NekDouble> errArr(nPts);
Vmath::Vsub(nPts, phys, 1, sol, 1, errArr, 1);
printSolution(E, "out.dat", x, y, z, errArr);
printSolution(E, "sol.dat", x, y, z, phys);
//--------------------------------------------
// Calculate L_inf error
cout << "L infinity error: " << E->Linf(phys,sol) << endl;
......@@ -547,16 +639,20 @@ NekDouble Tet_sol(NekDouble x, NekDouble y, NekDouble z,
int l,k,m;
NekDouble sol = 0.0;
//return x*x*y*y;
return x*x*y*y;
for(k = 0; k < order1; ++k)
{
for(l = 0; l < order2-k; ++l)
{
for(m = 0; m < order3-k-l; ++m)
{
sol += pow(x,k)*pow(y,l)*pow(z,m);
sol += pow(x,k)*pow(y,max(0,l-1))*pow(z,m);
}
}
}
return sol;
}
......@@ -586,7 +682,7 @@ NekDouble Prism_sol(NekDouble x, NekDouble y, NekDouble z,
{
for(m = 0; m < order3-k; ++m)
{
sol += pow(x,k)*pow(y,l)*pow(z,m);
sol += pow(z,m);
}
}
}
......
......@@ -93,28 +93,28 @@ namespace Nektar
// Edge 0
for (int i = 2; i <= P; ++i)
{
m_map [cnt ] = triple(i, 0, 0);
m_map [cnt ] = triple (i, 0, 0);
m_rmap[cnt++] = GetTetMode(i, 0, 0);
}
// Edge 1
for (int i = 2; i <= Q; ++i)
{
m_map [cnt ] = triple(1, i, 0);
m_map [cnt ] = triple (1, i, 0);
m_rmap[cnt++] = GetTetMode(0, i, 0);
}
// Edge 2
for (int i = 2; i <= P; ++i)
{
m_map [cnt ] = triple(i, 1, 0);
m_map [cnt ] = triple (i, 1, 0);
m_rmap[cnt++] = GetTetMode(i, 0, 0);
}
// Edge 3
for (int i = 2; i <= Q; ++i)
{
m_map [cnt ] = triple(0, i, 0);
m_map [cnt ] = triple (0, i, 0);
m_rmap[cnt++] = GetTetMode(0, i, 0);
}
......@@ -147,11 +147,11 @@ namespace Nektar
}
// Face 0 - TODO check this
for (int i = 2; i <= P; ++i)
for (int j = 2; j <= Q; ++j)
{
for (int j = 2; j <= Q; ++j)
for (int i = 2; i <= P; ++i)
{
m_map [cnt ] = triple(i, j, 0);
m_map [cnt ] = triple (i, j, 0);
m_rmap[cnt++] = GetTetMode(i, 0, 0);
}
}
......@@ -161,7 +161,7 @@ namespace Nektar
{
for (int j = 1; j <= R-i; ++j)
{
m_map [cnt ] = triple(i, 0, j);
m_map [cnt ] = triple (i, 0, j);
m_rmap[cnt++] = GetTetMode(i, 0, j);
}
}
......@@ -171,7 +171,7 @@ namespace Nektar
{
for (int j = 1; j <= R-i; ++j)
{
m_map [cnt ] = triple(1, i, j);
m_map [cnt ] = triple (1, i, j);
m_rmap[cnt++] = GetTetMode(0, i, j);
}
}
......@@ -181,7 +181,7 @@ namespace Nektar
{
for (int j = 1; j <= R-i; ++j)
{
m_map [cnt ] = triple(i, 1, j);
m_map [cnt ] = triple (i, 1, j);
m_rmap[cnt++] = GetTetMode(i, 0, j);
}
}
......@@ -191,7 +191,7 @@ namespace Nektar
{
for (int j = 1; j <= R-i; ++j)
{
m_map [cnt ] = triple(0, i, j);
m_map [cnt ] = triple (0, i, j);
m_rmap[cnt++] = GetTetMode(0, i, j);
}
}
......@@ -203,31 +203,33 @@ namespace Nektar
{
for (int k = 1; k <= R-i-j; ++k)
{
m_map [cnt ] = triple(i, j+1, k);
m_rmap[cnt++] = GetTetMode(i, j, k);
// need to go to j+1-th mode in the 'b' direction to
// select correct modified_a mode
m_map [cnt ] = triple (i, j+1, k);
m_rmap[cnt++] = GetTetMode(i, j, k);
}
}
}
#if 0
int nPts = GetTotPoints();
Array<OneD, NekDouble> modei(nPts);
Array<OneD, NekDouble> modej(nPts);
for (int i = 0; i < m_ncoeffs; ++i)
{
v_FillMode(i, modei);
for (int j = i+1; j < m_ncoeffs; ++j)
{
v_FillMode(j, modej);
Vmath::Vsub(nPts, modei, 1, modej, 1, modej, 1);
if (fabs(Vmath::Vmax(nPts, modej, 1)) < 1e-4)
{
cout << "wut" << i << " " << j << endl;
}
}
}
#endif
// const int P2 = m_base[0]->GetNumModes();
// const int Q2 = m_base[1]->GetNumModes();
// const int R2 = m_base[2]->GetNumModes();
// int idx = 0;
// for (int i = 2; i < P2-2; ++i)
// {
// for (int j = 1; j < Q2-i-1; ++j)
// {
// for (int k = 1; k < R2-i-j; ++k)
// {
// idx++;
// }
// }
// }
//cout << "IDX = " << idx << endl;
cout << cnt << endl;
}
StdPyrExp::StdPyrExp(const StdPyrExp &T)
......@@ -584,7 +586,7 @@ namespace Nektar
Array<OneD, const NekDouble> by = m_base[1]->GetBdata();
Array<OneD, const NekDouble> bz = m_base[2]->GetBdata();
int Q = m_base[2]->GetNumModes();
int Q = m_base[2]->GetNumModes()-1;
for (int k = 0; k < Qz; ++k)
{
......@@ -605,10 +607,11 @@ namespace Nektar
by[j + Qy*q]*
bz[k + Qz*m_rmap[cnt]];
if (r == 0 && p >= 2 && q > Q-p)
if (r == 0 && p >= 2 && q >= 2 && q > Q-p)
{
tmp *= bz[k + Qz*GetTetMode(Q-q,0,0)];
tmp *= bz[k + Qz*GetTetMode(q-1,0,0)];
}
sum += tmp;
}
......@@ -712,7 +715,7 @@ namespace Nektar
const Array<OneD, const NekDouble> &by = m_base[1]->GetBdata();
const Array<OneD, const NekDouble> &bz = m_base[2]->GetBdata();
int Q = m_base[1]->GetNumModes();
int Q = m_base[1]->GetNumModes()-1;
// Initial pyramid implementation. We need to iterate over vertices,
// edge int, face int and then interior.
......@@ -738,9 +741,9 @@ namespace Nektar
by[j + Qy*q] *
bz[k + Qz*m_rmap[cnt]];
if (r == 0 && p >= 2 && q > Q-p)
if (r == 0 && p >= 2 && q >= 2 && q > Q-p)
{
tmp *= bz[k + Qz*GetTetMode(Q-q,0,0)];
tmp *= bz[k + Qz*GetTetMode(q-1,0,0)];
}
g_pqr[s] += tmp;
......
......@@ -93,6 +93,12 @@ namespace Nektar
return m_map;
}
STD_REGIONS_EXPORT vector<int> &GetRMap()
{
return m_rmap;
}
STD_REGIONS_EXPORT int GetTetMode(int I, int J, int K);
protected:
//---------------------------------------
// Differentiation/integration Methods
......@@ -242,7 +248,6 @@ namespace Nektar
// Private helper functions
//---------------------------------------
STD_REGIONS_EXPORT int GetMode(int I, int J, int K);
STD_REGIONS_EXPORT int GetTetMode(int I, int J, int K);
STD_REGIONS_EXPORT void MultiplyByQuadratureMetric(
const Array<OneD, const NekDouble>& inarray,
Array<OneD, NekDouble> &outarray);
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment