Add checking procedure for not valid 3D meshes.

AMDiS/src/parallel/ElementObjectData.h

@@ -345,18 +345,24 @@ namespace AMDiS {
     /// Returns to an element object data the appropriate vertex DOF.
     DegreeOfFreedom getVertexLocalMap(ElementObjectData &data)
     {
+      TEST_EXIT_DBG(vertexLocalMap.count(data))("Should not happen!\n");
+
       return vertexLocalMap[data];
     }
 
     /// Returns to an element object data the appropriate edge.
     DofEdge getEdgeLocalMap(ElementObjectData &data)
     {
+      TEST_EXIT_DBG(edgeLocalMap.count(data))("Should not happen!\n");
+
       return edgeLocalMap[data];
     }
 
     /// Returns to an element object data the appropriate face.
     DofFace getFaceLocalMap(ElementObjectData &data)
     {
+      TEST_EXIT_DBG(faceLocalMap.count(data))("Should not happen!\n");
+
       return faceLocalMap[data];
     }
 

AMDiS/src/parallel/MeshDistributor.cc

@@ -23,6 +23,7 @@
 #include "parallel/ParallelDebug.h"
 #include "parallel/StdMpi.h"
 #include "parallel/ParMetisPartitioner.h"
+#include "parallel/MpiHelper.h"
 #include "io/ElementFileWriter.h"
 #include "io/MacroInfo.h"
 #include "io/VtkWriter.h"
@@ -228,6 +229,12 @@ namespace AMDiS {
     updateLocalGlobalNumbering();
 
 
+    // === In 3D we have to make some test, if the resulting mesh is valid. If   ===
+    // === it is not valid, there is no possiblity yet to fix this problem, just ===
+    // === exit with an error message.                                           ===
+    
+    check3dValidMesh();
+
     // === If in debug mode, make some tests. ===
 
 #if (DEBUG != 0)
@@ -245,7 +252,7 @@ namespace AMDiS {
 #endif
 
 
-    // === Create periodic dof mapping, if there are periodic boundaries. ===
+    // === Create periodic DOF mapping, if there are periodic boundaries. ===
 
     createPeriodicMap();
 
@@ -469,6 +476,84 @@ namespace AMDiS {
   }
 
 
+  void MeshDistributor::check3dValidMesh()
+  {
+    FUNCNAME("MeshDistributor::check3dValidMesh()");
+
+    if (mesh->getDim() != 3)
+      return;
+
+    std::set<DofEdge> allEdges;
+    std::set<int> rankMacroEls;
+
+    TraverseStack stack;
+    ElInfo *elInfo = stack.traverseFirst(mesh, 0, Mesh::CALL_EL_LEVEL);
+    while (elInfo) {
+      for (int i = 0; i < 4; i++) {
+	ElementObjectData elData(elInfo->getElement()->getIndex(), i);
+	allEdges.insert(elObjects.getEdgeLocalMap(elData));
+      }
+
+      rankMacroEls.insert(elInfo->getElement()->getIndex());
+
+      elInfo = stack.traverseNext(elInfo);
+    }
+
+    bool meshIsValid = true;
+
+    for (std::set<DofEdge>::iterator it = allEdges.begin(); it != allEdges.end(); ++it) {
+      vector<ElementObjectData>& edgeEls = elObjects.getElements(*it);
+
+      TEST_EXIT_DBG(edgeEls.size() > 0)
+	("No edge %d/%d in elObjDB!\n", it->first, it->second);
+
+      std::set<int> el0, el1;
+      for (unsigned int i = 0; i < edgeEls.size(); i++)
+	if (rankMacroEls.count(edgeEls[i].elIndex))
+	  if (el0.empty())
+	    el0.insert(edgeEls[i].elIndex);
+	  else
+	    el1.insert(edgeEls[i].elIndex);
+            
+      TEST_EXIT_DBG(!el0.empty())("Should not happen!\n");
+
+      if (el1.empty())
+	continue;
+      
+      bool found = false;      
+      do {
+	found = false;
+	for (std::set<int>::iterator elIt = el0.begin(); elIt != el0.end(); ++elIt) {
+	  for (int i = 0; i < 4; i++) {
+	    int neighEl = macroElementNeighbours[*elIt][i];
+	    if (neighEl != -1 && el1.count(neighEl)) {
+	      el0.insert(neighEl);
+	      el1.erase(neighEl);
+	      found = true;
+	    }
+	  }
+	  
+	  if (found)
+	    break;
+	}
+      } while (found);
+      
+      if (!el1.empty()) {
+	meshIsValid = false;
+	break;
+      }
+    }
+
+    if (!meshIsValid) {
+      MSG("Error: mesh is not a valid 3D mesh on this rank!\n");
+    }
+
+    int foundError = !meshIsValid;
+    mpi::globalAdd(foundError);
+    TEST_EXIT(foundError == 0)("Error found on at least on rank!\n");
+  }
+
+
   void MeshDistributor::setRankDofs()
   {
     for (unsigned int i = 0; i < probStat.size(); i++) {
@@ -603,9 +688,7 @@ namespace AMDiS {
       recvAllValues = 0;
       mpiComm.Allreduce(&sendValue, &recvAllValues, 1, MPI_INT, MPI_SUM);
 
-#if (DEBUG != 0)
       MSG("Mesh changed on %d ranks!\n", recvAllValues);
-#endif
     } while (recvAllValues != 0);
 
 #if (DEBUG != 0)

AMDiS/src/parallel/MeshDistributor.h

@@ -274,6 +274,8 @@ namespace AMDiS {
      */
     void synchVector(SystemVector &vec);
 
+    void check3dValidMesh();
+
   protected:
     /** \brief
      * Determines the interior boundaries, i.e. boundaries between ranks, and stores

AMDiS/src/parallel/ParallelDebug.cc

@@ -10,14 +10,14 @@
 // See also license.opensource.txt in the distribution.
 
 
-#include "ParallelDebug.h"
-#include "MeshDistributor.h"
+#include "parallel/ParallelDebug.h"
+#include "parallel/MeshDistributor.h"
+#include "parallel/MpiHelper.h"
 #include "ProblemVec.h"
 #include "DOFVector.h"
 #include "FixVec.h"
 #include "StdMpi.h"
 #include "Debug.h"
-#include "MpiHelper.h"
 #include "io/VtkWriter.h"
 
 namespace AMDiS {