Repartitioning for higher order elements works now.

AMDiS/src/Element.cc

@@ -34,12 +34,11 @@ namespace AMDiS {
     newCoord = NULL;
     elementData = NULL;
     mark = 0;
-    dofValid = true;
 
     if (mesh)
       createNewDofPtrs();
     else
-      mesh = NULL;    
+      dof = NULL;    
   }
 
 
@@ -302,6 +301,8 @@ namespace AMDiS {
   void Element::changeDofs1(const DOFAdmin* admin, std::vector<int>& newDofIndex,
 			    int n0, int nd0, int nd, int pos)
   {
+    FUNCNAME("Element::changeDofs1()");
+
     DegreeOfFreedom *ldof = dof[n0 + pos] + nd0;
     for (int j = 0; j < nd; j++) {
       int k = ldof[j];
@@ -396,9 +397,10 @@ namespace AMDiS {
     int nodes = mesh->getNumberOfNodes();
     int j = 0;
     SerUtil::serialize(out, nodes);
+    int dofValid = (dof != NULL ? 1 : 0);
     SerUtil::serialize(out, dofValid);
    
-    if (dofValid) {
+    if (dof != NULL) {
       for (int pos = 0; pos <= dim; pos++) {
 	GeoIndex position = INDEX_OF_DIM(pos, dim);
 	int ndof = 0;
@@ -493,10 +495,11 @@ namespace AMDiS {
     // read dofs
     int nodes;
     SerUtil::deserialize(in, nodes);
+    int dofValid = 0;
     SerUtil::deserialize(in, dofValid);
 
-    TEST_EXIT_DBG(nodes > 0)("Should not happen!\n");      
-    dof = new DegreeOfFreedom*[nodes]; 
+    TEST_EXIT_DBG(nodes > 0)("Should not happen!\n");    
+    dof = (dofValid ? new DegreeOfFreedom*[nodes] : NULL);
      
     if (dofValid) {
       for (int i = 0; i < nodes; i++) {

AMDiS/src/Element.h

@@ -119,7 +119,7 @@ namespace AMDiS {
     /// Returns \ref dof[i][j] which is the j-th DOF of the i-th node of Element.
     inline DegreeOfFreedom getDof(int i, int j) const 
     { 
-      TEST_EXIT_DBG(dofValid)("DOFs are not valid in element %d!\n", index);
+      TEST_EXIT_DBG(dof != NULL)("DOFs are not valid in element %d!\n", index);
 
       return dof[i][j];
     }
@@ -127,7 +127,7 @@ namespace AMDiS {
     /// Returns \ref dof[i] which is a pointer to the DOFs of the i-th node.
     inline const DegreeOfFreedom* getDof(int i) const 
     {
-      TEST_EXIT_DBG(dofValid)("DOFs are not valid in element %d!\n", index);
+      TEST_EXIT_DBG(dof != NULL)("DOFs are not valid in element %d!\n", index);
 
       return dof[i];
     }
@@ -135,7 +135,7 @@ namespace AMDiS {
     /// Returns a pointer to the DOFs of this Element
     inline const DegreeOfFreedom** getDof() const 
     {
-      TEST_EXIT_DBG(dofValid)("DOFs are not valid in element %d!\n", index);
+      TEST_EXIT_DBG(dof != NULL)("DOFs are not valid in element %d!\n", index);
 
       return const_cast<const DegreeOfFreedom**>(dof);
     }
@@ -146,11 +146,6 @@ namespace AMDiS {
       return mesh; 
     }
 
-    inline void setDofValid(bool b)
-    {
-      dofValid = b;
-    }
-
     /// Returns \ref elementData's error estimation, if Element is a leaf element
     /// and has leaf data. 
     inline double getEstimation(int row) const
@@ -276,6 +271,15 @@ namespace AMDiS {
       dof[pos] = p;
     }
 
+    /// Delets the main DOF pointer. The caller must ensure that the DOFs are
+    /// either freed correctly elsewhere or that they are still used by other
+    /// elements.
+    void delDofPtr()
+    {
+      delete [] dof;
+      dof = NULL;
+    }
+
     /// Checks whether Element is a leaf element and whether it has leaf data.
     /// If the checks don't fail, leaf data's error estimation is set to est.
     inline void setEstimation(double est, int row)
@@ -336,10 +340,8 @@ namespace AMDiS {
      * \{ 
      */
 
-    /** \brief
-     * orient the vertices of edges/faces.
-     * Used by Estimator for the jumps => same quadrature nodes from both sides!
-     */
+    /// orient the vertices of edges/faces.
+    /// Used by Estimator for the jumps => same quadrature nodes from both sides!
     virtual const FixVec<int,WORLD>& 
       sortFaceIndices(int face, FixVec<int, WORLD> *vec) const = 0;
 
@@ -466,7 +468,8 @@ namespace AMDiS {
     }
 
     /// Coarsens Element's leaf data
-    inline void coarsenElementData(Element* child1, Element* child2, int elType = 0) 
+    inline void coarsenElementData(Element* child1, Element* child2, 
+				   int elType = 0) 
     {
       ElementData *childData;
       childData = child1->getElementData();
@@ -549,10 +552,8 @@ namespace AMDiS {
     void changeDofs2(int n0, int nd0, int nd, int pos);
 
   protected:
-    /** \brief
-     * Pointers to the two children of interior elements of the tree. Pointers
-     * to NULL for leaf elements.
-     */
+    /// Pointers to the two children of interior elements of the tree. Pointers
+    /// to NULL for leaf elements.
     Element *child[2];
 
     /** \brief
@@ -595,19 +596,9 @@ namespace AMDiS {
     /// Pointer to Element's leaf data
     ElementData* elementData;
 
-    /** \brief
-     * Is true, if the DOF pointers are valid. In sequential computations this
-     * is always the case. In parallel computations all macro elements are stored
-     * in memory though not all of them are part of the mesh (because they are owned
-     * by another rank). In this case, there are no valid DOFs on the element. 
-     */
-    bool dofValid;
-
-    /** \brief
-     * This map is used for deletion of all DOFs of all elements of a mesh. Once
-     * a DOF-vector (all DOFS at a node, edge, etc.) is deleted, its address is
-     * added to this map to note not to delete it a second time.
-     */
+    /// This map is used for deletion of all DOFs of all elements of a mesh. Once
+    /// a DOF-vector (all DOFS at a node, edge, etc.) is deleted, its address is
+    /// added to this map to note not to delete it a second time.
     static std::map<DegreeOfFreedom*, bool> deletedDOFs;
 
     friend class Mesh;

AMDiS/src/Lagrange.cc

@@ -1158,7 +1158,8 @@ namespace AMDiS {
     /*  newest vertex of child[0] and child[1]                                  */
     /****************************************************************************/
 
-    DegreeOfFreedom cdof = el->getChild(0)->getDof(node + 2, n0);  /*      newest vertex is DIM */
+    // newest vertex is DIM 
+    DegreeOfFreedom cdof = el->getChild(0)->getDof(node + 2, n0);
     (*drv)[cdof] = (*drv)[pdof[5]];
 
     node = drv->getFeSpace()->getMesh()->getNode(EDGE);        

AMDiS/src/Mesh.cc

@@ -269,7 +269,7 @@ namespace AMDiS {
   }
 
 
-  void Mesh::removeMacroElements(std::set<MacroElement*>& macros,
+  void Mesh::removeMacroElements(std::set<MacroElement*>& delMacros,
 				 vector<const FiniteElemSpace*>& feSpaces) 
   {
     FUNCNAME("Mesh::removeMacroElement()");
@@ -319,7 +319,7 @@ namespace AMDiS {
       // If the current mesh macro element should not be deleted, i.e., it is not
       // a member of the list of macro elements to be deleted, is is inserted to
       // the new macro element list.
-      if (macros.find(*elIter) == macros.end())
+      if (delMacros.find(*elIter) == delMacros.end())
 	newMacroElements.push_back(*elIter);     
     }
 
@@ -332,8 +332,8 @@ namespace AMDiS {
     // === neighbours of some other macro elements. Furtheremore, delete the  ===
     // === whole element hierarchie structure of the macro element.           ===
     
-    for (std::set<MacroElement*>::iterator macroIt = macros.begin();
-	 macroIt != macros.end(); ++macroIt) {
+    for (std::set<MacroElement*>::iterator macroIt = delMacros.begin();
+	 macroIt != delMacros.end(); ++macroIt) {
 
       // Go through all neighbours of the macro element and remove this macro
       // element to be neighbour of some other macro element.
@@ -343,26 +343,24 @@ namespace AMDiS {
 	    if ((*macroIt)->getNeighbour(i)->getNeighbour(j) == *macroIt)
 	      (*macroIt)->getNeighbour(i)->setNeighbour(j, NULL);
 
+      Element *mel = (*macroIt)->getElement();
       // Delete element hierarchie
-      if (!(*macroIt)->getElement()->isLeaf()) {
-	delete (*macroIt)->getElement()->getChild(0);
-	delete (*macroIt)->getElement()->getChild(1);
+      if (!(mel->isLeaf())) {
+	delete mel->getChild(0);
+	delete mel->getChild(1);
 
-	(*macroIt)->getElement()->child[0] = NULL;
-	(*macroIt)->getElement()->child[1] = NULL;
+	mel->child[0] = NULL;
+	mel->child[1] = NULL;
 
-	(*macroIt)->getElement()->setElementData(elementDataPrototype->clone()); 
+	mel->setElementData(elementDataPrototype->clone()); 
       }
 
-      // We will delete at least some element DOFs in the next but will keep the
-      // element object still in memory. Therefore the DOF pointer must be set
-      // invalid.
-      (*macroIt)->getElement()->setDofValid(false);
+      mel->delDofPtr();
     }     
 
 
-    // === Check now all the DOFs that have no owner anymore and therefore have ===
-    // === to be removed.                                                       ===
+    // === Check now all the DOFs that have no owner anymore and therefore  ===
+    // === have to be removed.                                              ===
 
     for (DofElMap::iterator dofsIt = dofsOwner.begin(); 
 	 dofsIt != dofsOwner.end(); ++dofsIt) {
@@ -371,8 +369,8 @@ namespace AMDiS {
 
       for (std::set<MacroElement*>::iterator elIter = dofsIt->second.begin();
 	   elIter != dofsIt->second.end(); ++elIter) {
-	std::set<MacroElement*>::iterator mIt = macros.find(*elIter);
-	if (mIt == macros.end()) {
+	std::set<MacroElement*>::iterator mIt = delMacros.find(*elIter);
+	if (mIt == delMacros.end()) {
 	  deleteDof = false;
 	  break;
 	}
@@ -525,7 +523,8 @@ namespace AMDiS {
       int n  = localAdmin->getNumberOfDofs(position);
       int n0 = localAdmin->getNumberOfPreDofs(position);
       
-      TEST_EXIT_DBG(n + n0 <= ndof)("n=%d, n0=%d too large: ndof=%d\n", n, n0, ndof);
+      TEST_EXIT_DBG(n + n0 <= ndof)
+	("n = %d, n0 = %d too large: ndof = %d\n", n, n0, ndof);
       
       for (int j = 0; j < n; j++)
 	dof[n0 + j] = const_cast<DOFAdmin*>(localAdmin)->getDOFIndex();
@@ -630,7 +629,6 @@ namespace AMDiS {
     }
 
     delete [] dof;
-    dof = NULL;
   }
 
 

AMDiS/src/Mesh.h

@@ -798,10 +798,9 @@ namespace AMDiS {
      */
     DimVec<int> nDof;
 
-    /** \brief
-     * Number of nodes on a single element where DOFs are located; needed for 
-     * the (de-) allocation of the dof-vector on the element (\ref Element::dof);
-     */
+    /// Number of nodes on a single element where DOFs are located. Needed for 
+    /// the (de-) allocation of the dof-vector on the element (\ref Element::dof).
+    /// Here "node" is equivalent to the number of basis functions on the element.
     int nNodeEl;
 
     /** \brief

AMDiS/src/io/MacroReader.cc

@@ -301,9 +301,7 @@ namespace AMDiS {
 	  setDof(mesh->getNode(CENTER), mesh->getDof(CENTER));
     }
 
-    /****************************************************************************/
-    /* domain size                                                              */
-    /****************************************************************************/
+    // === Domain size ===
 
     WorldVector<double> x_min, x_max;
 

AMDiS/src/parallel/MeshDistributor.cc

@@ -271,6 +271,7 @@ namespace AMDiS {
 
 	  macroElementNeighbours[(*it)->getIndex()][i] = -1;
 	  macroElementNeighbours[neighIndex][(*it)->getOppVertex(i)] = -1;
+
 	}
       }
     }
@@ -1138,11 +1139,13 @@ namespace AMDiS {
     // === Create new element weights. ===
 
     elemWeights.clear();
-    TraverseStack stack;
-    ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
-    while (elInfo) {
-      elemWeights[elInfo->getMacroElement()->getIndex()]++;      
-      elInfo = stack.traverseNext(elInfo);
+    {
+      TraverseStack stack;
+      ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
+      while (elInfo) {
+	elemWeights[elInfo->getMacroElement()->getIndex()]++;      
+	elInfo = stack.traverseNext(elInfo);
+      }
     }
 
     // === Run mesh partitioner to calculate a new mesh partitioning.  ===
@@ -1207,11 +1210,13 @@ namespace AMDiS {
 	mel->setNeighbour(i, NULL);
 
       // Create new DOFs for the macro element.
+      mel->getElement()->createNewDofPtrs();
+
       for (int i = 0; i < mesh->getGeo(VERTEX); i++)
 	mel->getElement()->setDof(i, mesh->getDof(VERTEX));
 
-      // The macro element now has valid DOF pointers.
-      mel->getElement()->setDofValid(true);
+      for (int i = 0; i < mesh->getGeo(EDGE); i++)
+	mel->getElement()->setDof(mesh->getGeo(VERTEX) + i, mesh->getDof(EDGE));
 
       // Push the macro element to all macro elements in mesh.
       mesh->getMacroElements().push_back(mel);
@@ -1852,7 +1857,7 @@ namespace AMDiS {
   void MeshDistributor::removeMacroElements()
   {
     FUNCNAME("MeshDistributor::removeMacroElements()");
-
+    
     std::set<MacroElement*> macrosToRemove;
     for (deque<MacroElement*>::iterator it = mesh->firstMacroElement();
 	 it != mesh->endOfMacroElements(); ++it)

AMDiS/src/parallel/MeshManipulation.cc

@@ -228,9 +228,9 @@ namespace AMDiS {
 
     elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_EVERY_EL_PREORDER);
     while (elInfo) {
-      for (int i = 0; i < mesh->getGeo(VERTEX); i++)
+      for (int i = 0; i < mesh->getNumberOfNodes(); i++)
 	if (mapDelDofs.count(elInfo->getElement()->getDof(i)) == 1)
-	  elInfo->getElement()->setDof(i, const_cast<DegreeOfFreedom*>(mapDelDofs[elInfo->getElement()->getDof(i)]));	
+	  elInfo->getElement()->setDof(i, const_cast<DegreeOfFreedom*>(mapDelDofs[elInfo->getElement()->getDof(i)]));
 
       elInfo = stack.traverseNext(elInfo);
     }