* Parallel assemblage

* BoundaryManager.hh deleted, was not included in any file

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

AMDiS/src/BoundaryManager.cc

 #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);
 

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;
   };
 
 }

AMDiS/src/BoundaryManager.hh

-#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);
-      }
-    }
-  }
-
-}

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);
+	}
+      }	         
+    }      
+  }
 }

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

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>

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
    */

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;
   }
 
 

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

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);
       }
 

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;