Commit 5986cdc1 authored by Spencer Sherwin's avatar Spencer Sherwin
Browse files

Removal of m_contCoeffs and m_contNcoeffs and the associated routines this calls


git-svn-id: https://gforge.sci.utah.edu/svn/nektar/trunk@4065 305cdda6-5ce1-45b3-a98d-dfc68c8b3305
parent 646cf42e
......@@ -40,93 +40,92 @@ int main(int argc, char *argv[])
int FFT;
vSession->LoadParameter("HomModesZ", nzpoints);
vSession->LoadParameter("LZ", lz);
vSession->LoadParameter("USEFFT", FFT);
bool useFFT = false;
bool deal = false;
if(FFT==1){useFFT = true;}
const LibUtilities::PointsKey PkeyZ(nzpoints,LibUtilities::eFourierEvenlySpaced);
vSession->LoadParameter("LZ", lz);
vSession->LoadParameter("USEFFT", FFT);
bool useFFT = false;
bool deal = false;
if(FFT==1){useFFT = true;}
const LibUtilities::PointsKey PkeyZ(nzpoints,LibUtilities::eFourierEvenlySpaced);
const LibUtilities::BasisKey BkeyZ(LibUtilities::eFourier,nzpoints,PkeyZ);
Exp_u = MemoryManager<MultiRegions::ContField3DHomogeneous1D>::AllocateSharedPtr(vSession,BkeyZ,lz,useFFT,deal,graph2D,vSession->GetVariable(0));
Exp_v = MemoryManager<MultiRegions::ContField3DHomogeneous1D>::AllocateSharedPtr(vSession,BkeyZ,lz,useFFT,deal,graph2D,vSession->GetVariable(1));
Exp_w = MemoryManager<MultiRegions::ContField3DHomogeneous1D>::AllocateSharedPtr(vSession,BkeyZ,lz,useFFT,deal,graph2D,vSession->GetVariable(2));
//----------------------------------------------
Exp_u = MemoryManager<MultiRegions::ContField3DHomogeneous1D>::AllocateSharedPtr(vSession,BkeyZ,lz,useFFT,deal,graph2D,vSession->GetVariable(0));
Exp_v = MemoryManager<MultiRegions::ContField3DHomogeneous1D>::AllocateSharedPtr(vSession,BkeyZ,lz,useFFT,deal,graph2D,vSession->GetVariable(1));
Exp_w = MemoryManager<MultiRegions::ContField3DHomogeneous1D>::AllocateSharedPtr(vSession,BkeyZ,lz,useFFT,deal,graph2D,vSession->GetVariable(2));
//----------------------------------------------
//----------------------------------------------
// Print summary of solution details
flags.set(eUseContCoeff, false);
flags.set(eUseGlobal, false);
const SpatialDomains::ExpansionMap &expansions = graph2D->GetExpansions();
LibUtilities::BasisKey bkey0 = expansions.begin()->second->m_basisKeyVector[0];
cout << "Calculating Derivatives (Homogeneous in z-plane):" << endl;
cout << " Lz : " << lz << endl;
cout << "Calculating Derivatives (Homogeneous in z-plane):" << endl;
cout << " Lz : " << lz << endl;
cout << " N.modes : " << bkey0.GetNumModes() << endl;
cout << " N.Z homo modes : " << BkeyZ.GetNumModes() << endl;
cout << " N.Z homo modes : " << BkeyZ.GetNumModes() << endl;
cout << endl;
//----------------------------------------------
//----------------------------------------------
// Set up coordinates of mesh for Forcing function evaluation
int dim = Exp_u->GetCoordim(0);
int nq = Exp_u->GetTotPoints();
Array<OneD,NekDouble> xc0,xc1,xc2;
int dim = Exp_u->GetCoordim(0);
int nq = Exp_u->GetTotPoints();
Array<OneD,NekDouble> xc0,xc1,xc2;
xc0 = Array<OneD,NekDouble>(nq,0.0);
xc1 = Array<OneD,NekDouble>(nq,0.0);
xc2 = Array<OneD,NekDouble>(nq,0.0);
Exp_u->GetCoords(xc0,xc1,xc2);
//----------------------------------------------
Array<OneD,NekDouble> dudx,dvdy,dwdz;
Array<OneD,NekDouble> dump;
dump = Array<OneD,NekDouble>(nq,0.0);
dudx = Array<OneD,NekDouble>(nq,0.0);
dvdy = Array<OneD,NekDouble>(nq,0.0);
dwdz = Array<OneD,NekDouble>(nq,0.0);
Array<OneD,NekDouble> dump;
dump = Array<OneD,NekDouble>(nq,0.0);
dudx = Array<OneD,NekDouble>(nq,0.0);
dvdy = Array<OneD,NekDouble>(nq,0.0);
dwdz = Array<OneD,NekDouble>(nq,0.0);
//----------------------------------------------
// Define initial fields
LibUtilities::EquationSharedPtr ffunc_u = vSession->GetFunction("InitialCondition", 0);
LibUtilities::EquationSharedPtr ffunc_v = vSession->GetFunction("InitialCondition", 1);
LibUtilities::EquationSharedPtr ffunc_w = vSession->GetFunction("InitialCondition", 2);
LibUtilities::EquationSharedPtr exac_u = vSession->GetFunction("ExactSolution", 0);
LibUtilities::EquationSharedPtr exac_v = vSession->GetFunction("ExactSolution", 1);
LibUtilities::EquationSharedPtr exac_w = vSession->GetFunction("ExactSolution", 2);
LibUtilities::EquationSharedPtr ffunc_v = vSession->GetFunction("InitialCondition", 1);
LibUtilities::EquationSharedPtr ffunc_w = vSession->GetFunction("InitialCondition", 2);
LibUtilities::EquationSharedPtr exac_u = vSession->GetFunction("ExactSolution", 0);
LibUtilities::EquationSharedPtr exac_v = vSession->GetFunction("ExactSolution", 1);
LibUtilities::EquationSharedPtr exac_w = vSession->GetFunction("ExactSolution", 2);
ffunc_u->Evaluate(xc0,xc1,xc2,Exp_u->UpdatePhys());
ffunc_v->Evaluate(xc0,xc1,xc2,Exp_v->UpdatePhys());
ffunc_w->Evaluate(xc0,xc1,xc2,Exp_w->UpdatePhys());
exac_u->Evaluate(xc0,xc1,xc2,dudx);
exac_v->Evaluate(xc0,xc1,xc2,dvdy);
exac_w->Evaluate(xc0,xc1,xc2,dwdz);
//----------------------------------------------
//Taking derivative and printing the error
Exp_u->PhysDeriv(Exp_u->GetPhys(),Exp_u->UpdatePhys(),dump,dump,false);
cout << "L infinity error: " << Exp_u->Linf(dudx) << endl;
cout << "L 2 error : " << Exp_u->L2 (dudx) << endl;
Exp_v->PhysDeriv(Exp_v->GetPhys(),dump,Exp_v->UpdatePhys(),dump,false);
cout << "L infinity error: " << Exp_v->Linf(dvdy) << endl;
cout << "L 2 error : " << Exp_v->L2 (dvdy) << endl;
Exp_w->PhysDeriv(Exp_w->GetPhys(),dump,dump,Exp_w->UpdatePhys(),false);
cout << "L infinity error: " << Exp_w->Linf(dwdz) << endl;
cout << "L 2 error : " << Exp_w->L2 (dwdz) << endl;
return 0;
//Taking derivative and printing the error
Exp_u->PhysDeriv(Exp_u->GetPhys(),Exp_u->UpdatePhys(),dump,dump);
cout << "L infinity error: " << Exp_u->Linf(dudx) << endl;
cout << "L 2 error : " << Exp_u->L2 (dudx) << endl;
Exp_v->PhysDeriv(Exp_v->GetPhys(),dump,Exp_v->UpdatePhys(),dump);
cout << "L infinity error: " << Exp_v->Linf(dvdy) << endl;
cout << "L 2 error : " << Exp_v->L2 (dvdy) << endl;
Exp_w->PhysDeriv(Exp_w->GetPhys(),dump,dump,Exp_w->UpdatePhys());
cout << "L infinity error: " << Exp_w->Linf(dwdz) << endl;
cout << "L 2 error : " << Exp_w->L2 (dwdz) << endl;
return 0;
}
......@@ -58,7 +58,7 @@ int main(int argc, char *argv[])
//----------------------------------------------
// Print summary of solution details
flags.set(eUseContCoeff, false);
flags.set(eUseGlobal, false);
const SpatialDomains::ExpansionMap &expansions = graph2D->GetExpansions();
......@@ -111,24 +111,24 @@ int main(int argc, char *argv[])
//----------------------------------------------
//Taking derivative and printing the error
cout << "Deriv u" << endl;
Exp_u->PhysDeriv(Exp_u->GetPhys(),Exp_u->UpdatePhys(),dump,dump,false);
cout << "Deriv u done" << endl;
cout << "L infinity error: " << Exp_u->Linf(dudx) << endl;
cout << "L 2 error : " << Exp_u->L2 (dudx) << endl;
Exp_v->PhysDeriv(Exp_v->GetPhys(),dump,Exp_v->UpdatePhys(),dump,false);
cout << "L infinity error: " << Exp_v->Linf(dvdy) << endl;
cout << "L 2 error : " << Exp_v->L2 (dvdy) << endl;
Exp_w->PhysDeriv(Exp_w->GetPhys(),dump,dump,Exp_w->UpdatePhys(),false);
cout << "L infinity error: " << Exp_w->Linf(dwdz) << endl;
cout << "L 2 error : " << Exp_w->L2 (dwdz) << endl;
return 0;
//Taking derivative and printing the error
cout << "Deriv u" << endl;
Exp_u->PhysDeriv(Exp_u->GetPhys(),Exp_u->UpdatePhys(),dump,dump);
cout << "Deriv u done" << endl;
cout << "L infinity error: " << Exp_u->Linf(dudx) << endl;
cout << "L 2 error : " << Exp_u->L2 (dudx) << endl;
Exp_v->PhysDeriv(Exp_v->GetPhys(),dump,Exp_v->UpdatePhys(),dump);
cout << "L infinity error: " << Exp_v->Linf(dvdy) << endl;
cout << "L 2 error : " << Exp_v->L2 (dvdy) << endl;
Exp_w->PhysDeriv(Exp_w->GetPhys(),dump,dump,Exp_w->UpdatePhys());
cout << "L infinity error: " << Exp_w->Linf(dwdz) << endl;
cout << "L 2 error : " << Exp_w->L2 (dwdz) << endl;
return 0;
}
......@@ -66,7 +66,7 @@ int main(int argc, char *argv[])
//----------------------------------------------
// Print summary of solution details
flags.set(eUseContCoeff, false);
flags.set(eUseGlobal, false);
const SpatialDomains::ExpansionMap &expansions = graph1D->GetExpansions();
......@@ -122,17 +122,17 @@ int main(int argc, char *argv[])
//Taking derivative and printing the error
Exp_u->PhysDeriv(Exp_u->GetPhys(),Exp_u->UpdatePhys(),dump,dump,false);
Exp_u->PhysDeriv(Exp_u->GetPhys(),Exp_u->UpdatePhys(),dump,dump);
cout << "L infinity error: " << Exp_u->Linf(dudx) << endl;
cout << "L 2 error : " << Exp_u->L2 (dudx) << endl;
Exp_v->PhysDeriv(Exp_v->GetPhys(),dump,Exp_v->UpdatePhys(),dump,false);
Exp_v->PhysDeriv(Exp_v->GetPhys(),dump,Exp_v->UpdatePhys(),dump);
cout << "L infinity error: " << Exp_v->Linf(dvdy) << endl;
cout << "L 2 error : " << Exp_v->L2 (dvdy) << endl;
Exp_w->PhysDeriv(Exp_w->GetPhys(),dump,dump,Exp_w->UpdatePhys(),false);
Exp_w->PhysDeriv(Exp_w->GetPhys(),dump,dump,Exp_w->UpdatePhys());
cout << "L infinity error: " << Exp_w->Linf(dwdz) << endl;
cout << "L 2 error : " << Exp_w->L2 (dwdz) << endl;
......
......@@ -103,12 +103,12 @@ int main(int argc, char *argv[])
//---------------------------------------------
// Project onto Expansion
Exp->FwdTrans(Fce->GetPhys(), Exp->UpdateCoeffs(), true);
Exp->FwdTrans(Fce->GetPhys(), Exp->UpdateCoeffs(), MultiRegions::eGlobal);
//---------------------------------------------
//-------------------------------------------
// Backward Transform Solution to get projected values
Exp->BwdTrans(Exp->GetCoeffs(), Exp->UpdatePhys(), true);
Exp->BwdTrans(Exp->GetCoeffs(), Exp->UpdatePhys(), MultiRegions::eGlobal);
//-------------------------------------------
//----------------------------------------------
......
......@@ -121,13 +121,13 @@ int main(int argc, char *argv[])
//----------------------------------------------
// Helmholtz solution taking physical forcing
Exp->LinearAdvectionReactionSolve(Vel, Fce->GetPhys(), Exp->UpdateContCoeffs(), lambda, true);
Exp->LinearAdvectionReactionSolve(Vel, Fce->GetPhys(), Exp->UpdateCoeffs(), lambda, MultiRegions::eGlobal);
//----------------------------------------------
Timing("Linear Advection Solve ..");
//----------------------------------------------
// Backward Transform Solution to get solved values
Exp->BwdTrans(Exp->GetContCoeffs(), Exp->UpdatePhys(), true);
Exp->BwdTrans(Exp->GetCoeffs(), Exp->UpdatePhys(), MultiRegions::eGlobal);
//----------------------------------------------
//----------------------------------------------
......
......@@ -50,7 +50,7 @@ int main(int argc, char *argv[])
//----------------------------------------------
// Print summary of solution details
flags.set(eUseContCoeff, true);
flags.set(eUseGlobal, true);
factors[StdRegions::eFactorLambda] = vSession->GetParameter("Lambda");
const SpatialDomains::ExpansionMap &expansions = graph2D->GetExpansions();
LibUtilities::BasisKey bkey0 = expansions.begin()->second->m_basisKeyVector[0];
......@@ -128,14 +128,14 @@ int main(int argc, char *argv[])
//----------------------------------------------
// Helmholtz solution taking physical forcing
Exp->HelmSolve(Fce->GetPhys(), Exp->UpdateContCoeffs(), flags, factors, varcoeffs);
Exp->HelmSolve(Fce->GetPhys(), Exp->UpdateCoeffs(), flags, factors, varcoeffs);
//----------------------------------------------
Timing("Helmholtz Solve ..");
#ifdef TIMING
for(i = 0; i < 1000; ++i)
{
Exp->HelmSolve(Fce->GetPhys(), Exp->UpdateContCoeffs(), flags, factors, varcoeffs);
Exp->HelmSolve(Fce->GetPhys(), Exp->UpdateCoeffs(), flags, factors, varcoeffs);
}
Timing("1000 Helmholtz Solves:... ");
......@@ -143,7 +143,7 @@ int main(int argc, char *argv[])
//----------------------------------------------
// Backward Transform Solution to get solved values
Exp->BwdTrans(Exp->GetContCoeffs(), Exp->UpdatePhys(), true);
Exp->BwdTrans(Exp->GetCoeffs(), Exp->UpdatePhys(), MultiRegions::eGlobal);
//----------------------------------------------
//-----------------------------------------------
......@@ -157,7 +157,7 @@ int main(int argc, char *argv[])
= Exp->GetFieldDefinitions();
std::vector<std::vector<NekDouble> > FieldData(FieldDef.size());
Exp->GlobalToLocal(Exp->GetContCoeffs(),Exp->UpdateCoeffs());
Exp->GlobalToLocal(Exp->GetCoeffs(),Exp->UpdateCoeffs());
for(i = 0; i < FieldDef.size(); ++i)
{
FieldDef[i]->m_fields.push_back("u");
......
......@@ -39,7 +39,7 @@ int main(int argc, char *argv[])
//----------------------------------------------
// Print summary of solution details
flags.set(eUseContCoeff, true);
flags.set(eUseGlobal, true);
factors[StdRegions::eFactorLambda] = vSession->GetParameter("Lambda");
const SpatialDomains::ExpansionMap &expansions = graph3D->GetExpansions();
LibUtilities::BasisKey bkey0 = expansions.begin()->second->m_basisKeyVector[0];
......@@ -93,12 +93,12 @@ int main(int argc, char *argv[])
//----------------------------------------------
// Helmholtz solution taking physical forcing
Exp->HelmSolve(Fce->GetPhys(), Exp->UpdateContCoeffs(), flags, factors);
Exp->HelmSolve(Fce->GetPhys(), Exp->UpdateCoeffs(), flags, factors);
//----------------------------------------------
//----------------------------------------------
// Backward Transform Solution to get solved values at
Exp->BwdTrans(Exp->GetContCoeffs(), Exp->UpdatePhys(), true);
Exp->BwdTrans(Exp->GetCoeffs(), Exp->UpdatePhys(), MultiRegions::eGlobal);
//----------------------------------------------
//----------------------------------------------
......@@ -118,7 +118,7 @@ int main(int argc, char *argv[])
= Exp->GetFieldDefinitions();
std::vector<std::vector<NekDouble> > FieldData(FieldDef.size());
Exp->GlobalToLocal(Exp->GetContCoeffs(),Exp->UpdateCoeffs());
Exp->GlobalToLocal(Exp->GetCoeffs(),Exp->UpdateCoeffs());
for(i = 0; i < FieldDef.size(); ++i)
{
FieldDef[i]->m_fields.push_back("u");
......
......@@ -51,11 +51,11 @@ int main(int argc, char *argv[])
//----------------------------------------------
// Print summary of solution details
flags.set(eUseContCoeff, true);
flags.set(eUseGlobal, true);
factors[StdRegions::eFactorLambda] = vSession->GetParameter("Lambda");
const SpatialDomains::ExpansionMap &expansions = graph2D->GetExpansions();
LibUtilities::BasisKey bkey0
= expansions.begin()->second->m_basisKeyVector[0];
= expansions.begin()->second->m_basisKeyVector[0];
cout << "Solving 3D Helmholtz (Homogeneous in z-direction):" << endl;
cout << " Lambda : " << factors[StdRegions::eFactorLambda] << endl;
......@@ -94,16 +94,14 @@ int main(int argc, char *argv[])
//----------------------------------------------
// Helmholtz solution taking physical forcing
//Exp->HelmSolve(Fce->GetPhys(), Exp->UpdateCoeffs(), lambda);
Exp->HelmSolve(Fce->GetPhys(), Exp->UpdateContCoeffs(), flags, factors);
Exp->HelmSolve(Fce->GetPhys(), Exp->UpdateCoeffs(), flags, factors);
//----------------------------------------------
//----------------------------------------------
// Backward Transform Solution to get solved values at
//Exp->BwdTrans(Exp->GetCoeffs(), Exp->UpdatePhys());
Exp->BwdTrans(Exp->GetContCoeffs(), Exp->UpdatePhys(),true);
Exp->BwdTrans(Exp->GetCoeffs(), Exp->UpdatePhys(),MultiRegions::eGlobal);
//----------------------------------------------
//----------------------------------------------
// See if there is an exact solution, if so
// evaluate and plot errors
......
......@@ -47,40 +47,40 @@ int main(int argc, char *argv[])
vSession->LoadParameter("HomModesY", nypoints);
vSession->LoadParameter("HomModesZ", nzpoints);
vSession->LoadParameter("LY", ly);
vSession->LoadParameter("LZ", lz);
vSession->LoadParameter("USEFFT", FFT);
bool useFFT = false;
bool deal = false;
if(FFT==1){useFFT = true;}
vSession->LoadParameter("LY", ly);
vSession->LoadParameter("LZ", lz);
vSession->LoadParameter("USEFFT", FFT);
bool useFFT = false;
bool deal = false;
if(FFT==1){useFFT = true;}
const LibUtilities::PointsKey PkeyY(nypoints,LibUtilities::eFourierEvenlySpaced);
const LibUtilities::PointsKey PkeyY(nypoints,LibUtilities::eFourierEvenlySpaced);
const LibUtilities::BasisKey BkeyY(LibUtilities::eFourier,nypoints,PkeyY);
const LibUtilities::PointsKey PkeyZ(nzpoints,LibUtilities::eFourierEvenlySpaced);
const LibUtilities::PointsKey PkeyZ(nzpoints,LibUtilities::eFourierEvenlySpaced);
const LibUtilities::BasisKey BkeyZ(LibUtilities::eFourier,nzpoints,PkeyZ);
Exp = MemoryManager<MultiRegions::ContField3DHomogeneous2D>::AllocateSharedPtr(vSession,BkeyY,BkeyZ,ly,lz,useFFT,deal,graph1D,vSession->GetVariable(0));
Exp = MemoryManager<MultiRegions::ContField3DHomogeneous2D>::AllocateSharedPtr(vSession,BkeyY,BkeyZ,ly,lz,useFFT,deal,graph1D,vSession->GetVariable(0));
//----------------------------------------------
//----------------------------------------------
// Print summary of solution details
flags.set(eUseContCoeff, true);
factors[StdRegions::eFactorLambda] = vSession->GetParameter("Lambda");
flags.set(eUseGlobal, true);
factors[StdRegions::eFactorLambda] = vSession->GetParameter("Lambda");
const SpatialDomains::ExpansionMap &expansions = graph1D->GetExpansions();
LibUtilities::BasisKey bkey0 = expansions.begin()->second->m_basisKeyVector[0];
cout << "Solving 3D Helmholtz (Homogeneous in yz-plane):" << endl;
cout << "Solving 3D Helmholtz (Homogeneous in yz-plane):" << endl;
cout << " Lambda : " << factors[StdRegions::eFactorLambda] << endl;
cout << " Ly : " << ly << endl;
cout << " Lz : " << lz << endl;
cout << " Ly : " << ly << endl;
cout << " Lz : " << lz << endl;
cout << " N.modes : " << bkey0.GetNumModes() << endl;
cout << " N.Y homo modes : " << BkeyY.GetNumModes() << endl;
cout << " N.Z homo modes : " << BkeyZ.GetNumModes() << endl;
cout << " N.Z homo modes : " << BkeyZ.GetNumModes() << endl;
cout << endl;
//----------------------------------------------
......@@ -115,13 +115,13 @@ int main(int argc, char *argv[])
//----------------------------------------------
// Helmholtz solution taking physical forcing
//Exp->HelmSolve(Fce->GetPhys(), Exp->UpdateCoeffs(), lambda);
Exp->HelmSolve(Fce->GetPhys(), Exp->UpdateContCoeffs(), flags, factors);
Exp->HelmSolve(Fce->GetPhys(), Exp->UpdateCoeffs(), flags, factors);
//----------------------------------------------
//----------------------------------------------
// Backward Transform Solution to get solved values at
//Exp->BwdTrans(Exp->GetCoeffs(), Exp->UpdatePhys());
Exp->BwdTrans(Exp->GetContCoeffs(), Exp->UpdatePhys(),true);
Exp->BwdTrans(Exp->GetCoeffs(), Exp->UpdatePhys(),MultiRegions::eGlobal);
//----------------------------------------------
//----------------------------------------------
......
......@@ -77,12 +77,12 @@ namespace Nektar
/// Enumeration of flags for passing a list of options.
enum FlagType
{
eUseContCoeff
eUseGlobal
};
/// String map for FlagType enumeration.
const char* const FlagTypeMap[] = {
"UseContCoeff"
"UseGlobal"
};
/// Defines a list of flags.
......
......@@ -361,13 +361,24 @@ namespace Nektar
const Array<OneD, const NekDouble>& loc,
Array<OneD, NekDouble>& global) const
{
Array<OneD, const NekDouble> local;
if(global.data() == loc.data())
{
local = Array<OneD, NekDouble>(local.num_elements(),local.data());
}
else
{
local = loc; // create reference
}
if(m_signChange)
{
Vmath::Scatr(m_numLocalCoeffs, m_localToGlobalSign.get(), loc.get(), m_localToGlobalMap.get(), global.get());
Vmath::Scatr(m_numLocalCoeffs, m_localToGlobalSign.get(), local.get(), m_localToGlobalMap.get(), global.get());
}
else
{
Vmath::Scatr(m_numLocalCoeffs, loc.get(), m_localToGlobalMap.get(), global.get());
Vmath::Scatr(m_numLocalCoeffs, local.get(), m_localToGlobalMap.get(), global.get());
}
}
......@@ -382,13 +393,24 @@ namespace Nektar
const Array<OneD, const NekDouble>& global,
Array<OneD, NekDouble>& loc) const
{
Array<OneD, const NekDouble> glo;
if(global.data() == loc.data())
{
glo = Array<OneD, NekDouble>(global.num_elements(),global.data());
}
else
{
glo = global; // create reference
}
if(m_signChange)
{
Vmath::Gathr(m_numLocalCoeffs, m_localToGlobalSign.get(), global.get(), m_localToGlobalMap.get(), loc.get());
Vmath::Gathr(m_numLocalCoeffs, m_localToGlobalSign.get(), glo.get(), m_localToGlobalMap.get(), loc.get());
}
else
{
Vmath::Gathr(m_numLocalCoeffs, global.get(), m_localToGlobalMap.get(), loc.get());
Vmath::Gathr(m_numLocalCoeffs, glo.get(), m_localToGlobalMap.get(), loc.get());
}
}
......@@ -403,17 +425,26 @@ namespace Nektar
const Array<OneD, const NekDouble> &loc,
Array<OneD, NekDouble> &global) const
{
ASSERTL1(loc.get() != global.get(),"Local and Global Arrays cannot be the same");
Array<OneD, const NekDouble> local;
if(global.data() == loc.data())
{
local = Array<OneD, NekDouble>(local.num_elements(),local.data());
}
else
{
local = loc; // create reference
}
//ASSERTL1(loc.get() != global.get(),"Local and Global Arrays cannot be the same");
Vmath::Zero(m_numGlobalCoeffs, global.get(), 1);
if(m_signChange)
{
Vmath::Assmb(m_numLocalCoeffs, m_localToGlobalSign.get(), loc.get(), m_localToGlobalMap.get(), global.get());
Vmath::Assmb(m_numLocalCoeffs, m_localToGlobalSign.get(), local.get(), m_localToGlobalMap.get(), global.get());
}
else
{
Vmath::Assmb(m_numLocalCoeffs, loc.get(), m_localToGlobalMap.get(), global.get());
Vmath::Assmb(m_numLocalCoeffs, local.get(), m_localToGlobalMap.get(), global.get());
}
UniversalAssemble(global);
}
......
......@@ -85,8 +85,6 @@ namespace Nektar
ContField1D::ContField1D():
DisContField1D(),
m_locToGloMap(),
m_contNcoeffs(0),
m_contCoeffs(),
m_globalLinSysManager(
boost::bind(&ContField1D::GenGlobalLinSys, this, _1),
std::string("GlobalLinSys"))
......@@ -104,7 +102,7 @@ namespace Nektar
* constructs the mapping array (contained in #m_locToGloMap) for the
* transformation between local elemental level and global level, it
* calculates the total number global expansio