Fixed a very bad bug for 3D mesh partitioning.

AMDiS/src/Debug.cc

@@ -821,6 +821,27 @@ namespace AMDiS {
     }
 
 
+    void testDofsByCoords(DOFVector<WorldVector<double> > &coords,
+			  DofContainer &dofs0, 
+			  DofContainer &dofs1)
+    {
+      FUNCNAME("debug::testDofsByCoords()");
+
+      TEST_EXIT(dofs0.size() == dofs1.size())
+	("The dof containers have different sizes %d %d!\n",
+	 dofs0.size(), dofs1.size());
+
+      for (unsigned int i = 0; i < dofs0.size(); i++) {
+	WorldVector<double> tmp = coords[*(dofs0[i])];
+	tmp -= coords[*(dofs1[i])];
+
+	TEST_EXIT(norm(tmp) < 1e-13)
+	  ("DOFs %d and %d (i = %d) have different coords!\n", 
+	   *(dofs0[i]), *(dofs1[i]), i);
+      }
+      
+    }
+
   } // namespace debug
   
 } // namespace AMDiS

AMDiS/src/Debug.h

@@ -211,9 +211,36 @@ namespace AMDiS {
 		  const DegreeOfFreedom* dof3,
 		  DofContainer &vec);    
 
+    /** \brief
+     * Takes to vectors of DOF indices and tests if they pairwise equal. To test
+     * for equality, coordinates are checked. This makes it possible to test 
+     * internal algorithms that manipulate the mesh having two DOFs at the same
+     * geometrical position.
+     *
+     * This function does not return any value but abort the execution if it
+     * find two DOFs not to be equal.
+     *
+     * \param[in]   feSpace   Finite element space which is used to create
+     *                        coordinates.
+     * \param[in]   dofs0     First DOF container.
+     * \paran[in]   dofs1     Second DOF container.
+     */
     void testDofsByCoords(const FiniteElemSpace *feSpace,
-			  DofContainer &dofs0, DofContainer &dofs1);
+			  DofContainer &dofs0, 
+			  DofContainer &dofs1);
 
+    /** \brief
+     * Works in the same way as described above, but the calling function
+     * must provide a DOFVector which contains all DOF coordinates. This is
+     * more efficient if the function is called multiple times.
+     *
+     * \param[in]   coords    DOFVector of DOF coordinates.
+     * \param[in]   dofs0     First DOF container.
+     * \paran[in]   dofs1     Second DOF container.
+     */
+    void testDofsByCoords(DOFVector<WorldVector<double> > &coords,
+			  DofContainer &dofs0, 
+			  DofContainer &dofs1); 
   }
 }
 

AMDiS/src/RCNeighbourList.h

@@ -187,29 +187,23 @@ namespace AMDiS {
       /// This is the i-th entry in the RCNeighbourList 
       int ith;
 
-      /** \brief
-       * Only used in the coarsening module: flag is true if the coarsening 
-       * edge of the Element is the coarsening edge of the patch, otherwise 
-       * flag is false;
-       */
+      /// Only used in the coarsening module: flag is true if the coarsening 
+      /// edge of the Element is the coarsening edge of the patch, otherwise 
+      /// flag is false;
       bool flag;
 
-      /** \brief
-       * neigh[0/1] neighbour of element to the right/left in the orientation 
-       * of the edge, or a NULL pointer in the case of a boundary face (only 3d)
-       */
+      /// neigh[0/1] neighbour of element to the right/left in the orientation 
+      /// of the edge, or a NULL pointer in the case of a boundary face (only 3d)
       RCListElement* neigh[2];
 
       /// opp vertex[0/1] the opposite vertex of neigh[0/1] (only 3d)
       int oppVertex[2];
 
-      /** \brief
-       * The element type; is set during looping around the 
-       * refinement/coarsening edge; if neighbour information is produced by the
-       * traversal routines, information about the type of an element can not be
-       * accessed via el->getType() and thus has to be stored in the 
-       * RCListElement vector (only 3d)
-       */
+      /// The element type; is set during looping around the 
+      /// refinement/coarsening edge; if neighbour information is produced by the
+      /// traversal routines, information about the type of an element can not be
+      /// accessed via el->getType() and thus has to be stored in the 
+      /// RCListElement vector (only 3d)
       unsigned char elType;
     };
 

