* Better paralization of buildAfterCoarsen

AMDiS/src/ProblemVec.cc

@@ -747,31 +747,11 @@ namespace AMDiS {
 	assembledMatrix_[i][j] = true;
       }
 
-      // fill boundary conditions
-      if (rhs_->getDOFVector(i)->getBoundaryManager())
-	rhs_->getDOFVector(i)->getBoundaryManager()->initVector(rhs_->getDOFVector(i));     
-      
-      if (solution_->getDOFVector(i)->getBoundaryManager())
-      	solution_->getDOFVector(i)->getBoundaryManager()->initVector(solution_->getDOFVector(i));
-
-      TraverseStack stack;
-      ElInfo *elInfo = stack.traverseFirst(componentMeshes[i], -1, assembleFlag);
-      while (elInfo) {
-	if (rhs_->getDOFVector(i)->getBoundaryManager())
-	  rhs_->getDOFVector(i)->getBoundaryManager()->
-	    fillBoundaryConditions(elInfo, rhs_->getDOFVector(i));
-
-	if (solution_->getDOFVector(i)->getBoundaryManager())
-	  solution_->getDOFVector(i)->getBoundaryManager()->
-	    fillBoundaryConditions(elInfo, solution_->getDOFVector(i));
-
-	elInfo = stack.traverseNext(elInfo);
-      }
-
-      if (rhs_->getDOFVector(i)->getBoundaryManager())
-	rhs_->getDOFVector(i)->getBoundaryManager()->exitVector(rhs_->getDOFVector(i));
-      if (solution_->getDOFVector(i)->getBoundaryManager())
-	solution_->getDOFVector(i)->getBoundaryManager()->exitVector(solution_->getDOFVector(i));    
+      // And now assemble boundary conditions on the vectors
+      assembleBoundaryConditions(rhs_->getDOFVector(i),
+				 solution_->getDOFVector(i),
+				 componentMeshes[i],
+				 assembleFlag);
     }
 
 #ifdef _OPENMP
@@ -1116,6 +1096,86 @@ namespace AMDiS {
     }
   }
 
