From aa1ec307be91c457b46a681a47254b5d88bf8fee Mon Sep 17 00:00:00 2001
From: Thomas Witkowski <thomas.witkowski@gmx.de>
Date: Mon, 13 Feb 2012 06:56:22 +0000
Subject: [PATCH] Fixed some problems for mesh redistribution with mixed finite
 elements.

---
 AMDiS/src/MeshStructure.cc            | 21 +++++++++++---
 AMDiS/src/parallel/MeshDistributor.cc | 40 +++++++++++++++++----------
 2 files changed, 42 insertions(+), 19 deletions(-)

diff --git a/AMDiS/src/MeshStructure.cc b/AMDiS/src/MeshStructure.cc
index 8e03083a..7bb02f3a 100644
--- a/AMDiS/src/MeshStructure.cc
+++ b/AMDiS/src/MeshStructure.cc
@@ -428,8 +428,13 @@ namespace AMDiS {
   
     const FiniteElemSpace *feSpace = vec->getFeSpace();
     Mesh *mesh = feSpace->getMesh();
+    int nVertexPreDofs = feSpace->getAdmin()->getNumberOfPreDofs(VERTEX);
     bool feSpaceHasNonVertexDofs = (feSpace->getBasisFcts()->getDegree() > 1);
     values.clear();
+#if (DEBUG != 0)    
+    // In debug mode we add the macro element index to the value code.
+    values.push_back(static_cast<double>(macroElIndex));
+#endif
 
     ElementDofIterator elDofIter(feSpace);
     TraverseStack stack;
@@ -439,14 +444,14 @@ namespace AMDiS {
       // For the macro element the mesh structure code stores all vertex values.
       if (elInfo->getLevel() == 0)
 	for (int i = 0; i < mesh->getGeo(VERTEX); i++)
-	  values.push_back((*vec)[elInfo->getElement()->getDof(i, 0)]);
+	  values.push_back((*vec)[elInfo->getElement()->getDof(i, nVertexPreDofs)]);
 
       if (!elInfo->getElement()->isLeaf()) {
 	// If no leaf element store the value of the "new" DOF that is created
 	// by bisectioning of this element.
 
 	DegreeOfFreedom dof0 = 
-	  elInfo->getElement()->getChild(0)->getDof(mesh->getDim(), 0);
+	  elInfo->getElement()->getChild(0)->getDof(mesh->getDim(), nVertexPreDofs);
 	values.push_back((*vec)[dof0]);
       } else {
 	// If leaf element store all non vertex values of this element, thus
@@ -477,7 +482,14 @@ namespace AMDiS {
     const FiniteElemSpace *feSpace = vec->getFeSpace();
     Mesh *mesh = feSpace->getMesh();
     bool feSpaceHasNonVertexDofs = (feSpace->getBasisFcts()->getDegree() > 1);
+    int nVertexPreDofs = feSpace->getAdmin()->getNumberOfPreDofs(VERTEX);
     unsigned int counter = 0;
+#if (DEBUG != 0)    
+    TEST_EXIT(static_cast<int>(values[0]) == macroElIndex)
+      ("Value structure code was created for macro element index %d, but should be set to macro element index %d\n",
+       static_cast<int>(values[0]), macroElIndex);
+    counter++;
+#endif
 
     TEST_EXIT_DBG(static_cast<int>(values.size()) >= mesh->getGeo(VERTEX))
       ("Should not happen!\n");
@@ -490,14 +502,15 @@ namespace AMDiS {
       // For the macro element all vertex nodes are set first.
       if (elInfo->getLevel() == 0)
 	for (int i = 0; i < mesh->getGeo(VERTEX); i++)
-	  (*vec)[elInfo->getElement()->getDof(i, 0)] = values[counter++];      
+	  (*vec)[elInfo->getElement()->getDof(i, nVertexPreDofs)] = 
+	    values[counter++];
 
       if (!elInfo->getElement()->isLeaf()) {
 	// If no leaf element set the value of the "new" DOF that is created
 	// by bisectioning of this element.
 	TEST_EXIT_DBG(counter < values.size())("Should not happen!\n");
 	
-	(*vec)[elInfo->getElement()->getChild(0)->getDof(mesh->getDim(), 0)] =
+	(*vec)[elInfo->getElement()->getChild(0)->getDof(mesh->getDim(), nVertexPreDofs)] =
 	  values[counter++];      
       } else {
 	// On leaf elements set all non vertex values (thus DOFs of higher order
diff --git a/AMDiS/src/parallel/MeshDistributor.cc b/AMDiS/src/parallel/MeshDistributor.cc
index 2d67e07a..dfe708d1 100644
--- a/AMDiS/src/parallel/MeshDistributor.cc
+++ b/AMDiS/src/parallel/MeshDistributor.cc
@@ -1154,9 +1154,10 @@ namespace AMDiS {
 
     // In the case the partitioner does not create a new mesh partition, return
     // without and changes.
-    if (!partitioner->meshChanged())
+    if (!partitioner->meshChanged()) {
+      MSG("Mesh partition does not create a new partition!\n");
       return;
-
+    }
 
     // === Create map that maps macro element indices to pointers to the ===
     // === macro elements.                                               ===
@@ -1264,24 +1265,14 @@ namespace AMDiS {
     for (map<int, vector<int> >::iterator it = partitioner->getRecvElements().begin();
 	 it != partitioner->getRecvElements().end(); ++it) {
       MeshCodeVec &recvCodes = stdMpi.getRecvData()[it->first];
-      vector<vector<double> > &recvValues = stdMpi2.getRecvData()[it->first];
       int i = 0;
-      int j = 0;
-      
+
       TEST_EXIT_DBG(recvCodes.size() == it->second.size())
 	("Should not happen!\n");
 
       for (vector<int>::iterator elIt = it->second.begin();
-	   elIt != it->second.end(); ++elIt) {
-	recvCodes[i].fitMeshToStructure(mesh, refineManager, false, *elIt);
-
-	for (unsigned int k = 0; k < interchangeVectors.size(); k++)
-	  recvCodes[i].setMeshStructureValues(*elIt, 
-					      interchangeVectors[k], 
-					      recvValues[j++]);
-	
-	i++;
-      }
+	   elIt != it->second.end(); ++elIt)
+	recvCodes[i++].fitMeshToStructure(mesh, refineManager, false, *elIt);
     }
 
 
@@ -1343,6 +1334,25 @@ namespace AMDiS {
     updateInteriorBoundaryInfo();
     updateLocalGlobalNumbering();
 
+    
+    // === After mesh adaption, set values of interchange vectors. ===
+    
+    for (map<int, vector<int> >::iterator it = partitioner->getRecvElements().begin();
+	 it != partitioner->getRecvElements().end(); ++it) {
+      MeshCodeVec &recvCodes = stdMpi.getRecvData()[it->first];
+      vector<vector<double> > &recvValues = stdMpi2.getRecvData()[it->first];
+      int i = 0, j = 0;
+      for (vector<int>::iterator elIt = it->second.begin();
+	   elIt != it->second.end(); ++elIt) {
+	for (unsigned int k = 0; k < interchangeVectors.size(); k++)
+	  recvCodes[i].setMeshStructureValues(*elIt, 
+					      interchangeVectors[k], 
+					      recvValues[j++]);	
+	i++;
+      }
+    }
+
+
     // === Update periodic mapping, if there are periodic boundaries. ===
 
     createPeriodicMap();
-- 
GitLab