Periodic BC implemented by Siqi Ling added to this branch

AMDiS/src/Mesh.cc

@@ -1128,24 +1128,13 @@ namespace AMDiS {
     // === Write periodic associations. ===
     int mapSize = periodicAssociations.size();
     SerUtil::serialize(out, mapSize);
-    for (map<BoundaryType, VertexVector*>::iterator it = periodicAssociations.begin();
+    for (map<BoundaryType, std::vector<VertexVector*> >::iterator it = periodicAssociations.begin();
 	 it != periodicAssociations.end(); ++it) {
       BoundaryType b = it->first;
 
-      // Check which DOFAdmin is used for the current VertexVector we want to serialize.
-      int ithAdmin = -1;
-      for (int i = 0; i < static_cast<int>(admin.size()); i++) {
-	if (admin[i] == it->second->getAdmin()) {
-	  ithAdmin = i;
-	  break;
-	}
-      }
-      TEST_EXIT(ithAdmin >= 0)
-	("No DOFAdmin found for serialization of periodic associations!\n");
-
       SerUtil::serialize(out, b);
-      SerUtil::serialize(out, ithAdmin);
-      it->second->serialize(out);
+      for (int i = 0; i < size; i++)
+        (it->second)[i]->serialize(out);
     }
 
     serializedDOFs.clear();
@@ -1278,7 +1267,7 @@ namespace AMDiS {
 
       VertexVector *tmpvec = new VertexVector(admin[ithAdmin], "");
       tmpvec->deserialize(in);
-      periodicAssociations[b] = tmpvec;
+      periodicAssociations[b].push_back(tmpvec);
     }
   }
 
@@ -1514,38 +1503,55 @@ namespace AMDiS {
 #endif
 
 
-  bool Mesh::associated(DegreeOfFreedom dof1, DegreeOfFreedom dof2)
+   VertexVector& Mesh::getPeriodicAssociations(BoundaryType b, const DOFAdmin* admin)
+   {
+     FUNCNAME_DBG("Mesh::getPeriodicAssociations()");
+
+     TEST_EXIT_DBG(periodicAssociations.count(b) == 1)
+       ("There are no periodic assoications for boundary type %d!\n", b);
+
+     return admin ? (*(periodicAssociations[b][getAdminIndex(admin)])) : (*(periodicAssociations[b][0]));
+  }
+
+
+  bool Mesh::associated(DegreeOfFreedom dof1, DegreeOfFreedom dof2, int iadmin)
   {
-    map<BoundaryType, VertexVector*>::iterator it;
-    map<BoundaryType, VertexVector*>::iterator end = periodicAssociations.end();
+    FUNCNAME_DBG("Mesh::associated()");
+    TEST_EXIT_DBG(iadmin < admin.size())("Wrong DOF admin index\n");
+
+    map<BoundaryType, std::vector<VertexVector*> >::iterator it;
+    map<BoundaryType, std::vector<VertexVector*> >::iterator end = periodicAssociations.end();
 
     for (it = periodicAssociations.begin(); it != end; ++it)
-      if ((*(it->second))[dof1] == dof2)
-	return true;
+      if ((*(it->second)[iadmin])[dof1] == dof2)
+        return true;
 
     return false;
   }
 
 
-  bool Mesh::indirectlyAssociated(DegreeOfFreedom dof1, DegreeOfFreedom dof2)
+  bool Mesh::indirectlyAssociated(DegreeOfFreedom dof1, DegreeOfFreedom dof2, int iadmin)
   {
+    FUNCNAME_DBG("Mesh::indirectlyAssociated()");
+    TEST_EXIT_DBG(iadmin < admin.size())("Wrong DOF admin index\n");
+
     vector<DegreeOfFreedom> associatedToDOF1;
-    map<BoundaryType, VertexVector*>::iterator it;
-    map<BoundaryType, VertexVector*>::iterator end = periodicAssociations.end();
+    map<BoundaryType, std::vector<VertexVector*> >::iterator it;
+    map<BoundaryType, std::vector<VertexVector*> >::iterator end = periodicAssociations.end();
     DegreeOfFreedom dof, assDOF;
 
     associatedToDOF1.push_back(dof1);
     for (it = periodicAssociations.begin(); it != end; ++it) {
       int size = static_cast<int>(associatedToDOF1.size());
       for (int i = 0; i < size; i++) {
-	dof = associatedToDOF1[i];
-	assDOF = (*(it->second))[dof];
-	if (assDOF == dof2) {
-	  return true;
-	} else {
-	  if (assDOF != dof)
-	    associatedToDOF1.push_back(assDOF);
-	}
+        dof = associatedToDOF1[i];
+        assDOF = (*(it->second)[iadmin])[dof];
+        if (assDOF == dof2) {
+          return true;
+        } else {
+          if (assDOF != dof)
+            associatedToDOF1.push_back(assDOF);
+        }
       }
     }
 

AMDiS/src/Mesh.h

@@ -192,6 +192,18 @@ namespace AMDiS {
       return *(admin[i]);
     }
 
+    /// Returns the index of a given DOFMadmin pointer in the admin vector
+    inline int getAdminIndex(const DOFAdmin* dofAdmin) const
+    {
+      size_t i = 0;
+      for (; i < admin.size(); i++)
+        if (admin[i] == dofAdmin)
+          break;
+
+      TEST_EXIT_DBG(i < admin.size())("Wrong DOF admin index: %d\n", i);
+      return i;
+    }
+
     /// Creates a DOFAdmin with name lname. nDof specifies how many DOFs
     /// are needed at the different positions (see \ref DOFAdmin::nrDOF).
     /// A pointer to the created DOFAdmin is returned.
@@ -549,25 +561,24 @@ namespace AMDiS {
     }
 
     ///
-    inline std::map<BoundaryType, VertexVector*>& getPeriodicAssociations()
+    inline std::map<BoundaryType, std::vector<VertexVector*> >& getPeriodicAssociations()
     {
       return periodicAssociations;
     }
 
     /// Returns the periodic association for a specific boundary type.
-    inline VertexVector& getPeriodicAssociations(BoundaryType b)
-    {
-      FUNCNAME_DBG("Mesh::getPeriodicAssociations()");
-
-      TEST_EXIT_DBG(periodicAssociations.count(b) == 1)
-	("There are no periodic assoications for boundary type %d!\n", b);
-
-      return (*(periodicAssociations[b]));
-    }
+    VertexVector& getPeriodicAssociations(BoundaryType b, const DOFAdmin* admin = NULL);
 
     inline void setPeriodicAssociations(BoundaryType b, VertexVector* vec)
     {
-      periodicAssociations[b] = vec;
+      std::map<BoundaryType, std::vector<VertexVector*> >::iterator lb = periodicAssociations.lower_bound(b);
+      if (lb != periodicAssociations.end() && !periodicAssociations.key_comp()(b, lb->first)) {
+        lb->second.push_back(vec);
+      } else {
+        std::vector<VertexVector*> v;
+        v.push_back(vec);
+        periodicAssociations.insert(lb, std::make_pair(b, v));
+      }
     }
 
 
@@ -579,10 +590,20 @@ namespace AMDiS {
     }
 
     ///
-    bool associated(DegreeOfFreedom dof1, DegreeOfFreedom dof2);
+    bool associated(DegreeOfFreedom dof1, DegreeOfFreedom dof2, int iadmin = 0);
+
+    inline bool associated(DegreeOfFreedom dof1, DegreeOfFreedom dof2, const DOFAdmin* dofAdmin)
+    {
+      return associated(dof1, dof2, getAdminIndex(dofAdmin));
+    }
 
     ///
-    bool indirectlyAssociated(DegreeOfFreedom dof1, DegreeOfFreedom dof2);
+    bool indirectlyAssociated(DegreeOfFreedom dof1, DegreeOfFreedom dof2, int iadmin = 0);
+
+    inline bool indirectlyAssociated(DegreeOfFreedom dof1, DegreeOfFreedom dof2, const DOFAdmin* dofAdmin)
+    {
+      return indirectlyAssociated(dof1, dof2, getAdminIndex(dofAdmin));
+    }
 
     /// Returns \macroFileInfo
     inline MacroInfo* getMacroFileInfo()
@@ -854,7 +875,7 @@ namespace AMDiS {
     bool initialized;
 
     /// Map of managed periodic vertex associations.
-    std::map<BoundaryType, VertexVector*> periodicAssociations;
+    std::map<BoundaryType, std::vector<VertexVector*> > periodicAssociations;
 
     /// If the mesh has been created by reading a macro file, here the information
     /// are stored about the content of the file.

AMDiS/src/PeriodicBC.cc

@@ -106,13 +106,14 @@ namespace AMDiS {
   }
 
 
-  PeriodicBC::PeriodicBC(BoundaryType type, const FiniteElemSpace *rowSpace)
+  PeriodicBC::PeriodicBC(BoundaryType type, const FiniteElemSpace *rowSpace, bool diagonal)
     : BoundaryCondition(type, rowSpace, NULL),
-      masterMatrix(NULL)
+      masterMatrix(NULL),
+      isDiagonal(diagonal)
   {
     if (rowFeSpace->getMesh()->getDim() > 1)
       periodicDOFMapping =
-	PeriodicDOFMapping::providePeriodicDOFMapping(rowFeSpace->getBasisFcts());
+        PeriodicDOFMapping::providePeriodicDOFMapping(rowFeSpace->getBasisFcts());
     else
       periodicDOFMapping = NULL;
   }
@@ -129,10 +130,10 @@ namespace AMDiS {
     if (!masterMatrix) {
       masterMatrix = matrix;
       Mesh *mesh = matrix->getRowFeSpace()->getMesh();
-      associated = mesh->getPeriodicAssociations()[boundaryType];
+      associated = &(mesh->getPeriodicAssociations(boundaryType, rowFeSpace->getAdmin()));
 
       TEST_EXIT(associated)
-	("No associations for periodic boundary condition %d!\n", boundaryType);
+        ("No associations for periodic boundary condition %d!\n", boundaryType);
     }
   }
 
@@ -153,6 +154,7 @@ namespace AMDiS {
       return;
 
     DOFAdmin *admin = rowFeSpace->getAdmin();
+    int iadmin = rowFeSpace->getMesh()->getAdminIndex(admin);
     FixVec<int, WORLD> elFace(dim, NO_INIT);
     FixVec<int, WORLD> neighFace(dim, NO_INIT);
     DimVec<int> vertexPermutation(dim, NO_INIT);
@@ -186,7 +188,7 @@ namespace AMDiS {
 
 	    int j = 0;
 	    for (; j < dim + 1; j++)
-	      if (neigh->getDof(j, 0) == periodicDOF)
+          if (neigh->getDof(j, iadmin) == periodicDOF)
 		break;
 
 	    vertexPermutation[i] = j;
@@ -213,42 +215,71 @@ namespace AMDiS {
     FUNCNAME("PeriodicBC::exitMatrix()");
 
     TEST_EXIT(matrix)("No matrix\n");
-
     TEST_EXIT(associated)("No associated vector!\n");
 
     if (matrix == masterMatrix)
       masterMatrix = NULL;
 
     using namespace mtl;
+    typedef DOFMatrix::base_matrix_type Matrix;
+    typedef mtl::traits::range_generator<mtl::tag::row, Matrix>::type c_type;
+    typedef mtl::traits::range_generator<mtl::tag::nz, c_type>::type  ic_type;
 
     DOFAdmin* admin = rowFeSpace->getAdmin();
-    std::vector<int> dofMap(admin->getUsedSize());
-    for (int i = 0; i < admin->getUsedSize(); i++)
-      dofMap[i] = (*associated)[i];
+    int adminSize = admin->getUsedSize();
+    Matrix &A= matrix->getBaseMatrix();
+
+    mtl::traits::col<Matrix>::type c(A);
+    mtl::traits::value<Matrix>::type v(A);
+
+    std::vector<std::vector<std::pair<DegreeOfFreedom, double> > > row_values;
+    row_values.resize(adminSize);
 
-    // Compute reorder matrix (newRow and newCol yields transposed!!!)
-    matrix::traits::reorder<>::type R= matrix::reorder(dofMap);
-    DOFMatrix::base_matrix_type &A= matrix->getBaseMatrix(), C;
+    for (DegreeOfFreedom i = 0; i < adminSize; i++) {
+      if (i < (*associated)[i]) {
+        c_type cursor(begin<mtl::tag::row>(A) + i);
 
-    C = R * A * trans(R) + A;
-    A = 0.5 * C;
+        for (ic_type icursor(begin<mtl::tag::nz>(cursor)), icend(end<mtl::tag::nz>(cursor));
+            icursor != icend; ++icursor) {
+          row_values[i].push_back(std::make_pair(c(*icursor), v(*icursor)));
+          v(*icursor, 0.0);
+        }
+      }
+    }
+    matrix::inserter<Matrix, update_plus<double> > ins(A);
+    if (isDiagonal) {
+      for (DegreeOfFreedom i = 0; i < adminSize; i++) {
+        if (i < (*associated)[i]) {
+          ins[i][i] << 1.0;
+          ins[i][(*associated)[i]] << -1.0;
+        }
+      }
+    }
+    for (DegreeOfFreedom i = 0; i < adminSize; i++) {
+      if (i < (*associated)[i]) {
+        for (size_t j = 0; j < row_values[i].size(); j++)
+          ins[(*associated)[i]][row_values[i][j].first] << row_values[i][j].second;
+      }
+    }
   }
 
 
   void PeriodicBC::exitVector(DOFVectorBase<double>* vector)
   {
+    FUNCNAME("PeriodicBC::exitVector()");
+    TEST_EXIT(rowFeSpace == vector->getFeSpace())("Should not happen.\n");
+
     DOFIterator<double> vecIt(vector, USED_DOFS);
     Mesh *mesh = vector->getFeSpace()->getMesh();
-    VertexVector *associated = mesh->getPeriodicAssociations()[boundaryType];
+    VertexVector *associated = &(mesh->getPeriodicAssociations(boundaryType, rowFeSpace->getAdmin()));
 
     for (vecIt.reset(); !vecIt.end(); ++vecIt) {
       DegreeOfFreedom dof = vecIt.getDOFIndex();
       DegreeOfFreedom newDOF = (*associated)[dof];
 
       if (dof < newDOF) {
-	double entry = ((*vector)[dof] + (*vector)[newDOF]) * 0.5;
-	(*vector)[dof] = entry;
-	(*vector)[newDOF] = entry;
+        (*vector)[newDOF] += (*vector)[dof];
+        (*vector)[dof] = 0.0;
       }
     }
   }

AMDiS/src/PeriodicBC.h

@@ -88,7 +88,7 @@ namespace AMDiS {
   {
   public:
     /// Constructor.
-    PeriodicBC(BoundaryType type, const FiniteElemSpace *rowFeSpace);
+    PeriodicBC(BoundaryType type, const FiniteElemSpace *rowFeSpace, bool diagonal);
 
     ~PeriodicBC();
 
@@ -116,6 +116,8 @@ namespace AMDiS {
     PeriodicDOFMapping *periodicDOFMapping;
 
     DOFMatrix *masterMatrix;
+
+    bool isDiagonal;
   };
 
 }

AMDiS/src/ProblemStat.cc

@@ -1549,7 +1549,7 @@ namespace AMDiS {
   void ProblemStatSeq::addPeriodicBC(BoundaryType type, int row, int col)
   {
     boundaryConditionSet = true;
-    PeriodicBC *periodic = new PeriodicBC(type, componentSpaces[row]);
+    PeriodicBC *periodic = new PeriodicBC(type, componentSpaces[row], row == col);
 
     if (systemMatrix && (*systemMatrix)[row][col])
       (*systemMatrix)[row][col]->getBoundaryManager()->addBoundaryCondition(periodic);

AMDiS/src/RefinementManager2d.cc

@@ -74,9 +74,9 @@ namespace AMDiS {
     DegreeOfFreedom *last_edge[2] = {NULL, NULL};
     int n_neigh_periodic;
 
-    DegreeOfFreedom newDOF = -1;
-    DegreeOfFreedom lastNewDOF = -1;
-    DegreeOfFreedom firstNewDOF = -1;
+    DegreeOfFreedom *newDOF = NULL;
+    DegreeOfFreedom *lastNewDOF = NULL;
+    DegreeOfFreedom *firstNewDOF = NULL;
 
     RCNeighbourList periodicList;
 
@@ -87,23 +87,25 @@ namespace AMDiS {
 
       newDOF = refinePatch(edge, periodicList, n_neigh_periodic, bound);
 
-      if (firstNewDOF == -1)
-	firstNewDOF = newDOF;
-
-      if (lastNewDOF != -1) {
-	for (std::map<int, VertexVector*>::iterator it = mesh->getPeriodicAssociations().begin();
-	     it != mesh->getPeriodicAssociations().end(); ++it) {
-	  if (it->second) {
-	    if (((*(it->second))[edge[0][0]] == last_edge[0][0] &&
-		 (*(it->second))[edge[1][0]] == last_edge[1][0]) ||
-		((*(it->second))[edge[0][0]] == last_edge[1][0] &&
-		 (*(it->second))[edge[1][0]] == last_edge[0][0])) {
-	      (*(it->second))[lastNewDOF] = newDOF;
-	      (*(it->second))[newDOF] = lastNewDOF;
-	    }
-	  }
-	}
-      }
+      if (firstNewDOF == NULL)
+        firstNewDOF = newDOF;
+
+      if (lastNewDOF != NULL) {
+        for (std::map<int, std::vector<VertexVector*> >::iterator it = mesh->getPeriodicAssociations().begin();
+            it != mesh->getPeriodicAssociations().end(); ++it) {
+          if (!it->second.empty()) {
+            for (size_t iadmin = 0; iadmin < mesh->getNumberOfDOFAdmin(); iadmin++) {
+              if (((*(it->second)[iadmin])[edge[0][iadmin]] == last_edge[0][iadmin] &&
+                   (*(it->second)[iadmin])[edge[1][iadmin]] == last_edge[1][iadmin]) ||
+                   ((*(it->second[iadmin]))[edge[0][iadmin]] == last_edge[1][iadmin] &&
+                   (*(it->second)[iadmin])[edge[1][iadmin]] == last_edge[0][iadmin])) {
+                (*(it->second)[iadmin])[lastNewDOF[iadmin]] = newDOF[iadmin];
+                (*(it->second)[iadmin])[newDOF[iadmin]] = lastNewDOF[iadmin];
+              }
+            } // for
+          } // if
+        } // for
+      } // if lastNewDOF
       lastNewDOF = newDOF;
 
       last_edge[0] = edge[0];
@@ -114,19 +116,21 @@ namespace AMDiS {
     }
 
     if (lastNewDOF != firstNewDOF) {
-      for (std::map<int, VertexVector*>::iterator it = mesh->getPeriodicAssociations().begin();
-	   it != mesh->getPeriodicAssociations().end(); ++it) {
-	if (it->second) {
-	  if (((*(it->second))[first_edge[0][0]] == last_edge[0][0] &&
-	       (*(it->second))[first_edge[1][0]] == last_edge[1][0]) ||
-	      ((*(it->second))[first_edge[0][0]] == last_edge[1][0] &&
-	       (*(it->second))[first_edge[1][0]] == last_edge[0][0])) {
-	    (*(it->second))[lastNewDOF] = firstNewDOF;
-	    (*(it->second))[firstNewDOF] = lastNewDOF;
-	  }
-	}
-      }
-    }
+      for (std::map<int, std::vector<VertexVector*> >::iterator it = mesh->getPeriodicAssociations().begin();
+           it != mesh->getPeriodicAssociations().end(); ++it) {
+        if (!it->second.empty()) {
+          for (size_t iadmin = 0; iadmin < mesh->getNumberOfDOFAdmin(); iadmin++) {
+            if (((*(it->second)[iadmin])[first_edge[0][iadmin]] == last_edge[0][iadmin] &&
+                 (*(it->second)[iadmin])[first_edge[1][iadmin]] == last_edge[1][iadmin]) ||
+                 ((*(it->second)[iadmin])[first_edge[0][iadmin]] == last_edge[1][iadmin] &&
+                 (*(it->second)[iadmin])[first_edge[1][iadmin]] == last_edge[0][iadmin])) {
+              (*(it->second)[iadmin])[lastNewDOF[iadmin]] = firstNewDOF[iadmin];
+              (*(it->second)[iadmin])[firstNewDOF[iadmin]] = lastNewDOF[iadmin];
+            }
+          } // for
+        } // if
+      } // for
+    } // if !=
 
     return elInfo;
   }
@@ -175,9 +179,9 @@ namespace AMDiS {
   }
 
 
-  DegreeOfFreedom RefinementManager2d::refinePatch(DegreeOfFreedom *edge[2],
-						   RCNeighbourList &refineList,
-						   int n_neigh, bool bound)
+  DegreeOfFreedom* RefinementManager2d::refinePatch(DegreeOfFreedom *edge[2],
+                                                    RCNeighbourList &refineList,
+                                                    int n_neigh, bool bound)
   {
     DegreeOfFreedom *dof[3] = {NULL, NULL, NULL};
     Triangle *el =
@@ -222,8 +226,8 @@ namespace AMDiS {
       DOFAdmin* admin = const_cast<DOFAdmin*>(&mesh->getDofAdmin(iadmin));
       std::list<DOFIndexedBase*>::iterator end = admin->endDOFIndexed();
       for (std::list<DOFIndexedBase*>::iterator it = admin->beginDOFIndexed();
-	   it != end; it++)
-	(*it)->refineInterpol(refineList, n_neigh);
+           it != end; it++)
+        (*it)->refineInterpol(refineList, n_neigh);
     }