Changed concept of local indices computation.

AMDiS/other/include/Makefile_AMDiS.mk

@@ -97,7 +97,7 @@ endif
 # ============================================================================
 
 ifeq ($(strip $(DEBUG)), 0)
-       CPPFLAGS += -O2
+       CPPFLAGS += -O3
 else
        CPPFLAGS += -g -O0
 endif

AMDiS/src/Assembler.cc

@@ -252,12 +252,11 @@ namespace AMDiS {
   {
     FUNCNAME("Assembler::matVecAssemble()");
 
-    Element *el = elInfo->getElement(); 
     double *uhOldLoc = new double[nRow];
 
-    operat->uhOld->getLocalVector(el, uhOldLoc);
+    operat->uhOld->getLocalVector(elInfo, uhOldLoc);
     
-    if (el != lastMatEl) {
+    if (elInfo->getElement() != lastMatEl) {
       set_to_zero(elementMatrix);
       calculateElementMatrix(elInfo, elementMatrix);
     }

AMDiS/src/BasisFunction.h

@@ -275,7 +275,7 @@ namespace AMDiS {
     virtual void coarseRestr(DOFVector<WorldVector<double> >*, RCNeighbourList*, int)
     {}
 
-    /// Returns local dof indices of the element for the given fe space.
+    /// Returns local dof indices of the element for the given DOF admin.
     virtual const DegreeOfFreedom *getLocalIndices(const Element *el,
 						   const DOFAdmin *admin,
 						   DegreeOfFreedom *dofPtr) const
@@ -283,17 +283,11 @@ namespace AMDiS {
       return NULL;
     }
 
-    inline void getLocalIndices(const Element *el,
-				const DOFAdmin *admin,
-				std::vector<DegreeOfFreedom> &indices) const
-    {
-      FUNCNAME("BasisFunction::getLocalIndices()");
-      
-      TEST_EXIT_DBG(static_cast<int>(indices.size()) >= nBasFcts)
-	("Index vector is too small!\n");
-
-      getLocalIndices(el, admin, &(indices[0]));
-    }
+    /// Calculates local dof indices of the element for the given DOF admin.
+    virtual void getLocalIndices(const Element *el,
+				 const DOFAdmin *admin,
+				 std::vector<DegreeOfFreedom> &indices) const
+    {}
 
     virtual void getLocalDofPtrVec(const Element *el, 
 				   const DOFAdmin *admin,

AMDiS/src/BoundaryManager.cc

@@ -15,7 +15,6 @@ namespace AMDiS {
   BoundaryManager::BoundaryManager(const FiniteElemSpace *feSpace)
   {
     localBounds.resize(omp_get_overall_max_threads());
-    dofIndices.resize(omp_get_overall_max_threads());
     allocatedMemoryLocalBounds = feSpace->getBasisFcts()->getNumber();
     for (int i = 0; i < static_cast<int>(localBounds.size()); i++)
       localBounds[i] = new BoundaryType[allocatedMemoryLocalBounds];
@@ -26,7 +25,6 @@ namespace AMDiS {
     localBCs = bm.localBCs;
     allocatedMemoryLocalBounds = bm.allocatedMemoryLocalBounds;
     localBounds.resize(bm.localBounds.size());
-    dofIndices.resize(bm.localBounds.size());
     for (int i = 0; i < static_cast<int>(localBounds.size()); i++)
       localBounds[i] = new BoundaryType[allocatedMemoryLocalBounds];
   }
@@ -49,54 +47,61 @@ namespace AMDiS {
     return result;
   }
 
+
   void BoundaryManager::fillBoundaryConditions(ElInfo *elInfo, 
 					       DOFVectorBase<double> *vec)
   {
-    if (localBCs.size() > 0) {
-      const FiniteElemSpace *feSpace = vec->getFESpace();
-      std::vector<DegreeOfFreedom> &dofVec = dofIndices[omp_get_thread_num()];
-      const BasisFunction *basisFcts = feSpace->getBasisFcts();
-      int nBasFcts = basisFcts->getNumber();
-      dofVec.resize(nBasFcts);
-
-      // get boundaries of all DOFs
-      BoundaryType *localBound = localBounds[omp_get_thread_num()];
-      basisFcts->getBound(elInfo, localBound);
-
-      // get dof indices
-      basisFcts->getLocalIndices(elInfo->getElement(), feSpace->getAdmin(), dofVec);
-
-      // apply non dirichlet boundary conditions
-      for (BoundaryIndexMap::iterator it = localBCs.begin(); it != localBCs.end(); ++it)
-	if ((*it).second && !(*it).second->isDirichlet())
-	  (*it).second->fillBoundaryCondition(vec, elInfo, &dofVec[0], 
-					      localBound, nBasFcts);
-
-      // apply dirichlet boundary conditions
-      for (BoundaryIndexMap::iterator it = localBCs.begin(); it != localBCs.end(); ++it)
-	if ((*it).second && (*it).second->isDirichlet())
-	  (*it).second->fillBoundaryCondition(vec, elInfo, &dofVec[0], 
-					      localBound, nBasFcts);
-    }
+    FUNCNAME("BoundaryManager::fillBoundaryConditions()");
+
+    if (localBCs.size() <= 0)
+      return;
+
+    const FiniteElemSpace *feSpace = vec->getFESpace();
+    const BasisFunction *basisFcts = feSpace->getBasisFcts();
+    int nBasFcts = basisFcts->getNumber();
+
+    // get boundaries of all DOFs
+    BoundaryType *localBound = localBounds[omp_get_thread_num()];
+    basisFcts->getBound(elInfo, localBound);
+    
+    // get dof indices
+    std::vector<DegreeOfFreedom> &dofVec = elInfo->getLocalIndices(feSpace);
+    TEST_EXIT_DBG(static_cast<int>(dofVec.size()) == nBasFcts)
+      ("Local index vector is too small!\n");
+    
+    // apply non dirichlet boundary conditions
+    for (BoundaryIndexMap::iterator it = localBCs.begin(); it != localBCs.end(); ++it)
+      if ((*it).second && !(*it).second->isDirichlet())
+	(*it).second->fillBoundaryCondition(vec, elInfo, &dofVec[0], 
+					    localBound, nBasFcts);
+    
+    // apply dirichlet boundary conditions
+    for (BoundaryIndexMap::iterator it = localBCs.begin(); it != localBCs.end(); ++it)
+      if ((*it).second && (*it).second->isDirichlet())
+	(*it).second->fillBoundaryCondition(vec, elInfo, &dofVec[0], 
+					    localBound, nBasFcts);    
   }
 
+
   void BoundaryManager::fillBoundaryConditions(ElInfo *elInfo, DOFMatrix *mat)
   {
+    FUNCNAME("BoundaryManager::fillBoundaryConditions()");
+
     if (localBCs.size() <= 0)
       return;
       
     const FiniteElemSpace *feSpace = mat->getRowFESpace();
-    std::vector<DegreeOfFreedom> &dofVec = dofIndices[omp_get_thread_num()];
     const BasisFunction *basisFcts = feSpace->getBasisFcts();
     int nBasFcts = basisFcts->getNumber();
-          dofVec.resize(nBasFcts);
 
     // get boundaries of all DOFs
     BoundaryType *localBound = localBounds[omp_get_thread_num()];
     basisFcts->getBound(elInfo, localBound);
     
     // get dof indices
-    basisFcts->getLocalIndices(elInfo->getElement(), feSpace->getAdmin(), dofVec);
+    std::vector<DegreeOfFreedom> &dofVec = elInfo->getLocalIndices(feSpace);
+    TEST_EXIT_DBG(static_cast<int>(dofVec.size()) == nBasFcts)
+      ("Local index vector is too small!\n");
     
     // apply non dirichlet boundary conditions
     for (BoundaryIndexMap::iterator it = localBCs.begin(); it != localBCs.end(); ++it)

AMDiS/src/BoundaryManager.h

@@ -125,9 +125,6 @@ namespace AMDiS {
     /// Temporary thread-safe variable for functions fillBoundaryconditions.
     std::vector<BoundaryType*> localBounds;
 
-    /// Temporary thread-safe variable for functions fillBoundaryconditions.
-    std::vector<std::vector<DegreeOfFreedom> > dofIndices;
-
     /** \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.

AMDiS/src/DOFMatrix.cc

@@ -168,9 +168,12 @@ namespace AMDiS {
  
     // === Get indices mapping from local to global matrix indices. ===
 
-    rowFESpace->getBasisFcts()->getLocalIndices(rowElInfo->getElement(),
-						rowFESpace->getAdmin(),
-						rowIndices);
+    std::vector<DegreeOfFreedom> &rowIndices2 = rowElInfo->getLocalIndices(rowFESpace);
+    std::vector<DegreeOfFreedom> &colIndices2 = rowIndices2;
+
+//     rowFESpace->getBasisFcts()->getLocalIndices(rowElInfo->getElement(),
+// 						rowFESpace->getAdmin(),
+// 						rowIndices);
     if (rowFESpace == colFESpace) {
       colIndices = rowIndices;
     } else {
@@ -188,21 +191,8 @@ namespace AMDiS {
 
     using namespace mtl;
 
-#if 0
-    std::cout << "----- PRINT MAT --------" << std::endl;
-    std::cout << elMat << std::endl;
-    std::cout << "rows: ";
-    for (int i = 0; i < rowIndices.size(); i++)
-      std::cout << rowIndices[i] << " ";
-    std::cout << std::endl;
-    std::cout << "cols: ";
-    for (int i = 0; i < colIndices.size(); i++)
-      std::cout << colIndices[i] << " ";
-    std::cout << std::endl;
-#endif
-         
     for (int i = 0; i < nRow; i++)  {
-      DegreeOfFreedom row = rowIndices[i];
+      DegreeOfFreedom row = rowIndices2[i];
 
       BoundaryCondition *condition = 
 	bound ? boundaryManager->getBoundaryCondition(bound[i]) : NULL;
@@ -210,20 +200,15 @@ namespace AMDiS {
       if (condition && condition->isDirichlet()) {
 	if (condition->applyBoundaryCondition()) {
 #ifdef HAVE_PARALLEL_DOMAIN_AMDIS
- 	  if ((*rankDofs)[rowIndices[i]]) 
+ 	  if ((*rankDofs)[row]) 
  	    applyDBCs.insert(static_cast<int>(row));
 #else
 	  applyDBCs.insert(static_cast<int>(row));
 #endif
 	}
       } else {
-	for (int j = 0; j < nCol; j++) {
-	  DegreeOfFreedom col = colIndices[j];
-	  double entry = elMat[i][j];
-  
-	  //	  std::cout << "ADD at " << row << " " << col << std::endl;
-	  ins[row][col] += entry;
-	}
+	for (int j = 0; j < nCol; j++)
+	  ins[row][colIndices2[j]] += elMat[i][j];	
       }
     }
   }

AMDiS/src/DOFVector.cc

@@ -108,13 +108,16 @@ namespace AMDiS {
     // traverse mesh
     std::vector<bool> visited(getUsedSize(), false);
     TraverseStack stack;
-    Flag fillFlag = Mesh::CALL_LEAF_EL | Mesh::FILL_GRD_LAMBDA | Mesh::FILL_COORDS;
+    stack.addFeSpace(feSpace);
+    Flag fillFlag = 
+      Mesh::CALL_LEAF_EL | Mesh::FILL_GRD_LAMBDA | 
+      Mesh::FILL_COORDS | Mesh::FILL_LOCAL_INDICES;
     ElInfo *elInfo = stack.traverseFirst(mesh, -1, fillFlag);
 
     while (elInfo) {
       const DegreeOfFreedom **dof = elInfo->getElement()->getDOF();
       const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
-      getLocalVector(elInfo->getElement(), localUh);
+      getLocalVector(elInfo, localUh);
 
       int localDOFNr = 0;
       for (int i = 0; i < nNodes; i++) { // for all nodes
@@ -177,9 +180,11 @@ namespace AMDiS {
     // traverse mesh
     Mesh *mesh = feSpace->getMesh();
     TraverseStack stack;
+    stack.addFeSpace(feSpace);
     ElInfo *elInfo = stack.traverseFirst(mesh, -1,
 					 Mesh::CALL_LEAF_EL | Mesh::FILL_DET | 
-					 Mesh::FILL_GRD_LAMBDA | Mesh::FILL_COORDS);
+					 Mesh::FILL_GRD_LAMBDA | Mesh::FILL_COORDS |
+					 Mesh::FILL_LOCAL_INDICES);
 
     double *localUh = new double[basFcts->getNumber()];
 
@@ -187,7 +192,7 @@ namespace AMDiS {
       double det = elInfo->getDet();
       const DegreeOfFreedom **dof = elInfo->getElement()->getDOF();
       const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
-      getLocalVector(elInfo->getElement(), localUh);
+      getLocalVector(elInfo, localUh);
       basFcts->evalGrdUh(bary, grdLambda, localUh, &grd);
 
       for (int i = 0; i < dim + 1; i++) {
@@ -249,7 +254,7 @@ namespace AMDiS {
     }
   
     double *localVec = localVectors[myRank];
-    getLocalVector(elInfo->getElement(), localVec);
+    getLocalVector(elInfo, localVec);
 
     DimVec<double> &grd1 = *grdTmp[myRank];
     int parts = Global::getGeo(PARTS, dim);
@@ -334,7 +339,7 @@ namespace AMDiS {
     }
   
     double *localVec = localVectors[myRank];
-    getLocalVector(largeElInfo->getElement(), localVec);
+    getLocalVector(largeElInfo, localVec);
 
     const BasisFunction *basFcts = feSpace->getBasisFcts();    
     mtl::dense2D<double> &m = smallElInfo->getSubElemCoordsMat(basFcts->getDegree());
@@ -391,8 +396,6 @@ namespace AMDiS {
     TEST_EXIT_DBG(!quadFast || quadFast->getBasisFunctions() == feSpace->getBasisFcts())
       ("invalid basis functions");
 
-    Element *el = elInfo->getElement(); 
-
     int myRank = omp_get_thread_num();
     int dow = Global::getGeo(WORLD);
     int nPoints = quadFast ? quadFast->getQuadrature()->getNumPoints() : quad->getNumPoints();
@@ -412,7 +415,7 @@ namespace AMDiS {
     }
   
     double *localVec = localVectors[myRank];
-    getLocalVector(el, localVec);
+    getLocalVector(elInfo, localVec);
 
     DimMat<double> D2Tmp(dim, DEFAULT_VALUE, 0.0);
     int parts = Global::getGeo(PARTS, dim);
@@ -776,7 +779,7 @@ namespace AMDiS {
     }
       
     double *localVec = localVectors[omp_get_thread_num()];
-    getLocalVector(largeElInfo->getElement(), localVec);
+    getLocalVector(largeElInfo, localVec);
 
     mtl::dense2D<double> &m = smallElInfo->getSubElemCoordsMat(basFcts->getDegree());
 

AMDiS/src/DOFVector.h

@@ -67,7 +67,9 @@ namespace AMDiS {
      * For the given element, this function returns an array of all DOFs of this
      * DOFVector that are defined on this element.
      */
-    virtual const T *getLocalVector(const Element *el, T* localVec) const;
+    const T *getLocalVector(const Element *el, T* localVec) const;
+
+    const T *getLocalVector(const ElInfo *elInfo, T* localVec) const;
 
     /** \brief
      * Evaluates the DOF vector at a set of quadrature points defined on the 

AMDiS/src/DOFVector.hh

@@ -147,9 +147,9 @@ namespace AMDiS {
   {
     FUNCNAME("DOFVector::addElementVector()");
 
-    std::vector<DegreeOfFreedom> indices(nBasFcts);
-    feSpace->getBasisFcts()->getLocalIndices(elInfo->getElement(), feSpace->getAdmin(),
-					     indices);
+    std::vector<DegreeOfFreedom> &indices = elInfo->getLocalIndices(feSpace);
+    TEST_EXIT_DBG(static_cast<int>(indices.size()) == nBasFcts)
+      ("Local index vector is too small!\n");
 
     for (DegreeOfFreedom i = 0; i < nBasFcts; i++) {
       BoundaryCondition *condition = 
@@ -454,9 +454,11 @@ namespace AMDiS {
     int nPoints = quadFast->getNumPoints();
     std::vector<T> uh_vec(nPoints);
     TraverseStack stack;
+    stack.addFeSpace(this->feSpace);
     ElInfo *elInfo = 
       stack.traverseFirst(mesh, -1, 
-			  Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_DET);
+			  Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | 
+			  Mesh::FILL_DET | Mesh::FILL_LOCAL_INDICES);
     while (elInfo) {
       double det = elInfo->getDet();
       double normT = 0.0;
@@ -490,9 +492,11 @@ namespace AMDiS {
     int nPoints = quadFast->getNumPoints();
     std::vector<T> uh_vec(nPoints);
     TraverseStack stack;
+    stack.addFeSpace(this->feSpace);
     ElInfo *elInfo = 
       stack.traverseFirst(mesh, -1, 
-			  Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_DET);
+			  Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | 
+			  Mesh::FILL_DET | Mesh::FILL_LOCAL_INDICES);
     while (elInfo) {
       double det = elInfo->getDet();
       double normT = 0.0;
@@ -526,9 +530,11 @@ namespace AMDiS {
     int nPoints = quadFast->getNumPoints();
     std::vector<T> uh_vec(nPoints);
     TraverseStack stack;
+    stack.addFeSpace(this->feSpace);
     ElInfo *elInfo = 
       stack.traverseFirst(mesh, -1, 
-			  Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_DET);
+			  Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | 
+			  Mesh::FILL_DET | Mesh::FILL_LOCAL_INDICES);
     while (elInfo) {
       double det = elInfo->getDet();
       double normT = 0.0;
@@ -986,6 +992,49 @@ namespace AMDiS {
   }
 
 
+  template<typename T>
+  const T *DOFVectorBase<T>::getLocalVector(const ElInfo *elInfo, T *d) const
+  {
+    FUNCNAME("DOFVectorBase<T>::getLocalVector()");
+
+    TEST_EXIT_DBG(feSpace->getMesh() == elInfo->getElement()->getMesh())
+      ("Element is defined on a different mesh than the DOF vector!\n");
+
+    static T* localVec = NULL;
+    static int localVecSize = 0;
+    T *result;
+  
+    if (d) {
+      result = d;
+    } else {
+#ifdef _OPENMP
+      ERROR_EXIT("Using static variable while using OpenMP parallelization!\n");
+#endif
+      if (localVec && nBasFcts > localVecSize)  {
+	delete [] localVec;
+	localVec = new T[nBasFcts]; 
+      }
+
+      if (!localVec)
+	localVec = new T[nBasFcts]; 
+
+      localVecSize = nBasFcts;
+      result = localVec;
+    }
+
+    std::vector<DegreeOfFreedom> &localIndices = 
+      const_cast<ElInfo*>(elInfo)->getLocalIndices(feSpace);
+
+    TEST_EXIT_DBG(static_cast<int>(localIndices.size()) == nBasFcts)
+      ("Local indices vector is too small!\n");
+
+    for (int i = 0; i < nBasFcts; i++)
+      result[i] = (*this)[localIndices[i]];
+
+    return result;
+  }
+
+
   template<typename T>
   const T *DOFVectorBase<T>::getVecAtQPs(const ElInfo *elInfo, 
 					 const Quadrature *quad,
@@ -1004,7 +1053,6 @@ namespace AMDiS {
     TEST_EXIT_DBG(!quadFast || quadFast->getBasisFunctions() == feSpace->getBasisFcts())
       ("invalid basis functions");
 
-    Element *el = elInfo->getElement();
     const Quadrature *quadrature = quadFast ? quadFast->getQuadrature() : quad;
     const BasisFunction *basFcts = feSpace->getBasisFcts();
     int nPoints = quadrature->getNumPoints();
@@ -1028,7 +1076,7 @@ namespace AMDiS {
     }
 
     T *localVec = localVectors[omp_get_thread_num()];
-    getLocalVector(el, localVec);
+    getLocalVector(elInfo, localVec);