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

Remove integral functions, replace with (presumably more efficient) version...

Remove integral functions, replace with (presumably more efficient) version using standard quadrature metric.
parent bc3238a0
......@@ -576,12 +576,14 @@ int main(int argc, char *argv[]){
E->BwdTrans(coeffs,phys);
//-------------------------------------------
#if 0
int nPts = E->GetTotPoints();
Array<OneD, NekDouble> errArr(nPts);
Vmath::Vsub(nPts, phys, 1, sol, 1, errArr, 1);
printSolution(E, "out.dat", x, y, z, errArr);
printSolution(E, "sol.dat", x, y, z, phys);
#endif
//--------------------------------------------
// Calculate L_inf error
cout << "L infinity error: " << E->Linf(phys,sol) << endl;
......@@ -629,8 +631,6 @@ NekDouble Tet_sol(NekDouble x, NekDouble y, NekDouble z,
int l,k,m;
NekDouble sol = 0.0;
return cos(x)*cos(y)*cos(z);
for(k = 0; k < order1; ++k)
{
for(l = 0; l < order2-k; ++l)
......@@ -665,7 +665,6 @@ NekDouble Prism_sol(NekDouble x, NekDouble y, NekDouble z,
int l,k,m;
NekDouble sol = 0;
return cos(x)*cos(y)*cos(z);
for(k = 0; k < order1; ++k)
{
for(l = 0; l < order2; ++l)
......
......@@ -308,6 +308,16 @@ namespace Nektar
return x->second;
}
NekDouble StdExpansion3D::v_Integral(
const Array<OneD, const NekDouble>& inarray)
{
const int nqtot = GetTotPoints();
Array<OneD, NekDouble> tmp(GetTotPoints());
MultiplyByStdQuadratureMetric(inarray, tmp);
return Vmath::Vsum(nqtot, tmp, 1);
}
void StdExpansion3D::v_NegateFaceNormal(const int face)
{
m_negatedNormals[face] = true;
......
......@@ -182,6 +182,8 @@ namespace Nektar
Array<OneD,NekDouble> &outarray,
const StdRegions::StdMatrixKey &mkey);
STD_REGIONS_EXPORT virtual NekDouble v_Integral(
const Array<OneD, const NekDouble>& inarray);
STD_REGIONS_EXPORT virtual void v_NegateFaceNormal(
const int face);
......@@ -208,88 +210,3 @@ namespace Nektar
} //end of namespace
#endif //STDEXP3D_H
/**
* $Log: StdExpansion3D.h,v $
* Revision 1.19 2008/09/16 13:37:03 pvos
* Restructured the LocalToGlobalMap classes
*
* Revision 1.18 2008/07/04 10:18:40 pvos
* Some updates
*
* Revision 1.17 2008/06/06 23:21:13 ehan
* Added doxygen documentation
*
* Revision 1.16 2008/05/30 00:33:49 delisi
* Renamed StdRegions::ShapeType to StdRegions::ExpansionType.
*
* Revision 1.15 2008/05/15 22:39:34 ehan
* Clean up the codes
*
* Revision 1.14 2008/04/06 06:04:15 bnelson
* Changed ConstArray to Array<const>
*
* Revision 1.13 2008/02/12 02:46:33 jfrazier
* Moved typedef to the top of the file.
*
* Revision 1.12 2008/02/12 01:30:07 ehan
* Added typedef StdExpansion3DSharedPtr.
*
* Revision 1.11 2007/11/08 14:27:53 ehan
* Fixed PhysTensorDerivative3D matrix and improved L1 error up to 1e-15.
*
* Revision 1.10 2007/10/29 20:31:04 ehan
* Fixed floating point approximation up to 1-e15 for PhysEvaluate.
*
* Revision 1.9 2007/07/20 02:16:54 bnelson
* Replaced boost::shared_ptr with Nektar::ptr
*
* Revision 1.8 2007/05/15 05:18:23 bnelson
* Updated to use the new Array object.
*
* Revision 1.7 2007/04/10 14:00:45 sherwin
* Update to include SharedArray in all 2D element (including Nodal tris). Have also remvoed all new and double from 2D shapes in StdRegions
*
* Revision 1.6 2007/03/29 19:35:09 bnelson
* Replaced boost::shared_array with SharedArray
*
* Revision 1.5 2007/03/20 16:58:43 sherwin
* Update to use Array<OneD, NekDouble> storage and NekDouble usage, compiling and executing up to Demos/StdRegions/Project1D
*
* Revision 1.4 2007/01/17 16:05:40 pvos
* updated doxygen documentation
*
* Revision 1.3 2006/07/02 17:16:18 sherwin
*
* Modifications to make MultiRegions work for a connected domain in 2D (Tris)
*
* Revision 1.2 2006/06/13 18:05:02 sherwin
* Modifications to make MultiRegions demo ProjectLoc2D execute properly.
*
* Revision 1.1 2006/05/04 18:58:31 kirby
* *** empty log message ***
*
* Revision 1.12 2006/04/25 20:23:33 jfrazier
* Various fixes to correct bugs, calls to ASSERT, etc.
*
* Revision 1.11 2006/03/12 14:20:44 sherwin
*
* First compiling version of SpatialDomains and associated modifications
*
* Revision 1.10 2006/03/06 17:12:45 sherwin
*
* Updated to properly execute all current StdRegions Demos.
*
* Revision 1.9 2006/03/04 20:26:54 bnelson
* Added comments after #endif.
*
* Revision 1.8 2006/03/01 08:25:03 sherwin
*
* First compiling version of StdRegions
*
* Revision 1.7 2006/02/27 23:47:23 sherwin
*
* Standard coding update upto compilation of StdHexExp.cpp
*
**/
......@@ -80,194 +80,6 @@ namespace Nektar
{
}
//////////////////////////////
// Integration Methods
//////////////////////////////
/**
* @param fx ?
* @param gy ?
* @param hz ?
* @param inarray ?
* @param outarray ?
*/
void StdHexExp::TripleTensorProduct(
const Array<OneD, const NekDouble>& fx,
const Array<OneD, const NekDouble>& gy,
const Array<OneD, const NekDouble>& hz,
const Array<OneD, const NekDouble>& inarray,
Array<OneD, NekDouble> & outarray )
{
// Using matrix operation, not sum-factorization.
// Regarding the 3D array, inarray[k][j][i], x is changing the
// fastest and z the slowest. Thus, the first x-vector of points
// refers to the first row of the first stack. The first y-vector
// refers to the first column of the first stack. The first z-
// vector refers to the vector of stacks intersecting the first row
// and first column. So in C++, i refers to column, j to row, and k
// to stack. Contrasting this with the usual C++ matrix convention,
// note that i does not refer to a C++ row, nor j to C++ column.
int nx = fx.num_elements();
int ny = gy.num_elements();
int nz = hz.num_elements();
// Multiply by integration constants...
// Hadamard multiplication refers to elementwise multiplication of
// two vectors.
// Hadamard each row with the first vector (x-vector); the index i
// is changing the fastest.
// For each j and k, iterate over each row in all of the stacks at
// once.
for (int jk = 0; jk < ny*nz; ++jk)
{
Vmath::Vmul(
nx, // Size of first weight vector
&inarray[0] + jk*nx, 1, // Offset and stride of each row-
// vector (x is changing fastest)
fx.get(), 1, // First weight vector (with stride
// of 1)
&outarray[0] + jk*nx, 1 // Output has same offset and
// stride as input
);
}
// Hadamard each column with the second vector (y-vector)
// For each stack in the 3D-array, do the following...
for (int k = 0; k < nz; ++k)
{
// Iterate over each column in the current stack
for (int i = 0; i < nx; ++i)
{
Vmath::Vmul(
ny, // Size of second weight vector
&outarray[0] + i + nx*ny*k, nx, // Offset and
// stride of each column-vector
gy.get(), 1, // second weight vector (with
// stride of 1)
&outarray[0] + i + nx*ny*k, nx // Output has same
// offset and stride as input
);
}
}
// Hadamard each stack-vector with the third vector (z-vector)
// Iterate over each element in the topmost stack
for (int ij = 0; ij < nx*ny; ++ij)
{
Vmath::Vmul(
nz, // Size of third weight vector
&outarray[0] + ij, nx*ny, // Offset and stride of each
// stack-vector
hz.get(), 1, // Third weight vector (with
// stride of 1)
&outarray[0] + ij, nx*ny // Output has same offset and
// stride as input
);
}
}
/**
* Inner-Product with respect to the weights: i.e., this is the triple
* sum of the product of the four inputs over the Hexahedron
* x-dimension is the row, it is the index that changes the fastest
* y-dimension is the column
* z-dimension is the stack, it is the index that changes the slowest
*/
NekDouble StdHexExp::TripleInnerProduct(
const Array<OneD, const NekDouble>& fxyz,
const Array<OneD, const NekDouble>& wx,
const Array<OneD, const NekDouble>& wy,
const Array<OneD, const NekDouble>& wz
)
{
int Qx = wx.num_elements();
int Qy = wy.num_elements();
int Qz = wz.num_elements();
if( fxyz.num_elements() != Qx*Qy*Qz ) {
cerr << "TripleTetrahedralInnerProduct expected "
<< fxyz.num_elements()
<< " quadrature points from the discretized input "
"function but got "
<< Qx*Qy*Qz << " instead." << endl;
}
// Sum-factorizing over the stacks
Array<OneD, NekDouble> A(Qx*Qy, 0.0);
for( int i = 0; i < Qx; ++i ) {
for( int j = 0; j < Qy; ++j ) {
for( int k = 0; k < Qz; ++k ) {
A[i + Qx*j] += fxyz[i + Qx*(j + Qy*k)] * wz[k];
}
}
}
// Sum-factorizing over the columns
Array<OneD, NekDouble> b(Qx, 0.0);
for( int i = 0; i < Qx; ++i ) {
for( int j = 0; j < Qy; ++j ) {
b[i] += A[i + Qx*j] * wy[j];
}
}
// Sum-factorizing over the rows
NekDouble c = 0;
for( int i = 0; i < Qx; ++i ) {
c += b[i] * wx[i];
}
return c;
}
/**
* @param inarray ?
* @param wx ?
* @param wy ?
* @param wz ?
*/
NekDouble StdHexExp::Integral3D(
const Array<OneD, const NekDouble>& inarray,
const Array<OneD, const NekDouble>& wx,
const Array<OneD, const NekDouble>& wy,
const Array<OneD, const NekDouble>& wz)
{
return TripleInnerProduct( inarray, wx, wy, wz );
}
/**
* @param inarray Definition of function to be returned at
* quadrature point of expansion.
* @returns
* \f$\int^1_{-1}\int^1_{-1}\int^1_{-1} u(\xi_1, \xi_2, \xi_3)
* J[i,j,k] d \xi_1 d \xi_2 d \xi_3 \f$ \n
* \f$ = \sum_{i=0}^{Q_1 - 1} \sum_{j=0}^{Q_2 - 1}
* \sum_{k=0}^{Q_3 - 1} u(\xi_{1i}, \xi_{2j},\xi_{3k})
* w_{i} w_{j} w_{k} \f$ \n
* where \f$inarray[i,j, k] = u(\xi_{1i},\xi_{2j}, \xi_{3k}) \f$ \n
* and \f$ J[i,j,k] \f$ is the Jacobian evaluated at the quadrature
* point.
*/
NekDouble StdHexExp::v_Integral(
const Array<OneD, const NekDouble>& inarray)
{
Array<OneD, const NekDouble> w0, w1, w2;
w0 = m_base[0]->GetW();
w1 = m_base[1]->GetW();
w2 = m_base[2]->GetW();
return Integral3D(inarray, w0, w1, w2);
}
bool StdHexExp::v_IsBoundaryInteriorExpansion()
{
return (m_base[0]->GetBasisType() == LibUtilities::eModified_A) &&
......
......@@ -65,28 +65,6 @@ namespace Nektar
protected:
//-------------------------------
// Integration Methods
//-------------------------------
STD_REGIONS_EXPORT void TripleTensorProduct(
const Array<OneD, const NekDouble>& fx,
const Array<OneD, const NekDouble>& gy,
const Array<OneD, const NekDouble>& hz,
const Array<OneD, const NekDouble>& inarray,
Array<OneD, NekDouble>& outarray );
STD_REGIONS_EXPORT NekDouble TripleInnerProduct(
const Array<OneD, const NekDouble>& fxyz,
const Array<OneD, const NekDouble>& wx,
const Array<OneD, const NekDouble>& wy,
const Array<OneD, const NekDouble>& wz);
STD_REGIONS_EXPORT NekDouble Integral3D(
const Array<OneD, const NekDouble>& inarray,
const Array<OneD, const NekDouble>& wx,
const Array<OneD, const NekDouble>& wy,
const Array<OneD, const NekDouble>& wz);
STD_REGIONS_EXPORT virtual NekDouble v_Integral(
const Array<OneD, const NekDouble>& inarray);
//----------------------------
// Differentiation Methods
//----------------------------
......@@ -227,25 +205,28 @@ namespace Nektar
//--------------------------
STD_REGIONS_EXPORT virtual void v_GetFaceToElementMap(
const int fid,
const Orientation faceOrient,
const Orientation faceOrient,
Array<OneD, unsigned int> &maparray,
Array<OneD, int> &signarray,
int nummodesA = -1,
int nummodesB = -1);
STD_REGIONS_EXPORT virtual int v_GetVertexMap(int localVertexId,
bool useCoeffPacking = false); STD_REGIONS_EXPORT virtual void v_GetEdgeInteriorMap(const int eid,
const Orientation edgeOrient,
STD_REGIONS_EXPORT virtual int v_GetVertexMap(
int localVertexId,
bool useCoeffPacking = false);
STD_REGIONS_EXPORT virtual void v_GetEdgeInteriorMap(
const int eid,
const Orientation edgeOrient,
Array<OneD, unsigned int> &maparray,
Array<OneD, int> &signarray);
Array<OneD, int> &signarray);
STD_REGIONS_EXPORT virtual void v_GetFaceInteriorMap(
const int fid,
const Orientation faceOrient,
const int fid,
const Orientation faceOrient,
Array<OneD, unsigned int> &maparray,
Array<OneD, int>& signarray);
STD_REGIONS_EXPORT virtual void v_GetInteriorMap(
Array<OneD, unsigned int>& outarray);
Array<OneD, unsigned int> &outarray);
STD_REGIONS_EXPORT virtual void v_GetBoundaryMap(
Array<OneD, unsigned int>& outarray);
Array<OneD, unsigned int> &outarray);
//---------------------------------------
......
......@@ -74,201 +74,6 @@ namespace Nektar
{
}
//---------------------------------------
// Integration Methods
//---------------------------------------
void StdPrismExp::TripleTensorProduct(const Array<OneD, const NekDouble>& fx,
const Array<OneD, const NekDouble>& gy,
const Array<OneD, const NekDouble>& hz,
const Array<OneD, const NekDouble>& inarray,
Array<OneD, NekDouble>& outarray)
{
// Using matrix operation, not sum-factorization. Regarding the
// 3D array, inarray[k][j][i], x is changing the fastest and z the
// slowest. Thus, the first x-vector of points refers to the
// first row of the first stack. The first y-vector refers to the
// first column of the first stack. The first z-vector refers to
// the vector of stacks intersecting the first row and first
// column. So in C++, i refers to column, j to row, and k to
// stack. Contrasting this with the usual C++ matrix convention,
// note that i does not refer to a C++ row, nor j to C++ column.
int nx = fx.num_elements();
int ny = gy.num_elements();
int nz = hz.num_elements();
// Multiply by integration constants... Hadamard multiplication
// refers to elementwise multiplication of two vectors. Hadamard
// each row with the first vector (x-vector); the index i is
// changing the fastest.
for (int jk = 0; jk < ny*nz; ++jk) // For each j and k, iterate over each row in all of the stacks at once
{
Vmath::Vmul(
nx, // Size of first weight vector
&inarray[0] + jk*nx, 1, // Offset and stride of each row-vector (x is changing fastest)
fx.get(), 1, // First weight vector (with stride of 1)
&outarray[0] + jk*nx, 1 // Output has same offset and stride as input
);
}
// Hadamard each column with the second vector (y-vector)
for (int k = 0; k < nz; ++k) // For each stack in the 3D-array, do the following...
{
for (int i = 0; i < nx; ++i) // Iterate over each column in the current stack
{
Vmath::Vmul(
ny, // Size of second weight vector
&outarray[0] + i + nx*ny*k, nx, // Offset and stride of each column-vector
gy.get(), 1, // second weight vector (with stride of 1)
&outarray[0] + i + nx*ny*k, nx // Output has same offset and stride as input
);
}
}
// Hadamard each stack-vector with the third vector (z-vector)
for (int ij = 0; ij < nx*ny; ++ij) // Iterate over each element in the topmost stack
{
Vmath::Vmul(
nz, // Size of third weight vector
&outarray[0] + ij, nx*ny, // Offset and stride of each stack-vector
hz.get(), 1, // Third weight vector (with stride of 1)
&outarray[0] + ij, nx*ny // Output has same offset and stride as input
);
}
}
// Inner-Product with respect to the weights: i.e., this is the triple
// sum of the product of the four inputs over the prism. x-dimension
// is the row, it is the index that changes the fastest; y-dimension
// is the column; z-dimension is the stack, it is the index that
// changes the slowest.
NekDouble StdPrismExp::TripleInnerProduct(const Array<OneD, const NekDouble>& fxyz,
const Array<OneD, const NekDouble>& wx,
const Array<OneD, const NekDouble>& wy,
const Array<OneD, const NekDouble>& wz)
{
int Qx = wx.num_elements();
int Qy = wy.num_elements();
int Qz = wz.num_elements();
if (fxyz.num_elements() != Qx*Qy*Qz)
{
cerr << "TripleInnerProduct expected " << fxyz.num_elements()
<< " quadrature points from the discretized input function but got "
<< Qx*Qy*Qz << " instead." << endl;
}
// Sum-factorizing over the stacks
Array<OneD, NekDouble> A(Qx*Qy, 0.0);
for (int i = 0; i < Qx; ++i)
{
for (int j = 0; j < Qy; ++j)
{
for (int k = 0; k < Qz; ++k)
{
A[i + Qx*j] += fxyz[i + Qx*(j + Qy*k)] * wz[k];
}
}
}
// Sum-factorizing over the columns
Array<OneD, NekDouble> b(Qx, 0.0);
for (int i = 0; i < Qx; ++i)
{
for (int j = 0; j < Qy; ++j)
{
b[i] += A[i + Qx*j] * wy[j];
}
}
// Sum-factorizing over the rows
NekDouble c = 0;
for (int i = 0; i < Qx; ++i)
{
c += b[i] * wx[i];
}
return c;
}
NekDouble StdPrismExp::Integral3D(
const Array<OneD, const NekDouble>& inarray,
const Array<OneD, const NekDouble>& wx,
const Array<OneD, const NekDouble>& wy,
const Array<OneD, const NekDouble>& wz)
{
return TripleInnerProduct(inarray, wx, wy, wz);
}
/**
* \brief Integrate the physical point list \a inarray over prismatic
* region and return the value.
*
* Inputs:\n
*
* - \a inarray: definition of function to be returned at quadrature
* point of expansion.
*
* Outputs:\n
*
* - returns \f$\int^1_{-1}\int^1_{-1}\int^1_{-1} u(\bar \eta_1,
* \xi_2, \xi_3) J[i,j,k] d \bar \eta_1 d \xi_2 d \xi_3 \f$ \n \f$ =
* \sum_{i=0}^{Q_1 - 1} \sum_{j=0}^{Q_2 - 1} \sum_{k=0}^{Q_3 - 1}
* u(\bar \eta_{1i}^{0,0}, \xi_{2j}^{0,0},\xi_{3k}^{1,0})w_{i}^{0,0}
* w_{j}^{0,0} \hat w_{k}^{1,0} \f$ \n where \f$ inarray[i,j, k] =
* u(\bar \eta_{1i}^{0,0}, \xi_{2j}^{0,0},\xi_{3k}^{1,0}) \f$, \n
* \f$\hat w_{i}^{1,0} = \frac {w_{j}^{1,0}} {2} \f$ \n and \f$
* J[i,j,k] \f$ is the Jacobian evaluated at the quadrature point.
*/
NekDouble StdPrismExp::v_Integral(
const Array<OneD, const NekDouble>& inarray)
{
// Using implementation from page 146 of Spencer Sherwin's book.
int Qz = m_base[2]->GetNumPoints();
// Get the point distributions:
// * x is assumed to be Gauss-Lobatto-Legendre (incl. -1 and 1)
// * y is assumed to be Gauss-Lobatto-Legendre (incl. -1 and 1)
Array<OneD, const NekDouble> z, wx, wy, wz;
wx = m_base[0]->GetW();
wy = m_base[1]->GetW();
m_base[2]->GetZW(z,wz);
Array<OneD, NekDouble> wz_hat = Array<OneD, NekDouble>(Qz, 0.0);
// Convert wz into wz_hat, which includes the 1/2 scale factor.
// Nothing else need be done if the point distribution is Jacobi
// (1,0) since (1 - xi_z) is already factored into the weights.
// Note by coincidence, xi_y = eta_y, xi_z = eta_z (xi_z = z
// according to our notation)
switch(m_base[2]->GetPointsType())
{
// Common case
case LibUtilities::eGaussRadauMAlpha1Beta0: // (1,0) Jacobi Inner product
Vmath::Smul(Qz, 0.5, (NekDouble *)wz.get(), 1, wz_hat.get(), 1);
break;
// Corner cases
case LibUtilities::eGaussLobattoLegendre:
case LibUtilities::eGaussRadauMLegendre:
for (int k = 0; k < Qz; ++k)
{
wz_hat[k] = 0.5*(1.0 - z[k]) * wz[k];
}
break;
default:
ASSERTL0(false, "Unsupported quadrature points type.");
break;
}