From 3b97cf2bb9d8acf118a548eb8d072d5deef5f6a2 Mon Sep 17 00:00:00 2001
From: Thomas Witkowski <thomas.witkowski@gmx.de>
Date: Tue, 23 Sep 2008 15:10:22 +0000
Subject: [PATCH] * Parallel assemblage * BoundaryManager.hh deleted, was not
included in any file
---
AMDiS/src/BasisFunction.h | 23 +++---
AMDiS/src/BoundaryManager.cc | 38 +++++++--
AMDiS/src/BoundaryManager.h | 22 +++--
AMDiS/src/BoundaryManager.hh | 130 ------------------------------
AMDiS/src/DOFMatrix.cc | 60 ++++++++++++--
AMDiS/src/DOFMatrix.h | 17 ++++
AMDiS/src/DOFVector.hh | 2 +-
AMDiS/src/Global.h | 7 --
AMDiS/src/Lagrange.cc | 27 +------
AMDiS/src/Lagrange.h | 2 +-
AMDiS/src/MultiGridSolver.cc | 14 ++--
AMDiS/src/NonLinUpdater.cc | 45 +++++++----
AMDiS/src/ProblemScal.cc | 9 ++-
AMDiS/src/ProblemVec.cc | 113 +++++++++++++++++++-------
AMDiS/src/SecondOrderAssembler.cc | 2 +-
AMDiS/src/Traverse.h | 5 +-
16 files changed, 268 insertions(+), 248 deletions(-)
delete mode 100644 AMDiS/src/BoundaryManager.hh
diff --git a/AMDiS/src/BasisFunction.h b/AMDiS/src/BasisFunction.h
index b40ee744..e11450a3 100644
--- a/AMDiS/src/BasisFunction.h
+++ b/AMDiS/src/BasisFunction.h
@@ -177,18 +177,17 @@ namespace AMDiS {
DegreeOfFreedom*) const = 0;
/** \brief
- * Pointer to a function which fills a vector with the boundary types of the
- * basis functions;
- * getBound(el info, bound) returns a pointer to this vector of length
- * \ref nBasFcts where the i-th entry is the boundary type of the i-th basis
- * function; bound may be a pointer to a vector which has to be filled
- * (compare the dof argument of \ref getDOFIindices);
- * such a function needs boundary information; thus, all routines using this
- * function on the elements need the FILL_BOUND flag during mesh traversal;
- */
- virtual const BoundaryType* getBound(const ElInfo*, BoundaryType *) const {
- return NULL;
- };
+ * The second argument 'bound' has to be a pointer to a vector which has
+ * to be filled. Its length is \ref nBasFcts (the number of basis functions
+ * in the used finite element space). After calling this function, the i-th
+ * entry of the array is the boundary type of the i-th basis function of this
+ * element.
+ *
+ * This function needs boundary information within the ElInfo object; thus,
+ * all routines using this function on the elements need the FILL_BOUND
+ * flag during mesh traversal;
+ */
+ virtual void getBound(const ElInfo*, BoundaryType *) const {};
/** \brief
* Returns \ref degree of BasisFunction
diff --git a/AMDiS/src/BoundaryManager.cc b/AMDiS/src/BoundaryManager.cc
index 85fa9016..8ae905e8 100644
--- a/AMDiS/src/BoundaryManager.cc
+++ b/AMDiS/src/BoundaryManager.cc
@@ -1,19 +1,45 @@
#include "FiniteElemSpace.h"
-//#include "BoundaryCondition.h"
#include "BoundaryManager.h"
#include "DOFIndexed.h"
#include "DOFVector.h"
#include "Traverse.h"
#include "BasisFunction.h"
#include "ElInfo.h"
+#include "OpenMP.h"
namespace AMDiS {
+ BoundaryManager::BoundaryManager(const FiniteElemSpace *feSpace)
+ {
+ localBounds.resize(omp_get_max_threads());
+ allocatedMemoryLocalBounds = feSpace->getBasisFcts()->getNumber();
+ for (int i = 0; i < static_cast<int>(localBounds.size()); i++) {
+ localBounds[i] = GET_MEMORY(BoundaryType, allocatedMemoryLocalBounds);
+ }
+ }
+
+ BoundaryManager::BoundaryManager(BoundaryManager &bm)
+ {
+ localBCs = bm.localBCs;
+ allocatedMemoryLocalBounds = bm.allocatedMemoryLocalBounds;
+ localBounds.resize(bm.localBounds.size());
+ for (int i = 0; i < static_cast<int>(localBounds.size()); i++) {
+ localBounds[i] = GET_MEMORY(BoundaryType, allocatedMemoryLocalBounds);
+ }
+ }
+
+ BoundaryManager::~BoundaryManager()
+ {
+ for (int i = 0; i < static_cast<int>(localBounds.size()); i++) {
+ FREE_MEMORY(localBounds[i], BoundaryType, allocatedMemoryLocalBounds);
+ }
+ }
+
double BoundaryManager::boundResidual(ElInfo *elInfo,
DOFMatrix *matrix,
const DOFVectorBase<double> *dv)
{
- double result = 0;
+ double result = 0.0;
std::map<BoundaryType, BoundaryCondition*>::iterator it;
for (it = localBCs.begin(); it != localBCs.end(); ++it) {
if ((*it).second)
@@ -37,7 +63,7 @@ namespace AMDiS {
if (localBCs.size() > 0) {
// get boundaries of all DOFs
- BoundaryType *localBound = GET_MEMORY(BoundaryType, nBasFcts);
+ BoundaryType *localBound = localBounds[omp_get_thread_num()];
basisFcts->getBound(elInfo, localBound);
// get dof indices
@@ -60,8 +86,6 @@ namespace AMDiS {
}
}
}
-
- FREE_MEMORY(localBound, BoundaryType, nBasFcts);
}
}
@@ -79,7 +103,9 @@ namespace AMDiS {
if (localBCs.size() > 0) {
// get boundaries of all DOFs
- const BoundaryType *localBound = basisFcts->getBound(elInfo, NULL);
+ BoundaryType *localBound = localBounds[omp_get_thread_num()];
+ basisFcts->getBound(elInfo, localBound);
+
// get dof indices
basisFcts->getLocalIndicesVec(elInfo->getElement(), admin, &dofIndices);
diff --git a/AMDiS/src/BoundaryManager.h b/AMDiS/src/BoundaryManager.h
index a593b89b..f90ac30d 100644
--- a/AMDiS/src/BoundaryManager.h
+++ b/AMDiS/src/BoundaryManager.h
@@ -50,6 +50,12 @@ namespace AMDiS {
public:
MEMORY_MANAGED(BoundaryManager);
+ BoundaryManager(const FiniteElemSpace *feSpace);
+
+ BoundaryManager(BoundaryManager &bm);
+
+ ~BoundaryManager();
+
/** \brief
* Adds a local boundary condition to the list of managed conditions.
*/
@@ -92,11 +98,6 @@ namespace AMDiS {
return localBCs[type];
};
- // /** \brief
- // * Fills boundary types in DOFVectorBase bound.
- // */
- // void fillGlobalBoundVector(DOFVectorBase<BoundaryType> *bound);
-
const std::map<BoundaryType, BoundaryCondition*>& getBoundaryConditionMap() {
return localBCs;
};
@@ -110,6 +111,17 @@ namespace AMDiS {
* Map of managed local boundary conditions.
*/
std::map<BoundaryType, BoundaryCondition*> localBCs;
+
+ /** \brief
+ * Temporary thread-safe variable for functions fillBoundaryconditions.
+ */
+ std::vector<BoundaryType*> localBounds;
+
+ /** \brief
+ * Stores the number of byte that were allocated in the constructor for
+ * each localBounds value. Is used to free the memory in the destructor.
+ */
+ int allocatedMemoryLocalBounds;
};
}
diff --git a/AMDiS/src/BoundaryManager.hh b/AMDiS/src/BoundaryManager.hh
deleted file mode 100644
index 51333d35..00000000
--- a/AMDiS/src/BoundaryManager.hh
+++ /dev/null
@@ -1,130 +0,0 @@
-#include "FiniteElemSpace.h"
-#include "DOFIndexed.h"
-#include "DOFVector.h"
-#include "Traverse.h"
-#include "BasisFunction.h"
-#include "ElInfo.h"
-#include "BoundaryCondition.h"
-#include <list>
-
-namespace AMDiS {
-
- template<typename T>
- double BoundaryManager<T>::boundResidual(ElInfo *elInfo, Estimator<T> *estimator)
- {
- double result = 0;
- typename std::list<LocalBC<T>*>::iterator it;
- for(it = localBCs.begin(); it != localBCs.end(); ++it) {
- result += (*it)->boundResidual(elInfo, estimator);
- }
- return result;
- }
-
- template<typename T>
- void BoundaryManager<T>::fillGlobalBoundVector(DOFVector<BoundaryType> *globalBound)
- {
- int i;
- TraverseStack stack;
- const FiniteElemSpace *feSpace = globalBound->getFESpace();
-
- const BasisFunction *basisFcts = feSpace->getBasisFcts();
- int nBasFcts = basisFcts->getNumber();
- Mesh *mesh = feSpace->getMesh();
- DOFAdmin *admin = feSpace->getAdmin();
-
- const BoundaryType *localBound = NULL;
- const DegreeOfFreedom *dofIndices = NULL;
-
- ElInfo *elInfo =
- stack.traverseFirst(mesh, -1,
- Mesh::CALL_LEAF_EL |
- Mesh::FILL_BOUND |
- Mesh::FILL_COORDS);
-
- // for all elements ...
- while(elInfo) {
- // get boundaries of all DOFs
- localBound = basisFcts->getBound(elInfo, NULL);
-
- // get dof indices
- dofIndices = basisFcts->getLocalIndices(elInfo->getElement(),
- admin, NULL);
-
- // for all DOFs
- for(i=0; i < nBasFcts; i++) {
- (*globalBound)[dofIndices[i]] = localBound[i];
- }
-
- elInfo = stack.traverseNext(elInfo);
- }
- }
-
- template<typename T>
- void BoundaryManager<T>::fillLocalBCs(ElInfo *elInfo,
- DOFVector<T> *vec)
- {
- int i;
-
- // ===== fill local conditions ==============================================
- const FiniteElemSpace *feSpace = vec->getFESpace();
- const BoundaryType *localBound = NULL;
- const DegreeOfFreedom *dofIndices = NULL;
- Mesh *mesh = feSpace->getMesh();
- const BasisFunction *basisFcts = feSpace->getBasisFcts();
- int nBasFcts = basisFcts->getNumber();
- LocalBC<T> *localBC;
- DOFAdmin *admin = feSpace->getAdmin();
-
- typename std::list<LocalBC<T>*>::iterator it;
-
- if(localBCs.size() > 0) {
-
- // get boundaries of all DOFs
- localBound = basisFcts->getBound(elInfo, NULL);
-
- // get dof indices
- dofIndices = basisFcts->getLocalIndices(elInfo->getElement(),
- admin, NULL);
-
- // apply boundary conditions
- for(it = localBCs.begin(); it != localBCs.end(); ++it) {
- (*it)->fillLocalBC(vec, elInfo, dofIndices, localBound, nBasFcts);
- }
- }
- }
-
- template<typename T>
- void BoundaryManager<T>::fillLocalBCs(ElInfo *elInfo,
- DOFMatrix *mat)
- {
- int i;
-
- // ===== fill local conditions ==============================================
- const FiniteElemSpace *feSpace = mat->getRowFESpace();
- const BoundaryType *localBound = NULL;
- const DegreeOfFreedom *dofIndices = NULL;
- Mesh *mesh = feSpace->getMesh();
- const BasisFunction *basisFcts = feSpace->getBasisFcts();
- int nBasFcts = basisFcts->getNumber();
- LocalBC<T> *localBC;
- DOFAdmin *admin = feSpace->getAdmin();
-
- typename std::list<LocalBC<T>*>::iterator it;
-
- if(localBCs.size() > 0) {
-
- // get boundaries of all DOFs
- localBound = basisFcts->getBound(elInfo, NULL);
-
- // get dof indices
- dofIndices = basisFcts->getLocalIndices(elInfo->getElement(),
- admin, NULL);
-
- // apply boundary conditions
- for(it = localBCs.begin(); it != localBCs.end(); ++it) {
- (*it)->fillLocalBC(mat, elInfo, dofIndices, localBound, nBasFcts);
- }
- }
- }
-
-}
diff --git a/AMDiS/src/DOFMatrix.cc b/AMDiS/src/DOFMatrix.cc
index d1358583..2398e57f 100644
--- a/AMDiS/src/DOFMatrix.cc
+++ b/AMDiS/src/DOFMatrix.cc
@@ -42,9 +42,11 @@ namespace AMDiS {
if (rowFESpace && rowFESpace->getAdmin())
(const_cast<DOFAdmin*>(rowFESpace->getAdmin()))->addDOFIndexed(this);
- boundaryManager = NEW BoundaryManager;
+ boundaryManager = NEW BoundaryManager(rowFESpace_);
elementMatrix = NEW ElementMatrix(rowFESpace->getBasisFcts()->getNumber(),
colFESpace->getBasisFcts()->getNumber());
+
+ applyDBCs.clear();
}
DOFMatrix::DOFMatrix(const DOFMatrix& rhs)
@@ -215,7 +217,7 @@ namespace AMDiS {
operatorFactor = rhs.operatorFactor;
matrix = rhs.matrix;
if (rhs.boundaryManager) {
- boundaryManager = new BoundaryManager(*rhs.boundaryManager);
+ boundaryManager = NEW BoundaryManager(*rhs.boundaryManager);
} else {
boundaryManager = NULL;
}
@@ -248,14 +250,15 @@ namespace AMDiS {
bound ? boundaryManager->getBoundaryCondition(bound[i]) : NULL;
if (condition && condition->isDirichlet()) {
- MatrixRow *matrixRow = &(matrix[row]);
+ /* MatrixRow *matrixRow = &(matrix[row]);
if (coupleMatrix) {
matrixRow->resize(0);
} else {
matrixRow->resize(1);
((*matrixRow)[0]).col = row;
- ((*matrixRow)[0]).entry = 1.0;
- }
+ ((*matrixRow)[0]).entry = 1.0;
+ } */
+ applyDBCs.insert(static_cast<int>(row));
} else {
for (int j = 0; j < nCol; j++) { // for all columns
addSparseDOFEntry(sign, row, elMat.colIndices[j],
@@ -653,6 +656,21 @@ namespace AMDiS {
}
}
+ void DOFMatrix::removeRowsWithDBC(std::set<int> *rows)
+ {
+ for (std::set<int>::iterator it = rows->begin();
+ it != rows->end();
+ ++it) {
+ if (coupleMatrix) {
+ matrix[*it].resize(0);
+ } else {
+ matrix[*it].resize(1);
+ matrix[*it][0].col = *it;
+ matrix[*it][0].entry = 1.0;
+ }
+ }
+ }
+
void DOFMatrix::createPictureFile(const char* filename, int dim)
{
png_structp png_ptr = png_create_write_struct(PNG_LIBPNG_VER_STRING, NULL, NULL, NULL);
@@ -814,4 +832,36 @@ namespace AMDiS {
}
}
+ void addDOFMatrix(DOFMatrix *result, const DOFMatrix *a)
+ {
+ DOFMatrix::Iterator resultIterator(result, USED_DOFS);
+ DOFMatrix::Iterator aIterator(const_cast<DOFMatrix*>(a), USED_DOFS);
+
+ for (resultIterator.reset(), aIterator.reset();
+ !aIterator.end();
+ ++resultIterator, ++aIterator) {
+ std::vector<MatEntry>::iterator resultRowIt;
+ std::vector<MatEntry>::const_iterator aRowIt;
+
+ for (aRowIt = aIterator->begin();
+ aRowIt != aIterator->end();
+ ++aRowIt) {
+ bool added = false;
+ for (resultRowIt = resultIterator->begin();
+ resultRowIt != resultIterator->end();
+ ++resultRowIt) {
+ if (aRowIt->col == resultRowIt->col) {
+ resultRowIt->entry += aRowIt->entry;
+
+ added = true;
+ break;
+ }
+ }
+
+ if (!added) {
+ resultIterator->push_back(*aRowIt);
+ }
+ }
+ }
+ }
}
diff --git a/AMDiS/src/DOFMatrix.h b/AMDiS/src/DOFMatrix.h
index 77b51fe2..61386b3f 100644
--- a/AMDiS/src/DOFMatrix.h
+++ b/AMDiS/src/DOFMatrix.h
@@ -27,6 +27,7 @@
// ============================================================================
#include <vector>
+#include <set>
#include <memory>
#include <list>
#include "Global.h"
@@ -581,6 +582,8 @@ namespace AMDiS {
void addRow(std::vector<MatEntry> row);
+ void removeRowsWithDBC(std::set<int> *rows);
+
/** \brief
* Prints \ref matrix to stdout
*/
@@ -619,6 +622,10 @@ namespace AMDiS {
return boundaryManager;
};
+ std::set<int>* getApplyDBCs() {
+ return &applyDBCs;
+ }
+
inline void setBoundaryManager(BoundaryManager *bm) {
boundaryManager = bm;
};
@@ -720,6 +727,8 @@ namespace AMDiS {
*/
ElementMatrix *elementMatrix;
+ std::set<int> applyDBCs;
+
friend class DOFAdmin;
friend class DOFVector<double>;
friend class DOFVector<unsigned char>;
@@ -739,8 +748,16 @@ namespace AMDiS {
double max(std::vector<MatEntry> *row);
+ /** \brief
+ * Addes two matrices. The result matrix is overwritten.
+ */
void addDOFMatrix(DOFMatrix *result, const DOFMatrix *a, const DOFMatrix *b);
+ /** \brief
+ * Addes matrix a to matrix result.
+ */
+ void addDOFMatrix(DOFMatrix *result, const DOFMatrix *a);
+
}
#endif // AMDIS_DOFMATRIX_H
diff --git a/AMDiS/src/DOFVector.hh b/AMDiS/src/DOFVector.hh
index 3fe022e4..e58c6449 100644
--- a/AMDiS/src/DOFVector.hh
+++ b/AMDiS/src/DOFVector.hh
@@ -70,7 +70,7 @@ namespace AMDiS {
(feSpace->getAdmin())->addDOFIndexed(this);
}
- this->boundaryManager = NEW BoundaryManager;
+ this->boundaryManager = NEW BoundaryManager(f);
}
template<typename T>
diff --git a/AMDiS/src/Global.h b/AMDiS/src/Global.h
index b9aacee3..589bf359 100644
--- a/AMDiS/src/Global.h
+++ b/AMDiS/src/Global.h
@@ -88,13 +88,6 @@ namespace AMDiS {
*/
#define DIM_OF_INDEX(ind, dim) ((static_cast<int>(ind) == 0) ? dim : static_cast<int>(ind) - 1)
-
- /** \brief
- * Returns boundary->getBound() if boundary is not NULL. Returns INTERIOR
- * otherwise.
- */
- //#define GET_BOUND(boundary) ((boundary) ? (boundary)->getBound() : INTERIOR)
-
/** \brief
* Calculates factorial of i
*/
diff --git a/AMDiS/src/Lagrange.cc b/AMDiS/src/Lagrange.cc
index 81fdbd01..a86aa0da 100644
--- a/AMDiS/src/Lagrange.cc
+++ b/AMDiS/src/Lagrange.cc
@@ -765,27 +765,8 @@ namespace AMDiS {
return NULL;
}
- const BoundaryType* Lagrange::getBound(const ElInfo* el_info,
- BoundaryType* bound) const
+ void Lagrange::getBound(const ElInfo* el_info, BoundaryType* bound) const
{
- static BoundaryType *bound_vec = NULL;
- static int bound_vec_size = 0;
- BoundaryType *result;
-
- if (bound) {
- result = bound;
- } else {
- // realloc memory if neccessary
- if (bound_vec_size < nBasFcts) {
- if (bound_vec)
- FREE_MEMORY(bound_vec, BoundaryType, bound_vec_size);
- bound_vec = GET_MEMORY(BoundaryType, nBasFcts);
- bound_vec_size = nBasFcts;
- }
-
- result = bound_vec;
- }
-
el_info->testFlag(Mesh::FILL_BOUND);
// boundaries
@@ -801,7 +782,7 @@ namespace AMDiS {
boundaryType = el_info->getBoundary(j);
int kto = (*nDOF)[INDEX_OF_DIM(i, dim)];
for (int k = 0; k < kto; k++) {
- result[index++] = boundaryType;
+ bound[index++] = boundaryType;
}
}
offset -= Global::getGeo(INDEX_OF_DIM(i + 1, dim), dim);
@@ -809,12 +790,10 @@ namespace AMDiS {
// interior nodes in the center
for (int i = 0; i < (*nDOF)[CENTER]; i++) {
- result[index++] = INTERIOR;
+ bound[index++] = INTERIOR;
}
TEST_EXIT_DBG(index == nBasFcts)("found not enough boundarys\n");
-
- return result;
}
diff --git a/AMDiS/src/Lagrange.h b/AMDiS/src/Lagrange.h
index 5b6061e7..36a0165e 100644
--- a/AMDiS/src/Lagrange.h
+++ b/AMDiS/src/Lagrange.h
@@ -100,7 +100,7 @@ namespace AMDiS {
/** \brief
* Implements BasisFunction::getBound
*/
- const BoundaryType* getBound(const ElInfo*, BoundaryType *) const;
+ void getBound(const ElInfo*, BoundaryType *) const;
/** \brief
* Calculates the local vertex indices which are involved in evaluating
diff --git a/AMDiS/src/MultiGridSolver.cc b/AMDiS/src/MultiGridSolver.cc
index b9a69e33..689f67c2 100644
--- a/AMDiS/src/MultiGridSolver.cc
+++ b/AMDiS/src/MultiGridSolver.cc
@@ -157,7 +157,7 @@ namespace AMDiS {
std::vector<double*>::iterator factorIt;
std::vector<double*>::iterator factorBegin = systemMatrix_->getOperatorFactorBegin();
ElementMatrix *elementMatrix = NULL;
- const BoundaryType *bound;
+ BoundaryType *bound = GET_MEMORY(BoundaryType, numDOFs);
TraverseStack stack;
@@ -190,7 +190,7 @@ namespace AMDiS {
element = elInfo->getElement();
level = elInfo->getLevel();
basFcts->getLocalIndices(element, admin, dofs);
- bound = basFcts->getBound(elInfo, NULL);
+ basFcts->getBound(elInfo, bound);
elementMatrix = (*opBegin)->getAssembler(omp_get_thread_num())->initElementMatrix(elementMatrix, elInfo);
bool levelUsed = isMGLevel_[level];
@@ -274,6 +274,8 @@ namespace AMDiS {
} // if-else element is leaf
elInfo = stack.traverseNext(elInfo);
} // while elInfo
+
+ FREE_MEMORY(bound, BoundaryType, numDOFs);
if (!galerkin_) {
// Fill boundary conditions
@@ -1048,16 +1050,17 @@ namespace AMDiS {
elInfo = stack.traverseNext(elInfo);
}
+ BoundaryType *bound = GET_MEMORY(BoundaryType, numDOFs);
+
for (int col = 0; col < numComponents_; col++) {
DOFMatrix *dofMatrix = (*systemMatrix_)[row][col];
- if(dofMatrix) {
+ if (dofMatrix) {
std::vector<Operator*>::iterator it;
std::vector<Operator*>::iterator opBegin = dofMatrix->getOperatorsBegin();
std::vector<Operator*>::iterator opEnd = dofMatrix->getOperatorsEnd();
std::vector<double*>::iterator factorIt;
std::vector<double*>::iterator factorBegin = dofMatrix->getOperatorFactorBegin();
ElementMatrix *elementMatrix = NULL;
- const BoundaryType *bound;
// assemble level matrices and create dof sets
elInfo = stack.traverseFirst(mesh, -1,
@@ -1071,7 +1074,7 @@ namespace AMDiS {
element = elInfo->getElement();
level = elInfo->getLevel();
basFcts->getLocalIndices(element, admin, dofs);
- bound = basFcts->getBound(elInfo, NULL);
+ basFcts->getBound(elInfo, bound);
elementMatrix = (*opBegin)->getAssembler(omp_get_thread_num())->initElementMatrix(elementMatrix, elInfo);
bool levelUsed = isMGLevel_[level];
@@ -1163,6 +1166,7 @@ namespace AMDiS {
}
}
+ FREE_MEMORY(bound, BoundaryType, numDOFs);
FREE_MEMORY(dofs, DegreeOfFreedom, numDOFs);
}
diff --git a/AMDiS/src/NonLinUpdater.cc b/AMDiS/src/NonLinUpdater.cc
index e83a2fbb..9f866f11 100644
--- a/AMDiS/src/NonLinUpdater.cc
+++ b/AMDiS/src/NonLinUpdater.cc
@@ -18,8 +18,8 @@ namespace AMDiS {
FUNCNAME("NonLinUpdaterScal::update()");
TraverseStack stack;
- ElInfo *el_info;
- Flag fill_flag;
+ ElInfo *el_info;
+ Flag fill_flag;
TEST_EXIT(F || df)("neither F nor df are set\n");
@@ -43,9 +43,12 @@ namespace AMDiS {
Mesh::FILL_DET|
Mesh::FILL_GRD_LAMBDA;
+ BoundaryType *bound = GET_MEMORY(BoundaryType, feSpace->getBasisFcts()->getNumber());
+
el_info = stack.traverseFirst(feSpace->getMesh(), -1, fill_flag);
while (el_info) {
- const BoundaryType *bound = basFcts->getBound(el_info, NULL);
+ basFcts->getBound(el_info, bound);
+
if (updateDF) {
df->assemble(1.0, el_info, bound);
}
@@ -56,6 +59,8 @@ namespace AMDiS {
el_info = stack.traverseNext(el_info);
}
+
+ FREE_MEMORY(bound, BoundaryType, feSpace->getBasisFcts()->getNumber());
}
void NonLinUpdaterVec::update(bool updateDF,
@@ -63,19 +68,19 @@ namespace AMDiS {
{
FUNCNAME("NonLinUpdaterVec::update()");
- int i, j, size = df->getNumRows();
-
TraverseStack stack;
- ElInfo *el_info;
- Flag fill_flag;
+ ElInfo *el_info;
+ Flag fill_flag;
TEST_EXIT(F || df)("neither F nor df are set\n");
- if ((!updateDF) && (F==NULL)) {
+ if (!updateDF && (F == NULL)) {
WARNING("No RHS vector and no update for system matrix! Updater does nothing!\n");
return;
}
+ int size = df->getNumRows();
+
const FiniteElemSpace *feSpace =
F ?
F->getDOFVector(0)->getFESpace() :
@@ -83,9 +88,10 @@ namespace AMDiS {
if (updateDF) {
TEST_EXIT(df)("df not set but update tried!\n");
- for(i = 0; i < size; i++) {
- for(j = 0; j < size; j++) {
- if((*df)[i][j]) {
+
+ for (int i = 0; i < size; i++) {
+ for (int j = 0; j < size; j++) {
+ if ((*df)[i][j]) {
(*df)[i][j]->clear();
}
}
@@ -95,7 +101,7 @@ namespace AMDiS {
BasisFunction *basFcts = const_cast<BasisFunction*>(feSpace->getBasisFcts());
if (F) {
- for(i = 0; i < size; i++) {
+ for (int i = 0; i < size; i++) {
F->getDOFVector(i)->set(0.0);
}
}
@@ -107,13 +113,16 @@ namespace AMDiS {
Mesh::FILL_DET|
Mesh::FILL_GRD_LAMBDA;
+ BoundaryType *bound = GET_MEMORY(BoundaryType, basFcts->getNumber());
+
el_info = stack.traverseFirst(feSpace->getMesh(), -1, fill_flag);
while (el_info) {
- const BoundaryType *bound = basFcts->getBound(el_info, NULL);
+ basFcts->getBound(el_info, bound);
+
if (updateDF) {
- for(i = 0; i < size; i++) {
- for(j = 0; j < size; j++) {
- if((*df)[i][j]) {
+ for (int i = 0; i < size; i++) {
+ for (int j = 0; j < size; j++) {
+ if ((*df)[i][j]) {
(*df)[i][j]->assemble(1.0, el_info, bound);
}
}
@@ -121,13 +130,15 @@ namespace AMDiS {
}
if (F) {
- for(i = 0; i < size; i++) {
+ for(int i = 0; i < size; i++) {
F->getDOFVector(i)->assemble(1.0, el_info, bound);
}
}
el_info = stack.traverseNext(el_info);
}
+
+ FREE_MEMORY(bound, BoundaryType, basFcts->getNumber());
};
}
diff --git a/AMDiS/src/ProblemScal.cc b/AMDiS/src/ProblemScal.cc
index 8a46a287..3d65398e 100644
--- a/AMDiS/src/ProblemScal.cc
+++ b/AMDiS/src/ProblemScal.cc
@@ -629,10 +629,11 @@ namespace AMDiS {
int ProblemScal::buildAfterCoarsenFct(ElInfo *elInfo)
{
- const BoundaryType *bound;
+ BoundaryType *bound;
if (traversePtr_->getBoundUsed()) {
- bound = traversePtr_->getFESpace()->getBasisFcts()->getBound(elInfo, NULL);
+ bound = GET_MEMORY(BoundaryType, traversePtr_->getFESpace()->getBasisFcts()->getNumber());
+ traversePtr_->getFESpace()->getBasisFcts()->getBound(elInfo, bound);
} else {
bound = NULL;
}
@@ -640,6 +641,10 @@ namespace AMDiS {
traversePtr_->getSystemMatrix()->assemble(1.0, elInfo, bound);
traversePtr_->getRHS()->assemble(1.0, elInfo, bound);
+ if (traversePtr_->getBoundUsed()) {
+ FREE_MEMORY(bound, BoundaryTyppe, traversePtr_->getFESpace()->getBasisFcts()->getNumber());
+ }
+
return 0;
}
diff --git a/AMDiS/src/ProblemVec.cc b/AMDiS/src/ProblemVec.cc
index f11a7f1c..6a32fc54 100644
--- a/AMDiS/src/ProblemVec.cc
+++ b/AMDiS/src/ProblemVec.cc
@@ -22,6 +22,7 @@
#include "PeriodicBC.h"
#include "Lagrange.h"
#include "Flag.h"
+#include "TraverseParallel.h"
namespace AMDiS {
@@ -694,11 +695,7 @@ namespace AMDiS {
}
}
- int i;
-#ifdef _OPENMP
-#pragma omp parallel for
-#endif
- for (i = 0; i < nComponents; i++) {
+ for (int i = 0; i < nComponents; i++) {
for (int j = 0; j < nComponents; j++) {
// Only if this variable is true, the current matrix will be assembled.
bool assembleMatrix = true;
@@ -942,39 +939,95 @@ namespace AMDiS {
{
Mesh *mesh = feSpace->getMesh();
const BasisFunction *basisFcts = feSpace->getBasisFcts();
+ TraverseParallelStack stack;
- BoundaryType *bound = NULL;
- if (useGetBound_) {
- bound = GET_MEMORY(BoundaryType, basisFcts->getNumber());
- }
+#pragma omp parallel
+ {
+ BoundaryType *bound = useGetBound_ ? GET_MEMORY(BoundaryType, basisFcts->getNumber()) : NULL;
+
+ // Create for every thread its private matrix and vector, on that
+ // the thread will assemble its part of the mesh.
+ DOFMatrix *tmpMatrix = NULL;
+ DOFVector<double> *tmpVector = NULL;
- TraverseStack stack;
- ElInfo *elInfo = stack.traverseFirst(mesh, -1, assembleFlag);
-
- while (elInfo) {
- if (useGetBound_) {
- basisFcts->getBound(elInfo, bound);
- }
-
if (matrix) {
- matrix->assemble(1.0, elInfo, bound);
+ tmpMatrix = NEW DOFMatrix(matrix->getRowFESpace(),
+ matrix->getColFESpace(),
+ "tmp");
+
+ // Copy the global matrix to the private matrix, because we need the
+ // operators defined on the global matrix in the private one. Only the
+ // values have to be set to zero.
+ *tmpMatrix = *matrix;
+ tmpMatrix->clear();
+ }
+
+ if (vector) {
+ tmpVector = NEW DOFVector<double>(vector->getFESpace(), "tmp");
+
+ // Copy the global vector to the private vector, because we need the
+ // operatirs defined on the global vector in the private one. But set
+ // the values to zero of the private vector after copying.
+ *tmpVector = *vector;
+ tmpVector->set(0.0);
+ }
+
+ // Because we are using the parallel traverse stack, each thread will
+ // traverse only a part of the mesh.
+ ElInfo *elInfo = stack.traverseFirst(mesh, -1, assembleFlag);
+ while (elInfo) {
+ if (useGetBound_) {
+ basisFcts->getBound(elInfo, bound);
+ }
- if (matrix->getBoundaryManager()) {
- matrix->getBoundaryManager()->
- fillBoundaryConditions(elInfo, matrix);
- }
+ if (matrix) {
+ tmpMatrix->assemble(1.0, elInfo, bound);
+
+ // Take the matrix boundary manager from the public matrix,
+ // but assemble the boundary conditions on the thread private matrix.
+ if (matrix->getBoundaryManager()) {
+ matrix->getBoundaryManager()->
+ fillBoundaryConditions(elInfo, tmpMatrix);
+ }
+ }
+
+ if (vector) {
+ tmpVector->assemble(1.0, elInfo, bound);
+ }
+
+ elInfo = stack.traverseNext(elInfo);
}
-
+
+ // After mesh traverse, all thread have to added their private matrices and
+ // vectors to the global public matrix and public vector. Therefore, this is
+ // a critical section, which is allowed to be executed by on thread only at
+ // the same time.
+
+ if (matrix) {
+#pragma omp critical
+ {
+ addDOFMatrix(matrix, tmpMatrix);
+
+ // Remove rows corresponding to DOFs on a Dirichlet boundary.
+ matrix->removeRowsWithDBC(tmpMatrix->getApplyDBCs());
+ }
+
+ DELETE tmpMatrix;
+ }
+
if (vector) {
- vector->assemble(1.0, elInfo, bound);
+#pragma omp critical
+ *vector += *tmpVector;
+
+ DELETE tmpVector;
}
-
- elInfo = stack.traverseNext(elInfo);
- }
- if (useGetBound_) {
- FREE_MEMORY(bound, BoundaryType, basisFcts->getNumber());
- }
+ if (useGetBound_) {
+ FREE_MEMORY(bound, BoundaryType, basisFcts->getNumber());
+ }
+
+ } // pragma omp parallel
+
}
void ProblemVec::assembleOnDifMeshes(FiniteElemSpace *rowFeSpace, FiniteElemSpace *colFeSpace,
diff --git a/AMDiS/src/SecondOrderAssembler.cc b/AMDiS/src/SecondOrderAssembler.cc
index 538ef7ce..bd1dc76e 100644
--- a/AMDiS/src/SecondOrderAssembler.cc
+++ b/AMDiS/src/SecondOrderAssembler.cc
@@ -110,7 +110,7 @@ namespace AMDiS {
const double *values;
int myRank = omp_get_thread_num();
- DimMat<double> **LALt = tmpLALt[myRank];
+ DimMat<double> **LALt = tmpLALt[0];
DimMat<double> &tmpMat = *LALt[0];
tmpMat.set(0.0);
diff --git a/AMDiS/src/Traverse.h b/AMDiS/src/Traverse.h
index 32a2fc26..01090252 100644
--- a/AMDiS/src/Traverse.h
+++ b/AMDiS/src/Traverse.h
@@ -147,18 +147,19 @@ namespace AMDiS {
VectorOfFixVecs<DimVec<double> > *coords);
/** \brief
- *
+ * Is used for parallel mesh traverse.
*/
inline void setMyThreadId(int myThreadId) {
myThreadId_ = myThreadId;
}
/** \brief
- *
+ * Is used for parallel mesh traverse.
*/
inline void setMaxThreads(int maxThreads) {
maxThreads_ = maxThreads;
}
+
private:
/** \brief
* Enlargement of the stack
--
GitLab