Skip to content

GitLab

Commit 9bac580d authored by Spencer Sherwin's avatar Spencer Sherwin
Browse files

Address comments from Chris

parent 2307d9db
Hide whitespace changes
Inline Side-by-side
docs/doxygen/Doxyfile.in
View file @ 9bac580d
......@@ -757,6 +757,7 @@ INPUT = @CMAKE_SOURCE_DIR@/docs/doxygen/ \
@CMAKE_SOURCE_DIR@/library/LibUtilities/ \
@CMAKE_SOURCE_DIR@/library/StdRegions/ \
@CMAKE_SOURCE_DIR@/library/SpatialDomains/ \
@CMAKE_SOURCE_DIR@/library/Collections/ \
@CMAKE_SOURCE_DIR@/library/LocalRegions/ \
@CMAKE_SOURCE_DIR@/library/MultiRegions/ \
@CMAKE_SOURCE_DIR@/library/GlobalMapping/ \
......
library/Collections/BwdTrans.cpp
View file @ 9bac580d
......@@ -1052,7 +1052,7 @@ class BwdTrans_SumFac_Pyr : public Operator
int mode = 0;
int mode1 = 0;
int cnt = 0;
for (i = mode = mode1 = 0; i < m_nmodes0; ++i)
for (i = 0; i < m_nmodes0; ++i)
{
for (j = 0; j < m_nmodes1; ++j, ++cnt)
{
......
library/Collections/IProductWRTDerivBase.cpp
View file @ 9bac580d
......@@ -697,7 +697,38 @@ class IProductWRTDerivBase_SumFac_Tri : public Operator
{
}
virtual void operator()(
/**
* This method calculates:
*
* \f[ (d\phi/dx,in[0]) + (d\phi/dy,in[1]) \f]
*
* which can be represented in terms of local cartesian
* derivaties as:
*
* \f[ ((d\phi/d\xi_0\, d\xi_0/dx +
* d\phi/d\xi_1\, d\xi_1/dx),in[0]) + \f]
*
* \f[ ((d\phi/d\xi_0\, d\xi_0/dy +
* d\phi/d\xi_1\, d\xi_1/dy),in[1]) + \f]
*
* where we note that
*
* \f[ d\phi/d\xi_0 = d\phi/d\eta_0\, d\eta_0/d\xi_0 =
* d\phi/d\eta_0 2/(1-\eta_1) \f]
*
* \f[ d\phi/d\xi_1 = d\phi/d\eta_1\, d\eta_1/d\xi_1 +
* d\phi/d\eta_1\, d\eta_1/d\xi_1 = d\phi/d\eta_0 (1+\eta_0)/(1-\eta_1)
* + d\phi/d\eta_1 \f]
*
* and so the full inner products are
*
* \f[ (d\phi/dx,in[0]) + (dphi/dy,in[1]) =
* (d\phi/d\eta_0, ((2/(1-\eta_1) (d\xi_0/dx in[0] + d\xi_0/dy in[1])
* + (1-\eta_0)/(1-\eta_1) (d\xi_1/dx in[0]+d\xi_1/dy in[1]))
* + (d\phi/d\eta_1, (d\xi_1/dx in[0] + d\xi_1/dy in[1])) \f]
*
*/
virtual void operator()(
const Array<OneD, const NekDouble> &entry0,
Array<OneD, NekDouble> &entry1,
Array<OneD, NekDouble> &entry2,
......@@ -719,26 +750,6 @@ class IProductWRTDerivBase_SumFac_Tri : public Operator
tmp[0] = wsp; tmp[1] = wsp + nmax;
wsp1 = wsp + 2*nmax;
// calculate (dphi/dx,in[0]) = ((dphi/dxi_0 dxi_0/dx +
// dphi/dxi_1 dxi_1/dx),in[0])
// + (dphi/dy,in[1]) = ((dphi/dxi_0 dxi_0/dy +
// dphi/dxi_1 dxi_1/dy),in[1])
//
// Note dphi/dxi_0 =
// dphi/deta_0 deta_0/dxi_0 = dphi/deta_0 2/(1-eta_1)
//
// dphi/dxi_1 =
// dphi/deta_1 deta_1/dxi_1 + dphi/deta_1 deta_1/dxi_1 =
// dphi/deta_0 (1+eta_0)/(1-eta_1) + dphi/deta_1
//
// and so the full inner products are
//
// (dphi/dx,in[0]) + (dphi/dy,in[1])
// = (dphi/deta_0, ((2/(1-eta_1) (dxi_0/dx in[0]+dxi_0/dy in[1])
// + (1_eta_0)/(1-eta_1) (dxi_1/dx in[0]+dxi_1/dy in[1]))
// + (dphi/deta_1, (dxi_1/dx in[0] + dxi_1/dy in[1]))
for(int i = 0; i < 2; ++i)
{
Vmath::Vmul (ntot,m_derivFac[i],1, in[0],1, tmp[i],1);
......@@ -1017,6 +1028,60 @@ class IProductWRTDerivBase_SumFac_Tet : public Operator
public:
OPERATOR_CREATE(IProductWRTDerivBase_SumFac_Tet)
/**
* This method calculates:
*
* \f[ (d\phi/dx,in[0]) + (d\phi/dy,in[1]) + (d\phi/dz,in[2]) \f]
*
* which can be represented in terms of local cartesian
* derivaties as:
*
* \f[ ((d\phi/d\xi_0\, d\xi_0/dx +
* d\phi/d\xi_1\, d\xi_1/dx +
* d\phi/d\xi_2\, d\xi_2/dx),in[0]) + \f]
*
* \f[ ((d\phi/d\xi_0\, d\xi_0/dy +
* d\phi/d\xi_1\, d\xi_1/dy +
* d\phi/d\xi_2\, d\xi_2/dy),in[1]) + \f]
*
* \f[ ((d\phi/d\xi_0\, d\xi_0/dz +
* d\phi/d\xi_1\, d\xi_1/dz +
* d\phi/d\xi_2\, d\xi_2/dz),in[2]) \, \f]
*
* where we note that
*
* \f[ d\phi/d\xi_0 = d\phi/d\eta_0 4/((1-\eta_1)(1-\eta_2)) /f]
*
* \f[ d\phi/d\xi_1 = d\phi/d\eta_0 2(1+\eta_0)/((1-\eta_1)(1-\eta_2))
* + d\phi/d\eta_1 2/(1-\eta_2) \f]
*
* \f[ d\phi/d\xi_2 = d\phi/d\eta_0 2(1+\eta_0)/((1-\eta_1)(1-\eta_2))
* + d\phi/d\eta_1 (1+\eta_1)/(1-\eta_2) + d\phi/d\eta_2 \f]
*
* and so the full inner products are
*
* \f[ (d\phi/dx,in[0]) + (d\phi/dy,in[1]) + (d\phi/dz,in[2]) = \f]
*
* \f[ (d\phi/d\eta_0, fac0 (tmp0 + fac1(tmp1 + tmp2)))
* + (d\phi/d\eta_1, fac2 (tmp1 + fac3 tmp2))
* + (d\phi/d\eta_2, tmp2) \f]
*
* where
*
* \f[ \begin{array}{lcl}
* tmp0 &=& (d\xi_0/dx in[0] + d\xi_0/dy in[1] + d\xi_0/dz in[2]) \\
* tmp1 &=& (d\xi_1/dx in[0] + d\xi_1/dy in[1] + d\xi_1/dz in[2]) \\
* tmp2 &=& (d\xi_2/dx in[0] + d\xi_2/dy in[1] + d\xi_2/dz in[2])
* \end{array} \f]
*
* \f[ \begin{array}{lcl}
* fac0 &= & 4/((1-\eta_1)(1-\eta_2)) \\
* fac1 &= & (1+\eta_0)/2 \\
* fac2 &= & 2/(1-\eta_2) \\
* fac3 &= & (1+\eta_1)/2 \end{array} \f]
*
*/
virtual void operator()(
const Array<OneD, const NekDouble> &entry0,
Array<OneD, NekDouble> &entry1,
......@@ -1042,44 +1107,6 @@ class IProductWRTDerivBase_SumFac_Tet : public Operator
tmp[i] = wsp + i*nmax;
}
// calculate (dphi/dx,in[0]) = ((dphi/dxi_0 dxi_0/dx +
// dphi/dxi_1 dxi_1/dx +
// dphi/dxi_2 dxi_2/dx),in[0])
// + (dphi/dy,in[1]) = ((dphi/dxi_0 dxi_0/dy +
// dphi/dxi_1 dxi_1/dy +
// dphi/dxi_2 dxi_2/dy),in[1])
// + (dphi/dz,in[2]) = ((dphi/dxi_0 dxi_0/dz +
// dphi/dxi_1 dxi_1/dz +
// dphi/dxi_2 dxi_2/dz),in[1])
//
// Note dphi/dxi_0 =
// dphi/deta_0 4/((1-eta_1)(1-eta2))
//
// dphi/dxi_1 =
// dphi/deta_0 2(1+eta_0)/((1-eta_1)(1-eta_2)) +
// dphi/deta_1 2/(1-eta_2)
//
// dphi/dxi_2 =
// dphi/deta_0 2(1+eta_0)/((1-eta_1)(1-eta_2)) +
// dphi/deta_1 (1+eta_1)/(1-eta_2) + dphi/deta_2
//
// and so the full inner products are
//
// (dphi/dx,in[0]) + (dphi/dy,in[1]) + (dphi/dz,in[2])
// = (dphi/deta_0, fac0 (tmp0 + fac1(tmp1 + tmp2)))
// + (dphi/deta_1, fac2 (tmp1 + fac3 tmp2))
// + (dphi/deta_2, tmp2)
//
// tmp0 = (dxi_0/dx in[0] + dxi_0/dy in[1] + dxi_0/dz in[2])
// tmp1 = (dxi_1/dx in[0] + dxi_1/dy in[1] + dxi_1/dz in[2])
// tmp2 = (dxi_2/dx in[0] + dxi_2/dy in[1] + dxi_2/dz in[2])
// fac0 = 4/((1-eta_1)(1-eta2))
// fac1 = (1+eta_0)/2
// fac2 = 2/(1-eta_2)
// fac3 = (1+eta_1)/2
for(int i = 0; i < 3; ++i)
{
Vmath::Vmul (ntot,m_derivFac[i],1, in[0],1,tmp[i],1);
......@@ -1259,6 +1286,52 @@ class IProductWRTDerivBase_SumFac_Prism : public Operator
{
}
/**
* This method calculates:
*
* \f[ (d\phi/dx,in[0]) + (d\phi/dy,in[1]) + (d\phi/dz,in[2]) \f]
*
* which can be represented in terms of local cartesian
* derivaties as:
*
* \f[ ((d\phi/d\xi_0\, d\xi_0/dx +
* d\phi/d\xi_1\, d\xi_1/dx +
* d\phi/d\xi_2\, d\xi_2/dx),in[0]) + \f]
*
* \f[ ((d\phi/d\xi_0\, d\xi_0/dy +
* d\phi/d\xi_1\, d\xi_1/dy +
* d\phi/d\xi_2\, d\xi_2/dy),in[1]) + \f]
*
* \f[ ((d\phi/d\xi_0\, d\xi_0/dz +
* d\phi/d\xi_1\, d\xi_1/dz +
* d\phi/d\xi_2\, d\xi_2/dz),in[2]) \, \f]
*
* where we note that
*
* \f[ d\phi/d\xi_0 =
* d\phi/d\eta_0 d\eta_0/d\xi_0 = d\phi/d\eta_0 2/(1-\eta_2) \f]
*
* \f[ d\phi/d\xi_2 =
* d\phi/d\eta_0 d\eta_0/d\xi_2 + d\phi/d\eta_2 d\eta_2/d\xi_2 =
* d\phi/d\eta_0 (1+\eta_0)/(1-\eta_2) + d\phi/d\eta_2 \f]
*
*
* and so the full inner products are
*
* \f[ (d\phi/dx,in[0]) + (d\phi/dy,in[1]) + (d\phi/dz,in[2]) = \f]
*
* \f[ (d\phi/d\eta_0, ((2/(1-\eta_2) (d\xi_0/dx in[0] + d\xi_0/dy in[1]
* + d\xi_0/dz in[2])
* + (1-\eta_0)/(1-\eta_2) (d\xi_2/dx in[0] + d\xi_2/dy in[1]
* + d\xi_2/dz in[2] )) + \f]
*
* \f[ (d\phi/d\eta_1, (d\xi_1/dx in[0] + d\xi_1/dy in[1]
* + d\xi_1/dz in[2])) + \f]
*
* \f[ (d\phi/d\eta_2, (d\xi_2/dx in[0] + d\xi_2/dy in[1]
* + d\xi_2/dz in[2])) \f]
*
*/
virtual void operator()(
const Array<OneD, const NekDouble> &entry0,
Array<OneD, NekDouble> &entry1,
......@@ -1284,32 +1357,6 @@ class IProductWRTDerivBase_SumFac_Prism : public Operator
tmp[i] = wsp + i*nmax;
}
// calculate (dphi/dx,in[0]) = ((dphi/dxi_0 dxi_0/dx +
// dphi/dxi_1 dxi_1/dx),in[0])
// + (dphi/dy,in[1]) = ((dphi/dxi_0 dxi_0/dy +
// dphi/dxi_1 dxi_1/dy),in[1])
// + (dphi/dz,in[2]) = ((dphi/dxi_0 dxi_0/dz +
// dphi/dxi_1 dxi_1/dz),in[2])
//
// Note dphi/dxi_0 =
// dphi/deta_0 deta_0/dxi_0 = dphi/deta_0 2/(1-eta_2)
//
// dphi/dxi_2 =
// dphi/deta_0 deta_0/dxi_2 + dphi/deta_2 deta_2/dxi_2 =
// dphi/deta_0 (1+eta_0)/(1-eta_2) + dphi/deta_2
//
// and so the full inner products are
//
// (dphi/dx,in[0]) + (dphi/dy,in[1]) + (dphi/dz,in[2])
// = (dphi/deta_0, ((2/(1-eta_2) (dxi_0/dx in[0] + dxi_0/dy in[1]
// + dxi_0/dz in[2])
// + (1_eta_0)/(1-eta_2) (dxi_2/dx in[0] + dxi_2/dy in[1]
// + dxi_2/dz in[2] ))
// + (dphi/deta_1, (dxi_1/dx in[0] + dxi_1/dy in[1]
// + dxi_1/dz in[2]))
// + (dphi/deta_2, (dxi_2/dx in[0] + dxi_2/dy in[1]
// + dxi_2/dz in[2]))
for(int i = 0; i < 3; ++i)
{
Vmath::Vmul (ntot,m_derivFac[i],1, in[0],1,
......@@ -1467,7 +1514,61 @@ class IProductWRTDerivBase_SumFac_Pyr : public Operator
{
}
virtual void operator()(
/**
* This method calculates:
*
* \f[ (d\phi/dx,in[0]) + (d\phi/dy,in[1]) + (d\phi/dz,in[2]) \f]
*
* which can be represented in terms of local cartesian
* derivaties as:
*
* \f[ ((d\phi/d\xi_0\, d\xi_0/dx +
* d\phi/d\xi_1\, d\xi_1/dx +
* d\phi/d\xi_2\, d\xi_2/dx),in[0]) + \f]
*
* \f[ ((d\phi/d\xi_0\, d\xi_0/dy +
* d\phi/d\xi_1\, d\xi_1/dy +
* d\phi/d\xi_2\, d\xi_2/dy),in[1]) + \f]
*
* \f[ ((d\phi/d\xi_0\, d\xi_0/dz +
* d\phi/d\xi_1\, d\xi_1/dz +
* d\phi/d\xi_2\, d\xi_2/dz),in[2]) \, \f]
*
* where we note that
*
* \f[ d\phi/d\xi_0 =
* d\phi/d\eta_0\, d\eta_0/d\xi_0 =
* d\phi/d\eta_0\, 2/(1-\eta_2). \f]
*
* \f[ d\phi/d\xi_1 =
* d\phi/d\eta_1\, d\eta_1/d\xi_1 =
* d\phi/d\eta_1\, 2/(1-\eta_2) \f]
*
* \f[ d\phi/d\xi_2 =
* d\phi/d\eta_0\, d\eta_0/d\xi_2 +
* d\phi/d\eta_1\, d\eta_1/d\xi_2 +
* d\phi/d\eta_2\, d\eta_2/d\xi_2 =
* d\phi/d\eta_0 (1+\eta_0)/(1-\eta_2) +
* d\phi/d\eta_1 (1+\eta_1)/(1-\eta_2) + d\phi/d\eta_2 \f]
*
* and so the full inner products are
*
* \f[ (d\phi/dx,in[0]) + (d\phi/dy,in[1]) + (d\phi/dz,in[2]) = \f]
*
* \f[ (d\phi/d\eta_0, ((2/(1-\eta_2) (d\xi_0/dx in[0] +
* d\xi_0/dy in[1] +
* (1-\eta_0)/(1-\eta_2) (d\xi_2/dx in[0] + d\xi_2/dy in[1]
* + d\xi_2/dz in[2] )) + \f]
* \f[ (d\phi/d\eta_1, ((2/(1-\eta_2) (d\xi_1/dx in[0] +
* d\xi_0/dy in[1] + d\xi_0/dz in[2]) +
* (1-\eta_1)/(1-\eta_2) (d\xi_2/dx in[0] + d\xi_2/dy in[1] +
* d\xi_2/dz in[2] )) \f]
*
* \f[ (d\phi/d\eta_2, (d\xi_2/dx in[0] + d\xi_2/dy in[1] +
* d\xi_2/dz in[2])) \f]
*/
virtual void operator()(
const Array<OneD, const NekDouble> &entry0,
Array<OneD, NekDouble> &entry1,
Array<OneD, NekDouble> &entry2,
......@@ -1491,40 +1592,6 @@ class IProductWRTDerivBase_SumFac_Pyr : public Operator
{
tmp[i] = wsp + i*nmax;
}
// calculate (dphi/dx,in[0]) = ((dphi/dxi_0 dxi_0/dx +
// dphi/dxi_1 dxi_1/dx),in[0])
// + (dphi/dy,in[1]) = ((dphi/dxi_0 dxi_0/dy +
// dphi/dxi_1 dxi_1/dy),in[1])
// + (dphi/dz,in[2]) = ((dphi/dxi_0 dxi_0/dz +
// dphi/dxi_1 dxi_1/dz),in[2])
//
// Note dphi/dxi_0 =
// dphi/deta_0 deta_0/dxi_0 = dphi/deta_0 2/(1-eta_2)
//
// Note dphi/dxi_1 =
// dphi/deta_1 deta_1/dxi_1 = dphi/deta_1 2/(1-eta_2)
//
// dphi/dxi_2 =
// dphi/deta_0 deta_0/dxi_2 + dphi/deta_1 deta_1/dxi_2 +
// + dphi/deta_2 deta_2/dxi_2 =
// dphi/deta_0 (1+eta_0)/(1-eta_2) +
// dphi/deta_1 (1+eta_1)/(1-eta_2) + dphi/deta_2
//
// and so the full inner products are
//
// (dphi/dx,in[0]) + (dphi/dy,in[1]) + (dphi/dz,in[2])
// = (dphi/deta_0, ((2/(1-eta_2) (dxi_0/dx in[0] + dxi_0/dy in[1]
// + dxi_0/dz in[2])
// + (1_eta_0)/(1-eta_2) (dxi_2/dx in[0] + dxi_2/dy in[1]
// + dxi_2/dz in[2] ))
// + (dphi/deta_1, ((2/(1-eta_2) (dxi_1/dx in[0] + dxi_0/dy in[1]
// + dxi_0/dz in[2])
// + (1_eta_1)/(1-eta_2) (dxi_2/dx in[0] + dxi_2/dy in[1]
// + dxi_2/dz in[2] ))
// + (dphi/deta_2, (dxi_2/dx in[0] + dxi_2/dy in[1]
// + dxi_2/dz in[2]))
for(int i = 0; i < 3; ++i)
{
......
library/Collections/PhysDeriv.cpp
View file @ 9bac580d
......@@ -938,7 +938,7 @@ class PhysDeriv_SumFac_Hex : public Operator
Array<OneD, Array<OneD, NekDouble> > out(3);
out[0] = output0; out[1] = output1; out[2] = output2;
for(int i = 0; i < m_dim; ++i)
for(int i = 0; i < 3; ++i)
{
Diff[i] = wsp + i*ntot;
}
......@@ -967,10 +967,10 @@ class PhysDeriv_SumFac_Hex : public Operator
// calculate full derivative
for(int i = 0; i < m_coordim; ++i)
{
Vmath::Vmul(ntot,m_derivFac[i*m_dim],1,Diff[0],1,out[i],1);
for(int j = 1; j < m_dim; ++j)
Vmath::Vmul(ntot,m_derivFac[i*3],1,Diff[0],1,out[i],1);
for(int j = 1; j < 3; ++j)
{
Vmath::Vvtvp (ntot, m_derivFac[i*m_dim+j], 1,
Vmath::Vvtvp (ntot, m_derivFac[i*3+j], 1,
Diff[j], 1,
out[i], 1,
out[i], 1);
......@@ -989,7 +989,7 @@ class PhysDeriv_SumFac_Hex : public Operator
Array<OneD, NekDouble> tmp0,tmp1,tmp2;
Array<OneD, Array<OneD, NekDouble> > Diff(3);
for(int i = 0; i < m_dim; ++i)
for(int i = 0; i < 3; ++i)
{
Diff[i] = wsp + i*ntot;
}
......@@ -1016,10 +1016,10 @@ class PhysDeriv_SumFac_Hex : public Operator
}
// calculate full derivative
Vmath::Vmul(ntot,m_derivFac[dir*m_dim],1,Diff[0],1,output,1);
for(int j = 1; j < m_dim; ++j)
Vmath::Vmul(ntot,m_derivFac[dir*3],1,Diff[0],1,output,1);
for(int j = 1; j < 3; ++j)
{
Vmath::Vvtvp (ntot, m_derivFac[dir*m_dim+j], 1,
Vmath::Vvtvp (ntot, m_derivFac[dir*3+j], 1,
Diff[j], 1,
output, 1,
output, 1);
......@@ -1028,7 +1028,6 @@ class PhysDeriv_SumFac_Hex : public Operator
protected:
Array<TwoD, const NekDouble> m_derivFac;
int m_dim;
int m_coordim;
const int m_nquad0;
const int m_nquad1;
......@@ -1048,7 +1047,6 @@ class PhysDeriv_SumFac_Hex : public Operator
{
LibUtilities::PointsKeyVector PtsKey = m_stdExp->GetPointsKeys();
m_dim = PtsKey.size();
m_coordim = m_stdExp->GetCoordim();
m_derivFac = pGeomData->GetDerivFactors(pCollExp);
......@@ -1096,7 +1094,7 @@ class PhysDeriv_SumFac_Tet : public Operator
Array<OneD, Array<OneD, NekDouble> > out(3);
out[0] = output0; out[1] = output1; out[2] = output2;
for(int i = 0; i < m_dim; ++i)
for(int i = 0; i < 3; ++i)
{
Diff[i] = wsp + i*ntot;
}
......@@ -1161,10 +1159,10 @@ class PhysDeriv_SumFac_Tet : public Operator
// calculate full derivative
for(int i = 0; i < m_coordim; ++i)
{
Vmath::Vmul(ntot,m_derivFac[i*m_dim],1,Diff[0],1,out[i],1);
for(int j = 1; j < m_dim; ++j)
Vmath::Vmul(ntot,m_derivFac[i*3],1,Diff[0],1,out[i],1);
for(int j = 1; j < 3; ++j)
{
Vmath::Vvtvp (ntot, m_derivFac[i*m_dim+j], 1,
Vmath::Vvtvp (ntot, m_derivFac[i*3+j], 1,
Diff[j], 1, out[i], 1, out[i], 1);
}
}
......@@ -1181,7 +1179,7 @@ class PhysDeriv_SumFac_Tet : public Operator
Array<OneD, NekDouble> tmp0,tmp1,tmp2;
Array<OneD, Array<OneD, NekDouble> > Diff(3);
for(int i = 0; i < m_dim; ++i)
for(int i = 0; i < 3; ++i)
{
Diff[i] = wsp + i*ntot;
}
......@@ -1244,17 +1242,16 @@ class PhysDeriv_SumFac_Tet : public Operator
}
// calculate full derivative
Vmath::Vmul(ntot,m_derivFac[dir*m_dim],1,Diff[0],1,output,1);
for(int j = 1; j < m_dim; ++j)
Vmath::Vmul(ntot,m_derivFac[dir*3],1,Diff[0],1,output,1);
for(int j = 1; j < 3; ++j)
{
Vmath::Vvtvp (ntot, m_derivFac[dir*m_dim+j], 1,
Vmath::Vvtvp (ntot, m_derivFac[dir*3+j], 1,
Diff[j], 1, output, 1, output, 1);
}
}
protected:
Array<TwoD, const NekDouble> m_derivFac;
int m_dim;
int m_coordim;
const int m_nquad0;
const int m_nquad1;
......@@ -1278,7 +1275,6 @@ class PhysDeriv_SumFac_Tet : public Operator
{
LibUtilities::PointsKeyVector PtsKey = m_stdExp->GetPointsKeys();
m_dim = PtsKey.size();
m_coordim = m_stdExp->GetCoordim();
m_derivFac = pGeomData->GetDerivFactors(pCollExp);
......@@ -1358,7 +1354,7 @@ class PhysDeriv_SumFac_Prism : public Operator
Array<OneD, Array<OneD, NekDouble> > out(3);
out[0] = output0; out[1] = output1; out[2] = output2;
for(int i = 0; i < m_dim; ++i)
for(int i = 0; i < 3; ++i)
{
Diff[i] = wsp + i*ntot;
}
......@@ -1401,10 +1397,10 @@ class PhysDeriv_SumFac_Prism : public Operator
// calculate full derivative
for(int i = 0; i < m_coordim; ++i)
{
Vmath::Vmul(ntot,m_derivFac[i*m_dim],1,Diff[0],1,out[i],1);
for(int j = 1; j < m_dim; ++j)
Vmath::Vmul(ntot,m_derivFac[i*3],1,Diff[0],1,out[i],1);
for(int j = 1; j < 3; ++j)
{
Vmath::Vvtvp (ntot, m_derivFac[i*m_dim+j], 1,
Vmath::Vvtvp (ntot, m_derivFac[i*3+j], 1,
Diff[j], 1, out[i], 1, out[i], 1);
}
}
......@@ -1421,7 +1417,7 @@ class PhysDeriv_SumFac_Prism : public Operator
Array<OneD, NekDouble> tmp0,tmp1,tmp2;
Array<OneD, Array<OneD, NekDouble> > Diff(3);
for(int i = 0; i < m_dim; ++i)
for(int i = 0; i < 3; ++i)
{
Diff[i] = wsp + i*ntot;
}
......@@ -1462,17 +1458,16 @@ class PhysDeriv_SumFac_Prism : public Operator
}
// calculate full derivative
Vmath::Vmul(ntot,m_derivFac[dir*m_dim],1,Diff[0],1,output,1);
for(int j = 1; j < m_dim; ++j)
Vmath::Vmul(ntot,m_derivFac[dir*3],1,Diff[0],1,output,1);
for(int j = 1; j < 3; ++j)
{
Vmath::Vvtvp (ntot, m_derivFac[dir*m_dim+j], 1,
Vmath::Vvtvp (ntot, m_derivFac[dir*3+j], 1,
Diff[j], 1, output, 1, output, 1);
}
}
protected:
Array<TwoD, const NekDouble> m_derivFac;
int m_dim;
int m_coordim;
const int m_nquad0;
const int m_nquad1;
......@@ -1494,7 +1489,6 @@ class PhysDeriv_SumFac_Prism : public Operator
{
LibUtilities::PointsKeyVector PtsKey = m_stdExp->GetPointsKeys();
m_dim = PtsKey.size();
m_coordim = m_stdExp->GetCoordim();
m_derivFac = pGeomData->GetDerivFactors(pCollExp);
......@@ -1563,7 +1557,7 @@ class PhysDeriv_SumFac_Pyr : public Operator
Array<OneD, Array<OneD, NekDouble> > out(3);
out[0] = output0; out[1] = output1; out[2] = output2;
for(int i = 0; i < m_dim; ++i)
for(int i = 0; i < 3; ++i)
{
Diff[i] = wsp + i*ntot;
}
......@@ -1613,10 +1607,10 @@ class PhysDeriv_SumFac_Pyr : public Operator
// calculate full derivative
for(int i = 0; i < m_coordim; ++i)
{
Vmath::Vmul(ntot,m_derivFac[i*m_dim],1,Diff[0],1,out[i],1);
for(int j = 1; j < m_dim; ++j)
Vmath::Vmul(ntot,m_derivFac[i*3],1,Diff[0],1,out[i],1);
for(int j = 1; j < 3; ++j)
{
Vmath::Vvtvp (ntot, m_derivFac[i*m_dim+j], 1,
Vmath::Vvtvp (ntot, m_derivFac[i*3+j], 1,
Diff[j], 1, out[i], 1, out[i], 1);
}
}
......@@ -1633,7 +1627,7 @@ class PhysDeriv_SumFac_Pyr : public Operator
Array<OneD, NekDouble> tmp0,tmp1,tmp2;
Array<OneD, Array<OneD, NekDouble> > Diff(3);
for(int i = 0; i < m_dim; ++i)
for(int i = 0; i < 3; ++i)
{
Diff[i] = wsp + i*ntot;
}
......@@ -1680,17 +1674,16 @@ class PhysDeriv_SumFac_Pyr : public Operator
}
// calculate full derivative
Vmath::Vmul(ntot,m_derivFac[dir*m_dim],1,Diff[0],1,output,1);
for(int j = 1; j < m_dim; ++j)
Vmath::Vmul(ntot,m_derivFac[dir*3],1,Diff[0],1,output,1);
for(int j = 1; j < 3; ++j)
{
Vmath::Vvtvp (ntot, m_derivFac[dir*m_dim+j], 1,
Vmath::Vvtvp (ntot, m_derivFac[dir*3+j], 1,
Diff[j], 1, output, 1, output, 1);
}
}
protected: