Commit 9060fa9f authored by Spencer Sherwin's avatar Spencer Sherwin

Updates relating to static condensation implementation


git-svn-id: https://gforge.sci.utah.edu/svn/nektar/trunk@932 305cdda6-5ce1-45b3-a98d-dfc68c8b3305
parent 1e12abb3
......@@ -72,7 +72,7 @@ MACRO(SET_COMMON_PROPERTIES name)
IF ( NEKTAR_FULL_DEBUG )
SET(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} -DNEKTAR_DEBUG -DNEKTAR_FULLDEBUG")
ELSE( NEKTAR_FULL_DEBUG )
SET(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} -DNEKTAR_DEBUG")
SET(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} -DNEKTAR_DEBUG -fpermissive")
ENDIF( NEKTAR_FULL_DEBUG )
SET(CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS_RELEASE} -DNEKTAR_RELEASE")
......
......@@ -143,7 +143,7 @@ namespace Nektar
{
value = true;
}
return value;
}
......@@ -191,7 +191,6 @@ namespace Nektar
CreateFuncType m_globalCreateFunc;
CreateFuncContainer m_keySpecificCreateFuncs;
};
template <typename KeyType, typename ValueT, typename opLessCreator> typename NekManager<KeyType, ValueT, opLessCreator>::ValueContainerPool NekManager<KeyType, ValueT, opLessCreator>::m_ValueContainerPool;
}
}
......
......@@ -352,7 +352,7 @@ namespace Nektar
static ConstArray<OneD, DataType> CreateWithOffset(const ConstArray<OneD, DataType>& rhs, unsigned int offset)
{
ConstArray<OneD, DataType> result(rhs);
result.m_offset = offset;
result.m_offset += offset;
result.m_size = rhs.m_size - offset;
return result;
}
......@@ -573,7 +573,7 @@ namespace Nektar
static Array<OneD, DataType> CreateWithOffset(const Array<OneD, DataType>& rhs, unsigned int offset)
{
Array<OneD, DataType> result(rhs);
result.m_offset = offset;
result.m_offset += offset;
result.m_size = rhs.m_size - offset;
return result;
}
......
......@@ -260,7 +260,7 @@ namespace Nektar
/** \brief return basis definition array m_bdata */
inline const ConstArray<OneD, NekDouble>& GetBdata() const
{
return m_bdata;
return m_bdata;
}
/** \brief return basis definition array m_dbdata */
......
......@@ -23,6 +23,17 @@
#include <errno.h>
#include <stdlib.h>
// Microsoft fix for math functions
#ifndef _MSC_VER
#define _hypot hypot
#define _jn jn
#define _j0 j0
#define _j1 j1
#define _yn yn
#define _y0 y0
#define _y1 y1
#endif
namespace Nektar
{
namespace LibUtilities
......
......@@ -251,6 +251,8 @@ namespace Nektar
}
unsigned int GetRows() const { return A.GetRows(); }
unsigned int GetColumns() const { return A.GetColumns(); }
private:
void FactorMatrix(NekMatrix<DataType,DiagonalMatrixTag> &theA)
{
......
......@@ -58,6 +58,10 @@ namespace Nektar
typedef LinearSystem<DNekScalMat> DNekScalLinSys;
typedef boost::shared_ptr<DNekScalLinSys> DNekScalLinSysSharedPtr;
typedef NekMatrix<DNekMat, FullMatrixTag, BlockMatrixTag> DNekBlkMat;
typedef boost::shared_ptr<DNekBlkMat> DNekBlkMatSharedPtr;
typedef NekMatrix<DNekScalMat, FullMatrixTag, BlockMatrixTag> DNekScalBlkMat;
typedef boost::shared_ptr<DNekScalBlkMat> DNekScalBlkMatSharedPtr;
......@@ -67,6 +71,9 @@ namespace Nektar
/**
$Log: NekTypeDefs.hpp,v $
Revision 1.8 2007/07/26 08:40:49 sherwin
Update to use generalised i/o hooks in Helmholtz1D
Revision 1.7 2007/07/22 23:03:28 bnelson
Backed out Nektar::ptr.
......
......@@ -40,7 +40,11 @@ namespace Nektar
{
/// \brief Specifies if the pointer passed to a NekMatrix or NekVector is copied into an internal
/// representation or used directly.
enum PointerWrapper { eWrapper, eCopy };
enum PointerWrapper
{
eWrapper,
eCopy
};
}
#endif //NEKTAR_LIB_UTILITIES_LINEAR_ALGEBRA_POINTER_WRAPPER_H
......
......@@ -45,30 +45,30 @@ namespace Nektar
SegExp::SegExp(const LibUtilities::BasisKey &Ba,
const SpatialDomains::SegGeomSharedPtr &geom):
StdRegions::StdSegExp(Ba),
m_matrixManager(std::string("StdExp"))
m_matrixManager(std::string("StdExp")),
m_staticCondMatrixManager(std::string("StdExpStdCondMat"))
{
m_geom = geom;
for(int i = 0; i < StdRegions::SIZE_MatrixType; ++i)
{
m_matrixManager.RegisterCreator(MatrixKey((StdRegions::MatrixType) i,
StdRegions::eNoShapeType,*this),
boost::bind(&SegExp::CreateMatrix, this, _1));
m_matrixManager.RegisterCreator(MatrixKey((StdRegions::MatrixType) i, StdRegions::eNoShapeType,*this), boost::bind(&SegExp::CreateMatrix, this, _1));
m_staticCondMatrixManager.RegisterCreator(MatrixKey((StdRegions::MatrixType) i, StdRegions::eNoShapeType,*this), boost::bind(&SegExp::CreateStaticCondMatrix, this, _1));
}
GenMetricInfo();
}
SegExp::SegExp(const LibUtilities::BasisKey &Ba):
StdRegions::StdSegExp(Ba),
m_matrixManager(std::string("StdExp"))
StdRegions::StdSegExp(Ba),
m_matrixManager(std::string("StdExp")),
m_staticCondMatrixManager(std::string("StdExpStdCondMat"))
{
for(int i = 0; i < StdRegions::SIZE_MatrixType; ++i)
{
m_matrixManager.RegisterCreator(MatrixKey((StdRegions::MatrixType) i,
StdRegions::eNoShapeType,*this),
boost::bind(&SegExp::CreateMatrix, this, _1));
m_matrixManager.RegisterCreator(MatrixKey((StdRegions::MatrixType) i,StdRegions::eNoShapeType,*this), boost::bind(&SegExp::CreateMatrix, this, _1));
m_staticCondMatrixManager.RegisterCreator(MatrixKey((StdRegions::MatrixType) i,StdRegions::eNoShapeType,*this), boost::bind(&SegExp::CreateStaticCondMatrix, this, _1));
}
// Set up unit geometric factors.
......@@ -82,8 +82,9 @@ namespace Nektar
// copy constructor
SegExp::SegExp(const SegExp &S):
StdRegions::StdSegExp(S),
m_matrixManager(std::string("StdExp"))
StdRegions::StdSegExp(S),
m_matrixManager(std::string("StdExp")),
m_staticCondMatrixManager(std::string("StdExpStdCondMat"))
{
m_geom = S.m_geom;
m_metricinfo = S.m_metricinfo;
......@@ -142,7 +143,7 @@ namespace Nektar
// interpolate Jacobian
ndata = Array<OneD,NekDouble>(nq);
odata = Xgfac->GetJac();
Interp1D(CBasis0->GetBasisKey(),odata,
m_base[0]->GetBasisKey(), ndata);
......@@ -418,10 +419,10 @@ namespace Nektar
{
IProductWRTBase(inarray,outarray);
// get Mass matrix inverse
MatrixKey masskey(StdRegions::eInvMass, DetShapeType(),*this);
DNekScalMatSharedPtr matsys = m_matrixManager[masskey];
// get Mass matrix inverse
MatrixKey masskey(StdRegions::eInvMass, DetShapeType(),*this);
DNekScalMatSharedPtr matsys = m_matrixManager[masskey];
// copy inarray in case inarray == outarray
DNekVec in (m_ncoeffs,outarray);
DNekVec out(m_ncoeffs,outarray,eWrapper);
......@@ -621,6 +622,26 @@ namespace Nektar
}
}
DNekBlkMatSharedPtr SegExp::GetStdStaticCondMatrix(const StdRegions::StdMatrixKey &mkey)
{
// Need to check if matrix exists in stdStaticCondMatrixManager.
// If not then make a local expansion with standard metric info
// and generate matrix. Otherwise direct call is OK.
if(!StdStaticCondMatManagerAlreadyCreated(mkey))
{
LibUtilities::BasisKey bkey = m_base[0]->GetBasisKey();
SegExpSharedPtr tmp = MemoryManager<SegExp>::AllocateSharedPtr(bkey);
return tmp->StdSegExp::GetStdStaticCondMatrix(mkey);
}
else
{
return StdSegExp::GetStdStaticCondMatrix(mkey);
}
}
NekDouble SegExp::PhysEvaluate(const ConstArray<OneD,NekDouble>& coord)
{
Array<OneD,NekDouble> Lcoord = Array<OneD,NekDouble>(1);
......@@ -694,33 +715,23 @@ namespace Nektar
break;
case StdRegions::eHelmholtz:
{
NekDouble factor = mkey.GetScaleFactor();
MatrixKey masskey(StdRegions::eMass,
mkey.GetShapeType(), *this);
DNekScalMatSharedPtr mass = this->m_matrixManager[masskey];
#if 1
DNekScalMat &MassMat = *(this->m_matrixManager[masskey]);
MatrixKey lapkey(StdRegions::eLaplacian,
mkey.GetShapeType(), *this);
DNekScalMatSharedPtr lap = this->m_matrixManager[lapkey];
#else
MatrixKey lapkey(StdRegions::eLaplacian,
mkey.GetShapeType(), *this);
DNekMatSharedPtr mat = GetStdMatrix(*lapkey.GetStdMatKey());
NekDouble gfac = m_metricinfo->GetGmat()[0][0];
NekDouble jac = m_metricinfo->GetJac()[0];
NekDouble fac = gfac*gfac*jac;
DNekScalMatSharedPtr lap = MemoryManager<DNekScalMat>::AllocateSharedPtr(fac,mat);
#endif
int rows = lap->GetRows();
int cols = lap->GetColumns();
DNekScalMat &LapMat = *(this->m_matrixManager[lapkey]);
int rows = LapMat.GetRows();
int cols = LapMat.GetColumns();
DNekMatSharedPtr helm = MemoryManager<DNekMat>::AllocateSharedPtr(rows,cols);
(*helm) = (*lap) + 1.0/factor*(*mass);
// Even better: helm = lap + 1.0/factor*mass;
NekDouble one = 1.0;
(*helm) = factor*LapMat + MassMat;
returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(factor,helm);
returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,helm);
}
break;
default:
......@@ -731,11 +742,122 @@ namespace Nektar
return returnval;
}
DNekScalBlkMatSharedPtr SegExp::CreateStaticCondMatrix(const MatrixKey &mkey)
{
DNekScalBlkMatSharedPtr returnval;
ASSERTL2(m_metricinfo->GetGtype == SpatialDomains::eNoGeomType,"Geometric information is not set up");
// set up block matrix system
int nbdry = NumBndryCoeffs();
int nint = m_ncoeffs - nbdry;
unsigned int exp_size[] = {nbdry,nint};
int nblks = 2;
returnval = MemoryManager<DNekScalBlkMat>::AllocateSharedPtr(nblks,nblks,exp_size,exp_size); //Really need a constructor which takes Arrays
NekDouble factor = 1.0;
switch(mkey.GetMatrixType())
{
case StdRegions::eHelmholtz: // special case since Helmholtz not defined in StdRegions
// use Deformed case for both regular and deformed geometries
factor = 1.0;
goto UseLocRegionsMatrix;
break;
default:
if(m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
{
factor = 1.0;
goto UseLocRegionsMatrix;
}
else
{
DNekScalMatSharedPtr mat = GetLocMatrix(mkey);
factor = mat->Scale();
goto UseStdRegionsMatrix;
}
break;
UseStdRegionsMatrix:
{
NekDouble invfactor = 1.0/factor;
NekDouble one = 1.0;
DNekBlkMatSharedPtr mat = GetStdStaticCondMatrix(*(mkey.GetStdMatKey()));
DNekScalMatSharedPtr Atmp;
DNekMatSharedPtr Asubmat;
returnval->SetBlock(0,0,Atmp = MemoryManager<DNekScalMat>::AllocateSharedPtr(factor,Asubmat = mat->GetBlock(0,0)));
returnval->SetBlock(0,1,Atmp = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,Asubmat = mat->GetBlock(0,1)));
returnval->SetBlock(1,0,Atmp = MemoryManager<DNekScalMat>::AllocateSharedPtr(factor,Asubmat = mat->GetBlock(1,0)));
returnval->SetBlock(1,1,Atmp = MemoryManager<DNekScalMat>::AllocateSharedPtr(invfactor,Asubmat = mat->GetBlock(1,1)));
}
break;
UseLocRegionsMatrix:
{
int i,j;
NekDouble invfactor = 1.0/factor;
NekDouble one = 1.0;
DNekScalMat &mat = *GetLocMatrix(mkey);
DNekMatSharedPtr A = MemoryManager<DNekMat>::AllocateSharedPtr(nbdry,nbdry);
DNekMatSharedPtr B = MemoryManager<DNekMat>::AllocateSharedPtr(nbdry,nint);
DNekMatSharedPtr C = MemoryManager<DNekMat>::AllocateSharedPtr(nint,nbdry);
DNekMatSharedPtr D = MemoryManager<DNekMat>::AllocateSharedPtr(nint,nint);
for(i = 0; i < nbdry; ++i)
{
for(j = 0; j < nbdry; ++j)
{
(*A)(i,j) = mat(i,j);
}
for(j = 0; j < nint; ++j)
{
(*B)(i,j) = mat(i,nbdry+j);
}
}
for(i = 0; i < nint; ++i)
{
for(j = 0; j < nbdry; ++j)
{
(*C)(i,j) = mat(nbdry+i,j);
}
for(j = 0; j < nint; ++j)
{
(*D)(i,j) = mat(nbdry+i,nbdry+j);
}
}
// Calculate static condensed system
if(nint)
{
D->Invert();
(*B) = (*B)*(*D);
(*A) = (*A) - (*B)*(*C);
}
DNekScalMatSharedPtr Atmp;
returnval->SetBlock(0,0,Atmp = MemoryManager<DNekScalMat>::AllocateSharedPtr(factor,A));
returnval->SetBlock(0,1,Atmp = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,B));
returnval->SetBlock(1,0,Atmp = MemoryManager<DNekScalMat>::AllocateSharedPtr(factor,C));
returnval->SetBlock(1,1,Atmp = MemoryManager<DNekScalMat>::AllocateSharedPtr(invfactor,D));
}
}
return returnval;
}
} // end of namespace
}//end of namespace
//
// $Log: SegExp.cpp,v $
// Revision 1.28 2007/08/10 03:38:08 jfrazier
// Updated with new rev of NekManager.
//
// Revision 1.27 2007/07/30 21:00:06 sherwin
// Temporary fix in CreateMatrix
//
......
......@@ -142,10 +142,14 @@ namespace Nektar
DNekMatSharedPtr GetStdMatrix(const StdRegions::StdMatrixKey &mkey);
DNekScalMatSharedPtr CreateMatrix(const MatrixKey &mkey);
DNekBlkMatSharedPtr GetStdStaticCondMatrix(const StdRegions::StdMatrixKey &mkey);
DNekScalBlkMatSharedPtr CreateStaticCondMatrix(const MatrixKey &mkey);
SpatialDomains::SegGeomSharedPtr m_geom;
SpatialDomains::GeomFactorsSharedPtr m_metricinfo;
LibUtilities::NekManager<MatrixKey, DNekScalMat, MatrixKey::opLess> m_matrixManager;
LibUtilities::NekManager<MatrixKey, DNekScalBlkMat, MatrixKey::opLess> m_staticCondMatrixManager;
/// \brief Inner product of \a inarray over region with respect to
/// the expansion basis \a base and return in \a outarray
......@@ -342,10 +346,15 @@ namespace Nektar
return StdExpansion::L2();
}
virtual DNekScalMatSharedPtr v_GetLocMatrix(MatrixKey &mkey)
virtual DNekScalMatSharedPtr v_GetLocMatrix(const MatrixKey &mkey)
{
return m_matrixManager[mkey];
}
virtual DNekScalBlkMatSharedPtr v_GetLocStaticCondMatrix(const MatrixKey &mkey)
{
return m_staticCondMatrixManager[mkey];
}
};
// type defines for use of SegExp in a boost vector
......@@ -359,6 +368,9 @@ namespace Nektar
//
// $Log: SegExp.h,v $
// Revision 1.23 2007/08/10 03:38:15 jfrazier
// Updated with new rev of NekManager.
//
// Revision 1.22 2007/07/28 05:09:33 sherwin
// Fixed version with updated MemoryManager
//
......
......@@ -59,7 +59,7 @@ namespace Nektar
}
ContExpList1D::ContExpList1D(const LibUtilities::BasisKey &Ba,
const SpatialDomains::MeshGraph1D &graph1D):
const SpatialDomains::MeshGraph1D &graph1D):
ExpList1D(Ba,graph1D)
{
ASSERTL1((Ba.GetBasisType() == LibUtilities::eModified_A)
......@@ -72,7 +72,7 @@ namespace Nektar
// setup mapping array
m_locToGloMap = MemoryManager<LocalToGlobalMap1D>::AllocateSharedPtr(m_ncoeffs,*m_exp,graph1D.GetDomain());
m_contNcoeffs = m_locToGloMap->GetTotGloLen();
m_contNcoeffs = m_locToGloMap->GetTotGloDofs();
m_contCoeffs = Array<OneD,NekDouble>(m_contNcoeffs,0.0);
}
......@@ -95,7 +95,7 @@ namespace Nektar
// setup mapping array
m_locToGloMap = MemoryManager<LocalToGlobalMap1D>::AllocateSharedPtr(m_ncoeffs,*m_exp,domain);
m_contNcoeffs = m_locToGloMap->GetTotGloLen();
m_contNcoeffs = m_locToGloMap->GetTotGloDofs();
m_contCoeffs = Array<OneD,NekDouble>(m_contNcoeffs,0.0);
}
......@@ -132,7 +132,7 @@ namespace Nektar
if(matrixIter == m_globalMat->end())
{
mass_matrix = GenGlobalLinSys(key);
mass_matrix = GenGlobalLinSys(key,m_locToGloMap);
(*m_globalMat)[key] = mass_matrix;
}
else
......@@ -140,8 +140,7 @@ namespace Nektar
mass_matrix = matrixIter->second;
}
DNekVec v(m_contNcoeffs,m_contCoeffs,eWrapper);
mass_matrix->GetLinSys()->Solve(v,v);
mass_matrix->Solve(m_contCoeffs,m_contCoeffs,m_locToGloMap);
m_transState = eContinuous;
m_physState = false;
}
......@@ -157,7 +156,7 @@ namespace Nektar
if(matrixIter == m_globalMat->end())
{
helm_matrix = GenGlobalLinSys(key);
helm_matrix = GenGlobalLinSys(key,m_locToGloMap);
(*m_globalMat)[key] = helm_matrix;
}
else
......@@ -166,8 +165,7 @@ namespace Nektar
}
DNekVec v(m_contNcoeffs,m_contCoeffs,eWrapper);
helm_matrix->GetLinSys()->Solve(v,v);
helm_matrix->Solve(m_contCoeffs,m_contCoeffs,m_locToGloMap);
m_transState = eContinuous;
m_physState = false;
}
......@@ -183,52 +181,14 @@ namespace Nektar
ExpList1D::BwdTrans(In);
}
GlobalLinSysSharedPtr ContExpList1D::GenGlobalLinSys(const GlobalLinSysKey &mkey)
{
int i,j,n,gid1,gid2,loc_lda,cnt;
DNekScalMatSharedPtr loc_mat;
StdRegions::StdExpansionVectorIter def;
DNekLinSysSharedPtr linsys;
GlobalLinSysSharedPtr returnlinsys;
NekDouble zero = 0.0;
DNekMatSharedPtr Gmat = MemoryManager<DNekMat>::AllocateSharedPtr(m_contNcoeffs,m_contNcoeffs,zero);
// fill global matrix
for(n = cnt = 0; n < (*m_exp).size(); ++n)
{
LocalRegions::MatrixKey matkey(mkey.GetLinSysType(),
(*m_exp)[n]->DetShapeType(),
*(*m_exp)[n],mkey.GetScaleFactor());
loc_mat = (*m_exp)[n]->GetLocMatrix(matkey);
loc_lda = loc_mat->GetColumns();
for(i = 0; i < loc_lda; ++i)
{
gid1 = m_locToGloMap->GetMap(cnt + i);
for(j = 0; j < loc_lda; ++j)
{
gid2 = m_locToGloMap->GetMap(cnt + j);
(*Gmat)(gid1,gid2) += (*loc_mat)(i,j);
}
}
cnt += (*m_exp)[n]->GetNcoeffs();
}
linsys = MemoryManager<DNekLinSys>::AllocateSharedPtr(Gmat);
returnlinsys = MemoryManager<GlobalLinSys>::AllocateSharedPtr(mkey,linsys);
return returnlinsys;
}
} //end of namespace
} //end of namespace
/**
* $Log: ContExpList1D.cpp,v $
* Revision 1.23 2007/09/25 14:25:29 pvos
* Update for helmholtz1D with different expansion orders
*
* Revision 1.22 2007/08/11 23:43:25 sherwin
* Expansion bases reader part for Helmholtz1D
*
......
......@@ -40,17 +40,20 @@
#include <MultiRegions/ExpList1D.h>
#include <MultiRegions/LocalToGlobalMap1D.h>
#include <MultiRegions/GlobalLinSys.h>
namespace Nektar
{
namespace MultiRegions
{
class ContExpList1D: public ExpList1D
{
public:
ContExpList1D();
ContExpList1D(const LibUtilities::BasisKey &Ba,
const SpatialDomains::MeshGraph1D &graph1D);
class ContExpList1D:
public ExpList1D
{
public:
ContExpList1D();
//ContExpList1D(const LibUtilities::BasisKey &Ba,
//const SpatialDomains::MeshGraph1D &graph1D);
ContExpList1D(const LibUtilities::BasisKey &Ba,
const SpatialDomains::MeshGraph1D &graph1D);
ContExpList1D(SpatialDomains::MeshGraph1D &graph1D);
ContExpList1D(const ContExpList1D &In);
~ContExpList1D();
......@@ -70,7 +73,7 @@ namespace Nektar
{
m_locToGloMap->ContToLocal(inarray,outarray);
}
inline void LocalToCont()
{
m_locToGloMap->LocalToCont(m_coeffs,m_contCoeffs);
......@@ -96,20 +99,16 @@ namespace Nektar
void BwdTrans(const ExpList &In);
void GeneralMatrixOp(const StdRegions::MatrixType mtype,
const ConstArray<OneD,NekDouble> &inarray,
Array<OneD, NekDouble> &outarray,
NekDouble lambda);
protected:
int m_contNcoeffs;
Array<OneD, NekDouble> m_contCoeffs;