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

Add Laplacian metrics to pyramid

parent 9cf23142
......@@ -451,7 +451,8 @@ namespace Nektar
break;
case StdRegions::eLaplacian:
{
if(m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
if (true)
//if(m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
{
NekDouble one = 1.0;
DNekMatSharedPtr mat = GenMatrix(mkey);
......@@ -641,5 +642,227 @@ namespace Nektar
return returnval;
}
void PyrExp::v_ComputeLaplacianMetric()
{
if (m_metrics.count(MetricQuadrature) == 0)
{
ComputeQuadratureMetric();
}
int i, j;
const unsigned int nqtot = GetTotPoints();
const unsigned int dim = 3;
const MetricType m[3][3] = {
{ MetricLaplacian00, MetricLaplacian01, MetricLaplacian02 },
{ MetricLaplacian01, MetricLaplacian11, MetricLaplacian12 },
{ MetricLaplacian02, MetricLaplacian12, MetricLaplacian22 }
};
for (unsigned int i = 0; i < dim; ++i)
{
for (unsigned int j = i; j < dim; ++j)
{
m_metrics[m[i][j]] = Array<OneD, NekDouble>(nqtot);
}
}
// Define shorthand synonyms for m_metrics storage
Array<OneD,NekDouble> g0 (m_metrics[m[0][0]]);
Array<OneD,NekDouble> g1 (m_metrics[m[1][1]]);
Array<OneD,NekDouble> g2 (m_metrics[m[2][2]]);
Array<OneD,NekDouble> g3 (m_metrics[m[0][1]]);
Array<OneD,NekDouble> g4 (m_metrics[m[0][2]]);
Array<OneD,NekDouble> g5 (m_metrics[m[1][2]]);
// Allocate temporary storage
Array<OneD,NekDouble> alloc(9*nqtot,0.0);
Array<OneD,NekDouble> h0 (nqtot, alloc);
Array<OneD,NekDouble> h1 (nqtot, alloc+ 1*nqtot);
Array<OneD,NekDouble> h2 (nqtot, alloc+ 2*nqtot);
Array<OneD,NekDouble> wsp1 (nqtot, alloc+ 3*nqtot);
Array<OneD,NekDouble> wsp2 (nqtot, alloc+ 4*nqtot);
Array<OneD,NekDouble> wsp3 (nqtot, alloc+ 5*nqtot);
Array<OneD,NekDouble> wsp4 (nqtot, alloc+ 6*nqtot);
Array<OneD,NekDouble> wsp5 (nqtot, alloc+ 7*nqtot);
Array<OneD,NekDouble> wsp6 (nqtot, alloc+ 8*nqtot);
const Array<TwoD, const NekDouble>& df =
m_metricinfo->GetDerivFactors(GetPointsKeys());
const Array<OneD, const NekDouble>& z0 = m_base[0]->GetZ();
const Array<OneD, const NekDouble>& z1 = m_base[1]->GetZ();
const Array<OneD, const NekDouble>& z2 = m_base[2]->GetZ();
const unsigned int nquad0 = m_base[0]->GetNumPoints();
const unsigned int nquad1 = m_base[1]->GetNumPoints();
const unsigned int nquad2 = m_base[2]->GetNumPoints();
// Populate collapsed coordinate arrays h0, h1 and h2.
for(j = 0; j < nquad2; ++j)
{
for(i = 0; i < nquad1; ++i)
{
Vmath::Fill(nquad0, 2.0/(1.0-z2[j]), &h0[0]+i*nquad0 + j*nquad0*nquad1,1);
Vmath::Fill(nquad0, 1.0/(1.0-z2[j]), &h1[0]+i*nquad0 + j*nquad0*nquad1,1);
Vmath::Fill(nquad0, (1.0+z1[i])/(1.0-z2[j]), &h2[0]+i*nquad0 + j*nquad0*nquad1,1);
}
}
for(i = 0; i < nquad0; i++)
{
Blas::Dscal(nquad1*nquad2, 1+z0[i], &h1[0]+i, nquad0);
}
// Step 3. Construct combined metric terms for physical space to
// collapsed coordinate system.
// Order of construction optimised to minimise temporary storage
if(m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
{
// f_{1k}
Vmath::Vvtvvtp(nqtot, &df[0][0], 1, &h0[0], 1, &df[2][0], 1, &h1[0], 1, &wsp1[0], 1);
Vmath::Vvtvvtp(nqtot, &df[3][0], 1, &h0[0], 1, &df[5][0], 1, &h1[0], 1, &wsp2[0], 1);
Vmath::Vvtvvtp(nqtot, &df[6][0], 1, &h0[0], 1, &df[8][0], 1, &h1[0], 1, &wsp3[0], 1);
// g0
Vmath::Vvtvvtp(nqtot, &wsp1[0], 1, &wsp1[0], 1, &wsp2[0], 1, &wsp2[0], 1, &g0[0], 1);
Vmath::Vvtvp (nqtot, &wsp3[0], 1, &wsp3[0], 1, &g0[0], 1, &g0[0], 1);
// g4
Vmath::Vvtvvtp(nqtot, &df[2][0], 1, &wsp1[0], 1, &df[5][0], 1, &wsp2[0], 1, &g4[0], 1);
Vmath::Vvtvp (nqtot, &df[8][0], 1, &wsp3[0], 1, &g4[0], 1, &g4[0], 1);
// f_{2k}
Vmath::Vvtvvtp(nqtot, &df[1][0], 1, &h0[0], 1, &df[2][0], 1, &h2[0], 1, &wsp4[0], 1);
Vmath::Vvtvvtp(nqtot, &df[4][0], 1, &h0[0], 1, &df[5][0], 1, &h2[0], 1, &wsp5[0], 1);
Vmath::Vvtvvtp(nqtot, &df[7][0], 1, &h0[0], 1, &df[8][0], 1, &h2[0], 1, &wsp6[0], 1);
// g1
Vmath::Vvtvvtp(nqtot, &wsp4[0], 1, &wsp4[0], 1, &wsp5[0], 1, &wsp5[0], 1, &g1[0], 1);
Vmath::Vvtvp (nqtot, &wsp6[0], 1, &wsp6[0], 1, &g1[0], 1, &g1[0], 1);
// g3
Vmath::Vvtvvtp(nqtot, &wsp1[0], 1, &wsp4[0], 1, &wsp2[0], 1, &wsp5[0], 1, &g3[0], 1);
Vmath::Vvtvp (nqtot, &wsp3[0], 1, &wsp6[0], 1, &g3[0], 1, &g3[0], 1);
// g5
Vmath::Vvtvvtp(nqtot, &df[2][0], 1, &wsp4[0], 1, &df[5][0], 1, &wsp5[0], 1, &g5[0], 1);
Vmath::Vvtvp (nqtot, &df[8][0], 1, &wsp6[0], 1, &g5[0], 1, &g5[0], 1);
// g2
Vmath::Vvtvvtp(nqtot, &df[2][0], 1, &df[2][0], 1, &df[5][0], 1, &df[5][0], 1, &g2[0], 1);
Vmath::Vvtvp (nqtot, &df[8][0], 1, &df[8][0], 1, &g2[0], 1, &g2[0], 1);
}
else
{
// f_{1k}
Vmath::Svtsvtp(nqtot, df[0][0], &h0[0], 1, df[2][0], &h1[0], 1, &wsp1[0], 1);
Vmath::Svtsvtp(nqtot, df[3][0], &h0[0], 1, df[5][0], &h1[0], 1, &wsp2[0], 1);
Vmath::Svtsvtp(nqtot, df[6][0], &h0[0], 1, df[8][0], &h1[0], 1, &wsp3[0], 1);
// g0
Vmath::Vvtvvtp(nqtot, &wsp1[0], 1, &wsp1[0], 1, &wsp2[0], 1, &wsp2[0], 1, &g0[0], 1);
Vmath::Vvtvp (nqtot, &wsp3[0], 1, &wsp3[0], 1, &g0[0], 1, &g0[0], 1);
// g4
Vmath::Svtsvtp(nqtot, df[2][0], &wsp1[0], 1, df[5][0], &wsp2[0], 1, &g4[0], 1);
Vmath::Svtvp (nqtot, df[8][0], &wsp3[0], 1, &g4[0], 1, &g4[0], 1);
// f_{2k}
Vmath::Svtsvtp(nqtot, df[1][0], &h0[0], 1, df[2][0], &h2[0], 1, &wsp4[0], 1);
Vmath::Svtsvtp(nqtot, df[4][0], &h0[0], 1, df[5][0], &h2[0], 1, &wsp5[0], 1);
Vmath::Svtsvtp(nqtot, df[7][0], &h0[0], 1, df[8][0], &h2[0], 1, &wsp6[0], 1);
// g1
Vmath::Vvtvvtp(nqtot, &wsp4[0], 1, &wsp4[0], 1, &wsp5[0], 1, &wsp5[0], 1, &g1[0], 1);
Vmath::Vvtvp (nqtot, &wsp6[0], 1, &wsp6[0], 1, &g1[0], 1, &g1[0], 1);
// g3
Vmath::Vvtvvtp(nqtot, &wsp1[0], 1, &wsp4[0], 1, &wsp2[0], 1, &wsp5[0], 1, &g3[0], 1);
Vmath::Vvtvp (nqtot, &wsp3[0], 1, &wsp6[0], 1, &g3[0], 1, &g3[0], 1);
// g5
Vmath::Svtsvtp(nqtot, df[2][0], &wsp4[0], 1, df[5][0], &wsp5[0], 1, &g5[0], 1);
Vmath::Svtvp (nqtot, df[8][0], &wsp6[0], 1, &g5[0], 1, &g5[0], 1);
// g2
Vmath::Fill(nqtot, df[2][0]*df[2][0] + df[5][0]*df[5][0] + df[8][0]*df[8][0], &g2[0], 1);
}
for (unsigned int i = 0; i < dim; ++i)
{
for (unsigned int j = i; j < dim; ++j)
{
MultiplyByQuadratureMetric(m_metrics[m[i][j]],
m_metrics[m[i][j]]);
}
}
}
void PyrExp::v_LaplacianMatrixOp_MatFree_Kernel(
const Array<OneD, const NekDouble> &inarray,
Array<OneD, NekDouble> &outarray,
Array<OneD, NekDouble> &wsp)
{
// 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 nquad2 = m_base[2]->GetNumPoints();
int nqtot = nquad0*nquad1*nquad2;
ASSERTL1(wsp.num_elements() >= 6*nqtot,
"Insufficient workspace size.");
ASSERTL1(m_ncoeffs <= nqtot,
"Workspace not set up for ncoeffs > nqtot");
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();
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>& dbase2 = m_base[2]->GetDbdata();
const Array<OneD, const NekDouble>& metric00 = m_metrics[MetricLaplacian00];
const Array<OneD, const NekDouble>& metric01 = m_metrics[MetricLaplacian01];
const Array<OneD, const NekDouble>& metric02 = m_metrics[MetricLaplacian02];
const Array<OneD, const NekDouble>& metric11 = m_metrics[MetricLaplacian11];
const Array<OneD, const NekDouble>& metric12 = m_metrics[MetricLaplacian12];
const Array<OneD, const NekDouble>& metric22 = m_metrics[MetricLaplacian22];
// Allocate temporary storage
Array<OneD,NekDouble> wsp0 (2*nqtot, wsp);
Array<OneD,NekDouble> wsp1 ( nqtot, wsp+1*nqtot);
Array<OneD,NekDouble> wsp2 ( nqtot, wsp+2*nqtot);
Array<OneD,NekDouble> wsp3 ( nqtot, wsp+3*nqtot);
Array<OneD,NekDouble> wsp4 ( nqtot, wsp+4*nqtot);
Array<OneD,NekDouble> wsp5 ( nqtot, wsp+5*nqtot);
// LAPLACIAN MATRIX OPERATION
// wsp1 = du_dxi1 = D_xi1 * inarray = D_xi1 * u
// wsp2 = du_dxi2 = D_xi2 * inarray = D_xi2 * u
// wsp2 = du_dxi3 = D_xi3 * inarray = D_xi3 * u
StdExpansion3D::PhysTensorDeriv(inarray,wsp0,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,&wsp0[0],1,&metric01[0],1,&wsp1[0],1,&wsp3[0],1);
Vmath::Vvtvp (nqtot,&metric02[0],1,&wsp2[0],1,&wsp3[0],1,&wsp3[0],1);
Vmath::Vvtvvtp(nqtot,&metric01[0],1,&wsp0[0],1,&metric11[0],1,&wsp1[0],1,&wsp4[0],1);
Vmath::Vvtvp (nqtot,&metric12[0],1,&wsp2[0],1,&wsp4[0],1,&wsp4[0],1);
Vmath::Vvtvvtp(nqtot,&metric02[0],1,&wsp0[0],1,&metric12[0],1,&wsp1[0],1,&wsp5[0],1);
Vmath::Vvtvp (nqtot,&metric22[0],1,&wsp2[0],1,&wsp5[0],1,&wsp5[0],1);
// outarray = m = (D_xi1 * B)^T * k
// wsp1 = n = (D_xi2 * B)^T * l
IProductWRTBase_SumFacKernel(dbase0,base1,base2,wsp3,outarray,wsp0,false,true,true);
IProductWRTBase_SumFacKernel(base0,dbase1,base2,wsp4,wsp2, wsp0,true,false,true);
Vmath::Vadd(m_ncoeffs,wsp2.get(),1,outarray.get(),1,outarray.get(),1);
IProductWRTBase_SumFacKernel(base0,base1,dbase2,wsp5,wsp2, wsp0,true,true,false);
Vmath::Vadd(m_ncoeffs,wsp2.get(),1,outarray.get(),1,outarray.get(),1);
}
}//end of namespace
}//end of namespace
......@@ -138,10 +138,16 @@ namespace Nektar
const MatrixKey &mkey);
LOCAL_REGIONS_EXPORT DNekScalBlkMatSharedPtr CreateStaticCondMatrix(
const MatrixKey &mkey);
LOCAL_REGIONS_EXPORT virtual void v_ComputeLaplacianMetric();
private:
LibUtilities::NekManager<MatrixKey, DNekScalMat, MatrixKey::opLess> m_matrixManager;
LibUtilities::NekManager<MatrixKey, DNekScalBlkMat, MatrixKey::opLess> m_staticCondMatrixManager;
virtual void v_LaplacianMatrixOp_MatFree_Kernel(
const Array<OneD, const NekDouble> &inarray,
Array<OneD, NekDouble> &outarray,
Array<OneD, NekDouble> &wsp);
};
// type defines for use of PyrExp in a boost vector
......
......@@ -36,11 +36,11 @@
#include <StdRegions/StdPyrExp.h>
#include <LibUtilities/Foundations/ManagerAccess.h>
#include <iomanip>
namespace Nektar
{
namespace StdRegions
{
StdPyrExp::StdPyrExp() // Deafult construct of standard expansion directly called.
{
}
......
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