Work on parallel jump residuum. Is working, but source code must be improved.

AMDiS/src/AdaptInfo.h

@@ -389,24 +389,12 @@ namespace AMDiS {
     /// Sets \ref est_sum.
     inline void setEstSum(double e, int index) 
     {
-#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
-    double send_est_sum = e;
-    double est_sum = 0.0;
-    MPI::COMM_WORLD.Allreduce(&send_est_sum, &est_sum, 1, MPI_DOUBLE, MPI_SUM);
-    e = est_sum;
-#endif
       scalContents[index]->est_sum = e;
     }
 
     /// Sets \ref est_max.
     inline void setEstMax(double e, int index) 
     {
-#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
-    double send_est_max = e;
-    double est_max = 0.0;
-    MPI::COMM_WORLD.Allreduce(&send_est_max, &est_max, 1, MPI_DOUBLE, MPI_MAX);
-    e = est_max;
-#endif
       scalContents[index]->est_max = e;
     }
 

AMDiS/src/BoundaryObject.cc

@@ -128,6 +128,19 @@ namespace AMDiS {
   }
 
 
+  bool BoundaryObject::operator<(const BoundaryObject& other) const
+  {
+    if (elIndex == other.elIndex) {
+      if (subObj == other.subObj)
+	return ithObj < other.ithObj;
+
+      return subObj < other.subObj;
+    }
+     
+    return (elIndex < other.elIndex);
+  }
+
+
   bool AtomicBoundary::operator==(const AtomicBoundary& other) const
   {
     return (rankObj == other.rankObj && 

AMDiS/src/BoundaryObject.h

@@ -55,6 +55,8 @@ namespace AMDiS {
     
     bool operator!=(const BoundaryObject& other) const;
 
+    bool operator<(const BoundaryObject& other) const;
+
     /// The macro element to which the boundary element corresponds to.
     Element* el;
 
@@ -64,11 +66,9 @@ namespace AMDiS {
     /// Element type index, only used in 3d.
     int elType;
 
-    /** \brief
-     * Defines the geometrical object at the boundary. It must be "a part" of the
-     * macro element \ref el, i.e., either 1 (a vertex), 2 (an edge) or 3 
-     * (a face).
-     */
+    /// Defines the geometrical object at the boundary. It must be "a part" of
+    /// the macro element \ref el, i.e., either 1 (a vertex), 2 (an edge) or 3 
+    /// (a face).
     GeoIndex subObj;
 
     /** \brief

AMDiS/src/ElInfo.h

@@ -472,41 +472,29 @@ namespace AMDiS {
     /// \ref element is child of element parent
     Element *parent;
 
-    /** \brief 
-     * \ref element is an element of the binary tree located at MacroElement 
-     * macroElement
-     */
+    /// \ref element is an element of the binary tree located at MacroElement 
+    /// macroElement
     MacroElement *macroElement;
 
-    /** \brief
-     * Indicates wich elements will be called and wich information should be
-     * present while mesh traversal.
-     */
+    /// Indicates wich elements will be called and wich information should be
+    /// present while mesh traversal.
     Flag fillFlag;
 
-    /** \brief 
-     * Level of the element. The level is zero for macro elements and the level
-     * of the children is (level of the parent + 1). level_ is filled always by
-     * the traversal routines.
-     */
+    /// Level of the element. The level is zero for macro elements and the level
+    /// of the children is (level of the parent + 1). level_ is filled always by
+    /// the traversal routines.
     unsigned char level;
 
-    /** \brief
-     * Elements type index. This is used only for 3d, where type can be either 0, 1 or
-     * 2. In all other cases type is not used and the variable is set to 0.
-     * In 3d, it is filled automatically by the traversal routines.
-     */
+    /// Elements type index. This is used only for 3d, where type can be either 0, 1 or
+    /// 2. In all other cases type is not used and the variable is set to 0.
+    /// In 3d, it is filled automatically by the traversal routines.
     int elType;
 
-    /** \brief
-     * This ElInfo is the iChild-th child of the parent element.
-     */
+    /// This ElInfo is the iChild-th child of the parent element.
     int iChild;
 
-    /** \brief 
-     * \ref coord[i] is a WorldVector<double> storing the world coordinates of the
-     * i-th vertex of element \ref element.
-     */
+    /// \ref coord[i] is a WorldVector<double> storing the world coordinates of the
+    /// i-th vertex of element \ref element.
     FixVec<WorldVector<double>, VERTEX> coord;
 
     /** \brief 
@@ -518,34 +506,24 @@ namespace AMDiS {
      */
     FixVec<BoundaryType, BOUNDARY> boundary;
 
-    /** \brief
-     * Vector storing pointers to projections for each face, edge, vertex.
-     */
+    /// Vector storing pointers to projections for each face, edge, vertex.
     FixVec<Projection*, PROJECTION> projection;
 
-    /** \brief 
-     * oppCoord[i] coordinates of the i-th neighbour vertex opposite the
-     * common edge/face.
-     */
+    /// oppCoord[i] coordinates of the i-th neighbour vertex opposite the
+    /// common edge/face.
     FixVec<WorldVector<double>, NEIGH> oppCoord;
 
-    /** \brief 
-     * neighbour[i] pointer to the element at the edge/face with local index i.
-     * It is a pointer to NULL for boundary edges/faces.
-     */
+    /// neighbour[i] pointer to the element at the edge/face with local index i.
+    /// It is a pointer to NULL for boundary edges/faces.
     FixVec<Element*, NEIGH> neighbour;
 
-    /** \brief
-     * neighbourCoord[i][j] are the coordinate of the j-th vertex of the i-th
-     * neighbour element with the common edge/face.
-     */
+    /// neighbourCoord[i][j] are the coordinate of the j-th vertex of the i-th
+    /// neighbour element with the common edge/face.
     FixVec<FixVec<WorldVector<double>, VERTEX>, NEIGH> neighbourCoord;
 
-    /** \brief 
-     * oppVertex[i] is undefined if neighbour[i] is a pointer to NULL. 
-     * Otherwise it is the local index of the neighbour's vertex opposite the
-     * common edge/face.
-     */
+    /// oppVertex[i] is undefined if neighbour[i] is a pointer to NULL. 
+    /// Otherwise it is the local index of the neighbour's vertex opposite the
+    /// common edge/face.
     FixVec<int, NEIGH> oppVertex;
 
     /// Elements determinant.
@@ -569,18 +547,14 @@ namespace AMDiS {
     
     static std::vector<std::map<std::pair<int, unsigned long>, mtl::dense2D<double> > > subElemGradMatrices;
 
-    /** \brief 
-     * child_vertex[el_type][child][i] = father's local vertex index of new 
-     * vertex i. 4 stands for the newly generated vertex .
-     */
+    /// child_vertex[el_type][child][i] = father's local vertex index of new 
+    /// vertex i. 4 stands for the newly generated vertex .     
     static const int childVertex[3][2][4];
 
-    /** \brief 
-     * child_edge[el_type][child][i] = father's local edge index of new edge i.
-     * new edge 2 is half of old edge 0, new edges 4,5 are really new edges, and
-     * value is different: child_edge[][][4,5] = index of same edge in other
-     * child.
-     */
+    /// child_edge[el_type][child][i] = father's local edge index of new edge i.
+    /// new edge 2 is half of old edge 0, new edges 4,5 are really new edges, and
+    /// value is different: child_edge[][][4,5] = index of same edge in other
+    /// child.
     static const int childEdge[3][2][6];
 
     friend class ElInfo1d;

AMDiS/src/Marker.cc

@@ -140,7 +140,7 @@ namespace AMDiS {
     
     markRLimit = ESThetaP * epsP / nLeaves;
     markCLimit = ESThetaCP * epsP / nLeaves;
-    
+
     INFO(info, 2)("start mark_limits: %.3le %.3le   nt = %d\n", 
 		  markRLimit, markCLimit, nLeaves);
   }

AMDiS/src/ProblemStat.cc

@@ -1726,11 +1726,12 @@ namespace AMDiS {
 
     map<int, double> vec;    
     TraverseStack stack;
-    ElInfo *elInfo = stack.traverseFirst(this->getMesh(comp),  -1, 
-					 Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS);
+    ElInfo *elInfo = 
+      stack.traverseFirst(this->getMesh(comp),  -1, Mesh::CALL_LEAF_EL);
     
     while (elInfo) {		  
-      vec[elInfo->getElement()->getIndex()] = elInfo->getElement()->getEstimation(comp);
+      vec[elInfo->getElement()->getIndex()] = 
+	elInfo->getElement()->getEstimation(comp);
       elInfo = stack.traverseNext(elInfo);
     }
 

AMDiS/src/Quadrature.h

@@ -157,19 +157,15 @@ namespace AMDiS {
   
 
 
-    /** \brief
-     * Returns \ref lambda[a][b] which is the b-th coordinate entry of the a-th
-     * quadrature point
-     */
+    /// Returns \ref lambda[a][b] which is the b-th coordinate entry of the a-th
+    /// quadrature point
     inline double getLambda(int a, int b) const 
     {
       return (lambda ? (*lambda)[a][b] : 0.0);
     }
 
-    /** \brief
-     * Returns \ref lambda[a] which is a DimVec<double> containing the 
-     * coordiantes of the a-th quadrature point
-     */
+    /// Returns \ref lambda[a] which is a DimVec<double> containing the 
+    /// coordiantes of the a-th quadrature point
     inline const DimVec<double>& getLambda(int a) const 
     {
       return (*lambda)[a]; 
@@ -206,10 +202,8 @@ namespace AMDiS {
     double *w;
 
   protected:
-    /** \brief
-     * Initialisation of all static Quadrature objects which will be returned
-     * by \ref provideQuadrature()
-     */
+    /// Initialisation of all static Quadrature objects which will be returned
+    /// by \ref provideQuadrature()
     static void initStaticQuadratures();
 
     /** \name static quadratures, used weights, and barycentric coords
@@ -288,21 +282,15 @@ namespace AMDiS {
 
 
 
-  /** \brief
-   * Pre-compute the values of all basis functions at all quadrature nodes;  
-   */
+  /// Pre-compute the values of all basis functions at all quadrature nodes;  
   const Flag INIT_PHI=1;
 
-  /** \brief
-   * Pre-compute the gradients (with respect to the barycentric coordinates) of
-   * all basis functions at all quadrature nodes
-   */
+  /// Pre-compute the gradients (with respect to the barycentric coordinates) of
+  /// all basis functions at all quadrature nodes
   const Flag INIT_GRD_PHI=2; 
 
-  /** \brief
-   * pre-compute all 2nd derivatives (with respect to the barycentric 
-   * coordinates) of all basis functions at all quadrature nodes;
-   */
+  /// pre-compute all 2nd derivatives (with respect to the barycentric 
+  /// coordinates) of all basis functions at all quadrature nodes;
   const Flag INIT_D2_PHI=4;
 
 
@@ -325,10 +313,8 @@ namespace AMDiS {
   class FastQuadrature
   {
   protected:
-    /** \brief
-     * Constructs a FastQuadrature for the given Quadrature, BasisFunction, and
-     * flag.
-     */
+    /// Constructs a FastQuadrature for the given Quadrature, BasisFunction, and
+    /// flag.
     FastQuadrature(BasisFunction* basFcts, Quadrature* quad, Flag flag)
       : init_flag(flag),
 	phi(0, 0),
@@ -481,10 +467,8 @@ namespace AMDiS {
     }
 
   protected:
-    /** \brief
-     * Specifies which information should be pre-computed. Can be \ref INIT_PHI, 
-     * \ref INIT_GRD_PHI, or \ref INIT_D2_PHI
-     */
+    /// Specifies which information should be pre-computed. Can be \ref INIT_PHI, 
+    /// \ref INIT_GRD_PHI, or \ref INIT_D2_PHI
     Flag init_flag;
 
     /** \brief
@@ -513,16 +497,12 @@ namespace AMDiS {
      */
     MatrixOfFixVecs<DimMat<double> > *D2Phi;
 
-    /** \brief
-     * List of all used FastQuadratures
-     */
+    /// List of all used FastQuadratures
     static list<FastQuadrature*> fastQuadList;
 
-    /** \brief
-     * Maximal number of quadrature points for all yet initialised FastQuadrature
-     * objects. This value may change after a new initialisation of a
-     * FastQuadrature
-     */
+    /// Maximal number of quadrature points for all yet initialised FastQuadrature
+    /// objects. This value may change after a new initialisation of a
+    /// FastQuadrature
     static int max_points;
 
     /// This FastQuadrature stores values for Quadrature quadrature

AMDiS/src/est/ResidualEstimator.cc

@@ -173,7 +173,6 @@ namespace AMDiS {
 
     InteriorBoundary &intBoundary = 
       MeshDistributor::globalMeshDistributor->getIntBoundary();
-    RankToBoundMap& otherBound = intBoundary.getOther();
 
     ElInfo *elInfo = mesh->createNewElInfo();
     elInfo->setFillFlag(Mesh::FILL_COORDS);
@@ -181,56 +180,66 @@ namespace AMDiS {
     StdMpi<vector<double> > stdMpiDet(MeshDistributor::globalMeshDistributor->getMpiComm());
     StdMpi<vector<vector<WorldVector<double> > > > stdMpiGrdUh(MeshDistributor::globalMeshDistributor->getMpiComm());
 
-    for (RankToBoundMap::iterator it = otherBound.begin();
-	 it != otherBound.end(); ++it) {
+    RankToBoundMap allBounds = intBoundary.getOther();
+    allBounds.insert(intBoundary.getOwn().begin(), intBoundary.getOwn().end());
+
+    for (RankToBoundMap::iterator it = allBounds.begin();
+	 it != allBounds.end(); ++it) {
+
+      vector<BoundaryObject> subBound;
+
       for (unsigned int i = 0; i < it->second.size(); i++) {
 	BoundaryObject &bObj = it->second[i].rankObj;
 	if (bObj.subObj == VERTEX)
 	  continue;
 
-	vector<BoundaryObject> subBound;
 	bObj.el->getSubBoundary(bObj, subBound);
+      }
 
-	WorldVector<int> faceIndices;
-
-	for (unsigned int j = 0; j < subBound.size(); j++) {
-	  Element *el = subBound[j].el;	  
-	  int oppV = subBound[j].ithObj;
-
-	  elInfo->setElement(el);
-	  el->sortFaceIndices(oppV, faceIndices);
-
-	  elInfo->getCoord(oppV) = coords[el->getDof(oppV, 0)];
-	  for (int k = 0; k < dim; k++)
-	    elInfo->getCoord(faceIndices[k]) = coords[el->getDof(k, 0)];
-
-	  double detNeigh = abs(elInfo->calcGrdLambda(*lambdaNeigh));
-	  stdMpiDet.getSendData(it->first).push_back(detNeigh);
-
-	  for (int system = 0; system < nSystems; system++) {
-	    if (matrix[system] == NULL || secondOrderTerms[system] == false) 
-	      continue;
-
-	    uh[system]->getLocalVector(el, uhNeigh[system]);
-
-	    for (int iq = 0; iq < nPointsSurface; iq++) {
+      if (subBound.size() == 0)
+	continue;
 
-	      (*lambda)[oppV] = 0.0;
-	      for (int k = 0; k < dim; k++)
-		(*lambda)[faceIndices[k]] = surfaceQuad->getLambda(iq, k);
+      WorldVector<int> faceIndices;
 
-	      basFcts[system]->evalGrdUh(*lambda, *lambdaNeigh, uhNeigh[system], grdUhNeigh[iq]);
-	    }
+      for (unsigned int i = 0; i < subBound.size(); i++) {
+	Element *el = subBound[i].el;	  
+	int oppV = subBound[i].ithObj;
+	
+	elInfo->setElement(el);
+	el->sortFaceIndices(oppV, faceIndices);	
+	for (int k = 0; k <= dim; k++)
+	  elInfo->getCoord(k) = coords[el->getDof(k, 0)];
+	
+	double detNeigh = abs(elInfo->calcGrdLambda(*lambdaNeigh));
+	stdMpiDet.getSendData(it->first).push_back(detNeigh);
 
-	    stdMpiGrdUh.getSendData(it->first).push_back(grdUhNeigh);
+	for (int system = 0; system < nSystems; system++) {
+	  if (matrix[system] == NULL || secondOrderTerms[system] == false)
+	    continue;
+	  
+	  uh[system]->getLocalVector(el, uhNeigh[system]);
+	  
+	  for (int iq = 0; iq < nPointsSurface; iq++) {
+	    
+	    (*lambda)[oppV] = 0.0;
+	    for (int k = 0; k < dim; k++)
+	      (*lambda)[faceIndices[k]] = surfaceQuad->getLambda(iq, k);
+	    
+	    basFcts[system]->evalGrdUh(*lambda, *lambdaNeigh, uhNeigh[system], grdUhNeigh[iq]);
 	  }
+	  
+	  stdMpiGrdUh.getSendData(it->first).push_back(grdUhNeigh);	  
 	}
       }
+
+      stdMpiDet.recv(it->first);
+      stdMpiGrdUh.recv(it->first);      
     }
 
     stdMpiDet.updateSendDataSize();
     stdMpiGrdUh.updateSendDataSize();
 
+#if 0
     RankToBoundMap& ownBound = intBoundary.getOwn();
 
     for (RankToBoundMap::iterator it = ownBound.begin();
@@ -245,13 +254,14 @@ namespace AMDiS {
 	}
       }
     }
+#endif
 
     stdMpiDet.startCommunication();
     stdMpiGrdUh.startCommunication();
 
 
-    for (RankToBoundMap::iterator it = ownBound.begin();
-	 it != ownBound.end(); ++it) {
+    for (RankToBoundMap::iterator it = allBounds.begin();
+	 it != allBounds.end(); ++it) {
       vector<BoundaryObject> subBound;
 	
       for (unsigned int i = 0; i < it->second.size(); i++) {
@@ -270,6 +280,11 @@ namespace AMDiS {
       
       TEST_EXIT_DBG(subBound.size() == stdMpiGrdUh.getRecvData(it->first).size())
 	("Should not happen!\n");      
+
+      for (unsigned int i = 0; i < subBound.size(); i++) {
+	elBoundDet[subBound[i]] = stdMpiDet.getRecvData(it->first)[i];
+	elBoundGrdUhNeigh[subBound[i]] = stdMpiGrdUh.getRecvData(it->first)[i];
+      }
     }    
 
     MSG("INIT PARALLEL FINISCHED!\n");
@@ -482,20 +497,41 @@ namespace AMDiS {
     const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
     double det = elInfo->getDet();
     double h2 = h2_from_det(det, dim);
+    BoundaryObject blub;
 
     for (int face = 0; face < nNeighbours; face++) {  
       Element *neigh = const_cast<Element*>(elInfo->getNeighbour(face));
 
+      bool parallelMode = false;
+
+#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
+      if (neigh == NULL) {
+	BoundaryObject testObj(el, elInfo->getType(), EDGE, face);
+	
+	if (elBoundDet.count(testObj)) {
+	  parallelMode = true;
+	  blub = testObj;
+	} else {
+	  continue;
+	}
+      }
+
+      if (!parallelMode && !(neigh && neigh->getMark()))
+	continue;
+#else
       if (!(neigh && neigh->getMark()))
 	continue;
+#endif
 
-      int oppV = elInfo->getOppVertex(face);
-	
+      int oppV = elInfo->getOppVertex(face);	
       el->sortFaceIndices(face, faceIndEl);
-      neigh->sortFaceIndices(oppV, faceIndNeigh);	
-      neighInfo->setElement(const_cast<Element*>(neigh));
-      neighInfo->setFillFlag(Mesh::FILL_COORDS);     
-      neighInfo->getCoord(oppV) = elInfo->getOppCoord(face);
+
+      if (!parallelMode) {
+	neigh->sortFaceIndices(oppV, faceIndNeigh);	
+	neighInfo->setElement(const_cast<Element*>(neigh));
+	neighInfo->setFillFlag(Mesh::FILL_COORDS);     
+	neighInfo->getCoord(oppV) = elInfo->getOppCoord(face);
+      }
       
       // periodic leaf data ?
       ElementData *ldp = el->getElementData()->getElementData(PERIODIC);	
@@ -526,7 +562,7 @@ namespace AMDiS {
 	}
       }  // if (ldp)
       
-      if (!periodicCoords) {
+      if (!periodicCoords && !parallelMode) {
 	for (int i = 0; i < dim; i++) {
 	  int i1 = faceIndEl[i];
 	  int i2 = faceIndNeigh[i];
@@ -534,11 +570,19 @@ namespace AMDiS {
 	}
       }
 	
+      double detNeigh = 0.0;
       Parametric *parametric = mesh->getParametric();
-      if (parametric)
-	neighInfo = parametric->addParametricInfo(neighInfo);	  
+      if (!parallelMode) {
+	if (parametric)
+	  neighInfo = parametric->addParametricInfo(neighInfo);	  
 
-      double detNeigh = abs(neighInfo->calcGrdLambda(*lambdaNeigh));
+	detNeigh = abs(neighInfo->calcGrdLambda(*lambdaNeigh));
+      } else {
+	TEST_EXIT_DBG(!parametric)("No yet implemented!\n");
+
+	detNeigh = elBoundDet[blub];
+      }
+      
            
       for (int iq = 0; iq < nPointsSurface; iq++)
 	jump[iq].set(0.0);     
@@ -548,7 +592,8 @@ namespace AMDiS {
 	  continue;
 	      
 	uh[system]->getLocalVector(el, uhEl[system]);	
-	uh[system]->getLocalVector(neigh, uhNeigh[system]);
+	if (!parallelMode)
+	  uh[system]->getLocalVector(neigh, uhNeigh[system]);
 	  
 	for (int iq = 0; iq < nPointsSurface; iq++) {
 	  (*lambda)[face] = 0.0;
@@ -556,12 +601,18 @@ namespace AMDiS {
 	    (*lambda)[faceIndEl[i]] = surfaceQuad->getLambda(iq, i);
 	  
 	  basFcts[system]->evalGrdUh(*lambda, grdLambda, uhEl[system], grdUhEl[iq]);
+
+	  if (!parallelMode) {
+	    (*lambda)[oppV] = 0.0;
+	    for (int i = 0; i < dim; i++)
+	      (*lambda)[faceIndNeigh[i]] = surfaceQuad->getLambda(iq, i);
+
+	    basFcts[system]->evalGrdUh(*lambda, *lambdaNeigh, uhNeigh[system], grdUhNeigh[iq]);
+	  } else {
+	    grdUhNeigh[iq] = elBoundGrdUhNeigh[blub][iq];
+	  }
 	  
-	  (*lambda)[oppV] = 0.0;
-	  for (int i = 0; i < dim; i++)
-	    (*lambda)[faceIndNeigh[i]] = surfaceQuad->getLambda(iq, i);
-	  
-	  basFcts[system]->evalGrdUh(*lambda, *lambdaNeigh, uhNeigh[system], grdUhNeigh[iq]);
+
 	  
 	  grdUhEl[iq] -= grdUhNeigh[iq];
 	} // for iq				
@@ -599,11 +650,13 @@ namespace AMDiS {
 	val *= C1 * h2_from_det(d, dim) * d;
       else
 	val *= C1 * d;
-      
-      if (parametric)
-	neighInfo = parametric->removeParametricInfo(neighInfo);
-      
-      neigh->setEstimation(neigh->getEstimation(row) + val, row);
+
+      if (!parallelMode) {
+	if (parametric)
+	  neighInfo = parametric->removeParametricInfo(neighInfo);      
+	neigh->setEstimation(neigh->getEstimation(row) + val, row);
+      }
+
       result += val;
     } // for face
     

AMDiS/src/est/ResidualEstimator.h

@@ -188,6 +188,12 @@ namespace AMDiS {
     /// are used to ommit computations of the jump residual that is defined
     /// only on second order terms.
     std::vector<bool> secondOrderTerms;
+
+#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
+    map<BoundaryObject, double> elBoundDet;
+
+    map<BoundaryObject, std::vector<WorldVector<double> > > elBoundGrdUhNeigh;
+#endif
   };
 }