Commit 80ab2d10 authored by Chris Cantwell's avatar Chris Cantwell
Browse files

Removed redundant Laplacian routines in 2D and 3D elements.

parent 349dcbe4
......@@ -1669,113 +1669,6 @@ namespace Nektar
}
/**
* @param inarray Input coefficients.
* @param output Output coefficients.
* @param mkey Matrix key
*/
void HexExp::v_LaplacianMatrixOp_MatFree(
const Array<OneD, const NekDouble> &inarray,
Array<OneD,NekDouble> &outarray,
const StdRegions::StdMatrixKey &mkey)
{
if(mkey.GetNVarCoeff() == 0)
{
// This implementation is only valid when there are no
// coefficients associated to the Laplacian operator
int nqtot = GetTotPoints();
const Array<OneD, const NekDouble>& base0 = m_base[0]->GetBdata();
const Array<OneD, const NekDouble>& base1 = m_base[1]->GetBdata();
const Array<OneD, const NekDouble>& base2 = m_base[2]->GetBdata();
// Allocate temporary storage
Array<OneD,NekDouble> wsp0(7*nqtot);
Array<OneD,NekDouble> wsp1(wsp0+nqtot);
if(!(m_base[0]->Collocation() && m_base[1]->Collocation() &&
m_base[2]->Collocation()))
{
// LAPLACIAN MATRIX OPERATION
// wsp0 = u = B * u_hat
// wsp1 = du_dxi1 = D_xi1 * wsp0 = D_xi1 * u
// wsp2 = du_dxi2 = D_xi2 * wsp0 = D_xi2 * u
BwdTrans_SumFacKernel(base0,base1,base2,inarray,wsp0,wsp1,true,true,true);
LaplacianMatrixOp_MatFree_Kernel(wsp0,outarray,wsp1);
}
else
{
LaplacianMatrixOp_MatFree_Kernel(inarray,outarray,wsp1);
}
}
else
{
StdExpansion::LaplacianMatrixOp_MatFree_GenericImpl(inarray,outarray,mkey);
}
}
void HexExp::v_HelmholtzMatrixOp_MatFree(
const Array<OneD, const NekDouble> &inarray,
Array<OneD,NekDouble> &outarray,
const StdRegions::StdMatrixKey &mkey)
{
if(mkey.GetNVarCoeff() == 0)
{
int nquad0 = m_base[0]->GetNumPoints();
int nquad1 = m_base[1]->GetNumPoints();
int nquad2 = m_base[2]->GetNumPoints();
int nqtot = nquad0*nquad1*nquad2;
int nmodes0 = m_base[0]->GetNumModes();
int nmodes1 = m_base[1]->GetNumModes();
int nmodes2 = m_base[2]->GetNumModes();
int wspsize = max(nquad0*nmodes2*(nmodes1+nquad1),
nquad0*nquad1*(nquad2+nmodes0)+
nmodes0*nmodes1*nquad2);
NekDouble lambda = mkey.GetConstFactor(StdRegions::eFactorLambda);
const Array<OneD, const NekDouble>& base0 = m_base[0]->GetBdata ();
const Array<OneD, const NekDouble>& base1 = m_base[1]->GetBdata ();
const Array<OneD, const NekDouble>& base2 = m_base[2]->GetBdata ();
Array<OneD,NekDouble> wsp0(8*wspsize);
Array<OneD,NekDouble> wsp1(wsp0+1*wspsize);
Array<OneD,NekDouble> wsp2(wsp0+2*wspsize);
if(!(m_base[0]->Collocation() && m_base[1]->Collocation() &&
m_base[2]->Collocation()))
{
// MASS MATRIX OPERATION
// The following is being calculated:
// wsp0 = B * u_hat = u
// wsp1 = W * wsp0
// outarray = B^T * wsp1 = B^T * W * B * u_hat = M * u_hat
BwdTrans_SumFacKernel (base0,base1,base2,inarray,
wsp0,wsp2,true,true,true);
MultiplyByQuadratureMetric (wsp0,wsp1);
IProductWRTBase_SumFacKernel (base0,base1,base2,wsp1,
outarray,wsp2,true,true,true);
LaplacianMatrixOp_MatFree_Kernel(wsp0,wsp1,wsp2);
}
else
{
// specialised implementation for the classical spectral
// element method
MultiplyByQuadratureMetric (inarray,outarray);
LaplacianMatrixOp_MatFree_Kernel(inarray,wsp1,wsp2);
}
// outarray = lambda * outarray + wsp1
// = (lambda * M + L ) * u_hat
Vmath::Svtvp(m_ncoeffs,lambda,&outarray[0],1,&wsp1[0],1,
&outarray[0],1);
}
else
{
StdExpansion::HelmholtzMatrixOp_MatFree_GenericImpl(inarray,outarray,mkey);
}
}
//-----------------------------
......@@ -2247,7 +2140,7 @@ namespace Nektar
m_staticCondMatrixManager.DeleteObject(mkey);
}
void HexExp::LaplacianMatrixOp_MatFree_Kernel(
void HexExp::v_LaplacianMatrixOp_MatFree_Kernel(
const Array<OneD, const NekDouble> &inarray,
Array<OneD, NekDouble> &outarray,
Array<OneD, NekDouble> &wsp)
......
......@@ -247,16 +247,6 @@ namespace Nektar
Array<OneD,NekDouble> &outarray,
const StdRegions::StdMatrixKey &mkey);
LOCAL_REGIONS_EXPORT virtual void v_LaplacianMatrixOp_MatFree(
const Array<OneD, const NekDouble> &inarray,
Array<OneD,NekDouble> &outarray,
const StdRegions::StdMatrixKey &mkey);
LOCAL_REGIONS_EXPORT virtual void v_HelmholtzMatrixOp_MatFree(
const Array<OneD, const NekDouble> &inarray,
Array<OneD,NekDouble> &outarray,
const StdRegions::StdMatrixKey &mkey);
//---------------------------------------
// Matrix creation functions
......@@ -292,7 +282,7 @@ namespace Nektar
HexExp();
LOCAL_REGIONS_EXPORT void LaplacianMatrixOp_MatFree_Kernel(
virtual void v_LaplacianMatrixOp_MatFree_Kernel(
const Array<OneD, const NekDouble> &inarray,
Array<OneD, NekDouble> &outarray,
Array<OneD, NekDouble> &wsp);
......
......@@ -1300,48 +1300,6 @@ namespace Nektar
StdExpansion::LaplacianMatrixOp_MatFree(k1,k2,inarray,outarray,mkey);
}
void PrismExp::v_LaplacianMatrixOp_MatFree(
const Array<OneD, const NekDouble> &inarray,
Array<OneD, NekDouble> &outarray,
const StdRegions::StdMatrixKey &mkey)
{
if(mkey.GetNVarCoeff() == 0)
{
// This implementation is only valid when there are no
// coefficients associated to the Laplacian operator
// if(m_metricinfo->IsUsingLaplMetrics())
// {
// ASSERTL0(false,"Finish implementing HexExp for Lap metrics");
// // Get this from HexExp
// }
// else
// {
int nquad0 = m_base[0]->GetNumPoints();
int nquad1 = m_base[1]->GetNumPoints();
int nquad2 = m_base[2]->GetNumPoints();
int nmodes0 = m_base[0]->GetNumModes ();
int nmodes1 = m_base[1]->GetNumModes ();
int nqtot = nquad0*nquad1*nquad2;
const Array<OneD, const NekDouble>& base0 = m_base[0]->GetBdata ();
const Array<OneD, const NekDouble>& base1 = m_base[1]->GetBdata ();
const Array<OneD, const NekDouble>& base2 = m_base[2]->GetBdata ();
Array<OneD,NekDouble> wsp (nquad2*nmodes0*(nquad1+nmodes1));
Array<OneD,NekDouble> wsp1(nqtot);
// Backwards transform to obtain u = B * u_hat.
BwdTrans_SumFacKernel (base0,base1,base2,inarray,wsp1,wsp,true,true,true);
LaplacianMatrixOp_Kernel(wsp1, outarray, wsp);
// }
}
else
{
StdExpansion::LaplacianMatrixOp_MatFree_GenericImpl(
inarray,outarray,mkey);
}
}
void PrismExp::v_HelmholtzMatrixOp(
const Array<OneD, const NekDouble> &inarray,
Array<OneD, NekDouble> &outarray,
......@@ -1372,52 +1330,6 @@ namespace Nektar
}
}
void PrismExp::v_HelmholtzMatrixOp_MatFree(
const Array<OneD, const NekDouble> &inarray,
Array<OneD, NekDouble> &outarray,
const StdRegions::StdMatrixKey &mkey)
{
// if(m_metricinfo->IsUsingLaplMetrics())
// {
// ASSERTL0(false,"Finish implementing PrismExp Helmholtz for Lapl Metrics");
// }
// else
// {
int nquad0 = m_base[0]->GetNumPoints();
int nquad1 = m_base[1]->GetNumPoints();
int nquad2 = m_base[2]->GetNumPoints();
int nmodes0 = m_base[0]->GetNumModes ();
int nmodes1 = m_base[1]->GetNumModes ();
int nqtot = nquad0*nquad1*nquad2;
const Array<OneD, const NekDouble>& base0 = m_base[0]->GetBdata ();
const Array<OneD, const NekDouble>& base1 = m_base[1]->GetBdata ();
const Array<OneD, const NekDouble>& base2 = m_base[2]->GetBdata ();
Array<OneD,NekDouble> wsp (nquad2*nmodes0*(nquad1+nmodes1));
Array<OneD,NekDouble> wsp0(nqtot);
Array<OneD,NekDouble> wsp1(nqtot);
NekDouble lambda = mkey.GetConstFactor(StdRegions::eFactorLambda);
// MASS MATRIX OPERATION
// The following is being calculated:
// wsp0 = B * u_hat = u
// wsp1 = W * wsp0
// outarray = B^T * wsp1 = B^T * W * B * u_hat = M * u_hat
BwdTrans_SumFacKernel (base0,base1,base2,inarray,
wsp0,wsp,true,true,true);
MultiplyByQuadratureMetric (wsp0,wsp1);
IProductWRTBase_SumFacKernel(base0,base1,base2,wsp1,
outarray,wsp,true,true,true);
LaplacianMatrixOp_Kernel (wsp0,wsp1,wsp);
// outarray = lambda * outarray + wsp1
// = (lambda * M + L ) * u_hat
Vmath::Svtvp(m_ncoeffs,lambda,&outarray[0],1,&wsp1[0],1,
&outarray[0],1);
// }
}
//---------------------------------------
// Matrix creation functions
......@@ -1898,7 +1810,7 @@ namespace Nektar
*
* @see %TetExp::v_HelmholtzMatrixOp_MatFree
*/
void PrismExp::LaplacianMatrixOp_Kernel(
void PrismExp::v_LaplacianMatrixOp_MatFree_Kernel(
const Array<OneD, const NekDouble> &inarray,
Array<OneD, NekDouble> &outarray,
Array<OneD, NekDouble> &wsp)
......
......@@ -188,14 +188,6 @@ namespace Nektar
const Array<OneD, const NekDouble> &inarray,
Array<OneD,NekDouble> &outarray,
const StdRegions::StdMatrixKey &mkey);
LOCAL_REGIONS_EXPORT virtual void v_LaplacianMatrixOp_MatFree(
const Array<OneD, const NekDouble> &inarray,
Array<OneD,NekDouble> &outarray,
const StdRegions::StdMatrixKey &mkey);
LOCAL_REGIONS_EXPORT virtual void v_HelmholtzMatrixOp_MatFree(
const Array<OneD, const NekDouble> &inarray,
Array<OneD, NekDouble> &outarray,
const StdRegions::StdMatrixKey &mkey);
//---------------------------------------
// Matrix creation functions
......@@ -219,7 +211,7 @@ namespace Nektar
LibUtilities::NekManager<MatrixKey, DNekScalMat, MatrixKey::opLess> m_matrixManager;
LibUtilities::NekManager<MatrixKey, DNekScalBlkMat, MatrixKey::opLess> m_staticCondMatrixManager;
LOCAL_REGIONS_EXPORT void LaplacianMatrixOp_Kernel(
virtual void v_LaplacianMatrixOp_MatFree_Kernel(
const Array<OneD, const NekDouble> &inarray,
Array<OneD, NekDouble> &outarray,
Array<OneD, NekDouble> &wsp);
......
......@@ -2429,115 +2429,7 @@ namespace Nektar
}
void QuadExp::v_LaplacianMatrixOp_MatFree(
const Array<OneD, const NekDouble> &inarray,
Array<OneD,NekDouble> &outarray,
const StdRegions::StdMatrixKey &mkey)
{
if (mkey.GetNVarCoeff() == 0
&&!mkey.ConstFactorExists(StdRegions::eFactorSVVCutoffRatio))
{
// This implementation is only valid when there are no
// coefficients associated to the Laplacian operator
int nquad0 = m_base[0]->GetNumPoints();
int nquad1 = m_base[1]->GetNumPoints();
int nqtot = nquad0*nquad1;
int nmodes0 = m_base[0]->GetNumModes();
int nmodes1 = m_base[1]->GetNumModes();
int wspsize = max(max(max(nqtot,m_ncoeffs),nquad1*nmodes0),nquad0*nmodes1);
const Array<OneD, const NekDouble>& base0 = m_base[0]->GetBdata();
const Array<OneD, const NekDouble>& base1 = m_base[1]->GetBdata();
// Allocate temporary storage
Array<OneD,NekDouble> wsp0(4*wspsize); // size wspsize
Array<OneD,NekDouble> wsp1(wsp0+wspsize); // size 3*wspsize
if(!(m_base[0]->Collocation() && m_base[1]->Collocation()))
{
// LAPLACIAN MATRIX OPERATION
// wsp0 = u = B * u_hat
// wsp1 = du_dxi1 = D_xi1 * wsp0 = D_xi1 * u
// wsp2 = du_dxi2 = D_xi2 * wsp0 = D_xi2 * u
BwdTrans_SumFacKernel(base0,base1,inarray,wsp0,wsp1,true,true);
LaplacianMatrixOp_MatFree_Kernel(wsp0, outarray, wsp1);
}
else
{
LaplacianMatrixOp_MatFree_Kernel(inarray, outarray, wsp1);
}
}
else
{
StdExpansion::LaplacianMatrixOp_MatFree_GenericImpl(
inarray,outarray,mkey);
}
}
void QuadExp::v_HelmholtzMatrixOp_MatFree(
const Array<OneD, const NekDouble> &inarray,
Array<OneD,NekDouble> &outarray,
const StdRegions::StdMatrixKey &mkey)
{
if (mkey.GetNVarCoeff() == 0
&&!mkey.ConstFactorExists(StdRegions::eFactorSVVCutoffRatio))
{
int nquad0 = m_base[0]->GetNumPoints();
int nquad1 = m_base[1]->GetNumPoints();
int nqtot = nquad0*nquad1;
int nmodes0 = m_base[0]->GetNumModes();
int nmodes1 = m_base[1]->GetNumModes();
int wspsize = max(max(max(nqtot,m_ncoeffs),nquad1*nmodes0),
nquad0*nmodes1);
NekDouble lambda =
mkey.GetConstFactor(StdRegions::eFactorLambda);
const Array<OneD, const NekDouble>& base0 = m_base[0]->GetBdata();
const Array<OneD, const NekDouble>& base1 = m_base[1]->GetBdata();
// Allocate temporary storage
Array<OneD,NekDouble> wsp0(5*wspsize); // size wspsize
Array<OneD,NekDouble> wsp1(wsp0 + wspsize); // size wspsize
Array<OneD,NekDouble> wsp2(wsp0 + 2*wspsize);// size 3*wspsize
if (!(m_base[0]->Collocation() && m_base[1]->Collocation()))
{
// MASS MATRIX OPERATION
// The following is being calculated:
// wsp0 = B * u_hat = u
// wsp1 = W * wsp0
// outarray = B^T * wsp1 = B^T * W * B * u_hat = M * u_hat
BwdTrans_SumFacKernel (base0, base1, inarray,
wsp0, wsp2,true,true);
MultiplyByQuadratureMetric (wsp0, wsp1);
IProductWRTBase_SumFacKernel(base0, base1, wsp1, outarray,
wsp2, true, true);
LaplacianMatrixOp_MatFree_Kernel(wsp0, wsp1, wsp2);
}
else
{
MultiplyByQuadratureMetric(inarray,outarray);
LaplacianMatrixOp_MatFree_Kernel(inarray, wsp1, wsp2);
}
// outarray = lambda * outarray + wsp1
// = (lambda * M + L ) * u_hat
Vmath::Svtvp(m_ncoeffs, lambda, &outarray[0], 1,
&wsp1[0], 1, &outarray[0], 1);
}
else
{
StdExpansion::HelmholtzMatrixOp_MatFree_GenericImpl(
inarray,outarray,mkey);
}
}
void QuadExp::LaplacianMatrixOp_MatFree_Kernel(
void QuadExp::v_LaplacianMatrixOp_MatFree_Kernel(
const Array<OneD, const NekDouble> &inarray,
Array<OneD, NekDouble> &outarray,
Array<OneD, NekDouble> &wsp)
......
......@@ -264,15 +264,7 @@ namespace Nektar
const Array<OneD, const NekDouble> &inarray,
Array<OneD, NekDouble> &outarray,
const StdRegions::StdMatrixKey &mkey);
LOCAL_REGIONS_EXPORT virtual void v_LaplacianMatrixOp_MatFree(
const Array<OneD, const NekDouble> &inarray,
Array<OneD, NekDouble> &outarray,
const StdRegions::StdMatrixKey &mkey);
LOCAL_REGIONS_EXPORT virtual void v_HelmholtzMatrixOp_MatFree(
const Array<OneD, const NekDouble> &inarray,
Array<OneD, NekDouble> &outarray,
const StdRegions::StdMatrixKey &mkey);
LOCAL_REGIONS_EXPORT virtual void LaplacianMatrixOp_MatFree_Kernel(
LOCAL_REGIONS_EXPORT virtual void v_LaplacianMatrixOp_MatFree_Kernel(
const Array<OneD, const NekDouble> &inarray,
Array<OneD, NekDouble> &outarray,
Array<OneD, NekDouble> &wsp);
......
......@@ -1143,142 +1143,6 @@ namespace Nektar
}
/**
* \brief Evaluates the Helmholtz operator using a matrix-free approach.
*
* To construct the Helmholtz operator in a physical tetrahedron
* requires coordinate transforms from both the collapsed coordinate
* system to the standard region and from the standard region to the
* local region. This double application of the chain rule requires the
* calculation of two sets of geometric factors:
* @f[ h_{ij} = \frac{\partial \eta_i}{\partial \xi_j} @f]
* and
* @f[ g_{ij} = \frac{\partial \xi_i}{\partial x_j} @f]
*
* From the definition of the collapsed coordinates, the @f$h_{ij}@f$
* terms are (Sherwin & Karniadakis, p152)
* @f[
* \mathbf{H} = \left[\begin{array}{ccc}
* \frac{4}{(1-\eta_2)(1-\eta_3)} &
* \frac{2(1+\eta_1)}{(1-\eta_2)(1-\eta_3)} &
* \frac{2(1+\eta_1)}{(1-\eta_2)(1-\eta_3)} \\
* 0 &
* \frac{2}{1-eta_3} &
* \frac{1+\eta_2}{1-\eta_3} \\
* 0 &
* 0 &
* 1
* \end{array}\right]
* @f]
* This maps from the collapsed coordinate system to the standard
* tetrahedral region. The mapping to the local region is then given
* by the @f$g_{ij}@f$ computed in the GeomFactors3D class. The
* cumulative factors for mapping the collapsed coordinate system to
* the physical region are therefore given by
* @f$\mathbf{F} = \mathbf{GH^{\top}}@f$, i.e.
* @f[
* f_{ij} = \frac{\partial \eta_i}{\partial x_j}
* = \sum_k g_{ik} h_{kj}
* @f]
*
* Finally, the evaluation of the Helmholtz matrix operator requires
* the summation of these factors as follows. For the case of deformed
* elements, these coefficients are vectors, whereas for regular
* elements they are just scalars.
* @f[
* \begin{array}{l}
* p_0 = \sum_k f_{1k}^2 \\
* p_1 = \sum_k f_{2k}^2 \\
* p_2 = \sum_k f_{3k}^2 \\
* p_3 = \sum_k f_{1k}f_{2k} \\
* p_4 = \sum_k f_{1k}f_{3k} \\
* p_5 = \sum_k f_{2k}f_{3k}
* \end{array}
* @f]
* to give the Helmholtz operator:
* @f{align}
* \mathbf{L^e\hat{u}}
* = \mathbf{B^{\top}D_{\eta_1}^{\top}Wp_0D_{\eta_1}B\hat{u}}
* + \mathbf{B^{\top}D_{\eta_2}^{\top}Wp_1D_{\eta_2}B\hat{u}}
* + \mathbf{B^{\top}D_{\eta_3}^{\top}Wp_2D_{\eta_3}B\hat{u}}\\
* + \mathbf{B^{\top}D_{\eta_1}^{\top}Wp_3D_{\eta_2}B\hat{u}}
* + \mathbf{B^{\top}D_{\eta_1}^{\top}Wp_4D_{\eta_3}B\hat{u}}
* + \mathbf{B^{\top}D_{\eta_2}^{\top}Wp_5D_{\eta_3}B\hat{u}}
* @f}
* Therefore, we construct the operator as follows:
* -# Apply the mass matrix for the @f$\lambda@f$ term
* @f$ \mathbf{B^{\top}WB\hat{u}} @f$.
* and compute the derivatives @f$ \mathbf{D_{\xi_i}B} @f$.
* -# Compute the non-trivial @f$ \mathbf{H} @f$ matrix terms.
* -# Compute the intermediate factors @f$ \mathbf{G} @f$ and
* @f$ f_{ij} @f$ and then compute the combined terms @f$ p_i @f$.
* -# Apply quadrature weights and inner product with respect to the
* derivative bases.
* -# Combine to produce the complete operator.
*/
void TetExp::v_HelmholtzMatrixOp_MatFree(
const Array<OneD, const NekDouble> &inarray,
Array<OneD,NekDouble> &outarray,
const StdRegions::StdMatrixKey &mkey)
{
if(mkey.GetNVarCoeff() == 0)
{
int nquad0 = m_base[0]->GetNumPoints();
int nquad1 = m_base[1]->GetNumPoints();
int nquad2 = m_base[2]->GetNumPoints();
int nqtot = nquad0*nquad1*nquad2;
int nmodes0 = m_base[0]->GetNumModes();
int nmodes1 = m_base[1]->GetNumModes();
int nmodes2 = m_base[2]->GetNumModes();
int wspsize = max(nquad0*nmodes2*(nmodes1+nquad1),
nquad2*nmodes0*nmodes1*(nmodes1+1)/2+
nquad2*nquad1*nmodes0);
NekDouble lambda = mkey.GetConstFactor(StdRegions::eFactorLambda);
const Array<OneD, const NekDouble>& base0 = m_base[0]->GetBdata ();
const Array<OneD, const NekDouble>& base1 = m_base[1]->GetBdata ();
const Array<OneD, const NekDouble>& base2 = m_base[2]->GetBdata ();
Array<OneD,NekDouble> wsp (wspsize);
Array<OneD,NekDouble> wsp0(nqtot);
Array<OneD,NekDouble> wsp1(nqtot);
if(!(m_base[0]->Collocation() && m_base[1]->Collocation() &&
m_base[2]->Collocation()))
{
// MASS MATRIX OPERATION
// The following is being calculated:
// wsp0 = B * u_hat = u
// wsp1 = W * wsp0
// outarray = B^T * wsp1 = B^T * W * B * u_hat = M * u_hat
BwdTrans_SumFacKernel (base0,base1,base2,inarray,
wsp0,wsp,true,true,true);
MultiplyByQuadratureMetric (wsp0,wsp1);
IProductWRTBase_SumFacKernel (base0,base1,base2,wsp1,
outarray,wsp,true,true,true);
LaplacianMatrixOp_MatFree_Kernel(wsp0,wsp1,wsp);
}
else
{
// specialised implementation for the classical spectral
// element method
MultiplyByQuadratureMetric (inarray,outarray);
LaplacianMatrixOp_MatFree_Kernel(inarray,wsp1,wsp);
}
// outarray = lambda * outarray + wsp1
// = (lambda * M + L ) * u_hat
Vmath::Svtvp(m_ncoeffs,lambda,&outarray[0],1,&wsp1[0],1,
&outarray[0],1);
}
else
{
StdExpansion::HelmholtzMatrixOp_MatFree_GenericImpl(inarray,outarray,mkey);
}
}
void TetExp::v_LaplacianMatrixOp(
const Array<OneD, const NekDouble> &inarray,
Array<OneD,NekDouble> &outarray,
......@@ -1298,38 +1162,6 @@ namespace Nektar
mkey);
}
void TetExp::v_LaplacianMatrixOp_MatFree(
const Array<OneD, const NekDouble> &inarray,
Array<OneD,NekDouble> &outarray,
const StdRegions::StdMatrixKey &mkey)
{
if(mkey.GetNVarCoeff() == 0)
{
// This implementation is only valid when there are no
// coefficients associated to the Laplacian operator
int nqtot = GetTotPoints();
const Array<OneD, const NekDouble>& base0 = m_base[0]->GetBdata();
const Array<OneD, const NekDouble>& base1 = m_base[1]->GetBdata();
const Array<OneD, const NekDouble>& base2 = m_base[2]->GetBdata();
// Allocate temporary storage
Array<OneD,NekDouble> wsp0(7*nqtot);
Array<OneD,NekDouble> wsp1(wsp0+nqtot);
// LAPLACIAN MATRIX OPERATION