AMDiS/src/RefinementManager3d.cc

@@ -653,6 +653,7 @@ namespace AMDiS {
       edge[0] = const_cast<DegreeOfFreedom*>(el->getDof(1));
     }
 
+
 #ifdef HAVE_PARALLEL_DOMAIN_AMDIS
     Element *otherEl = NULL;
     int otherEdge = -1;
@@ -744,7 +745,7 @@ namespace AMDiS {
 
     // ============ Check for periodic boundary ============
 
-    DegreeOfFreedom *next_edge[2];
+    DegreeOfFreedom *next_edge[2] = {NULL, NULL};
     DegreeOfFreedom *first_edge[2] = {edge[0], edge[1]};
     DegreeOfFreedom *last_edge[2] = {NULL, NULL};
     int n_neigh_periodic;
@@ -759,7 +760,8 @@ namespace AMDiS {
 			       &n_neigh, &n_neigh_periodic, 
 			       periodicList);
 
-      DegreeOfFreedom newDof = refinePatch(edge, periodicList, n_neigh_periodic, bound);
+      DegreeOfFreedom newDof = 
+	refinePatch(edge, periodicList, n_neigh_periodic, bound);
 
       if (firstNewDof == -1)
 	firstNewDof = newDof;
@@ -828,27 +830,35 @@ namespace AMDiS {
       for (int i = 0; i < elInfo->getLevel(); i++) {
 	TEST_EXIT_DBG(checkIndex >= 1)("Should not happen!\n");
 	
+
 	int parentType = elInfos[checkIndex - 1]->getType();
 	int ithParentChild = 0;
 	if (elInfos[checkIndex - 1]->getElement()->getChild(1) ==
 	    elInfos[checkIndex]->getElement())
 	  ithParentChild = 1;
-	localEdgeNo = Tetrahedron::childEdge[parentType][ithParentChild][localEdgeNo];
-	
-	if (localEdgeNo >= 4) {
-	  localEdgeNo = -1;
-	  break;
-	}
+
+	// If local edge number if equal or greater to 4, its a new edge and
+	// thus cannot be part of the parent edge.
+ 	if (localEdgeNo >= 4) {
+ 	  localEdgeNo = -1;
+ 	  break;
+ 	}
 	
+	localEdgeNo = 
+	  Tetrahedron::childEdge[parentType][ithParentChild][localEdgeNo];
+
 	checkIndex--;
       }
-      
-      // If the refinement edge is part of an edge on the macro level, we must check
-      // if the refinement edge is part of an edge that must be fixed.
+
+     
+      // If the refinement edge is part of an edge on the macro level, we must 
+      // check if the refinement edge is part of an edge that must be fixed.
       if (localEdgeNo >= 0) {
-	TEST_EXIT_DBG(elInfos[checkIndex]->getLevel() == 0)("Should not happen!\n");
+	TEST_EXIT_DBG(elInfos[checkIndex]->getLevel() == 0)
+	  ("Should not happen!\n");
 	TEST_EXIT_DBG(elInfos[checkIndex]->getElement()->getIndex() ==
-		      elInfo->getMacroElement()->getIndex())("Should not happen!\n");
+		      elInfo->getMacroElement()->getIndex())
+	  ("Should not happen!\n");
 	TEST_EXIT_DBG(localEdgeNo <= 5)("Should not happen!\n");
 	
 	for (unsigned int i = 0; i < FixRefinementPatch::connectedEdges.size(); i++) {

AMDiS/src/io/ElementFileWriter.cc

@@ -138,13 +138,15 @@ namespace AMDiS {
     if (writeVtkFormat || writeVtkVectorFormat) {
       TEST_EXIT(mesh)("no mesh\n");
 
-      writeVtkValues(fn, vtkExt,-1,writeVtkVectorFormat==1);
+      writeVtkValues(fn, vtkExt, -1, writeVtkVectorFormat == 1);
       MSG("VTK file written to %s\n", (fn + vtkExt).c_str());		     
     }
 
-    if (writeParaViewAnimation) {
-      VtkWriter::updateAnimationFile(adaptInfo, (fn + vtkExt), &paraViewAnimationFrames, (filename + pvdExt));
-    }
+    if (writeParaViewAnimation)
+      VtkWriter::updateAnimationFile(adaptInfo, 
+				     (fn + vtkExt), 
+				     &paraViewAnimationFrames, 
+				     (filename + pvdExt));
   }
 
 

AMDiS/src/parallel/MatrixNnzStructure.cc

@@ -237,7 +237,7 @@ namespace AMDiS {
 	stdMpi.recv(*it);
       stdMpi.startCommunication();
       
-      
+     
       // === Evaluate the nnz structure this rank got from other ranks and add ===
       // === it to the PETSc nnz data structure.                               ===
       

AMDiS/src/parallel/MeshDistributor.cc

@@ -229,7 +229,7 @@ namespace AMDiS {
       Parameters::get("dbg->write part mesh", writePartMesh);
 
       if (writePartMesh > 0) {
-	debug::writeElementIndexMesh(mesh, debugOutputDir + "elementIndex.vtu");
+	debug::writeElementIndexMesh(mesh, debugOutputDir + "elementIndex");
 	ParallelDebug::writePartitioning(*this, debugOutputDir + "part.vtu");
       }
     }
@@ -295,9 +295,10 @@ namespace AMDiS {
     // === 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)
@@ -573,15 +574,19 @@ namespace AMDiS {
     if (mesh->getDim() != 3)
       return;
 
+    // === Create set of all edges and all macro element indices in rank's ===
+    // === subdomain.                                                      ===
+
     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++) {
+      for (int i = 0; i < 6; i++) {
 	ElementObjectData elData(elInfo->getElement()->getIndex(), i);
-	allEdges.insert(elObjDb.getEdgeLocalMap(elData));
+	DofEdge e = elObjDb.getEdgeLocalMap(elData);
+	allEdges.insert(e);
       }
 
       rankMacroEls.insert(elInfo->getElement()->getIndex());
@@ -590,15 +595,27 @@ namespace AMDiS {
     }
 
 
+    // === Reset fixed refinement edges and assume the mesh is valid. Search ===
+    // === for the opposite.                                                 ===
+
     FixRefinementPatch::connectedEdges.clear();
     bool valid3dMesh = true;
 
+    // Traverse all edges
     for (std::set<DofEdge>::iterator it = allEdges.begin(); it != allEdges.end(); ++it) {
+      // For each edge get all macro elements that contain this edge.
       vector<ElementObjectData>& edgeEls = elObjDb.getElements(*it);
 
       TEST_EXIT_DBG(edgeEls.size() > 0)
 	("No edge %d/%d in elObjDB!\n", it->first, it->second);
 
+
+      // We have now to check, if the current edge connects two macro elements,
+      // which are otherwise not connected. The basic idea to check this is very
+      // simple: We take the first macro element in rank's subdomain that contain
+      // this edge and add it to the set variable "el0". All other macro elements
+      // which share this edge are added to "el1". 
+
       std::set<int> el0, el1;
       map<int, int> edgeNoInEl;
 
@@ -612,9 +629,11 @@ namespace AMDiS {
 	  edgeNoInEl[edgeEls[i].elIndex] = edgeEls[i].ithObject;
 	}
       }
-            
+           
       TEST_EXIT_DBG(!el0.empty())("Should not happen!\n");
 
+      // If el1 is empty, there is only on macro element in the mesh which
+      // contains this edge. Hence, we can continue with checking another edge.
       if (el1.empty())
 	continue;
       
@@ -649,6 +668,7 @@ namespace AMDiS {
 	    pair<Element*, int> edge1 = 
 	      make_pair(elObjDb.getElementPtr(*elIt1), edgeNoInEl[*elIt1]);
 
+#if (DEBUG != 0)
 	    DofEdge dofEdge0 = edge0.first->getEdge(edge0.second);
 	    DofEdge dofEdge1 = edge1.first->getEdge(edge1.second);
 
@@ -656,12 +676,15 @@ namespace AMDiS {
 	    mesh->getDofIndexCoords(dofEdge0.first, feSpaces[0], c0);
 	    mesh->getDofIndexCoords(dofEdge0.second, feSpaces[0], c1);
 
-	    /*
-	    MSG_DBG("FOUND EDGE: %d %d with coords %f %f %f   %f %f %f\n",
-		    dofEdge0.first, dofEdge0.second, 
-		    c0[0], c0[1], c0[2],
-		    c1[0], c1[1], c1[2]);
-	    */
+	    MSG("FOUND EDGE %d/%d <-> %d/%d\n", 
+		edge0.first->getIndex(), edge0.second,
+		edge1.first->getIndex(), edge1.second);
+
+	    MSG("FOUND EDGE: %d %d with coords %f %f %f   %f %f %f\n",
+		dofEdge0.first, dofEdge0.second, 
+		c0[0], c0[1], c0[2],
+		c1[0], c1[1], c1[2]);
+#endif
 
 	    TEST_EXIT_DBG(dofEdge0 == dofEdge1)("Should noth happen!\n");
 
@@ -1022,9 +1045,9 @@ namespace AMDiS {
       repartitionMesh();
       nMeshChangesAfterLastRepartitioning = 0;
     } else {
-      MSG_DBG("Repartitioning not tried because tryRepartitioning = %d   repartitioningAllowed = %d   nMeshChange = %d   repartitionIthChange = %d\n",
-	      tryRepartition, repartitioningAllowed, 
-	      nMeshChangesAfterLastRepartitioning, repartitionIthChange);
+      MSG("Repartitioning not tried because tryRepartitioning = %d   repartitioningAllowed = %d   nMeshChange = %d   repartitionIthChange = %d\n",
+	  tryRepartition, repartitioningAllowed, 
+	  nMeshChangesAfterLastRepartitioning, repartitionIthChange);
     }
 
     // === Print imbalance factor. ===
@@ -1322,7 +1345,6 @@ namespace AMDiS {
     // === Run mesh partitioner to calculate a new mesh partitioning.  ===
 
     partitioner->setLocalGlobalDofMap(&(dofMap[feSpaces[0]].getMap()));
-
     bool partitioningSucceed = 
       partitioner->partition(elemWeights, ADAPTIVE_REPART);
     if (!partitioningSucceed) {
@@ -1331,6 +1353,7 @@ namespace AMDiS {
       return;
     }
 
+
     // In the case the partitioner does not create a new mesh partition, return
     // without and changes.
     if (!partitioner->meshChanged()) {
@@ -1339,6 +1362,7 @@ namespace AMDiS {
       return;
     }
 
+
     TEST_EXIT_DBG(!(partitioner->getSendElements().size() == mesh->getMacroElements().size() && 
 		    partitioner->getRecvElements().size() == 0))
       ("Partition is empty, should not happen!\n");
@@ -1548,6 +1572,7 @@ namespace AMDiS {
     // In 3D we have to make some test, if the resulting mesh is valid.
     check3dValidMesh();
 
+
     MPI::COMM_WORLD.Barrier();
     MSG("Mesh repartitioning needed %.5f seconds\n", MPI::Wtime() - timePoint);   
   }
@@ -1695,6 +1720,7 @@ namespace AMDiS {
   
     lastMeshChangeIndex = mesh->getChangeIndex();
 
+
 #if (DEBUG != 0)
     ParallelDebug::testDofContainerCommunication(*this);
 
@@ -1718,8 +1744,8 @@ namespace AMDiS {
       }
     }
 
-    debug::writeElementIndexMesh(mesh, debugOutputDir + "elementIndex-" + 
-				 lexical_cast<string>(mpiRank) + ".vtu");
+//     debug::writeElementIndexMesh(mesh, debugOutputDir + "elementIndex-" + 
+// 				 lexical_cast<string>(mpiRank) + ".vtu");
     ParallelDebug::writeDebugFile(feSpaces[feSpaces.size() - 1], dofMap,
 				  debugOutputDir + "mpi-dbg", "dat");
     debug::testSortedDofs(mesh, elMap);

AMDiS/src/parallel/MeshDistributor.h

@@ -262,6 +262,12 @@ namespace AMDiS {
     /// DOFVector, but instead sends all values of all DOFVectors all at once.
     void synchVector(SystemVector &vec, int level = 0);
 
+    /// In 3D, a subdomain may not be a valid AMDiS mesh if it contains two
+    /// parts which are only connected by an edge. In this case, the standard
+    /// refinement algorithm does not work correctly, as two elements connected
+    /// only on one edge are not neighours by definition. This functions checks
+    /// for this situation and fix the problem. For this, the mesh is search for
+    /// all edges connecting two elements that are otherwise not connected.
     void check3dValidMesh();
 
     void setBoundaryDofRequirement(Flag flag)

AMDiS/src/parallel/MeshManipulation.cc

@@ -78,6 +78,12 @@ namespace AMDiS {
     }
 
 
+#if (DEBUG != 0)
+    DOFVector<WorldVector<double> > coords(feSpace, "dofCorrds");
+    feSpace->getMesh()->getDofIndexCoords(feSpace, coords);
+#endif
+
+
     // === Traverse all new elements and connect them to the overall mesh     ===
     // === structure by deleting all double DOFs on macro element's vertices, ===
     // === edges and faces.                                                   ===
@@ -147,9 +153,10 @@ namespace AMDiS {
 
 #if (DEBUG != 0)
 	    if (feSpaces.size() == 1)
-	      debug::testDofsByCoords(feSpace, dofs0, dofs1);
+	      debug::testDofsByCoords(coords, dofs0, dofs1);
 	    else 
-	      TEST_EXIT_DBG(dofs0.size() == dofs1.size())("Should not happen!\n");
+	      TEST_EXIT_DBG(dofs0.size() == dofs1.size())
+		("Should not happen!\n");
 #endif
 
  	    for (unsigned int i = 0; i < dofs0.size(); i++) {
@@ -196,9 +203,10 @@ namespace AMDiS {
 
 #if (DEBUG != 0)
 	    if (feSpaces.size() == 1)
-	      debug::testDofsByCoords(feSpace, dofs0, dofs1);
+	      debug::testDofsByCoords(coords, dofs0, dofs1);
 	    else 
-	      TEST_EXIT_DBG(dofs0.size() == dofs1.size())("Should not happen!\n");
+	      TEST_EXIT_DBG(dofs0.size() == dofs1.size())
+		("Should not happen!\n");
 #endif
 
  	    for (unsigned int i = 0; i < dofs0.size(); i++) {
@@ -206,7 +214,8 @@ namespace AMDiS {
 	      dofPosIndex[dofs0[i]] = dofGeoIndex1[i];
 
 	      TEST_EXIT_DBG(dofGeoIndex0[i] == dofGeoIndex1[i])
-		("Should not happen: %d %d\n", dofGeoIndex0[i], dofGeoIndex1[i]);
+		("Should not happen: %d %d\n", 
+		 dofGeoIndex0[i], dofGeoIndex1[i]);
 	    }
 
 	    break;

AMDiS/src/parallel/PetscProblemStat.cc

@@ -54,6 +54,7 @@ namespace AMDiS {
 
     tmp = "";
     Parameters::get(nameStr + "->solver->petsc prefix", tmp);
+    MSG("SET PREFIX \"%s\"\n", tmp.c_str());
     petscSolver->setKspPrefix(tmp);    
   }
 

AMDiS/src/parallel/PetscSolver.h

@@ -113,6 +113,7 @@ namespace AMDiS {
 
     void setKspPrefix(std::string s)
     {
+      MSG("======>>>>>> SET TO %s\n", s.c_str());
       kspPrefix = s;
     }
 

AMDiS/src/parallel/PetscSolverFeti.cc

@@ -1339,7 +1339,10 @@ namespace AMDiS {
     createFetiData();
     
     // === Create matrices for the FETI-DP method. ===
-   
+
+    if (printTimings) {
+      MPI::COMM_WORLD.Barrier();
+    }
     double wtime = MPI::Wtime();
   
     subdomain->fillPetscMatrix(mat);
@@ -1616,8 +1619,13 @@ namespace AMDiS {
     // tmp_primal0 = f_Pi - M_PiB * inv(K_BB) * f_B
     VecAXPBY(tmp_primal0, 1.0, -1.0, subdomain->getRhsCoarseSpace());
 
+    double wtime2 = MPI::Wtime();
     // tmp_primal0 = inv(S_PiPi) (f_Pi - M_PiB * inv(K_BB) * f_B)
     KSPSolve(ksp_schur_primal, tmp_primal0, tmp_primal0);
+    if (printTimings) {
+      MSG("FETI-DP timing 09a: %.5f seconds (create rhs vector)\n", 
+	  MPI::Wtime() - wtime2);
+    }
 
     //
     MatMult(subdomain->getMatIntCoarse(), tmp_primal0, tmp_b0);

AMDiS/src/parallel/PetscSolverGlobalMatrix.cc

@@ -97,6 +97,7 @@ namespace AMDiS {
     KSPSetOperators(kspInterior, matIntInt, matIntInt, SAME_NONZERO_PATTERN); 
     KSPSetTolerances(kspInterior, 0.0, 1e-8, PETSC_DEFAULT, PETSC_DEFAULT);
     KSPSetType(kspInterior, KSPBCGS);
+    MSG("SET PREFIX TOOO: %s\n", kspPrefix.c_str());
     KSPSetOptionsPrefix(kspInterior, kspPrefix.c_str());
     KSPSetFromOptions(kspInterior);
 
@@ -125,7 +126,7 @@ namespace AMDiS {
     FUNCNAME("PetscSolverGlobalMatrix::fillPetscMatrixWithCoarseSpace()");
 
     TEST_EXIT_DBG(interiorMap)("Should not happen!\n");
-    
+
     vector<const FiniteElemSpace*> feSpaces = getFeSpaces(mat);
 
     int nRowsRankInterior = interiorMap->getRankDofs();
@@ -155,7 +156,8 @@ namespace AMDiS {
 		   0, nnzInterior.dnnz, 0, nnzInterior.onnz, 
 		   &matIntInt);
     }
-    
+
+   
     if (coarseSpaceMap) {
       int nRowsRankCoarse = coarseSpaceMap->getRankDofs();
       int nRowsOverallCoarse = coarseSpaceMap->getOverallDofs();
@@ -179,6 +181,7 @@ namespace AMDiS {
 		   &matIntCoarse);
     }
 
+
     // === Prepare traverse of sequentially created matrices. ===
 
     using mtl::tag::row; using mtl::tag::nz; using mtl::begin; using mtl::end;
@@ -256,13 +259,8 @@ namespace AMDiS {
 
 	    if (colsOther.size()) {
 	      if (subdomainLevel == 0) {
-		for (unsigned int k = 0; k < colsOther.size(); k++) {
-		  int local = interiorMap->getMatIndex(j, colsOther[k]) - interiorMap->getStartDofs();
-		  if (rowIndex == 0 && local == 9297)
-		    MSG("FOUND: %d %d %d\n", colsOther[k], j, interiorMap->getMatIndex(j, colsOther[k]));
-
+		for (unsigned int k = 0; k < colsOther.size(); k++)
 		  colsOther[k] = interiorMap->getMatIndex(j, colsOther[k]);
-		}
 	      } else {
 		for (unsigned int k = 0; k < colsOther.size(); k++)
 		  colsOther[k] = 
@@ -304,11 +302,13 @@ namespace AMDiS {
       }
     }
 
+
     // === Start global assembly procedure. ===
 
     MatAssemblyBegin(matIntInt, MAT_FINAL_ASSEMBLY);
     MatAssemblyEnd(matIntInt, MAT_FINAL_ASSEMBLY);
 
+
     if (coarseSpaceMap) {
       MatAssemblyBegin(matCoarseCoarse, MAT_FINAL_ASSEMBLY);
       MatAssemblyEnd(matCoarseCoarse, MAT_FINAL_ASSEMBLY);
@@ -320,6 +320,7 @@ namespace AMDiS {
       MatAssemblyEnd(matCoarseInt, MAT_FINAL_ASSEMBLY);
     }
 
+
     // === Remove Dirichlet BC DOFs. ===
 
     //    removeDirichletBcDofs(mat);