Merge branch 'issue/periodic_bc' into 'master'

Periodic BC implemented by Siqi Ling added to this branch

Rewritten periodic BC, to be compatible in sequential and parallel and to allow mixed finite element spaces.

See merge request !9

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);
     }
 
 
@@ -231,18 +235,18 @@ namespace AMDiS {
       // === If there should be no dof information on interior leaf elements ===
       // === remove dofs from edges and the centers of parents.              ===
       if (mesh->getNumberOfDofs(EDGE)) {
-	int node = mesh->getNode(EDGE);
+        int node = mesh->getNode(EDGE);
 
-	// === The only DOF that can be freed is that in the refinement edge; all ===
-	// === other DOFs are handed on the children.                             ===
+        // === The only DOF that can be freed is that in the refinement edge; all ===
+        // === other DOFs are handed on the children.                             ===
 
-	mesh->freeDof(const_cast<DegreeOfFreedom*>(el->getDof(node+2)), EDGE);
+        mesh->freeDof(const_cast<DegreeOfFreedom*>(el->getDof(node+2)), EDGE);
       }
       if (mesh->getNumberOfDofs(EDGE) || mesh->getNumberOfDofs(CENTER))
-	refineList.removeDOFParents(n_neigh);
+        refineList.removeDOFParents(n_neigh);
     }
 
-    return dof[0][0];
+    return dof[0];
   }
 
 

AMDiS/src/RefinementManager2d.h

@@ -58,12 +58,12 @@ namespace AMDiS {
      *  looping around the refinement edge, otherwise  0
      */
     void getRefinePatch(ElInfo **elInfo, DegreeOfFreedom *edge[2], int dir,
-			RCNeighbourList &refineList, int *n_neigh);
+                        RCNeighbourList &refineList, int *n_neigh);
 
     /// Refines all elements in the patch.
-    DegreeOfFreedom refinePatch(DegreeOfFreedom *edge[2],
-				RCNeighbourList &refineList,
-				int n_neigh, bool bound);
+    DegreeOfFreedom* refinePatch(DegreeOfFreedom *edge[2],
+                                 RCNeighbourList &refineList,
+                                 int n_neigh, bool bound);
 
     /// Implements RefinementManager::refineFunction.
     ElInfo* refineFunction(ElInfo* elInfo);

AMDiS/src/RefinementManager3d.cc

@@ -330,9 +330,9 @@ namespace AMDiS {
   }
 
 
-  DegreeOfFreedom RefinementManager3d::refinePatch(DegreeOfFreedom *edge[2],
-						   RCNeighbourList &refineList,
-						   int n_neigh, bool bound)
+  DegreeOfFreedom* RefinementManager3d::refinePatch(DegreeOfFreedom *edge[2],
+                                                    RCNeighbourList &refineList,
+                                                    int n_neigh, bool bound)
   {
     Tetrahedron *el =
       dynamic_cast<Tetrahedron*>(const_cast<Element*>(refineList.getElement(0)));
@@ -402,7 +402,7 @@ namespace AMDiS {
       mesh->incrementNumberOfFaces(2 * n_neigh);
     }
 
-    return dof[0][0];
+    return dof[0];
   }
 
 
@@ -770,8 +770,8 @@ namespace AMDiS {
     DegreeOfFreedom *last_edge[2] = {NULL, NULL};
     int n_neigh_periodic = 0;
 
-    DegreeOfFreedom lastNewDof = -1;
-    DegreeOfFreedom firstNewDof = -1;
+    DegreeOfFreedom *lastNewDof = NULL;
+    DegreeOfFreedom *firstNewDof = NULL;
 
     RCNeighbourList periodicList;
 
@@ -780,26 +780,28 @@ namespace AMDiS {
 			       &n_neigh, &n_neigh_periodic,
 			       periodicList);
 
-      DegreeOfFreedom 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;
-	    }
-	  }
-	}
-      }
+      DegreeOfFreedom *newDof =
+        refinePatch(edge, periodicList, n_neigh_periodic, bound);
+
+      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;