Commit c18a5128 authored by Giacomo Castiglioni's avatar Giacomo Castiglioni Committed by Spencer Sherwin

Feature/array mod

parent 49b288f4
......@@ -6,6 +6,8 @@ v5.1.0
**Library**
- Refactored time integration code using factory pattern (!1034)
- Fix to preprocessor logic for boost with Visual Studio >= 2015 (!1115)
- Fix type consistency and real comparison in SharedArray.hpp, replaced
num_elements with size() (!1127)
**CardiacEPSolver**
- Added additional parameter sets to Fenton-Karma model (!1119)
......
......@@ -437,7 +437,7 @@ class BwdTrans_SumFac_Quad : public Operator
}
else
{
ASSERTL1(wsp.num_elements() == m_wspSize,
ASSERTL1(wsp.size() == m_wspSize,
"Incorrect workspace size");
// Those two calls correpsond to the operation
......@@ -523,7 +523,7 @@ class BwdTrans_SumFac_Tri : public Operator
{
boost::ignore_unused(output1, output2);
ASSERTL1(wsp.num_elements() == m_wspSize,
ASSERTL1(wsp.size() == m_wspSize,
"Incorrect workspace size");
int ncoeffs = m_stdExp->GetNcoeffs();
......@@ -637,7 +637,7 @@ class BwdTrans_SumFac_Hex : public Operator
}
else
{
ASSERTL1(wsp.num_elements() == m_wspSize,
ASSERTL1(wsp.size() == m_wspSize,
"Incorrect workspace size");
// Assign second half of workspace for 2nd DGEMM operation.
......@@ -744,7 +744,7 @@ class BwdTrans_SumFac_Tet : public Operator
{
boost::ignore_unused(output1, output2);
ASSERTL1(wsp.num_elements() == m_wspSize,
ASSERTL1(wsp.size() == m_wspSize,
"Incorrect workspace size");
Array<OneD, NekDouble > tmp = wsp;
......@@ -921,7 +921,7 @@ class BwdTrans_SumFac_Prism : public Operator
{
boost::ignore_unused(output1, output2);
ASSERTL1(wsp.num_elements() == m_wspSize,
ASSERTL1(wsp.size() == m_wspSize,
"Incorrect workspace size");
// Assign second half of workspace for 2nd DGEMM operation.
......@@ -1065,7 +1065,7 @@ class BwdTrans_SumFac_Pyr : public Operator
{
boost::ignore_unused(output1, output2);
ASSERTL1(wsp.num_elements() == m_wspSize,
ASSERTL1(wsp.size() == m_wspSize,
"Incorrect workspace size");
// Assign second half of workspace for 2nd DGEMM operation.
......
......@@ -441,14 +441,14 @@ void PyrIProduct(bool sortTopVertex, int numElmt,
int cnt;
int mode, mode1;
ASSERTL1(wsp.num_elements() >= numElmt*(nquad1*nquad2*nmodes0 +
ASSERTL1(wsp.size() >= numElmt*(nquad1*nquad2*nmodes0 +
nquad2*max(nquad0*nquad1,nmodes0*nmodes1)),
"Insufficient workspace size");
Vmath::Vmul(numElmt*totpoints,jac,1,input,1,wsp,1);
Array<OneD, NekDouble> wsp1 = wsp + numElmt * nquad2
* (max(nquad0*nquad1,
* (max(nquad0*nquad1,
nmodes0*nmodes1));
// Perform iproduct with respect to the '0' direction
......@@ -489,7 +489,7 @@ void PyrIProduct(bool sortTopVertex, int numElmt,
mode += nmodes2-ijmax;
}
}
// fix for modified basis for top singular vertex component
// Already have evaluated (1+c)/2 (1-b)/2 (1-a)/2
if(sortTopVertex)
......
......@@ -72,10 +72,10 @@ class IProductWRTBase_StdMat : public Operator
{
boost::ignore_unused(output1, output2);
ASSERTL1(wsp.num_elements() == m_wspSize,
ASSERTL1(wsp.size() == m_wspSize,
"Incorrect workspace size");
Vmath::Vmul(m_jac.num_elements(),m_jac,1,input,1,wsp,1);
Vmath::Vmul(m_jac.size(),m_jac,1,input,1,wsp,1);
Blas::Dgemm('N', 'N', m_mat->GetRows(), m_numElmt,
m_mat->GetColumns(), 1.0, m_mat->GetRawPtr(),
......@@ -170,14 +170,14 @@ class IProductWRTBase_IterPerExp : public Operator
{
boost::ignore_unused(output1, output2);
ASSERTL1(wsp.num_elements() == m_wspSize,
ASSERTL1(wsp.size() == m_wspSize,
"Incorrect workspace size");
const int nCoeffs = m_stdExp->GetNcoeffs();
const int nPhys = m_stdExp->GetTotPoints();
Array<OneD, NekDouble> tmp;
Vmath::Vmul(m_jac.num_elements(),m_jac,1,input,1,wsp,1);
Vmath::Vmul(m_jac.size(),m_jac,1,input,1,wsp,1);
for (int i = 0; i < m_numElmt; ++i)
{
......@@ -463,7 +463,7 @@ class IProductWRTBase_SumFac_Quad : public Operator
{
boost::ignore_unused(output1, output2);
ASSERTL1(wsp.num_elements() == m_wspSize,
ASSERTL1(wsp.size() == m_wspSize,
"Incorrect workspace size");
QuadIProduct(m_colldir0,m_colldir1,m_numElmt,
......@@ -542,7 +542,7 @@ class IProductWRTBase_SumFac_Tri : public Operator
{
boost::ignore_unused(output1, output2);
ASSERTL1(wsp.num_elements() == m_wspSize,
ASSERTL1(wsp.size() == m_wspSize,
"Incorrect workspace size");
TriIProduct(m_sortTopVertex, m_numElmt, m_nquad0, m_nquad1,
......@@ -625,7 +625,7 @@ class IProductWRTBase_SumFac_Hex : public Operator
{
boost::ignore_unused(output1, output2);
ASSERTL1(wsp.num_elements() == m_wspSize,
ASSERTL1(wsp.size() == m_wspSize,
"Incorrect workspace size");
HexIProduct(m_colldir0,m_colldir1,m_colldir2, m_numElmt,
......@@ -714,7 +714,7 @@ class IProductWRTBase_SumFac_Tet : public Operator
{
boost::ignore_unused(output1, output2);
ASSERTL1(wsp.num_elements() == m_wspSize,
ASSERTL1(wsp.size() == m_wspSize,
"Incorrect workspace size");
TetIProduct(m_sortTopEdge, m_numElmt,
......@@ -809,7 +809,7 @@ class IProductWRTBase_SumFac_Prism : public Operator
{
boost::ignore_unused(output1, output2);
ASSERTL1(wsp.num_elements() == m_wspSize,
ASSERTL1(wsp.size() == m_wspSize,
"Incorrect workspace size");
PrismIProduct(m_sortTopVertex, m_numElmt,
......@@ -904,7 +904,7 @@ class IProductWRTBase_SumFac_Pyr : public Operator
{
boost::ignore_unused(output1, output2);
ASSERTL1(wsp.num_elements() == m_wspSize,
ASSERTL1(wsp.size() == m_wspSize,
"Incorrect workspace size");
PyrIProduct(m_sortTopVertex, m_numElmt,
......
......@@ -525,9 +525,9 @@ class PhysDeriv_SumFac_Seg : public Operator
{
const int nqcol = m_nquad0*m_numElmt;
ASSERTL1(wsp.num_elements() == m_wspSize,
ASSERTL1(wsp.size() == m_wspSize,
"Incorrect workspace size");
ASSERTL1(input.num_elements() >= nqcol,
ASSERTL1(input.size() >= nqcol,
"Incorrect input size");
Array<OneD, NekDouble> diff0(nqcol, wsp);
......@@ -558,9 +558,9 @@ class PhysDeriv_SumFac_Seg : public Operator
{
const int nqcol = m_nquad0*m_numElmt;
ASSERTL1(wsp.num_elements() == m_wspSize,
ASSERTL1(wsp.size() == m_wspSize,
"Incorrect workspace size");
ASSERTL1(input.num_elements() >= nqcol,
ASSERTL1(input.size() >= nqcol,
"Incorrect input size");
Array<OneD, NekDouble> diff0(nqcol, wsp);
......@@ -627,9 +627,9 @@ class PhysDeriv_SumFac_Quad : public Operator
const int nqtot = m_nquad0 * m_nquad1;
const int nqcol = nqtot*m_numElmt;
ASSERTL1(wsp.num_elements() == m_wspSize,
ASSERTL1(wsp.size() == m_wspSize,
"Incorrect workspace size");
ASSERTL1(input.num_elements() >= nqcol,
ASSERTL1(input.size() >= nqcol,
"Incorrect input size");
Array<OneD, NekDouble> diff0(nqcol, wsp );
......@@ -673,9 +673,9 @@ class PhysDeriv_SumFac_Quad : public Operator
const int nqtot = m_nquad0 * m_nquad1;
const int nqcol = nqtot*m_numElmt;
ASSERTL1(wsp.num_elements() == m_wspSize,
ASSERTL1(wsp.size() == m_wspSize,
"Incorrect workspace size");
ASSERTL1(input.num_elements() >= nqcol,
ASSERTL1(input.size() >= nqcol,
"Incorrect input size");
Array<OneD, NekDouble> diff0(nqcol, wsp );
......@@ -756,9 +756,9 @@ class PhysDeriv_SumFac_Tri : public Operator
const int nqtot = m_nquad0 * m_nquad1;
const int nqcol = nqtot*m_numElmt;
ASSERTL1(wsp.num_elements() == m_wspSize,
ASSERTL1(wsp.size() == m_wspSize,
"Incorrect workspace size");
ASSERTL1(input.num_elements() >= nqcol,
ASSERTL1(input.size() >= nqcol,
"Incorrect input size");
Array<OneD, NekDouble> diff0(nqcol, wsp );
......@@ -813,9 +813,9 @@ class PhysDeriv_SumFac_Tri : public Operator
const int nqtot = m_nquad0 * m_nquad1;
const int nqcol = nqtot*m_numElmt;
ASSERTL1(wsp.num_elements() == m_wspSize,
ASSERTL1(wsp.size() == m_wspSize,
"Incorrect workspace size");
ASSERTL1(input.num_elements() >= nqcol,
ASSERTL1(input.size() >= nqcol,
"Incorrect input size");
Array<OneD, NekDouble> diff0(nqcol, wsp );
......@@ -1598,11 +1598,11 @@ class PhysDeriv_SumFac_Pyr : public Operator
// dxi1 = 2/(1-eta_2) d Eta_1
Vmath::Vmul(nPhys,&m_fac0[0],1,Diff[1].get()+cnt,1,
Diff[1].get()+cnt,1);
// dxi2 = (1+eta0)/(1-eta_2) d Eta_0 + d/dEta2;
Vmath::Vvtvp(nPhys,&m_fac1[0],1,Diff[0].get()+cnt,1,
Diff[2].get()+cnt,1,Diff[2].get()+cnt,1);
// dxi2 += (1+eta1)/(1-eta_2) d Eta_1
// dxi2 += (1+eta1)/(1-eta_2) d Eta_1
Vmath::Vvtvp(nPhys,&m_fac2[0],1,Diff[1].get()+cnt,1,
Diff[2].get()+cnt,1,Diff[2].get()+cnt,1);
cnt += nPhys;
......@@ -1667,7 +1667,7 @@ class PhysDeriv_SumFac_Pyr : public Operator
// dxi1 = 2/(1-eta_2) d Eta_1
Vmath::Vmul(nPhys,&m_fac0[0],1,Diff[1].get()+cnt,1,
Diff[1].get()+cnt,1);
// dxi2 = (1+eta0)/(1-eta_2) d Eta_0 + d/dEta2;
Vmath::Vvtvp(nPhys,&m_fac1[0],1,Diff[0].get()+cnt,1,
Diff[2].get()+cnt,1,Diff[2].get()+cnt,1);
......
......@@ -81,7 +81,7 @@ long double integrationFunction(int nPts, PointsType type)
integral = 0;
break;
default:
integral = (long double)20.0/3.0;
integral = (long double)20.0/3.0;
}
return integral;
}
......@@ -95,7 +95,7 @@ long double integrandWeightFunction(long double x, PointsType type)
case eGaussGaussChebyshev:
case eGaussRadauMChebyshev:
case eGaussRadauPChebyshev:
case eGaussLobattoChebyshev:
case eGaussLobattoChebyshev:
weight = 1.0 / sqrt(1.0 - x*x);
break;
......@@ -103,10 +103,10 @@ long double integrandWeightFunction(long double x, PointsType type)
// weight = 1.0 + x; // ?
break;
case eGaussRadauMAlpha0Beta2:
case eGaussRadauMAlpha0Beta2:
// weight = (1.0 + x)*(1.0 + x); // ?
break;
default:
weight = 1.0;
......@@ -152,7 +152,7 @@ int main(int argc, char *argv[])
// Show the set up to the user
cout << "Points type: " << kPointsTypeStr[pointsType] << endl;
cout << "Number of points: " << nPts << endl;
// Display the example test function to the user
if( pointsType == eFourierEvenlySpaced )
{
......@@ -167,12 +167,12 @@ int main(int argc, char *argv[])
// Obtain a reference to a Points object via an appropriate PointsKey object
PointsKey key(nPts, pointsType);
PointsSharedPtr points = PointsManager()[key];
//std::shared_ptr<Points<NekDouble> > points = PointsManager()[key];
//const ptr<Points<NekDouble> > points = PointsManager()[key];
// Get the abscissas and their matching quadrature weights
// SharedArray<const NekDouble> z, w;
Array<OneD, const NekDouble> z, w;
......@@ -187,7 +187,7 @@ int main(int argc, char *argv[])
{
y[i] = func( z[i], nPts, pointsType );
}
// /////////////////////////////////////////////////////////
......@@ -198,12 +198,12 @@ int main(int argc, char *argv[])
int nNodes = 2*nPts; // Number of interpolating nodes
std::shared_ptr<Points<NekDouble> > nodes = PointsManager()[PointsKey(nNodes, pointsType)];
Array<OneD, const NekDouble> zNode = nodes->GetZ();
// Get the interpolation matrix I
// Note that I is 2N rows by N columns
const Points<NekDouble>::MatrixSharedPtrType Iptr = points->GetI(nNodes,zNode);
const NekMatrix<NekDouble> & I = *Iptr;
// Interpolate the data values in the y vector using the interpolation matrix I
Array<OneD, NekDouble> u(I.GetRows());
for(int i = 0; i < int(I.GetRows()); ++i)
......@@ -214,7 +214,7 @@ int main(int argc, char *argv[])
u[i] += I(i,j) * y[j];
}
}
// Display the original samples
cout << setprecision(3);
cout << "\nOriginal data: \nx = ";
......@@ -236,7 +236,7 @@ int main(int argc, char *argv[])
{
cout << setw(6) << u[i] << " ";
}
// Determine the exact solutions
cout << "\nexact = ";
for(int i = 0; i < nNodes; ++i)
......@@ -266,7 +266,7 @@ int main(int argc, char *argv[])
cout << setprecision(6);
cout << "\nLinf = " << setw(6) << Linf;
cout << "\nRMS = " << setw(6) << RMS << endl;
// Show the interpolation matrix
cout << "\nI = " << endl;
for(int i = 0; i < int(I.GetRows()); ++i)
......@@ -285,19 +285,19 @@ int main(int argc, char *argv[])
///////////////////////////////////////////////////////////
// Generate a list of projection nodes
Array<OneD, NekDouble> yNode(zNode.num_elements());
Array<OneD, NekDouble> yNode(zNode.size());
for(int i = 0; i < nNodes; ++i)
{
yNode[i] = func( zNode[i], nNodes, pointsType );
}
PointsKey key1(nNodes, pointsType);
// Note that I is 2N rows by N columns
const Points<NekDouble>::MatrixSharedPtrType GPptr = points->GetGalerkinProjection(key1);
const NekMatrix<NekDouble> & GP = *GPptr;
// Project the data values in the yNode vector using the projection matrix GP
for(int i = 0; i < int(GP.GetRows()); ++i)
{
......@@ -307,12 +307,12 @@ int main(int argc, char *argv[])
u[i] += GP(i,j) * yNode[j];
}
}
// Display the original samples
cout << setprecision(3);
cout << "\n\n\n **** Galerkin Project ****";
cout << "\n\nResults of Galkerin Project with " << kPointsTypeStr[pointsType] << ":";
cout << "\nOriginal data: \nx = ";
for(int i = 0; i < nPts; ++i)
{
......@@ -329,7 +329,7 @@ int main(int argc, char *argv[])
{
cout << setw(6) << u[i] << " ";
}
// Display the pointwise error
......@@ -354,7 +354,7 @@ int main(int argc, char *argv[])
cout << setprecision(6);
cout << "\nLinf = " << setw(6) << Linf;
cout << "\nRMS = " << setw(6) << RMS << endl;
// Show the projection matrix
cout << "\nI = " << endl;
for(int i = 0; i < int(GP.GetRows()); ++i)
......@@ -366,7 +366,7 @@ int main(int argc, char *argv[])
}
cout << "" << endl;
}
// /////////////////////////////////////////////////////////
// Derivation //
......@@ -377,8 +377,8 @@ int main(int argc, char *argv[])
// Note that, unlike I, D is N rows by N columns
const Points<NekDouble>::MatrixSharedPtrType Dptr = points->GetD();
const NekMatrix<NekDouble> & D = *Dptr;
// Differentiate the data values in the y vector using the derivative matrix D
Array<OneD, NekDouble> v(nPts);
for(int i = 0; i < int(D.GetRows()); ++i)
......@@ -390,7 +390,7 @@ int main(int argc, char *argv[])
}
}
// Display the derivative approximations
cout << "\n\n\n **** Differentiation ****" << setprecision(3);
cout << "\n\nResults of differentiation with " << kPointsTypeStr[pointsType] << ":";
......@@ -399,7 +399,7 @@ int main(int argc, char *argv[])
{
cout << setw(6) << v[i] << " ";
}
// Determine the exact solutions
cout << "\nexact = ";
for(int i = 0; i < nPts; ++i)
......@@ -431,7 +431,7 @@ int main(int argc, char *argv[])
cout << setprecision(6);
cout << "\nLinf = " << setw(6) << Linf;
cout << "\nRMS = " << setw(6) << RMS << endl;
// Show the derivation matrix
cout << "\nD = " << endl;
......@@ -446,13 +446,13 @@ int main(int argc, char *argv[])
}
// /////////////////////////////////////////////////////////
// Integration //
// /////////////////////////////////////////////////////////
// Integrate the function by approximating the integral as a weighted
// linear combination of function evaluations at the z[i] points: y[i] = f(z[i])
long double numericalIntegration = 0.0;
......@@ -461,7 +461,7 @@ int main(int argc, char *argv[])
numericalIntegration += w[i] * y[i] / integrandWeightFunction(z[i], pointsType);
}
// Display the integral approximation
cout << "\n\n\n **** Integration ****" << setprecision(6);
cout << "\n\nResults of integration with " << kPointsTypeStr[pointsType] << ":";
......
......@@ -156,7 +156,7 @@ int main(int argc, char *argv[])
}
ASSERTL1(util, "Unknown shape type!");
const int nPoints = r.num_elements();
const int nPoints = r.size();
const int dim = (shape == eTriangle || shape == eQuadrilateral) ? 2 : 3;
if (vm.count("integral"))
......
......@@ -185,8 +185,8 @@ void InputDat::ReadTecplotFEBlockZone(std::ifstream &datFile,
int nelmt = atoi(tag.substr(start + 2, end).c_str());
// set-up or extend m_pts array;
int norigpts = pts[0].num_elements();
int totfields = pts.num_elements();
int norigpts = pts[0].size();
int totfields = pts.size();
Array<OneD, Array<OneD, NekDouble> > origpts(totfields);
for (int i = 0; i < totfields; ++i)
{
......
......@@ -336,7 +336,7 @@ void OutputFileBase::ConvertExpToEquispaced(po::variables_map &vm)
for (int i = 0; i < numFields; ++i)
{
BndExpOld = expOld[i]->GetBndCondExpansions();
for (int j = 0; j < BndExpOld.num_elements(); ++j)
for (int j = 0; j < BndExpOld.size(); ++j)
{
BndExp = m_f->m_exp[i]->UpdateBndCondExpansion(j);
......@@ -361,10 +361,10 @@ void OutputFileBase::PrintErrorFromPts()
// We can just grab everything from points. This should be a
// reference, not a copy.
m_f->m_fieldPts->GetPts(fields);
for (int i = 0; i < fields.num_elements(); ++i)
for (int i = 0; i < fields.size(); ++i)
{
// calculate L2 and Linf value
int npts = fields[i].num_elements();
int npts = fields[i].size();
NekDouble l2err = 0.0;
NekDouble linferr = 0.0;
......
......@@ -129,7 +129,7 @@ void OutputVtk::OutputFromPts(po::variables_map &vm)
int numBlocks = 0;
for (i = 0; i < ptsConn.size(); ++i)
{
numBlocks += ptsConn[i].num_elements() / nvert;
numBlocks += ptsConn[i].size() / nvert;
}
// write out pieces of data.
......@@ -164,7 +164,7 @@ void OutputVtk::OutputFromPts(po::variables_map &vm)
int cnt = 1;
for (i = 0; i < ptsConn.size(); ++i)
{
for (j = 0; j < ptsConn[i].num_elements(); ++j)
for (j = 0; j < ptsConn[i].size(); ++j)
{
outfile << ptsConn[i][j] << " ";
if ((!(cnt % nvert)) && cnt)
......
......@@ -72,7 +72,7 @@ void ProcessDeform::Process(po::variables_map &vm)
Array<OneD, MultiRegions::ExpListSharedPtr> exp(m_f->m_exp.size());
for (int i = 0; i < exp.num_elements(); ++i)
for (int i = 0; i < exp.size(); ++i)
{
exp[i] = m_f->m_exp[i];
}
......
......@@ -316,8 +316,8 @@ void ProcessEquiSpacedOutput::Process(po::variables_map &vm)
e->GetSimplexEquiSpacedConnectivity(conn);
}
}
Array<OneD, int> newconn(conn.num_elements());
for (int j = 0; j < conn.num_elements(); ++