+  void ProblemVec::assembleBoundaryConditions(DOFVector<double> *rhs,
+					      DOFVector<double> *solution,
+					      Mesh *mesh,
+					      Flag assembleFlag)
+  {
+    /* ================ Initialization of vectors ==================== */
+
+    if (rhs->getBoundaryManager())
+      rhs->getBoundaryManager()->initVector(rhs);      
+    if (solution->getBoundaryManager())
+      solution->getBoundaryManager()->initVector(solution);
+    
+#ifdef _OPENMP
+    TraverseParallelStack stack;
+#else
+    TraverseStack stack;
+#endif
+
+    /* ================= Parallel Boundary Assemblage ================= */
+#ifdef _OPENMP
+#pragma omp parallel
+#endif
+    {
+      // Each thread assembles on its own dof-vectors.
+      DOFVector<double> *tmpRhsVec = NEW DOFVector<double>(rhs->getFESpace(), "tmpRhs");
+      DOFVector<double> *tmpSolVec = NEW DOFVector<double>(solution->getFESpace(), "tmpSol");
+      tmpRhsVec->set(0.0);
+      tmpSolVec->set(0.0);
+
+
+      // (Parallel) traverse of mesh.
+      ElInfo *elInfo = stack.traverseFirst(mesh, -1, assembleFlag);
+      while (elInfo) {
+	if (rhs->getBoundaryManager())
+	  rhs->getBoundaryManager()-> fillBoundaryConditions(elInfo, tmpRhsVec);
+	
+	if (solution->getBoundaryManager())
+	  solution->getBoundaryManager()->fillBoundaryConditions(elInfo, tmpSolVec);
+      
+	elInfo = stack.traverseNext(elInfo);
+      }
+
+
+      // After (parallel) mesh traverse, the result is applied to the final
+      // vectors. This section is not allowed to be executed by more than one
+      // thread at the same time.
+#ifdef _OPENMP
+#pragma omp critical
+#endif
+      {
+	DOFVector<double>::Iterator rhsIt(rhs, USED_DOFS);
+	DOFVector<double>::Iterator solIt(solution, USED_DOFS);
+	DOFVector<double>::Iterator tmpRhsIt(tmpRhsVec, USED_DOFS);
+	DOFVector<double>::Iterator tmpSolIt(tmpSolVec, USED_DOFS);
+	for (rhsIt.reset(), solIt.reset(), tmpRhsIt.reset(), tmpSolIt.reset();
+	     !rhsIt.end();
+	     ++rhsIt, ++solIt, ++tmpRhsIt, ++tmpSolIt) {	     
+	  if (*tmpRhsIt != 0.0)
+	    *rhsIt = *tmpRhsIt;
+	  
+	  if (*tmpSolIt != 0.0)
+	    *solIt = *tmpSolIt;	  
+	}
+      } // pragma omp critical
+
+
+      DELETE tmpRhsVec;
+      DELETE tmpSolVec;
+    } // pragma omp parallel
+     
+
+    /* ======================= Finalize vectors ================== */
+
+    if (rhs->getBoundaryManager())
+      rhs->getBoundaryManager()->exitVector(rhs);
+    if (solution->getBoundaryManager())
+      solution->getBoundaryManager()->exitVector(solution);
+  }
+
+
   void ProblemVec::writeResidualMesh(AdaptInfo *adaptInfo, const std::string name)
   {
     FUNCNAME("ProblemVec::writeResidualMesh()");

AMDiS/src/ProblemVec.h

@@ -300,6 +300,11 @@ namespace AMDiS {
     void assembleOnDifMeshes(FiniteElemSpace *rowFeSpace, FiniteElemSpace *colFeSpace,
 			     Flag assembleFlag,
 			     DOFMatrix *matrix, DOFVector<double> *vector);
+
+    void assembleBoundaryConditions(DOFVector<double> *rhs,
+				    DOFVector<double> *solution,
+				    Mesh *mesh,
+				    Flag assembleFlag);
   
     // ===== getting-methods ======================================================
 

AMDiS/src/RobinBC.cc

@@ -175,39 +175,36 @@ namespace AMDiS {
   }
 
   void RobinBC::fillBoundaryCondition(DOFVectorBase<double>* vector, 
-				      ElInfo*                   elInfo,
-				      const DegreeOfFreedom*    dofIndices,
-				      const BoundaryType*       localBound,
-				      int                       nBasFcts)
+				      ElInfo* elInfo,
+				      const DegreeOfFreedom* dofIndices,
+				      const BoundaryType* localBound,
+				      int nBasFcts)
   {
     FUNCNAME("RobinBC::fillBoundaryCondition()");
-    TEST_EXIT(vector->getFESpace() == rowFESpace)
-      ("invalid row fe space\n");
+    TEST_EXIT_DBG(vector->getFESpace() == rowFESpace)("invalid row fe space\n");
 
     int dim = elInfo->getMesh()->getDim();
-    int i;
 
-    if(neumannOperators) {
-      for(i=0; i < dim+1; i++) {
-	if(elInfo->getBoundary(i) == boundaryType) {
+    if (neumannOperators) {
+      for (int i = 0; i < dim + 1; i++) {
+	if (elInfo->getBoundary(i) == boundaryType) {
 	  vector->assemble(1.0, elInfo, localBound, (*neumannOperators)[i]);
 	}
       }
     }
   }
 
-  void RobinBC::fillBoundaryCondition(DOFMatrix*     matrix,
-				      ElInfo*                elInfo,
+  void RobinBC::fillBoundaryCondition(DOFMatrix* matrix,
+				      ElInfo* elInfo,
 				      const DegreeOfFreedom* dofIndices,
-				      const BoundaryType*    localBound,
-				      int                    nBasFcts) 
+				      const BoundaryType* localBound,
+				      int nBasFcts) 
   {
     int dim = elInfo->getMesh()->getDim();
-    int i;
 
-    if(robinOperators) {
-      for(i=0; i < dim+1; i++) {
-	if(elInfo->getBoundary(i) == boundaryType) {
+    if (robinOperators) {
+      for (int i = 0; i < dim + 1; i++) {
+	if (elInfo->getBoundary(i) == boundaryType) {
 	  matrix->assemble(-1.0, elInfo, localBound, (*robinOperators)[i]);
 	}
       }