Commit 99719e31 authored by Chris Cantwell's avatar Chris Cantwell
Browse files

Removed references to Laplacian metrics in GeomFactors

Tidied up QuadExp, TriExp, HexExp LaplacianMatrixOp/HelmholtzMatrixOp ops.
parent 7c96ad93
This diff is collapsed.
......@@ -1309,13 +1309,13 @@ namespace Nektar
{
// 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
{
// 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();
......@@ -1333,7 +1333,7 @@ namespace Nektar
// 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
{
......@@ -1377,12 +1377,12 @@ namespace Nektar
Array<OneD, NekDouble> &outarray,
const StdRegions::StdMatrixKey &mkey)
{
if(m_metricinfo->IsUsingLaplMetrics())
{
ASSERTL0(false,"Finish implementing PrismExp Helmholtz for Lapl Metrics");
}
else
{
// 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();
......@@ -1416,7 +1416,7 @@ namespace Nektar
// = (lambda * M + L ) * u_hat
Vmath::Svtvp(m_ncoeffs,lambda,&outarray[0],1,&wsp1[0],1,
&outarray[0],1);
}
// }
}
//---------------------------------------
......
......@@ -2440,11 +2440,6 @@ namespace Nektar
{
// This implementation is only valid when there are no
// coefficients associated to the Laplacian operator
if (m_metrics.count(MetricLaplacian00) == 0)
{
ComputeLaplacianMetric();
}
int nquad0 = m_base[0]->GetNumPoints();
int nquad1 = m_base[1]->GetNumPoints();
int nqtot = nquad0*nquad1;
......@@ -2454,16 +2449,10 @@ namespace Nektar
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>& dbase0 = m_base[0]->GetDbdata();
const Array<OneD, const NekDouble>& dbase1 = m_base[1]->GetDbdata();
const Array<OneD, const NekDouble>& metric00 = m_metrics[MetricLaplacian00];
const Array<OneD, const NekDouble>& metric01 = m_metrics[MetricLaplacian01];
const Array<OneD, const NekDouble>& metric11 = m_metrics[MetricLaplacian11];
// Allocate temporary storage
Array<OneD,NekDouble> wsp0(3*wspsize);
Array<OneD,NekDouble> wsp1(wsp0+wspsize);
Array<OneD,NekDouble> wsp2(wsp0+2*wspsize);
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()))
{
......@@ -2472,28 +2461,12 @@ namespace Nektar
// 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);
StdExpansion2D::PhysTensorDeriv(wsp0,wsp1,wsp2);
LaplacianMatrixOp_MatFree_Kernel(wsp0, outarray, wsp1);
}
else
{
StdExpansion2D::PhysTensorDeriv(inarray,wsp1,wsp2);
LaplacianMatrixOp_MatFree_Kernel(inarray, outarray, wsp1);
}
// wsp0 = k = g0 * wsp1 + g1 * wsp2 = g0 * du_dxi1 + g1 * du_dxi2
// wsp2 = l = g1 * wsp1 + g2 * wsp2 = g0 * du_dxi1 + g1 * du_dxi2
// where g0, g1 and g2 are the metric terms set up in the GeomFactors class
// especially for this purpose
Vmath::Vvtvvtp(nqtot,&metric00[0],1,&wsp1[0],1,&metric01[0],1,&wsp2[0],1,&wsp0[0],1);
Vmath::Vvtvvtp(nqtot,&metric01[0],1,&wsp1[0],1,&metric11[0],1,&wsp2[0],1,&wsp2[0],1);
// outarray = m = (D_xi1 * B)^T * k
// wsp1 = n = (D_xi2 * B)^T * l
IProductWRTBase_SumFacKernel(dbase0,base1,wsp0,outarray,wsp1,false,true);
IProductWRTBase_SumFacKernel(base0,dbase1,wsp2,wsp1, wsp0,true,false);
// outarray = outarray + wsp1
// = L * u_hat
Vmath::Vadd(m_ncoeffs,wsp1.get(),1,outarray.get(),1,outarray.get(),1);
}
else
{
......@@ -2509,7 +2482,8 @@ namespace Nektar
Array<OneD,NekDouble> &outarray,
const StdRegions::StdMatrixKey &mkey)
{
if (m_metricinfo->IsUsingLaplMetrics())
if (mkey.GetNVarCoeff() == 0
&&!mkey.ConstFactorExists(StdRegions::eFactorSVVCutoffRatio))
{
int nquad0 = m_base[0]->GetNumPoints();
int nquad1 = m_base[1]->GetNumPoints();
......@@ -2521,22 +2495,13 @@ namespace Nektar
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>& dbase0 =
m_base[0]->GetDbdata();
const Array<OneD, const NekDouble>& dbase1 =
m_base[1]->GetDbdata();
const Array<TwoD, const NekDouble>& metric =
m_metricinfo->GetLaplacianMetrics();
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);
Array<OneD,NekDouble> wsp1(wsp0 + wspsize);
Array<OneD,NekDouble> wsp2(wsp0 + 2*wspsize);
Array<OneD,NekDouble> wsp3(wsp0 + 3*wspsize);
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()))
{
......@@ -2546,200 +2511,85 @@ namespace Nektar
// wsp1 = W * wsp0
// outarray = B^T * wsp1 = B^T * W * B * u_hat = M * u_hat
BwdTrans_SumFacKernel (base0, base1, inarray,
wsp0, wsp1,true,true);
MultiplyByQuadratureMetric (wsp0, wsp2);
IProductWRTBase_SumFacKernel(base0, base1, wsp2, outarray,
wsp1, true, true);
wsp0, wsp2,true,true);
MultiplyByQuadratureMetric (wsp0, wsp1);
IProductWRTBase_SumFacKernel(base0, base1, wsp1, outarray,
wsp2, true, true);
// LAPLACIAN MATRIX OPERATION
// wsp1 = du_dxi1 = D_xi1 * wsp0 = D_xi1 * u
// wsp2 = du_dxi2 = D_xi2 * wsp0 = D_xi2 * u
StdExpansion2D::PhysTensorDeriv(wsp0,wsp1,wsp2);
LaplacianMatrixOp_MatFree_Kernel(wsp0, wsp1, wsp2);
}
else
{
// specialised implementation for the
// classical spectral element method
StdExpansion2D::PhysTensorDeriv(inarray,wsp1,wsp2);
MultiplyByQuadratureMetric(inarray,outarray);
LaplacianMatrixOp_MatFree_Kernel(inarray, wsp1, wsp2);
}
// wsp0 = k = g0 * wsp1 + g1 * wsp2 = g0 * du_dxi1 + g1 * du_dxi2
// wsp2 = l = g1 * wsp1 + g2 * wsp2 = g1 * du_dxi1 + g2 * du_dxi2
// where g0, g1 and g2 are the metric terms set up in the
// GeomFactors class especially for this purpose
if(!m_metricinfo->LaplacianMetricIsZero(1))
{
Vmath::Vvtvvtp(nqtot, &metric[0][0], 1, &wsp1[0], 1,
&metric[1][0], 1, &wsp2[0], 1, &wsp0[0], 1);
Vmath::Vvtvvtp(nqtot, &metric[1][0], 1, &wsp1[0], 1,
&metric[2][0], 1, &wsp2[0], 1, &wsp2[0], 1);
}
else
{
// special implementation in case g1 = 0
// (which should hold for undistorted quads)
// wsp0 = k = g0 * wsp1 = g0 * du_dxi1
// wsp2 = l = g2 * wsp2 = g2 * du_dxi2
Vmath::Vmul(nqtot, &metric[0][0], 1,
&wsp1[0], 1, &wsp0[0], 1);
Vmath::Vmul(nqtot, &metric[2][0], 1,
&wsp2[0], 1, &wsp2[0], 1);
}
// wsp1 = m = (D_xi1 * B)^T * k
// wsp0 = n = (D_xi2 * B)^T * l
IProductWRTBase_SumFacKernel(dbase0, base1,
wsp0, wsp1, wsp3, false, true);
IProductWRTBase_SumFacKernel(base0, dbase1,
wsp2, wsp0, wsp3, true, false);
// outarray = lambda * outarray + (wsp0 + wsp1)
// outarray = lambda * outarray + wsp1
// = (lambda * M + L ) * u_hat
Vmath::Vstvpp(m_ncoeffs, lambda, &outarray[0], 1,
&wsp1[0], 1, &wsp0[0], 1, &outarray[0], 1);
Vmath::Svtvp(m_ncoeffs, lambda, &outarray[0], 1,
&wsp1[0], 1, &outarray[0], 1);
}
else
{
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();
const Array<OneD, const NekDouble>& dbase0 =
m_base[0]->GetDbdata();
const Array<OneD, const NekDouble>& dbase1 =
m_base[1]->GetDbdata();
// Allocate temporary storage
Array<OneD,NekDouble> wsp0(6*wspsize);
Array<OneD,NekDouble> wsp1(wsp0 + wspsize);
Array<OneD,NekDouble> wsp2(wsp0 + 2*wspsize);
Array<OneD,NekDouble> wsp3(wsp0 + 3*wspsize);
Array<OneD,NekDouble> wsp4(wsp0 + 4*wspsize);
Array<OneD,NekDouble> wsp5(wsp0 + 5*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, wsp1, true, true);
MultiplyByQuadratureMetric (wsp0, wsp2);
IProductWRTBase_SumFacKernel(base0, base1, wsp2, outarray,
wsp1, true, true);
// LAPLACIAN MATRIX OPERATION
// wsp1 = du_dxi1 = D_xi1 * wsp0 = D_xi1 * u
// wsp2 = du_dxi2 = D_xi2 * wsp0 = D_xi2 * u
StdExpansion2D::PhysTensorDeriv(wsp0, wsp1, wsp2);
}
else
{
// specialised implementation for the
// classical spectral element method
StdExpansion2D::PhysTensorDeriv(inarray,wsp1,wsp2);
MultiplyByQuadratureMetric(inarray,outarray);
}
// wsp0 = k = g0 * wsp1 + g1 * wsp2 = g0 * du_dxi1 + g1 * du_dxi2
// wsp2 = l = g1 * wsp1 + g2 * wsp2 = g1 * du_dxi1 + g2 * du_dxi2
// where g0, g1 and g2 are the metric terms set up in the
// GeomFactors class especially for this purpose
const Array<TwoD, const NekDouble>& df =
m_metricinfo->GetDerivFactors();
int dim = m_geom->GetCoordim();
if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
{
// wsp3 = g0
Vmath::Vvtvvtp(nqtot,
&df[0][0], 1, &df[0][0], 1,
&df[2][0], 1, &df[2][0], 1, &wsp3[0], 1);
// wsp4 = g1;
Vmath::Vvtvvtp(nqtot,
&df[0][0], 1, &df[1][0], 1,
&df[2][0], 1, &df[3][0], 1, &wsp4[0], 1);
// wsp5 = g2;
Vmath::Vvtvvtp(nqtot,
&df[1][0], 1, &df[1][0], 1,
&df[3][0], 1, &df[3][0], 1, &wsp5[0], 1);
if (dim == 3)
{
Vmath::Vvtvp(nqtot,
&df[4][0], 1, &df[4][0], 1,
&wsp3[0], 1, &wsp3[0], 1);
Vmath::Vvtvp(nqtot,
&df[4][0], 1, &df[5][0], 1,
&wsp4[0], 1, &wsp4[0], 1);
Vmath::Vvtvp(nqtot,
&df[5][0], 1, &df[5][0], 1,
&wsp5[0], 1, &wsp5[0], 1);
}
Vmath::Vvtvvtp(nqtot,
&wsp3[0], 1, &wsp1[0], 1, &wsp4[0], 1,
&wsp2[0], 1, &wsp0[0], 1);
Vmath::Vvtvvtp(nqtot,
&wsp4[0], 1, &wsp1[0], 1, &wsp5[0], 1,
&wsp2[0], 1, &wsp2[0], 1);
}
else
{
NekDouble g0 = df[0][0]*df[0][0] + df[2][0]*df[2][0];
NekDouble g1 = df[0][0]*df[1][0] + df[2][0]*df[3][0];
NekDouble g2 = df[1][0]*df[1][0] + df[3][0]*df[3][0];
if (dim == 3)
{
g0 += df[4][0]*df[4][0];
g1 += df[4][0]*df[5][0];
g2 += df[5][0]*df[5][0];
}
if(fabs(g1) < NekConstants::kGeomFactorsTol)
{
Vmath::Smul(nqtot,g0,&wsp1[0],1,&wsp0[0],1);
Vmath::Smul(nqtot,g2,&wsp2[0],1,&wsp2[0],1);
}
else
{
Vmath::Svtsvtp(nqtot, g0, &wsp1[0], 1,
g1, &wsp2[0], 1, &wsp0[0], 1);
Vmath::Svtsvtp(nqtot, g1, &wsp1[0], 1,
g2, &wsp2[0], 1, &wsp2[0],1);
}
}
MultiplyByQuadratureMetric(wsp0, wsp0);
MultiplyByQuadratureMetric(wsp2, wsp2);
StdExpansion::HelmholtzMatrixOp_MatFree_GenericImpl(
inarray,outarray,mkey);
}
}
// wsp1 = m = (D_xi1 * B)^T * k
// wsp0 = n = (D_xi2 * B)^T * l
IProductWRTBase_SumFacKernel(dbase0, base1,
wsp0, wsp1, wsp3, false, true);
IProductWRTBase_SumFacKernel(base0, dbase1,
wsp2, wsp0, wsp3, true, false);
// outarray = lambda * outarray + (wsp0 + wsp1)
// = (lambda * M + L ) * u_hat
Vmath::Vstvpp(m_ncoeffs, lambda, &outarray[0], 1, &wsp1[0], 1,
&wsp0[0], 1, &outarray[0], 1);
void QuadExp::LaplacianMatrixOp_MatFree_Kernel(
const Array<OneD, const NekDouble> &inarray,
Array<OneD, NekDouble> &outarray,
Array<OneD, NekDouble> &wsp)
{
if (m_metrics.count(MetricLaplacian00) == 0)
{
ComputeLaplacianMetric();
}
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);
ASSERTL1(wsp.num_elements() >= 3*wspsize,
"Workspace is of insufficient size.");
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>& dbase0 = m_base[0]->GetDbdata();
const Array<OneD, const NekDouble>& dbase1 = m_base[1]->GetDbdata();
const Array<OneD, const NekDouble>& metric00 = m_metrics[MetricLaplacian00];
const Array<OneD, const NekDouble>& metric01 = m_metrics[MetricLaplacian01];
const Array<OneD, const NekDouble>& metric11 = m_metrics[MetricLaplacian11];
// Allocate temporary storage
Array<OneD,NekDouble> wsp0(wsp);
Array<OneD,NekDouble> wsp1(wsp+wspsize);
Array<OneD,NekDouble> wsp2(wsp+2*wspsize);
StdExpansion2D::PhysTensorDeriv(inarray,wsp1,wsp2);
// wsp0 = k = g0 * wsp1 + g1 * wsp2 = g0 * du_dxi1 + g1 * du_dxi2
// wsp2 = l = g1 * wsp1 + g2 * wsp2 = g0 * du_dxi1 + g1 * du_dxi2
// where g0, g1 and g2 are the metric terms set up in the GeomFactors class
// especially for this purpose
Vmath::Vvtvvtp(nqtot,&metric00[0],1,&wsp1[0],1,&metric01[0],1,&wsp2[0],1,&wsp0[0],1);
Vmath::Vvtvvtp(nqtot,&metric01[0],1,&wsp1[0],1,&metric11[0],1,&wsp2[0],1,&wsp2[0],1);
// outarray = m = (D_xi1 * B)^T * k
// wsp1 = n = (D_xi2 * B)^T * l
IProductWRTBase_SumFacKernel(dbase0,base1,wsp0,outarray,wsp1,false,true);
IProductWRTBase_SumFacKernel(base0,dbase1,wsp2,wsp1, wsp0,true,false);
// outarray = outarray + wsp1
// = L * u_hat
Vmath::Vadd(m_ncoeffs,wsp1.get(),1,outarray.get(),1,outarray.get(),1);
}
void QuadExp::v_ComputeLaplacianMetric()
{
if (m_metrics.count(MetricQuadrature) == 0)
......
......@@ -272,7 +272,11 @@ namespace Nektar
const Array<OneD, const NekDouble> &inarray,
Array<OneD, NekDouble> &outarray,
const StdRegions::StdMatrixKey &mkey);
LOCAL_REGIONS_EXPORT virtual void LaplacianMatrixOp_MatFree_Kernel(
const Array<OneD, const NekDouble> &inarray,
Array<OneD, NekDouble> &outarray,
Array<OneD, NekDouble> &wsp);
LOCAL_REGIONS_EXPORT virtual void v_ComputeLaplacianMetric();
private:
......
......@@ -1221,81 +1221,7 @@ namespace Nektar
Array<OneD,NekDouble> &outarray,
const StdRegions::StdMatrixKey &mkey)
{
if(m_metricinfo->IsUsingLaplMetrics())
{
ASSERTL0(false,"Finish implementing TetExp Helmholtz for Lapl Metrics");
/*
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.GetConstant(0);
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>& dbase0 = m_base[0]->GetDbdata();
const Array<OneD, const NekDouble>& dbase1 = m_base[1]->GetDbdata();
const Array<TwoD, const NekDouble>& metric = m_metricinfo->GetLaplacianMetrics();
// Allocate temporary storage
Array<OneD,NekDouble> wsp0(4*wspsize);
Array<OneD,NekDouble> wsp1(wsp0+wspsize);
Array<OneD,NekDouble> wsp2(wsp0+2*wspsize);
Array<OneD,NekDouble> wsp3(wsp0+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, wsp1,true,true);
MultiplyByQuadratureMetric (wsp0,wsp2);
IProductWRTBase_SumFacKernel(base0,base1,wsp2, outarray,wsp1,true,true);
// LAPLACIAN MATRIX OPERATION
// wsp1 = du_dxi1 = D_xi1 * wsp0 = D_xi1 * u
// wsp2 = du_dxi2 = D_xi2 * wsp0 = D_xi2 * u
StdExpansion2D::PhysTensorDeriv(wsp0,wsp1,wsp2);
}
else
{
// specialised implementation for the classical spectral element method
StdExpansion2D::PhysTensorDeriv(inarray,wsp1,wsp2);
MultiplyByQuadratureMetric(inarray,outarray);
}
// wsp0 = k = g0 * wsp1 + g1 * wsp2 = g0 * du_dxi1 + g1 * du_dxi2
// wsp2 = l = g1 * wsp1 + g2 * wsp2 = g1 * du_dxi1 + g2 * du_dxi2
// where g0, g1 and g2 are the metric terms set up in the GeomFactors class
// especially for this purpose
if(!m_metricinfo->LaplacianMetricIsZero(1))
{
Vmath::Vvtvvtp(nqtot,&metric[0][0],1,&wsp1[0],1,&metric[1][0],1,&wsp2[0],1,&wsp0[0],1);
Vmath::Vvtvvtp(nqtot,&metric[1][0],1,&wsp1[0],1,&metric[2][0],1,&wsp2[0],1,&wsp2[0],1);
}
else
{
// special implementation in case g1 = 0 (which should hold for undistorted quads)
// wsp0 = k = g0 * wsp1 = g0 * du_dxi1
// wsp2 = l = g2 * wsp2 = g2 * du_dxi2
Vmath::Vmul(nqtot,&metric[0][0],1,&wsp1[0],1,&wsp0[0],1);
Vmath::Vmul(nqtot,&metric[2][0],1,&wsp2[0],1,&wsp2[0],1);
}
// wsp1 = m = (D_xi1 * B)^T * k
// wsp0 = n = (D_xi2 * B)^T * l
IProductWRTBase_SumFacKernel(dbase0,base1,wsp0,wsp1,wsp3,false,true);
IProductWRTBase_SumFacKernel(base0,dbase1,wsp2,wsp0,wsp3,true,false);
// outarray = lambda * outarray + (wsp0 + wsp1)
// = (lambda * M + L ) * u_hat
Vmath::Vstvpp(m_ncoeffs,lambda,&outarray[0],1,&wsp1[0],1,&wsp0[0],1,&outarray[0],1);
*/ }
else
if(mkey.GetNVarCoeff() == 0)
{
int nquad0 = m_base[0]->GetNumPoints();
int nquad1 = m_base[1]->GetNumPoints();
......@@ -1345,6 +1271,11 @@ namespace Nektar
Vmath::Svtvp(m_ncoeffs,lambda,&outarray[0],1,&wsp1[0],1,
&outarray[0],1);
}
else
{
StdExpansion::HelmholtzMatrixOp_MatFree_GenericImpl(inarray,outarray,mkey);
}
}
......@@ -1372,50 +1303,26 @@ namespace Nektar
Array<OneD,NekDouble> &outarray,
const StdRegions::StdMatrixKey &mkey)
{
if( mkey.GetNVarCoeff() == 0 &&
!mkey.ConstFactorExists(StdRegions::eFactorSVVCutoffRatio))
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 TetExp 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 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);
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> wsp1(nqtot);
// Backwards transform to obtain u = B * u_hat.
if(!(m_base[0]->Collocation() && m_base[1]->Collocation() &&
m_base[2]->Collocation()))
{
BwdTrans_SumFacKernel (base0,base1,base2,inarray,
wsp1,wsp,true,true,true);
LaplacianMatrixOp_MatFree_Kernel(wsp1, outarray, wsp);
}
else
{
LaplacianMatrixOp_MatFree_Kernel(wsp1, outarray, wsp);
}
}
int nqtot = GetTotPoints();
const Array<OneD, const NekDouble>& base0 = m_base[0]->GetBdata();
const Array<OneD, const NekDouble>& base1 = m_base[1]->GetBdata();