Periodic boundary condition work for higher elements.

AMDiS/src/ElInfo2d.cc

@@ -77,31 +77,82 @@ namespace AMDiS {
 	  Element *nb = const_cast<Element*>(neighbour_[i]);
 
 	  int edgeNo = oppVertex_[i] = mel->getOppVertex(i);
-	  if (nb->getFirstChild() && (edgeNo != 2)) {   // make nb nearest el.
+
+
+	  if (nb->getFirstChild() && (edgeNo != 2)) {  
+
+	    // Search for the next neighbour. In many cases, the neighbour element 
+	    // may be refinemed in a way, such that there is no new vertex on the 
+	    // common edge. This situation is shown in the following picture: 
+	    //
+	    //               /|\
+	    //              / | \
+	    //             /  |  \
+	    //            /\  |   \
+	    //           /  \ |    \ 
+	    //          /    \|     \
+	    //          -------------
+	    //
+	    //            nb     el
+	    //
+	    // Note that we know (because of the last if statement), that the 
+	    // neighbour element has children and the common edge is not the 
+	    // refinement edge, which has always the number 2, of our element.
+
+
 	    if (edgeNo == 0) {
+	      // The situation is as follows:
+	      //
+	      //          -------
+	      //          \    /|\
+	      //           \  / | \
+	      //            \/  |  \
+	      //             \  |   \
+	      //              \ |    \ 
+	      //               \|     \
+	      //                -------
+	      //
+	      //            nb     el
+	      // That means, the edge 0 of the same level neighbour is the common
+	      // edge, i.e., the direct neighbour is the second child of the same
+	      // level neighbour.
+
 	      nb = neighbour_[i] = nb->getSecondChild();
 	    } else {
+	      // The situation is as shown in the picture above. So the next
+	      // neighbour is the first child of the same level neighbour element.
 	      nb = neighbour_[i] = nb->getFirstChild();
 	    }
 
+	    // In both cases the opp vertex number is 2, as one can see in the 
+	    // pictures above.
 	    oppVertex_[i] = 2;
 
 	    if (fill_opp_coords) {
 	      if (nb->isNewCoordSet()) {
 		oppCoord_[i] = *(nb->getNewCoord());
 	      } else {
-		oppCoord_[i] = (macroNeighbour->coord[0] + macroNeighbour->coord[1]) * 0.5;
-	      }	    	      
+		// In both cases, that are shown in the pictures above, the opp
+		// vertex of the neighbour edge is the midpoint of the vertex 0
+		// and vertex 1 of the same level neighbour element.
+		oppCoord_[i] = (macroNeighbour->coord[0] + 
+				macroNeighbour->coord[1]) * 0.5;
+	      }
 	      
 	      switch (i) {
 	      case 0:
-		if (*(macroNeighbour->getElement()->getDOF(2)) == *(element_->getDOF(2))) {
+		// The common edge is the edge 0 of this element.
+
+		switch (edgeNo) {
+		case 1:
 		  neighbourCoord_[i][0] = macroNeighbour->coord[2];
 		  neighbourCoord_[i][1] = macroNeighbour->coord[0];
-		} else if (*(macroNeighbour->getElement()->getDOF(2)) == *(element_->getDOF(1))) {
+		  break;
+		case 0:		  
 		  neighbourCoord_[i][0] = macroNeighbour->coord[1];
 		  neighbourCoord_[i][1] = macroNeighbour->coord[2];
-		} else {
+		  break;
+		default:
 		  ERROR_EXIT("Should not happen!\n");
 		}
 	
@@ -109,13 +160,17 @@ namespace AMDiS {
 		break;
 		
 	      case 1:
-		if (*(macroNeighbour->getElement()->getDOF(2)) == *(element_->getDOF(2))) {
+		// The commonedge is the edge 1 of this element.
+		switch (edgeNo) {
+		case 0:
 		  neighbourCoord_[i][0] = macroNeighbour->coord[1];
 		  neighbourCoord_[i][1] = macroNeighbour->coord[2];
-		} else if (*(macroNeighbour->getElement()->getDOF(2)) == *(element_->getDOF(0))) {
+		  break;
+		case 1:
 		  neighbourCoord_[i][0] = macroNeighbour->coord[2];
 		  neighbourCoord_[i][1] = macroNeighbour->coord[0];
-		} else {
+		  break;
+		default:
 		  ERROR_EXIT("Should not happen!\n");
 		}
 		
@@ -123,15 +178,10 @@ namespace AMDiS {
 		break;
 		
 	      case 2:
-		if (*(macroNeighbour->getElement()->getDOF(2)) == *(element_->getDOF(0))) {
-		  neighbourCoord_[i][0] = macroNeighbour->coord[2];
-		  neighbourCoord_[i][1] = macroNeighbour->coord[1];
-		} else if (*(macroNeighbour->getElement()->getDOF(2)) == *(element_->getDOF(1))) {
-		  neighbourCoord_[i][0] = macroNeighbour->coord[0];
-		  neighbourCoord_[i][1] = macroNeighbour->coord[2];		 
-		} else {
-		  ERROR_EXIT("Should not happen!\n");
-		}
+		// I've deleted here some code, be I think that this case is not
+		// possible. If an error occurs in this line, please check AMDiS
+		// revision <= 476 at the same position.
+		ERROR_EXIT("Should not happen!\n");
 		break;
 
 	      default:
@@ -146,17 +196,24 @@ namespace AMDiS {
 		  } else {
 		    std::cout << "  Neighbour " << j << ": not existing" << std::endl;
 		  }
-		  std::cout << "  OppVertex " << j << ": " << static_cast<int>(mel->getOppVertex(j)) << std::endl;
-		  std::cout << std::endl;
+		  std::cout << "  OppVertex " << j << ": " 
+			    << static_cast<int>(mel->getOppVertex(j)) 
+			    << std::endl << std::endl;
 		}
 		ERROR_EXIT("should not happen!\n");
 		break;
 	      }
 	    }
 	  } else {
+
+	    // In this case, we know that the common edge is the refinement edge.
+	    // This makes everything much more simpler, because we know that the
+	    // next neighbour is equal to the samel level neighbour. If the same
+	    // level neighbour would be refinement, also this element must to be 
+	    // refinement, because they share the refinement edge.
+
 	    if (fill_opp_coords) {
 	      oppCoord_[i] = macroNeighbour->coord[edgeNo];
-
 	      neighbourCoord_[i] = macroNeighbour->coord;	      
 	    }
 	  }

AMDiS/src/PeriodicBC.cc

@@ -205,32 +205,19 @@ namespace AMDiS {
       masterMatrix_ = NULL;
     }
 
-    matrix->print();
-
     using namespace mtl;
 
     DOFAdmin* admin = rowFESpace->getAdmin();
     std::vector<int> dofMap(admin->getUsedSize());
-    for (int i = 0; i < admin->getUsedSize(); i++) {
+    for (int i = 0; i < admin->getUsedSize(); i++)
       dofMap[i] = (*associated)[i];
-      std::cout << "map " << i << " to " << dofMap[i] << std::endl;
-    }
 
     // Compute reorder matrix (newRow and newCol yields transposed!!!)
     matrix::traits::reorder<>::type R= matrix::reorder(dofMap);
-    DOFMatrix::base_matrix_type &A= matrix->getBaseMatrix(), B, D, E, TR;
-
-    A*= 0.5;
-    // Half of entries with decreased row index + half of the strict lower origing
-    B= strict_upper(R) * A + strict_lower(A);
-    // Add half of entries with same row and decreased column index + half of increased columns
-    // Remark regarding column permutation: trans(strict_upper(trans(X))) == strict_lower(X)
-    // B+= bands(trans(R), 0, 1) * (A * strict_lower(R) + lower(A));  -> make this work one day!!!
-    TR= trans(R);
-    D= bands(TR, 0, 1);
-    E= A * strict_lower(R) + lower(A);
-    B+= D * E;
-    swap(A, B);  
+    DOFMatrix::base_matrix_type &A= matrix->getBaseMatrix(), C;
+
+    C = R * A * trans(R) + A;
+    A = 0.5 * C;
   }
 
   void PeriodicBC::exitVector(DOFVectorBase<double>* vector)

AMDiS/src/ProblemScal.cc

@@ -586,6 +586,8 @@ namespace AMDiS {
 
     systemMatrix->removeRowsWithDBC(systemMatrix->getApplyDBCs());
 
+    systemMatrix->finishInsertion();
+
     // TODO: ExitMatrix should be called after finishInsertion!
     if (systemMatrix->getBoundaryManager())
       systemMatrix->getBoundaryManager()->exitMatrix(systemMatrix);
@@ -594,8 +596,6 @@ namespace AMDiS {
     if (solution->getBoundaryManager())
       solution->getBoundaryManager()->exitVector(solution);
 
-    systemMatrix->finishInsertion();
-
     INFO(info, 8)("buildAfterCoarsen needed %.5f seconds\n", TIME_USED(first,clock()));
 
 #ifdef HAVE_PARALLEL_DOMAIN_AMDIS
@@ -606,8 +606,6 @@ namespace AMDiS {
 #endif
     
     createPrecon();
-
-    systemMatrix->print();
   }
 
   void ProblemScal::writeResidualMesh(AdaptInfo *adaptInfo, const std::string name)

AMDiS/src/TimedObject.h

@@ -24,10 +24,6 @@
 
 namespace AMDiS {
 
-  // ===========================================================================
-  // ===== class TimedObject ===================================================
-  // ===========================================================================
-
   /** \brief
    * This class can be used as base class for time dependent objects where
    * different objects refer to the same time. Herefore a pointer to
@@ -37,24 +33,25 @@ namespace AMDiS {
   class TimedObject
   {
   public:
-    /** \brief
-     * Constructor.
-     */
-    TimedObject() : timePtr(NULL) {};
-
-    /** \brief
-     * Sets the time pointer.
-     */
-    inline void setTimePtr(double *timePtr_) { timePtr = timePtr_; };
-
-    /** \brief
-     * Returns the time pointer.
-     */
-    inline double *getTimePtr() { return timePtr; };
+    /// Constructor.
+    TimedObject() 
+      : timePtr(NULL) 
+    {}
+
+    /// Sets the time pointer.
+    inline void setTimePtr(double *timePtr_) 
+    { 
+      timePtr = timePtr_; 
+    }
+
+    /// Returns the time pointer.
+    inline double *getTimePtr() 
+    { 
+      return timePtr; 
+    }
+
   protected:
-    /** \brief
-     * Pointer to the externally managed time value.
-     */
+    /// Pointer to the externally managed time value.
     double *timePtr;
   };