Commit 4e17d3d2 authored by Sergey Yakovlev's avatar Sergey Yakovlev
Browse files

Optimized trace term computation in Expansion3D::v_GenMatrix

Added weak derivative matrix code to TetExp::v_GenMatrix


git-svn-id: https://gforge.sci.utah.edu/svn/nektar/trunk@4055 305cdda6-5ce1-45b3-a98d-dfc68c8b3305
parent ea195f58
......@@ -300,7 +300,6 @@ namespace Nektar
ASSERTL1(map1.num_elements() == map2.num_elements(), "There is an error with the GetFaceToElementMap");
//if(dir1 != dir2) cout << "face " << face << endl;
for(j = 0; j < map1.num_elements(); ++j)//index in the standard orientation
for(k = 0; k < map2.num_elements(); ++k)//index in the actual orientation
{
......@@ -469,38 +468,82 @@ namespace Nektar
Array<OneD,unsigned int> fmap;
Array<OneD,int> sign;
// Set up face segment expansions from local geom info
for(i = 0; i < nfaces; ++i)
{
FaceExp[i] = GetFaceExp(i);
}
// for each degree of freedom of the lambda space
// calculate Umat entry
// Generate Lambda to U_lambda matrix
for(j = 0; j < nbndry; ++j)
{
// standard basis vectors e_j
Vmath::Zero(nbndry,&lambda[0],1);
Vmath::Zero(ncoeffs,&f[0],1);
lambda[j] = 1.0;
//cout << Lambda;
SetTraceToGeomOrientation(lambda);
//cout << Lambda << endl;
// Compute F = [I D_1 M^{-1} D_2 M^{-1}] C e_j
AddHDGHelmholtzTraceTerms(tau, lambda, FaceExp, mkey.GetVarCoeffs(), f);
// Compute U^e_j
Ulam = invHmat*F; // generate Ulam from lambda
// fill column of matrix
for(k = 0; k < ncoeffs; ++k)
{
Umat(k,j) = Ulam[k];
}
}
//alternative way to add boundary terms contribution
int bndry_cnt = 0;
for(i = 0; i < nfaces; ++i)
{
FaceExp[i] = GetFaceExp(i);//temporary, need to rewrite AddHDGHelmholtzFaceTerms
int nface = GetFaceNcoeffs(i);
Array<OneD, NekDouble> face_lambda(nface);
const Array<OneD, const Array<OneD, NekDouble> > normals
= GetFaceNormal(i);
//cout << endl << "face #" << i;
//cout << endl << "nquad_f " << FaceExp[i]->GetNumPoints(0)*FaceExp[i]->GetNumPoints(1);
//cout << endl << "normals[0] " << normals[0].num_elements();
//cout << endl << "normals[1] " << normals[1].num_elements();
//cout << endl << "normals[2] " << normals[2].num_elements();
for(j = 0; j < nface; ++j)
{
Vmath::Zero(nface,&face_lambda[0],1);
Vmath::Zero(ncoeffs,&f[0],1);
face_lambda[j] = 1.0;
SetFaceToGeomOrientation(i, face_lambda);
Vmath::Vcopy(nface, face_lambda, 1, FaceExp[i]->UpdateCoeffs(), 1);
FaceExp[i]->BwdTrans(FaceExp[i]->GetCoeffs(), FaceExp[i]->UpdatePhys());
AddHDGHelmholtzFaceTerms(tau, i, FaceExp, mkey.GetVarCoeffs(), f);
Ulam = invHmat*F; // generate Ulam from lambda
// fill column of matrix
for(k = 0; k < ncoeffs; ++k)
{
Umat(k,bndry_cnt) = Ulam[k];
}
++bndry_cnt;
}
}
//// Set up face expansions from local geom info
//for(i = 0; i < nfaces; ++i)
//{
// FaceExp[i] = GetFaceExp(i);
//}
//
//// for each degree of freedom of the lambda space
//// calculate Umat entry
//// Generate Lambda to U_lambda matrix
//for(j = 0; j < nbndry; ++j)
//{
// // standard basis vectors e_j
// Vmath::Zero(nbndry,&lambda[0],1);
// Vmath::Zero(ncoeffs,&f[0],1);
// lambda[j] = 1.0;
//
// //cout << Lambda;
// SetTraceToGeomOrientation(lambda);
// //cout << Lambda << endl;
//
// // Compute F = [I D_1 M^{-1} D_2 M^{-1}] C e_j
// AddHDGHelmholtzTraceTerms(tau, lambda, FaceExp, mkey.GetVarCoeffs(), f);
//
// // Compute U^e_j
// Ulam = invHmat*F; // generate Ulam from lambda
//
// // fill column of matrix
// for(k = 0; k < ncoeffs; ++k)
// {
// Umat(k,j) = Ulam[k];
// }
//}
}
break;
// Q_0, Q_1, Q_2 matrices (P23)
......
......@@ -1393,6 +1393,64 @@ namespace Nektar
}
}
break;
case StdRegions::eWeakDeriv0:
case StdRegions::eWeakDeriv1:
case StdRegions::eWeakDeriv2:
{
if(m_metricinfo->GetGtype() == SpatialDomains::eDeformed ||
mkey.GetNVarCoeff())
{
NekDouble one = 1.0;
DNekMatSharedPtr mat = GenMatrix(mkey);
returnval = MemoryManager<DNekScalMat>
::AllocateSharedPtr(one,mat);
}
else
{
NekDouble jac = (m_metricinfo->GetJac())[0];
Array<TwoD, const NekDouble> gmat
= m_metricinfo->GetGmat();
int dir;
switch(mkey.GetMatrixType())
{
case StdRegions::eWeakDeriv0:
dir = 0;
break;
case StdRegions::eWeakDeriv1:
dir = 1;
break;
case StdRegions::eWeakDeriv2:
dir = 2;
break;
}
MatrixKey deriv0key(StdRegions::eWeakDeriv0,
mkey.GetExpansionType(), *this);
MatrixKey deriv1key(StdRegions::eWeakDeriv1,
mkey.GetExpansionType(), *this);
MatrixKey deriv2key(StdRegions::eWeakDeriv2,
mkey.GetExpansionType(), *this);
DNekMat &deriv0 = *GetStdMatrix(deriv0key);
DNekMat &deriv1 = *GetStdMatrix(deriv1key);
DNekMat &deriv2 = *GetStdMatrix(deriv2key);
int rows = deriv0.GetRows();
int cols = deriv1.GetColumns();
DNekMatSharedPtr WeakDeriv = MemoryManager<DNekMat>
::AllocateSharedPtr(rows,cols);
(*WeakDeriv) = gmat[3*dir][0]*deriv0
+ gmat[3*dir+1][0]*deriv1
+ gmat[3*dir+2][0]*deriv2;
returnval = MemoryManager<DNekScalMat>
::AllocateSharedPtr(jac,WeakDeriv);
}
}
break;
case StdRegions::eLaplacian:
{
if(m_metricinfo->GetGtype() == SpatialDomains::eDeformed ||
......
Supports Markdown
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