From 06dbd114430271eacf34b063df4225ed2b938a65 Mon Sep 17 00:00:00 2001
From: Thomas Witkowski <thomas.witkowski@gmx.de>
Date: Fri, 20 Nov 2009 08:35:36 +0000
Subject: [PATCH] Made some operators a little bit faster.
---
AMDiS/src/DOFVector.hh | 12 +-
AMDiS/src/FirstOrderAssembler.cc | 18 +--
AMDiS/src/Operator.cc | 88 +++++++++---
AMDiS/src/Operator.h | 178 +++++++++++++++----------
AMDiS/src/ProblemVec.cc | 213 ++++++++++++++++++++++++++++--
AMDiS/src/ProblemVec.h | 13 ++
AMDiS/src/SecondOrderAssembler.cc | 15 +--
AMDiS/src/ZeroOrderAssembler.cc | 34 ++---
AMDiS/src/ZeroOrderAssembler.h | 3 +
9 files changed, 436 insertions(+), 138 deletions(-)
diff --git a/AMDiS/src/DOFVector.hh b/AMDiS/src/DOFVector.hh
index e78ea500..83d38ed7 100644
--- a/AMDiS/src/DOFVector.hh
+++ b/AMDiS/src/DOFVector.hh
@@ -981,13 +981,13 @@ namespace AMDiS {
{
FUNCNAME("DOFVector<T>::getVecAtQPs()");
- TEST_EXIT(quad || quadFast)("neither quad nor quadFast defined\n");
+ TEST_EXIT_DBG(quad || quadFast)("neither quad nor quadFast defined\n");
if (quad && quadFast)
- TEST_EXIT(quad == quadFast->getQuadrature())
+ TEST_EXIT_DBG(quad == quadFast->getQuadrature())
("quad != quadFast->quadrature\n");
- TEST_EXIT(!quadFast || quadFast->getBasisFunctions() == feSpace->getBasisFcts())
+ TEST_EXIT_DBG(!quadFast || quadFast->getBasisFunctions() == feSpace->getBasisFcts())
("invalid basis functions");
Element *el = elInfo->getElement();
@@ -1012,8 +1012,8 @@ namespace AMDiS {
result = localvec;
}
-
- T *localVec = new T[nBasFcts];
+
+ T *localVec = localVectors[omp_get_thread_num()];
getLocalVector(el, localVec);
for (int i = 0; i < nPoints; i++) {
@@ -1025,8 +1025,6 @@ namespace AMDiS {
((*(basFcts->getPhi(j)))(quad->getLambda(i))));
}
- delete [] localVec;
-
return const_cast<const T*>(result);
}
diff --git a/AMDiS/src/FirstOrderAssembler.cc b/AMDiS/src/FirstOrderAssembler.cc
index 6122c7be..b1b34018 100644
--- a/AMDiS/src/FirstOrderAssembler.cc
+++ b/AMDiS/src/FirstOrderAssembler.cc
@@ -338,8 +338,6 @@ namespace AMDiS {
void Quad01::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix& mat)
{
- VectorOfFixVecs<DimVec<double> > *grdPhi;
-
if (firstCall) {
#ifdef _OPENMP
#pragma omp critical
@@ -363,25 +361,23 @@ namespace AMDiS {
for (int i = 0; i < static_cast<int>(terms[myRank].size()); i++)
(static_cast<FirstOrderTerm*>((terms[myRank][i])))->getLb(elInfo, nPoints, Lb);
-
+
for (int iq = 0; iq < nPoints; iq++) {
Lb[iq] *= elInfo->getDet();
-
const double *psi = psiFast->getPhi(iq);
- grdPhi = phiFast->getGradient(iq);
- for (int i = 0; i < nRow; i++) {
- double factor = quadrature->getWeight(iq) * psi[i];
- for (int j = 0; j < nCol; j++)
- mat[i][j] += factor * (Lb[iq] * (*grdPhi)[j]);
+ for (int i = 0; i < nCol; i++) {
+ double factor =
+ quadrature->getWeight(iq) * (Lb[iq] * phiFast->getGradient(iq, i));
+ for (int j = 0; j < nRow; j++)
+ mat[j][i] += factor * psi[j];
}
}
}
Pre01::Pre01(Operator *op, Assembler *assembler, Quadrature *quad)
: FirstOrderAssembler(op, assembler, quad, true, GRD_PHI)
- {
- }
+ {}
void Pre01::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix& mat)
{
diff --git a/AMDiS/src/Operator.cc b/AMDiS/src/Operator.cc
index 9f5bc807..a5757ecf 100644
--- a/AMDiS/src/Operator.cc
+++ b/AMDiS/src/Operator.cc
@@ -697,6 +697,17 @@ namespace AMDiS {
auxFESpaces.push_back(dv->getFESpace());
}
+ VecAtQP_FOT::VecAtQP_FOT(DOFVectorBase<double> *dv,
+ AbstractFunction<double, double> *af,
+ int bIdx)
+ : FirstOrderTerm(af->getDegree()), vec(dv), f(af)
+ {
+ TEST_EXIT(dv)("No vector!\n");
+
+ bOne = bIdx;
+ auxFESpaces.push_back(dv->getFESpace());
+ }
+
VectorGradient_FOT::VectorGradient_FOT(DOFVectorBase<double> *dv,
AbstractFunction<WorldVector<double>, WorldVector<double> > *af)
: FirstOrderTerm(af->getDegree()), vec(dv), f(af)
@@ -824,6 +835,15 @@ namespace AMDiS {
auxFESpaces.push_back(dv2->getFESpace());
}
+ Vec2AtQP_FOT::Vec2AtQP_FOT(DOFVectorBase<double> *dv1, DOFVectorBase<double> *dv2,
+ int bIdx)
+ : FirstOrderTerm(8), vec1(dv1), vec2(dv2)
+ {
+ bOne = bIdx;
+ auxFESpaces.push_back(dv1->getFESpace());
+ auxFESpaces.push_back(dv2->getFESpace());
+ }
+
void Vec2AtQP_FOT::initElement(const ElInfo* elInfo,
SubAssembler* subAssembler,
Quadrature *quad)
@@ -836,16 +856,21 @@ namespace AMDiS {
VectorOfFixVecs<DimVec<double> >& Lb) const
{
const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
- for (int iq = 0; iq < nPoints; iq++) {
- if (b)
+
+ if (bOne > -1) {
+ for (int iq = 0; iq < nPoints; iq++)
+ lb_one(Lambda, Lb[iq], vec1AtQPs[iq] * vec2AtQPs[iq]);
+ } else if (b) {
+ for (int iq = 0; iq < nPoints; iq++)
lb(Lambda, *b, Lb[iq], vec1AtQPs[iq] * vec2AtQPs[iq]);
- else
+ } else {
+ for (int iq = 0; iq < nPoints; iq++)
l1(Lambda, Lb[iq], vec1AtQPs[iq] * vec2AtQPs[iq]);
}
}
void Vec2AtQP_FOT::eval(int nPoints,
- const double *uhAtQP,
+ const double *uhAtQP,
const WorldVector<double> *grdUhAtQP,
const WorldMatrix<double> *D2UhAtQP,
double *result,
@@ -858,14 +883,34 @@ namespace AMDiS {
Vec3FctAtQP_FOT::Vec3FctAtQP_FOT(DOFVectorBase<double> *dv1, DOFVectorBase<double> *dv2,DOFVectorBase<double> *dv3,
BinaryAbstractFunction<double, double, double> *f_,
- WorldVector<double> *b_)
- : FirstOrderTerm(10), vec1(dv1), vec2(dv2), vec3(dv3), f(f_), b(b_)
+ WorldVector<double> *bvec)
+ : FirstOrderTerm(10),
+ vec1(dv1),
+ vec2(dv2),
+ vec3(dv3),
+ f(f_),
+ b(bvec)
{
auxFESpaces.push_back(dv1->getFESpace());
auxFESpaces.push_back(dv2->getFESpace());
auxFESpaces.push_back(dv3->getFESpace());
}
+ Vec3FctAtQP_FOT::Vec3FctAtQP_FOT(DOFVectorBase<double> *dv1, DOFVectorBase<double> *dv2,DOFVectorBase<double> *dv3,
+ BinaryAbstractFunction<double, double, double> *f_,
+ int b)
+ : FirstOrderTerm(10),
+ vec1(dv1),
+ vec2(dv2),
+ vec3(dv3),
+ f(f_)
+ {
+ bOne = b;
+ auxFESpaces.push_back(dv1->getFESpace());
+ auxFESpaces.push_back(dv2->getFESpace());
+ auxFESpaces.push_back(dv3->getFESpace());
+ }
+
void Vec3FctAtQP_FOT::initElement(const ElInfo* elInfo,
SubAssembler* subAssembler,
Quadrature *quad)
@@ -876,19 +921,25 @@ namespace AMDiS {
}
- void Vec3FctAtQP_FOT::getLb(const ElInfo *elInfo, int nPoints, VectorOfFixVecs<DimVec<double> >& Lb) const
+ void Vec3FctAtQP_FOT::getLb(const ElInfo *elInfo, int nPoints,
+ VectorOfFixVecs<DimVec<double> >& Lb) const
{
const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
- for (int iq = 0; iq < nPoints; iq++) {
- if (b)
- lb(Lambda, *b, Lb[iq], (vec1AtQPs[iq])*(*f)(vec2AtQPs[iq],vec3AtQPs[iq]));
- else
- l1(Lambda, Lb[iq], (vec1AtQPs[iq])*(*f)(vec2AtQPs[iq],vec3AtQPs[iq]));
+ if (bOne > -1) {
+ for (int iq = 0; iq < nPoints; iq++)
+ lb_one(Lambda, Lb[iq], vec1AtQPs[iq] * (*f)(vec2AtQPs[iq], vec3AtQPs[iq]));
+ } else {
+ for (int iq = 0; iq < nPoints; iq++) {
+ if (b)
+ lb(Lambda, *b, Lb[iq], vec1AtQPs[iq] * (*f)(vec2AtQPs[iq], vec3AtQPs[iq]));
+ else
+ l1(Lambda, Lb[iq], vec1AtQPs[iq] * (*f)(vec2AtQPs[iq], vec3AtQPs[iq]));
+ }
}
}
void Vec3FctAtQP_FOT::eval(int nPoints,
- const double *uhAtQP,
+ const double *uhAtQP,
const WorldVector<double> *grdUhAtQP,
const WorldMatrix<double> *D2UhAtQP,
double *result,
@@ -1582,10 +1633,15 @@ namespace AMDiS {
VectorOfFixVecs<DimVec<double> >& Lb) const
{
const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
- for (int iq = 0; iq < nPoints; iq++) {
- if (b)
+
+ if (bOne > -1) {
+ for (int iq = 0; iq < nPoints; iq++)
+ lb_one(Lambda, Lb[iq], (*f)(vecAtQPs[iq]));
+ } else if (b) {
+ for (int iq = 0; iq < nPoints; iq++)
lb(Lambda, *b, Lb[iq], (*f)(vecAtQPs[iq]));
- else
+ } else {
+ for (int iq = 0; iq < nPoints; iq++)
l1(Lambda, Lb[iq], (*f)(vecAtQPs[iq]));
}
}
diff --git a/AMDiS/src/Operator.h b/AMDiS/src/Operator.h
index 1094f876..c274d643 100644
--- a/AMDiS/src/Operator.h
+++ b/AMDiS/src/Operator.h
@@ -54,7 +54,9 @@ namespace AMDiS {
OperatorTerm(int deg)
: properties(0),
degree(deg),
- auxFESpaces(0)
+ dimOfWorld(Global::getGeo(WORLD)),
+ auxFESpaces(0),
+ bOne(-1)
{}
/// Destructor.
@@ -65,14 +67,11 @@ namespace AMDiS {
* each OperatorTerm belonging to this SubAssembler. E.g., vectors
* and coordinates at quadrature points can be calculated here.
*/
- virtual void initElement(const ElInfo*,
- SubAssembler*,
- Quadrature *quad = NULL)
+ virtual void initElement(const ElInfo*, SubAssembler*, Quadrature *quad = NULL)
{}
- virtual void initElement(const ElInfo* largeElInfo,
- const ElInfo* smallElInfo,
- SubAssembler*,
+ virtual void initElement(const ElInfo* largeElInfo, const ElInfo* smallElInfo,
+ SubAssembler*,
Quadrature *quad = NULL)
{}
@@ -100,6 +99,12 @@ namespace AMDiS {
return degree;
}
+ /// Sets one component of the b vector to be one. See \ref bOne.
+ void setB(int b)
+ {
+ bOne = b;
+ }
+
/// Evaluation of the OperatorTerm at all quadrature points.
virtual void eval(int nPoints,
const double *uhAtQP,
@@ -166,12 +171,11 @@ namespace AMDiS {
double factor) const;
/// Evaluation of \f$ \Lambda \cdot b\f$.
- static inline void lb(const DimVec<WorldVector<double> >& Lambda,
- const WorldVector<double>& b,
- DimVec<double>& Lb,
- double factor)
+ inline void lb(const DimVec<WorldVector<double> >& Lambda,
+ const WorldVector<double>& b,
+ DimVec<double>& Lb,
+ double factor) const
{
- const int dimOfWorld = Global::getGeo(WORLD);
const int dim = Lambda.size() - 1;
for (int i = 0; i <= dim; i++) {
@@ -184,15 +188,25 @@ namespace AMDiS {
}
}
+ /// Evaluation of \f$ \Lambda \cdot b\f$.
+ inline void lb_one(const DimVec<WorldVector<double> >& Lambda,
+ DimVec<double>& Lb,
+ double factor) const
+ {
+ const int dim = Lambda.size();
+
+ for (int i = 0; i < dim; i++)
+ Lb[i] += Lambda[i][bOne] * factor;
+ }
+
/** \brief
* Evaluation of \f$ \Lambda \cdot b\f$ if b contains the value 1.0 in
* each component.
*/
- static inline void l1(const DimVec<WorldVector<double> >& Lambda,
- DimVec<double>& Lb,
- double factor)
+ inline void l1(const DimVec<WorldVector<double> >& Lambda,
+ DimVec<double>& Lb,
+ double factor) const
{
- const int dimOfWorld = Global::getGeo(WORLD);
const int dim = Lambda.size() - 1;
for (int i = 0; i <= dim; i++) {
@@ -212,12 +226,24 @@ namespace AMDiS {
/// Polynomial degree of the term. Used to detemine the degree of the quadrature.
int degree;
+ /// Stores the dimension of the world.
+ int dimOfWorld;
+
/// List off all fe spaces, the operator term makes use off.
std::vector<const FiniteElemSpace*> auxFESpaces;
/// Pointer to the Operator this OperatorTerm belongs to.
Operator* operat;
+ /** \brief
+ * In many cases, the vector b in the evaluation \f$ \Lambda \cdot b\f$ has zeros
+ * in all components expect one that is set to one. Using the function \ref lb is
+ * then unnecessary time consuming. Instead, this variable defines the component
+ * of the vector b to be one. The function \ref lb_one is used if this variable is
+ * not -1.
+ */
+ int bOne;
+
/// Flag for piecewise constant terms
static const Flag PW_CONST;
@@ -1466,17 +1492,27 @@ namespace AMDiS {
: FirstOrderTerm(0), b(wv)
{}
- /** \brief
- * Implements FirstOrderTerm::getLb().
- */
+ /// Constructor.
+ Vector_FOT(int bIdx)
+ : FirstOrderTerm(0)
+ {
+ bOne = bIdx;
+ }
+
+ /// Implements FirstOrderTerm::getLb().
inline void getLb(const ElInfo *elInfo,
int nPoints,
VectorOfFixVecs<DimVec<double> >&Lb) const
{
const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
- for (int iq = 0; iq < nPoints; iq++)
- lb(Lambda, b, Lb[iq], 1.0);
+ if (bOne > -1) {
+ for (int iq = 0; iq < nPoints; iq++)
+ lb_one(Lambda, Lb[iq], 1.0);
+ } else {
+ for (int iq = 0; iq < nPoints; iq++)
+ lb(Lambda, b, Lb[iq], 1.0);
+ }
}
/// Implements FirstOrderTerm::eval().
@@ -1507,10 +1543,15 @@ namespace AMDiS {
{
public:
/// Constructor.
- VecAtQP_FOT(DOFVectorBase<double> *dv,
+ VecAtQP_FOT(DOFVectorBase<double> *dv,
AbstractFunction<double, double> *af,
WorldVector<double> *wv);
+ /// Constructor.
+ VecAtQP_FOT(DOFVectorBase<double> *dv,
+ AbstractFunction<double, double> *af,
+ int bIdx);
+
/// Implementation of \ref OperatorTerm::initElement().
void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
Quadrature *quad = NULL);
@@ -1931,22 +1972,22 @@ namespace AMDiS {
class FctVecAtQP_FOT : public FirstOrderTerm
{
public:
-
FctVecAtQP_FOT(DOFVectorBase<double> *dv,
AbstractFunction<double, WorldVector<double> > *f_,
WorldVector<double> *b_);
-
- void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
+
+ void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
Quadrature *quad = NULL);
-
- void getLb(const ElInfo *elInfo, int nPoints, VectorOfFixVecs<DimVec<double> >& Lb) const;
-
- void eval(int nPoints,
- const double *uhAtQP,
- const WorldVector<double> *grdUhAtQP,
- const WorldMatrix<double> *D2UhAtQP,
- double *result,
- double factor) const;
+
+ void getLb(const ElInfo *elInfo, int nPoints,
+ VectorOfFixVecs<DimVec<double> >& Lb) const;
+
+ void eval(int nPoints,
+ const double *uhAtQP,
+ const WorldVector<double> *grdUhAtQP,
+ const WorldMatrix<double> *D2UhAtQP,
+ double *result,
+ double factor) const;
protected:
DOFVectorBase<double>* vec;
@@ -1960,22 +2001,25 @@ namespace AMDiS {
class Vec2AtQP_FOT : public FirstOrderTerm
{
public:
-
Vec2AtQP_FOT(DOFVectorBase<double> *dv1, DOFVectorBase<double> *dv2,
WorldVector<double> *b_);
- void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
+ Vec2AtQP_FOT(DOFVectorBase<double> *dv1, DOFVectorBase<double> *dv2,
+ int bIdx);
+
+ void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
Quadrature *quad = NULL);
-
- void getLb(const ElInfo *elInfo, int nPoints, VectorOfFixVecs<DimVec<double> >& Lb) const;
-
- void eval(int nPoints,
- const double *uhAtQP,
+
+ void getLb(const ElInfo *elInfo, int nPoints,
+ VectorOfFixVecs<DimVec<double> >& Lb) const;
+
+ void eval(int nPoints,
+ const double *uhAtQP,
const WorldVector<double> *grdUhAtQP,
const WorldMatrix<double> *D2UhAtQP,
double *result,
double factor) const;
-
+
protected:
DOFVectorBase<double>* vec1;
DOFVectorBase<double>* vec2;
@@ -1988,18 +2032,22 @@ namespace AMDiS {
class Vec3FctAtQP_FOT : public FirstOrderTerm
{
public:
+ Vec3FctAtQP_FOT(DOFVectorBase<double> *dv1, DOFVectorBase<double> *dv2, DOFVectorBase<double> *dv3,
+ BinaryAbstractFunction<double, double, double> *f,
+ WorldVector<double> *b);
- Vec3FctAtQP_FOT(DOFVectorBase<double> *dv1, DOFVectorBase<double> *dv2,DOFVectorBase<double> *dv3,
- BinaryAbstractFunction<double, double, double> *f_,
- WorldVector<double> *b_);
+ Vec3FctAtQP_FOT(DOFVectorBase<double> *dv1, DOFVectorBase<double> *dv2, DOFVectorBase<double> *dv3,
+ BinaryAbstractFunction<double, double, double> *f,
+ int b);
void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
Quadrature *quad = NULL);
- void getLb(const ElInfo *elInfo, int nPoints, VectorOfFixVecs<DimVec<double> >& Lb) const;
+ void getLb(const ElInfo *elInfo, int nPoints,
+ VectorOfFixVecs<DimVec<double> >& Lb) const;
void eval(int nPoints,
- const double *uhAtQP,
+ const double *uhAtQP,
const WorldVector<double> *grdUhAtQP,
const WorldMatrix<double> *D2UhAtQP,
double *result,
@@ -3513,42 +3561,30 @@ namespace AMDiS {
const FiniteElemSpace *rowFESpace,
const FiniteElemSpace *colFESpace = NULL);
- /** \brief
- * Destructor.
- */
+ /// Destructor.
virtual ~Operator() {}
- /** \brief
- * Sets \ref optimized.
- */
+ /// Sets \ref optimized.
inline void useOptimizedAssembler(bool opt)
{
optimized = opt;
}
- /** \brief
- * Returns \ref optimized.
- */
+ /// Returns \ref optimized.
inline bool isOptimized()
{
return optimized;
}
- /** \brief
- * Adds a SecondOrderTerm to the Operator
- */
+ /// Adds a SecondOrderTerm to the Operator
template <typename T>
void addSecondOrderTerm(T *term);
- /** \brief
- * Adds a FirstOrderTerm to the Operator
- */
+ /// Adds a FirstOrderTerm to the Operator
template <typename T>
- void addFirstOrderTerm(T *term,
- FirstOrderType type = GRD_PHI);
- /** \brief
- * Adds a ZeroOrderTerm to the Operator
- */
+ void addFirstOrderTerm(T *term, FirstOrderType type = GRD_PHI);
+
+ /// Adds a ZeroOrderTerm to the Operator
template <typename T>
void addZeroOrderTerm(T *term);
@@ -3604,14 +3640,10 @@ namespace AMDiS {
return auxFESpaces;
}
- /** \brief
- * Sets \ref uhOld.
- */
+ /// Sets \ref uhOld.
void setUhOld(const DOFVectorBase<double> *uhOld);
- /** \brief
- * Returns \ref uhOld.
- */
+ /// Returns \ref uhOld.
inline const DOFVectorBase<double> *getUhOld()
{
return uhOld;
diff --git a/AMDiS/src/ProblemVec.cc b/AMDiS/src/ProblemVec.cc
index a5244b25..65233cd7 100644
--- a/AMDiS/src/ProblemVec.cc
+++ b/AMDiS/src/ProblemVec.cc
@@ -619,8 +619,9 @@ namespace AMDiS {
{
FUNCNAME("ProblemVec::buildAfterCoarsen()");
- clock_t first = clock();
+ // buildAfterCoarsen_sebastianMode(adaptInfo, flag);
+ clock_t first = clock();
#ifdef _OPENMP
double wtime = omp_get_wtime();
#endif
@@ -648,13 +649,6 @@ namespace AMDiS {
for (int i = 0; i < nComponents; i++) {
-#if 0
- std::vector<int> nnzInRow;
- int overallNnz, averageNnz;
- componentSpaces[i]->getMesh()->computeMatrixFillin(componentSpaces[i],
- nnzInRow, overallNnz, averageNnz);
-#endif
-
MSG("%d DOFs for %s\n",
componentSpaces[i]->getAdmin()->getUsedSize(),
componentSpaces[i]->getName().c_str());
@@ -795,6 +789,201 @@ namespace AMDiS {
#endif
}
+ void ProblemVec::buildAfterCoarsen_sebastianMode(AdaptInfo *adaptInfo, Flag flag)
+ {
+ FUNCNAME("ProblemVec::buildAfterCoarsen()");
+
+ clock_t first = clock();
+#ifdef _OPENMP
+ double wtime = omp_get_wtime();
+#endif
+
+ for (int i = 0; i < static_cast<int>(meshes.size()); i++)
+ meshes[i]->dofCompress();
+
+ Flag assembleFlag =
+ flag |
+ (*systemMatrix)[0][0]->getAssembleFlag() |
+ rhs->getDOFVector(0)->getAssembleFlag() |
+ Mesh::CALL_LEAF_EL |
+ Mesh::FILL_COORDS |
+ Mesh::FILL_DET |
+ Mesh::FILL_GRD_LAMBDA |
+ Mesh::FILL_NEIGH;
+
+ if (useGetBound)
+ assembleFlag |= Mesh::FILL_BOUND;
+
+ traverseInfo.updateStatus();
+
+ // Used to calculate the overall number of non zero entries.
+ int nnz = 0;
+
+
+ /// === INITIALIZE ===
+
+ for (int i = 0; i < nComponents; i++) {
+ MSG("%d DOFs for %s\n",
+ componentSpaces[i]->getAdmin()->getUsedSize(),
+ componentSpaces[i]->getName().c_str());
+
+ rhs->getDOFVector(i)->set(0.0);
+
+ for (int j = 0; j < nComponents; j++) {
+ // Only if this variable is true, the current matrix will be assembled.
+ bool assembleMatrix = true;
+ // The DOFMatrix which should be assembled (or not, if assembleMatrix
+ // will be set to false).
+ DOFMatrix *matrix = (*systemMatrix)[i][j];
+
+ if (matrix)
+ matrix->calculateNnz();
+
+ // If the matrix was assembled before and it is marked to be assembled
+ // only once, it will not be assembled.
+ if (assembleMatrixOnlyOnce[i][j] && assembledMatrix[i][j]) {
+ assembleMatrix = false;
+ } else if (matrix) {
+ matrix->getBaseMatrix().
+ change_dim(componentSpaces[i]->getAdmin()->getUsedSize(),
+ componentSpaces[j]->getAdmin()->getUsedSize());
+
+ set_to_zero(matrix->getBaseMatrix());
+ }
+
+ // If there is no DOFMatrix, e.g., if it is completly 0, do not assemble.
+ if (!matrix || !assembleMatrix)
+ assembleMatrix = false;
+
+ // If the matrix should not be assembled, the rhs vector has to be considered.
+ // This will be only done, if i == j. So, if both is not true, we can jump
+ // to the next matrix.
+ if (!assembleMatrix && i != j) {
+ if (matrix)
+ nnz += matrix->getBaseMatrix().nnz();
+
+ continue;
+ }
+
+ if (assembleMatrix && matrix->getBoundaryManager())
+ matrix->getBoundaryManager()->initMatrix(matrix);
+
+ if (matrix && assembleMatrix)
+ matrix->startInsertion(matrix->getNnz());
+ }
+ }
+
+
+ // === TRAVERSE ===
+
+ Mesh *mesh = componentMeshes[0];
+ const FiniteElemSpace *feSpace = componentSpaces[0];
+ const BasisFunction *basisFcts = feSpace->getBasisFcts();
+ ElementMatrix elMat(basisFcts->getNumber(), basisFcts->getNumber());
+ ElementMatrix tmpElMat(elMat);
+ ElementVector elVec(basisFcts->getNumber());
+ ElementVector tmpElVec(elVec);
+ TraverseStack stack;
+ BoundaryType *bound =
+ useGetBound ? new BoundaryType[basisFcts->getNumber()] : NULL;
+ ElInfo *elInfo = stack.traverseFirst(mesh, -1, assembleFlag);
+ while (elInfo) {
+ if (useGetBound)
+ basisFcts->getBound(elInfo, bound);
+
+ for (std::map<Operator*, std::vector<OperatorPos> >::iterator opIt = operators.begin();
+ opIt != operators.end(); ++opIt) {
+ if (opIt->first->getNeedDualTraverse() == true)
+ continue;
+
+ if (opFlags[opIt->first].isSet(Operator::MATRIX_OPERATOR)) {
+ set_to_zero(elMat);
+ opIt->first->getElementMatrix(elInfo, elMat, 1.0);
+ }
+ if (opFlags[opIt->first].isSet(Operator::VECTOR_OPERATOR)) {
+ set_to_zero(elVec);
+ opIt->first->getElementVector(elInfo, elVec, 1.0);
+ }
+
+ for (std::vector<OperatorPos>::iterator posIt = opIt->second.begin();
+ posIt != opIt->second.end(); ++posIt) {
+
+ if (posIt->operatorType.isSet(Operator::MATRIX_OPERATOR)) {
+ if (*(posIt->factor) == 1.0) {
+ (*systemMatrix)[posIt->row][posIt->col]->addElementMatrix(elMat, bound, elInfo, NULL);
+ } else {
+ tmpElMat = *(posIt->factor) * elMat;
+ (*systemMatrix)[posIt->row][posIt->col]->addElementMatrix(tmpElMat, bound, elInfo, NULL);
+ }
+ }
+
+ if (posIt->operatorType.isSet(Operator::VECTOR_OPERATOR)) {
+ if (*(posIt->factor) == 1.0) {
+ rhs->getDOFVector(posIt->row)->addElementVector(1.0, elVec, bound, elInfo);
+ } else {
+ tmpElVec = *(posIt->factor) * elVec;
+ rhs->getDOFVector(posIt->row)->addElementVector(1.0, tmpElVec, bound, elInfo);
+ }
+ }
+ }
+ }
+
+ elInfo = stack.traverseNext(elInfo);
+ }
+
+ if (useGetBound)
+ delete [] bound;
+
+ // === FINELIZE ===
+
+ for (int i = 0; i < nComponents; i++) {
+ for (int j = 0; j < nComponents; j++) {
+ bool assembleMatrix = true;
+ DOFMatrix *matrix = (*systemMatrix)[i][j];
+
+ if (assembleMatrixOnlyOnce[i][j] && assembledMatrix[i][j])
+ assembleMatrix = false;
+ if (!matrix || !assembleMatrix)
+ assembleMatrix = false;
+ if (!assembleMatrix && i != j)
+ continue;
+
+ assembledMatrix[i][j] = true;
+
+ if (assembleMatrix) {
+ matrix->removeRowsWithDBC(matrix->getApplyDBCs());
+ matrix->finishInsertion();
+ }
+ if (assembleMatrix && matrix->getBoundaryManager())
+ matrix->getBoundaryManager()->exitMatrix(matrix);
+ if (matrix)
+ nnz += matrix->getBaseMatrix().nnz();
+ }
+
+ assembleBoundaryConditions(rhs->getDOFVector(i),
+ solution->getDOFVector(i),
+ componentMeshes[i],
+ assembleFlag);
+ }
+
+ solverMatrix.setMatrix(*systemMatrix);
+
+ createPrecon();
+
+ INFO(info, 8)("fillin of assembled matrix: %d\n", nnz);
+
+#ifdef _OPENMP
+ INFO(info, 8)("buildAfterCoarsen needed %.5f seconds system time / %.5f seconds wallclock time\n",
+ TIME_USED(first, clock()), omp_get_wtime() - wtime);
+#else
+ INFO(info, 8)("buildAfterCoarsen needed %.5f seconds\n",
+ TIME_USED(first, clock()));
+#endif
+
+ exit(0);
+
+ }
+
void ProblemVec::createPrecon()
{
std::string preconType("no");
@@ -884,6 +1073,10 @@ namespace AMDiS {
break;
}
}
+
+ OperatorPos opPos = {i, j, factor, estFactor, Operator::MATRIX_OPERATOR};
+ operators[op].push_back(opPos);
+ opFlags[op].setFlag(Operator::MATRIX_OPERATOR);
}
void ProblemVec::addVectorOperator(Operator *op, int i,
@@ -902,6 +1095,10 @@ namespace AMDiS {
break;
}
}
+
+ OperatorPos opPos = {i, -1, factor, estFactor, Operator::VECTOR_OPERATOR};
+ operators[op].push_back(opPos);
+ opFlags[op].setFlag(Operator::VECTOR_OPERATOR);
}
void ProblemVec::addDirichletBC(BoundaryType type, int row, int col,
diff --git a/AMDiS/src/ProblemVec.h b/AMDiS/src/ProblemVec.h
index 4f6266e4..c6740d38 100644
--- a/AMDiS/src/ProblemVec.h
+++ b/AMDiS/src/ProblemVec.h
@@ -37,6 +37,13 @@
namespace AMDiS {
+ struct OperatorPos
+ {
+ int row, col;
+ double *factor, *estFactor;
+ Flag operatorType;
+ };
+
class ProblemVec : public ProblemStatBase,
public StandardProblemIteration
{
@@ -171,6 +178,8 @@ namespace AMDiS {
bool assembleMatrix = true,
bool assembleVector = true);
+ void buildAfterCoarsen_sebastianMode(AdaptInfo *adaptInfo, Flag flag);
+
void createPrecon();
@@ -594,6 +603,10 @@ namespace AMDiS {
* defined by \ref exactSolutionFcts.
*/
bool computeExactError;
+
+ std::map<Operator*, std::vector<OperatorPos> > operators;
+
+ std::map<Operator*, Flag> opFlags;
};
}
diff --git a/AMDiS/src/SecondOrderAssembler.cc b/AMDiS/src/SecondOrderAssembler.cc
index dbe5ad65..0d6fd931 100644
--- a/AMDiS/src/SecondOrderAssembler.cc
+++ b/AMDiS/src/SecondOrderAssembler.cc
@@ -210,33 +210,32 @@ namespace AMDiS {
for (int i = 0; i < static_cast<int>(terms[myRank].size()); i++)
(static_cast<SecondOrderTerm*>(terms[myRank][i]))->getLALt(elInfo, nPoints, LALt);
+ VectorOfFixVecs< DimVec<double> > *grdPsi, *grdPhi;
+
if (symmetric) {
// === Symmetric assembling. ===
-
TEST_EXIT_DBG(nCol == nRow)("nCol != nRow, but symmetric assembling!\n");
for (int iq = 0; iq < nPoints; iq++) {
(*LALt[iq]) *= elInfo->getDet();
+ grdPsi = psiFast->getGradient(iq);
+ grdPhi = phiFast->getGradient(iq);
for (int i = 0; i < nCol; i++) {
- amv(quadrature->getWeight(iq),
- (*LALt[iq]), phiFast->getGradient(iq, i), dimVec);
+ amv(quadrature->getWeight(iq), (*LALt[iq]), (*grdPhi)[i], dimVec);
- mat[i][i] += (psiFast->getGradient(iq, i) * dimVec);
+ mat[i][i] += (*grdPsi)[i] * dimVec;
for (int j = i + 1; j < nRow; j++) {
- double tmp = (phiFast->getGradient(iq, j) * dimVec);
+ double tmp = (*grdPhi)[j] * dimVec;
mat[i][j] += tmp;
mat[j][i] += tmp;
}
}
}
} else {
- VectorOfFixVecs< DimVec<double> > *grdPsi, *grdPhi;
-
// === Non symmetric assembling. ===
for (int iq = 0; iq < nPoints; iq++) {
(*LALt[iq]) *= elInfo->getDet();
-
grdPsi = psiFast->getGradient(iq);
grdPhi = phiFast->getGradient(iq);
diff --git a/AMDiS/src/ZeroOrderAssembler.cc b/AMDiS/src/ZeroOrderAssembler.cc
index c7975367..fc6b4000 100644
--- a/AMDiS/src/ZeroOrderAssembler.cc
+++ b/AMDiS/src/ZeroOrderAssembler.cc
@@ -155,7 +155,9 @@ namespace AMDiS {
FastQuadZOA::FastQuadZOA(Operator *op, Assembler *assembler, Quadrature *quad)
: ZeroOrderAssembler(op, assembler, quad, true)
- {}
+ {
+ tmpC.resize(omp_get_overall_max_threads());
+ }
void FastQuadZOA::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix& mat)
{
@@ -175,24 +177,28 @@ namespace AMDiS {
}
}
- std::vector<double> c(nPoints, 0.0);
- std::vector<OperatorTerm*>::iterator termIt;
- for (termIt = terms[myRank].begin(); termIt != terms[myRank].end(); ++termIt) {
+ std::vector<double> &c = tmpC[myRank];
+ if (c.size() < static_cast<int>(nPoints))
+ c.resize(nPoints);
+ for (int i = 0; i < static_cast<int>(nPoints); i++)
+ c[i] = 0.0;
+
+ for (std::vector<OperatorTerm*>::iterator termIt = terms[myRank].begin();
+ termIt != terms[myRank].end(); ++termIt)
(static_cast<ZeroOrderTerm*>((*termIt)))->getC(elInfo, nPoints, c);
- }
if (symmetric) {
TEST_EXIT_DBG(nCol == nRow)("nCol != nRow, but symmetric assembling!\n");
for (int iq = 0; iq < nPoints; iq++) {
- c[iq] *= elInfo->getDet();
-
+ c[iq] *= elInfo->getDet() * quadrature->getWeight(iq);
const double *psi = psiFast->getPhi(iq);
const double *phi = phiFast->getPhi(iq);
+
for (int i = 0; i < nRow; i++) {
- mat[i][i] += quadrature->getWeight(iq) * c[iq] * psi[i] * phi[i];
+ mat[i][i] += c[iq] * psi[i] * phi[i];
for (int j = i + 1; j < nCol; j++) {
- double val = quadrature->getWeight(iq) * c[iq] * psi[i] * phi[j];
+ double val = c[iq] * psi[i] * phi[j];
mat[i][j] += val;
mat[j][i] += val;
}
@@ -200,15 +206,13 @@ namespace AMDiS {
}
} else { /* non symmetric assembling */
for (int iq = 0; iq < nPoints; iq++) {
- c[iq] *= elInfo->getDet();
+ c[iq] *= elInfo->getDet() * quadrature->getWeight(iq);
const double *psi = psiFast->getPhi(iq);
const double *phi = phiFast->getPhi(iq);
- for (int i = 0; i < nRow; i++) {
- for (int j = 0; j < nCol; j++) {
- mat[i][j] += quadrature->getWeight(iq) * c[iq] * psi[i] * phi[j];
- }
- }
+ for (int i = 0; i < nRow; i++)
+ for (int j = 0; j < nCol; j++)
+ mat[i][j] += c[iq] * psi[i] * phi[j];
}
}
}
diff --git a/AMDiS/src/ZeroOrderAssembler.h b/AMDiS/src/ZeroOrderAssembler.h
index 0b3e63fa..8b349c1a 100644
--- a/AMDiS/src/ZeroOrderAssembler.h
+++ b/AMDiS/src/ZeroOrderAssembler.h
@@ -85,6 +85,9 @@ namespace AMDiS {
/// Implements SubAssembler::calculateElementVector().
void calculateElementVector(const ElInfo *elInfo, ElementVector& vec);
+
+ protected:
+ std::vector<std::vector<double> > tmpC;
};
--
GitLab