Commit 9b510dbe authored by Thomas Witkowski's avatar Thomas Witkowski
Browse files

Added new operators.

parent b2a88ee2
...@@ -80,7 +80,8 @@ namespace AMDiS { ...@@ -80,7 +80,8 @@ namespace AMDiS {
// estimate before first adaption // estimate before first adaption
if (adaptInfo->getTime() <= adaptInfo->getStartTime()) if (adaptInfo->getTime() <= adaptInfo->getStartTime())
problemIteration_->oneIteration(adaptInfo, ESTIMATE); problemIteration_->oneIteration(adaptInfo, ESTIMATE);
// increment time // increment time
adaptInfo->setTime(adaptInfo->getTime() + adaptInfo->getTimestep()); adaptInfo->setTime(adaptInfo->getTime() + adaptInfo->getTimestep());
......
...@@ -72,11 +72,13 @@ namespace AMDiS { ...@@ -72,11 +72,13 @@ namespace AMDiS {
DOFMatrix::~DOFMatrix() DOFMatrix::~DOFMatrix()
{ {
FUNCNAME("DOFMatrix::~DOFMatrix()"); FUNCNAME("DOFMatrix::~DOFMatrix()");
if (rowFESpace && rowFESpace->getAdmin()) if (rowFESpace && rowFESpace->getAdmin())
(const_cast<DOFAdmin*>(rowFESpace->getAdmin()))->removeDOFIndexed(this); (const_cast<DOFAdmin*>(rowFESpace->getAdmin()))->removeDOFIndexed(this);
if (boundaryManager)
if (boundaryManager) delete boundaryManager; delete boundaryManager;
if (inserter) delete inserter; if (inserter)
delete inserter;
} }
void DOFMatrix::print() const void DOFMatrix::print() const
...@@ -135,11 +137,10 @@ namespace AMDiS { ...@@ -135,11 +137,10 @@ namespace AMDiS {
int non_symmetric = !symmetric(); int non_symmetric = !symmetric();
if (non_symmetric) { if (non_symmetric)
MSG("matrix `%s' not symmetric.\n", name.data()); MSG("matrix `%s' not symmetric.\n", name.data());
} else { else
MSG("matrix `%s' is symmetric.\n", name.data()); MSG("matrix `%s' is symmetric.\n", name.data());
}
} }
DOFMatrix& DOFMatrix::operator=(const DOFMatrix& rhs) DOFMatrix& DOFMatrix::operator=(const DOFMatrix& rhs)
...@@ -154,12 +155,11 @@ namespace AMDiS { ...@@ -154,12 +155,11 @@ namespace AMDiS {
if (rhs.inserter == 0 && inserter == 0) if (rhs.inserter == 0 && inserter == 0)
matrix = rhs.matrix; matrix = rhs.matrix;
if (rhs.boundaryManager) { if (rhs.boundaryManager)
boundaryManager = new BoundaryManager(*rhs.boundaryManager); boundaryManager = new BoundaryManager(*rhs.boundaryManager);
} else { else
boundaryManager = NULL; boundaryManager = NULL;
}
nRow = rhs.nRow; nRow = rhs.nRow;
nCol = rhs.nCol; nCol = rhs.nCol;
elementMatrix.change_dim(nRow, nCol); elementMatrix.change_dim(nRow, nCol);
...@@ -302,9 +302,8 @@ namespace AMDiS { ...@@ -302,9 +302,8 @@ namespace AMDiS {
{ {
FUNCNAME("DOFMatrix::assemble2()"); FUNCNAME("DOFMatrix::assemble2()");
if (!op && operators.size() == 0) { if (!op && operators.size() == 0)
return; return;
}
set_to_zero(elementMatrix); set_to_zero(elementMatrix);
...@@ -335,8 +334,7 @@ namespace AMDiS { ...@@ -335,8 +334,7 @@ namespace AMDiS {
{ {
// call the operatos cleanup procedures // call the operatos cleanup procedures
for (std::vector<Operator*>::iterator it = operators.begin(); for (std::vector<Operator*>::iterator it = operators.begin();
it != operators.end(); it != operators.end(); ++it)
++it)
(*it)->finishAssembling(); (*it)->finishAssembling();
} }
......
...@@ -382,6 +382,23 @@ namespace AMDiS { ...@@ -382,6 +382,23 @@ namespace AMDiS {
boundaryManager = bm; boundaryManager = bm;
} }
/// Calculates the average of non zero entries per row in matrix.
void calculateNnz()
{
nnzPerRow = 0;
if (num_rows(matrix) != 0)
nnzPerRow = int(double(matrix.nnz()) / num_rows(matrix) * 1.2);
if (nnzPerRow < 5)
nnzPerRow= 5;
}
/// Returns \ref nnzPerRow.
int getNnz()
{
return nnzPerRow;
}
private: private:
template <typename T> template <typename T>
void s_write(::std::ostream &out, const T& value) void s_write(::std::ostream &out, const T& value)
...@@ -536,6 +553,13 @@ namespace AMDiS { ...@@ -536,6 +553,13 @@ namespace AMDiS {
*/ */
std::set<int> applyDBCs; std::set<int> applyDBCs;
/* \brief
* Number of non zero entries per row (average). For instationary problems this
* information may be used in the next timestep to accelerate insertion of
* elemnts during assembling.
*/
int nnzPerRow;
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS #ifdef HAVE_PARALLEL_DOMAIN_AMDIS
std::map<DegreeOfFreedom, bool> isRankDOF; std::map<DegreeOfFreedom, bool> isRankDOF;
#endif #endif
......
...@@ -970,11 +970,9 @@ namespace AMDiS { ...@@ -970,11 +970,9 @@ namespace AMDiS {
#ifdef _OPENMP #ifdef _OPENMP
#pragma omp parallel for schedule(dynamic, 25000) default(shared) private(i) #pragma omp parallel for schedule(dynamic, 25000) default(shared) private(i)
#endif #endif
for (i = 0; i < maxI; i++) { for (i = 0; i < maxI; i++)
if (!admin->isDOFFree(i)) { if (!admin->isDOFFree(i))
y[i] = alpha * y[i] + x[i]; y[i] = alpha * y[i] + x[i];
}
}
} }
template<typename T> template<typename T>
...@@ -984,12 +982,9 @@ namespace AMDiS { ...@@ -984,12 +982,9 @@ namespace AMDiS {
{ {
typename DOFVector<T>::Iterator vIterator(dynamic_cast<DOFIndexed<T>*>(const_cast<DOFVector<T>*>(&v)), USED_DOFS); typename DOFVector<T>::Iterator vIterator(dynamic_cast<DOFIndexed<T>*>(const_cast<DOFVector<T>*>(&v)), USED_DOFS);
typename DOFVector<T>::Iterator rIterator(dynamic_cast<DOFIndexed<T>*>(&result), USED_DOFS); typename DOFVector<T>::Iterator rIterator(dynamic_cast<DOFIndexed<T>*>(&result), USED_DOFS);
for(vIterator.reset(), rIterator.reset(); for (vIterator.reset(), rIterator.reset();
!vIterator.end(); !vIterator.end(); ++vIterator, ++rIterator)
++vIterator, ++rIterator) *rIterator = scal * (*vIterator);
{
*rIterator = scal * (*vIterator);
};
return result; return result;
} }
...@@ -1001,12 +996,10 @@ namespace AMDiS { ...@@ -1001,12 +996,10 @@ namespace AMDiS {
{ {
typename DOFVector<T>::Iterator vIterator(dynamic_cast<DOFIndexed<T>*>(const_cast<DOFVector<T>*>(&v)), USED_DOFS); typename DOFVector<T>::Iterator vIterator(dynamic_cast<DOFIndexed<T>*>(const_cast<DOFVector<T>*>(&v)), USED_DOFS);
typename DOFVector<T>::Iterator rIterator(dynamic_cast<DOFIndexed<T>*>(&result), USED_DOFS); typename DOFVector<T>::Iterator rIterator(dynamic_cast<DOFIndexed<T>*>(&result), USED_DOFS);
for(vIterator.reset(), rIterator.reset(); for (vIterator.reset(), rIterator.reset();
!vIterator.end(); !vIterator.end(); ++vIterator, ++rIterator)
++vIterator, ++rIterator) *rIterator = (*vIterator) + scal;
{
*rIterator = (*vIterator) + scal;
};
return result; return result;
} }
...@@ -1018,20 +1011,18 @@ namespace AMDiS { ...@@ -1018,20 +1011,18 @@ namespace AMDiS {
typename DOFVector<T>::Iterator v1Iterator(dynamic_cast<DOFIndexed<T>*>(const_cast<DOFVector<T>*>(&v1)), USED_DOFS); typename DOFVector<T>::Iterator v1Iterator(dynamic_cast<DOFIndexed<T>*>(const_cast<DOFVector<T>*>(&v1)), USED_DOFS);
typename DOFVector<T>::Iterator v2Iterator(dynamic_cast<DOFIndexed<T>*>(const_cast<DOFVector<T>*>(&v2)), USED_DOFS); typename DOFVector<T>::Iterator v2Iterator(dynamic_cast<DOFIndexed<T>*>(const_cast<DOFVector<T>*>(&v2)), USED_DOFS);
typename DOFVector<T>::Iterator rIterator(dynamic_cast<DOFIndexed<T>*>(&result), USED_DOFS); typename DOFVector<T>::Iterator rIterator(dynamic_cast<DOFIndexed<T>*>(&result), USED_DOFS);
for(v1Iterator.reset(), v2Iterator.reset(), rIterator.reset(); for (v1Iterator.reset(), v2Iterator.reset(), rIterator.reset();
!v1Iterator.end(); !v1Iterator.end(); ++v1Iterator, ++v2Iterator, ++rIterator)
++v1Iterator, ++v2Iterator, ++rIterator) *rIterator = (*v1Iterator) + (*v2Iterator);
{
*rIterator = (*v1Iterator) + (*v2Iterator);
};
return result;
return result;
} }
template<typename T> template<typename T>
const T *DOFVectorBase<T>::getLocalVector(const Element *el, T *d) const const T *DOFVectorBase<T>::getLocalVector(const Element *el, T *d) const
{ {
static T* localVec = NULL; static T* localVec = NULL;
static int localVecSize = 0;
const DOFAdmin* admin = feSpace->getAdmin(); const DOFAdmin* admin = feSpace->getAdmin();
T *result; T *result;
...@@ -1039,8 +1030,14 @@ namespace AMDiS { ...@@ -1039,8 +1030,14 @@ namespace AMDiS {
if (d) { if (d) {
result = d; result = d;
} else { } else {
if(localVec) delete [] localVec; if (localVec && nBasFcts > localVecSize) {
localVec = new T[nBasFcts]; delete [] localVec;
localVec = new T[nBasFcts];
}
if (!localVec)
localVec = new T[nBasFcts];
localVecSize = nBasFcts;
result = localVec; result = localVec;
} }
......
...@@ -47,8 +47,8 @@ namespace AMDiS { ...@@ -47,8 +47,8 @@ namespace AMDiS {
&standardSubAssemblersGrdPhi); &standardSubAssemblersGrdPhi);
int myRank = omp_get_thread_num(); int myRank = omp_get_thread_num();
std::vector<OperatorTerm*> opTerms std::vector<OperatorTerm*> opTerms =
= (type == GRD_PSI) ? op->firstOrderGrdPsi[myRank] : op->firstOrderGrdPhi[myRank]; (type == GRD_PSI) ? op->firstOrderGrdPsi[myRank] : op->firstOrderGrdPhi[myRank];
// check if a assembler is needed at all // check if a assembler is needed at all
if (opTerms.size() == 0) if (opTerms.size() == 0)
...@@ -103,15 +103,15 @@ namespace AMDiS { ...@@ -103,15 +103,15 @@ namespace AMDiS {
Stand10::Stand10(Operator *op, Assembler *assembler, Quadrature *quad) Stand10::Stand10(Operator *op, Assembler *assembler, Quadrature *quad)
: FirstOrderAssembler(op, assembler, quad, false, GRD_PSI) : FirstOrderAssembler(op, assembler, quad, false, GRD_PSI)
{} {
psi = owner->getRowFESpace()->getBasisFcts();
phi = owner->getColFESpace()->getBasisFcts();
}
void Stand10::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix& mat) void Stand10::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix& mat)
{ {
DimVec<double> grdPsi(dim, DEFAULT_VALUE, 0.0); DimVec<double> grdPsi(dim, DEFAULT_VALUE, 0.0);
const BasisFunction *psi = owner->getRowFESpace()->getBasisFcts();
const BasisFunction *phi = owner->getColFESpace()->getBasisFcts();
int nPoints = quadrature->getNumPoints(); int nPoints = quadrature->getNumPoints();
int myRank = omp_get_thread_num(); int myRank = omp_get_thread_num();
VectorOfFixVecs<DimVec<double> > &Lb = tmpLb[myRank]; VectorOfFixVecs<DimVec<double> > &Lb = tmpLb[myRank];
...@@ -132,7 +132,7 @@ namespace AMDiS { ...@@ -132,7 +132,7 @@ namespace AMDiS {
for (int i = 0; i < nRow; i++) { for (int i = 0; i < nRow; i++) {
(*(psi->getGrdPhi(i)))(quadrature->getLambda(iq), grdPsi); (*(psi->getGrdPhi(i)))(quadrature->getLambda(iq), grdPsi);
for (int j = 0; j < nCol; j++) for (int j = 0; j < nCol; j++)
mat[i][j] += quadrature->getWeight(iq) * (Lb[iq] * grdPsi) * phival[j]; mat[i][j] += quadrature->getWeight(iq) * phival[j] * (Lb[iq] * grdPsi);
} }
} }
} }
...@@ -140,7 +140,6 @@ namespace AMDiS { ...@@ -140,7 +140,6 @@ namespace AMDiS {
void Stand10::calculateElementVector(const ElInfo *elInfo, ElementVector& vec) void Stand10::calculateElementVector(const ElInfo *elInfo, ElementVector& vec)
{ {
DimVec<double> grdPsi(dim, DEFAULT_VALUE, 0.0); DimVec<double> grdPsi(dim, DEFAULT_VALUE, 0.0);
const BasisFunction *psi = owner->getRowFESpace()->getBasisFcts();
int nPoints = quadrature->getNumPoints(); int nPoints = quadrature->getNumPoints();
int myRank = omp_get_thread_num(); int myRank = omp_get_thread_num();
VectorOfFixVecs<DimVec<double> > &Lb = tmpLb[myRank]; VectorOfFixVecs<DimVec<double> > &Lb = tmpLb[myRank];
...@@ -189,13 +188,11 @@ namespace AMDiS { ...@@ -189,13 +188,11 @@ namespace AMDiS {
VectorOfFixVecs<DimVec<double> > &Lb = tmpLb[myRank]; VectorOfFixVecs<DimVec<double> > &Lb = tmpLb[myRank];
Lb.resize(nPoints); Lb.resize(nPoints);
for (int iq = 0; iq < nPoints; iq++) { for (int iq = 0; iq < nPoints; iq++)
Lb[iq].set(0.0); Lb[iq].set(0.0);
}
for (int i = 0; i < static_cast<int>(terms[myRank].size()); i++) { for (int i = 0; i < static_cast<int>(terms[myRank].size()); i++)
(static_cast<FirstOrderTerm*>((terms[myRank][i])))->getLb(elInfo, nPoints, Lb); (static_cast<FirstOrderTerm*>((terms[myRank][i])))->getLb(elInfo, nPoints, Lb);
}
for (int iq = 0; iq < nPoints; iq++) { for (int iq = 0; iq < nPoints; iq++) {
Lb[iq] *= elInfo->getDet(); Lb[iq] *= elInfo->getDet();
...@@ -203,10 +200,9 @@ namespace AMDiS { ...@@ -203,10 +200,9 @@ namespace AMDiS {
const double *phi = phiFast->getPhi(iq); const double *phi = phiFast->getPhi(iq);
grdPsi = psiFast->getGradient(iq); grdPsi = psiFast->getGradient(iq);
for (int i = 0; i < nRow; i++) { for (int i = 0; i < nRow; i++)
for (int j = 0; j < nCol; j++) for (int j = 0; j < nCol; j++)
mat[i][j] += quadrature->getWeight(iq) * (Lb[iq] * (*grdPsi)[i]) * phi[j]; mat[i][j] += quadrature->getWeight(iq) * (Lb[iq] * (*grdPsi)[i]) * phi[j];
}
} }
} }
...@@ -232,21 +228,17 @@ namespace AMDiS { ...@@ -232,21 +228,17 @@ namespace AMDiS {
VectorOfFixVecs<DimVec<double> > &Lb = tmpLb[myRank]; VectorOfFixVecs<DimVec<double> > &Lb = tmpLb[myRank];
Lb.resize(nPoints); Lb.resize(nPoints);
for (int iq = 0; iq < nPoints; iq++) { for (int iq = 0; iq < nPoints; iq++)
Lb[iq].set(0.0); Lb[iq].set(0.0);
} for (int i = 0; i < static_cast<int>(terms[myRank].size()); i++)
for (int i = 0; i < static_cast<int>(terms[myRank].size()); i++) {
(static_cast<FirstOrderTerm*>((terms[myRank][i])))->getLb(elInfo, nPoints, Lb); (static_cast<FirstOrderTerm*>((terms[myRank][i])))->getLb(elInfo, nPoints, Lb);
}
for (int iq = 0; iq < nPoints; iq++) { for (int iq = 0; iq < nPoints; iq++) {
Lb[iq] *= elInfo->getDet(); Lb[iq] *= elInfo->getDet();
grdPsi = psiFast->getGradient(iq); grdPsi = psiFast->getGradient(iq);
for (int i = 0; i < nRow; i++) { for (int i = 0; i < nRow; i++)
vec[i] += quadrature->getWeight(iq) * (Lb[iq] * (*grdPsi)[i]); vec[i] += quadrature->getWeight(iq) * (Lb[iq] * (*grdPsi)[i]);
}
} }
} }
...@@ -305,13 +297,13 @@ namespace AMDiS { ...@@ -305,13 +297,13 @@ namespace AMDiS {
VectorOfFixVecs<DimVec<double> > VectorOfFixVecs<DimVec<double> >
grdPhi(assembler->getRowFESpace()->getMesh()->getDim(), nCol, NO_INIT); grdPhi(assembler->getRowFESpace()->getMesh()->getDim(), nCol, NO_INIT);
tmpGrdPhi.resize(omp_get_overall_max_threads(), grdPhi); tmpGrdPhi.resize(omp_get_overall_max_threads(), grdPhi);
psi = owner->getRowFESpace()->getBasisFcts();
phi = owner->getColFESpace()->getBasisFcts();
} }
void Stand01::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix& mat) void Stand01::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix& mat)
{ {
const BasisFunction *psi = owner->getRowFESpace()->getBasisFcts();
const BasisFunction *phi = owner->getColFESpace()->getBasisFcts();
int nPoints = quadrature->getNumPoints(); int nPoints = quadrature->getNumPoints();
int myRank = omp_get_thread_num(); int myRank = omp_get_thread_num();
VectorOfFixVecs<DimVec<double> >& Lb = tmpLb[myRank]; VectorOfFixVecs<DimVec<double> >& Lb = tmpLb[myRank];
...@@ -320,6 +312,7 @@ namespace AMDiS { ...@@ -320,6 +312,7 @@ namespace AMDiS {
for (int iq = 0; iq < nPoints; iq++) for (int iq = 0; iq < nPoints; iq++)
Lb[iq].set(0.0); Lb[iq].set(0.0);
for (int i = 0; i < static_cast<int>(terms[myRank].size()); i++) for (int i = 0; i < static_cast<int>(terms[myRank].size()); i++)
(static_cast<FirstOrderTerm*>((terms[myRank][i])))->getLb(elInfo, nPoints, Lb); (static_cast<FirstOrderTerm*>((terms[myRank][i])))->getLb(elInfo, nPoints, Lb);
...@@ -332,7 +325,7 @@ namespace AMDiS { ...@@ -332,7 +325,7 @@ namespace AMDiS {
for (int i = 0; i < nRow; i++) { for (int i = 0; i < nRow; i++) {
double psival = (*(psi->getPhi(i)))(quadrature->getLambda(iq)); double psival = (*(psi->getPhi(i)))(quadrature->getLambda(iq));
for (int j = 0; j < nCol; j++) for (int j = 0; j < nCol; j++)
mat[i][j] += quadrature->getWeight(iq) * ((Lb[iq] * psival) * grdPhi[j]); mat[i][j] += quadrature->getWeight(iq) * psival * (Lb[iq] * grdPhi[j]);
} }
} }
} }
......
...@@ -78,6 +78,9 @@ namespace AMDiS { ...@@ -78,6 +78,9 @@ namespace AMDiS {
/// Implements SubAssembler::calculateElementVector(). /// Implements SubAssembler::calculateElementVector().
void calculateElementVector(const ElInfo *, ElementVector&); void calculateElementVector(const ElInfo *, ElementVector&);
protected:
const BasisFunction *psi, *phi;
}; };
...@@ -104,6 +107,8 @@ namespace AMDiS { ...@@ -104,6 +107,8 @@ namespace AMDiS {
protected: protected:
std::vector<VectorOfFixVecs<DimVec<double> > > tmpGrdPhi; std::vector<VectorOfFixVecs<DimVec<double> > > tmpGrdPhi;
const BasisFunction *psi, *phi;
}; };
......
...@@ -233,17 +233,20 @@ namespace AMDiS { ...@@ -233,17 +233,20 @@ namespace AMDiS {
} }
/// Returns the \ref size of this VectorOfFixVecs /// Returns the \ref size of this VectorOfFixVecs
inline int getSize() const { inline int getSize() const
{
return size; return size;
} }
/// Returns \ref dim /// Returns \ref dim
inline int getDim() const { inline int getDim() const
{
return dim; return dim;
} }
/// Returns the size of the contained FixVecs /// Returns the size of the contained FixVecs
inline int getSizeOfFixVec() const { inline int getSizeOfFixVec() const
{
return vec[0]->getSize(); return vec[0]->getSize();
} }
...@@ -515,7 +518,8 @@ namespace AMDiS { ...@@ -515,7 +518,8 @@ namespace AMDiS {
/// creates and inits and double array /// creates and inits and double array
double *createAndInitArray(int size, ...); double *createAndInitArray(int size, ...);
inline WorldVector<double> operator*(const WorldVector<double>& v, double d) { inline WorldVector<double> operator*(const WorldVector<double>& v, double d)
{
WorldVector<double> result = v; WorldVector<double> result = v;
result *= d; result *= d;
return result; return result;
......
...@@ -8,9 +8,8 @@ namespace AMDiS { ...@@ -8,9 +8,8 @@ namespace AMDiS {
thisIt != this->end(); thisIt != this->end();
thisIt++) { thisIt++) {
*thisIt = 0; *thisIt = 0;
for (T* vIt = v.begin(); vIt != v.end(); vIt++, mIt++) { for (T* vIt = v.begin(); vIt != v.end(); vIt++, mIt++)
*thisIt += *vIt * *mIt; *thisIt += *vIt * *mIt;
}
} }
} }
...@@ -19,9 +18,8 @@ namespace AMDiS { ...@@ -19,9 +18,8 @@ namespace AMDiS {
{ {
double erg = 0.0; double erg = 0.0;
for (int i = 0; i < a.getSize() ; i++) { for (int i = 0; i < a.getSize() ; i++)
erg = erg + ((a[i] - b[i]) * (a[i] - b[i])); erg = erg + ((a[i] - b[i]) * (a[i] - b[i]));
}
return sqrt(erg); return sqrt(erg);
} }
...@@ -39,9 +37,8 @@ namespace AMDiS { ...@@ -39,9 +37,8 @@ namespace AMDiS {
thisIt++) { thisIt++) {
*thisIt = 0; *thisIt = 0;
for (T* vIt = v.begin(); vIt != v.end(); vIt++, mIt++) { for (T* vIt = v.begin(); vIt != v.end(); vIt++, mIt++)
*thisIt += *vIt * *mIt; *thisIt += *vIt * *mIt;
}
} }
} }
...@@ -70,9 +67,8 @@ namespace AMDiS { ...@@ -70,9 +67,8 @@ namespace AMDiS {