Merge with local changes.

AMDiS/src/Element.h

@@ -457,10 +457,8 @@ namespace AMDiS {
     /// assignment operator
     Element& operator=(const Element& el);
 
-    /** \brief
-     * Checks whether the face with vertices dof[0],..,dof[DIM-1] is
-     * part of mel's boundary. returns the opposite vertex if true, -1 else
-     */
+    /// Checks whether the face with vertices dof[0],..,dof[DIM-1] is
+    /// part of mel's boundary. returns the opposite vertex if true, -1 else
     int oppVertex(FixVec<DegreeOfFreedom*, DIMEN> pdof) const;
 
     /// Refines Element's leaf data
@@ -565,38 +563,30 @@ namespace AMDiS {
     /// to NULL for leaf elements.
     Element *child[2];
 
-    /** \brief
-     * Vector of pointers to DOFs. These pointers must be available for elements
-     * vertices (for the geometric description of the mesh). There my be pointers
-     * for the edges, for faces and for the center of an element. They are 
-     * ordered the following way: The first N_VERTICES entries correspond to the
-     * DOFs at the vertices of the element. The next ones are those at the edges,
-     * if present, then those at the faces, if present, and then those at the 
-     * barycenter, if present.
-     */
+    /// Vector of pointers to DOFs. These pointers must be available for elements
+    /// vertices (for the geometric description of the mesh). There my be pointers
+    /// for the edges, for faces and for the center of an element. They are 
+    /// ordered the following way: The first N_VERTICES entries correspond to the
+    /// DOFs at the vertices of the element. The next ones are those at the edges,
+    /// if present, then those at the faces, if present, and then those at the 
+    /// barycenter, if present.
     DegreeOfFreedom **dof;
 
-    /** \brief
-     * Unique global index of the element. these indices are not strictly ordered
-     * and may be larger than the number of elements in the binary tree (the list
-     * of indices may have holes after coarsening).
-     */
+    /// Unique global index of the element. these indices are not strictly ordered
+    /// and may be larger than the number of elements in the binary tree (the list
+    /// of indices may have holes after coarsening).
     int index;
 
-    /** \brief
-     * Marker for refinement and coarsening. if mark is positive for a leaf
-     * element, this element is refined mark times. if mark is negative for
-     * a leaf element, this element is coarsened -mark times.
-     */
+    /// Marker for refinement and coarsening. if mark is positive for a leaf
+    /// element, this element is refined mark times. if mark is negative for
+    /// a leaf element, this element is coarsened -mark times.
     int mark;
  
-    /** \brief
-     * If the element has a boundary edge on a curved boundary, this is a pointer
-     * to the coordinates of the new vertex that is created due to the refinement
-     * of the element, otherwise it is a NULL pointer. Thus coordinate 
-     * information can be also produced by the traversal routines in the case of 
-     * curved boundary.
-     */
+    /// If the element has a boundary edge on a curved boundary, this is a pointer
+    /// to the coordinates of the new vertex that is created due to the refinement
+    /// of the element, otherwise it is a NULL pointer. Thus coordinate 
+    /// information can be also produced by the traversal routines in the case of 
+    /// curved boundary.
     WorldVector<double> *newCoord;
 
     /// Pointer to the Mesh this element belongs to

AMDiS/src/parallel/PetscSolverFeti.cc

@@ -699,6 +699,22 @@ namespace AMDiS {
     }
   }
 
+  bool PetscSolverFeti::testWirebasketEdge(BoundaryObject &edge, const FiniteElemSpace *feSpace)
+  {
+    Element *el = edge.el;
+    int i0 = el->getVertexOfEdge(edge.ithObj, 0);
+    int i1 = el->getVertexOfEdge(edge.ithObj, 1);
+    DegreeOfFreedom d0 = el->getDof(i0, 0);
+    DegreeOfFreedom d1 = el->getDof(i1, 0);
+    WorldVector<double> c0, c1;
+    el->getMesh()->getDofIndexCoords(d0, feSpace, c0);
+    el->getMesh()->getDofIndexCoords(d1, feSpace, c1);
+    bool xe = fabs(c0[0] - c1[0]) < 1e-8;
+    bool ye = fabs(c0[1] - c1[1]) < 1e-8;
+    bool ze = fabs(c0[2] - c1[2]) < 1e-8;
+    int counter = static_cast<int>(xe) + static_cast<int>(ye) + static_cast<int>(ze);
+    return (counter == 2);
+  }
 
   void PetscSolverFeti::createMatAugmentedLagrange(vector<const FiniteElemSpace*> &feSpaces)
   {
@@ -709,13 +725,16 @@ namespace AMDiS {
 
     double wtime = MPI::Wtime();
 
-    nRankEdges = 0;
     nOverallEdges = 0;
     InteriorBoundary &intBound = meshDistributor->getIntBoundary();
+    std::set<BoundaryObject> allEdges;
     for (InteriorBoundary::iterator it(intBound.getOwn()); !it.end(); ++it)
-      if (it->rankObj.subObj == EDGE)
-	nRankEdges++;
-    
+      if (it->rankObj.subObj == EDGE && 
+	  testWirebasketEdge(it->rankObj, feSpaces[0]) && 
+	  allEdges.count(it->rankObj) == 0)
+	allEdges.insert(it->rankObj);
+
+    nRankEdges = allEdges.size();    
     int rStartEdges = 0;
     mpi::getDofNumbering(mpiCommGlobal, nRankEdges, rStartEdges, nOverallEdges);
 
@@ -733,30 +752,30 @@ namespace AMDiS {
     MatSetOption(mat_augmented_lagrange, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
 
     int rowCounter = rStartEdges;
-    for (InteriorBoundary::iterator it(intBound.getOwn()); !it.end(); ++it) {
-      if (it->rankObj.subObj == EDGE) {
-	for (int i = 0; i < feSpaces.size(); i++) {
-	  DofContainer edgeDofs;
-	  it->rankObj.el->getAllDofs(feSpaces[i], it->rankObj, edgeDofs);
-
-	  MSG("SIZE = %d\n", edgeDofs.size());
-
-	  for (DofContainer::iterator it = edgeDofs.begin();
-	       it != edgeDofs.end(); it++) {
-	    TEST_EXIT_DBG(isPrimal(feSpaces[i], **it) == false)
-	      ("Should not be primal!\n");
-
-	    if (dirichletRows[feSpaces[i]].count(**it))
-	      continue;
+    for (std::set<BoundaryObject>::iterator edgeIt = allEdges.begin(); 
+	 edgeIt != allEdges.end(); ++edgeIt) {
 
-	    int col = lagrangeMap.getMatIndex(i, **it);
-	    double value = 1.0;
-	    MatSetValue(mat_augmented_lagrange, rowCounter, col, value, INSERT_VALUES);
-	  }
-
-	  rowCounter++;
+      for (int i = 0; i < feSpaces.size(); i++) {
+	DofContainer edgeDofs;
+	edgeIt->el->getAllDofs(feSpaces[i], *edgeIt, edgeDofs);
+	
+	MSG("SIZE = %d\n", edgeDofs.size());
+	
+	for (DofContainer::iterator it = edgeDofs.begin();
+	     it != edgeDofs.end(); ++it) {
+	  TEST_EXIT_DBG(isPrimal(feSpaces[i], **it) == false)
+	    ("Should not be primal!\n");
+	  
+	  if (dirichletRows[feSpaces[i]].count(**it))
+	    continue;
+	  
+	  int col = lagrangeMap.getMatIndex(i, **it);
+	  double value = 1.0;
+	  MatSetValue(mat_augmented_lagrange, rowCounter, col, value, INSERT_VALUES);
 	}
-      }
+	
+	rowCounter++;
+      }      
     }
 
     MatAssemblyBegin(mat_augmented_lagrange, MAT_FINAL_ASSEMBLY);
@@ -1013,7 +1032,38 @@ namespace AMDiS {
       MatShellSetOperation(tmp, MATOP_MULT, 
 			   (void(*)(void))petscMultMatSchurPrimalAugmented);
 
+      schurPrimalAugmentedData.testMode = false;
       MatComputeExplicitOperator(tmp, &mat_schur_primal);
+      schurPrimalAugmentedData.testMode = true;
+
+      {
+	Vec tvec0, tvec1;
+	VecCreateMPI(PETSC_COMM_WORLD, 
+		     primalDofMap.getRankDofs() + nRankEdges,
+		     primalDofMap.getOverallDofs() + nOverallEdges, 
+		     &tvec0);
+	VecDuplicate(tvec0, &tvec1);
+
+	if (meshDistributor->getMpiRank() == 0)
+	  VecSetValue(tvec0, 12, 0.0, INSERT_VALUES);
+
+	VecAssemblyBegin(tvec0);
+	VecAssemblyEnd(tvec1);
+
+	MatMult(tmp, tvec0, tvec1);
+	
+	//	VecView(tvec1, PETSC_VIEWER_STDOUT_WORLD);
+
+	PetscReal n, a, b;
+	VecNorm(tvec1, NORM_2, &n);
+	VecMax(tvec1, PETSC_NULL, &a);
+	VecMin(tvec1, PETSC_NULL, &b);
+
+	MSG("DIE WERTE SIND %e %e %e\n", n, a, b);
+      }
+
+
+      //      MatView(mat_schur_primal, PETSC_VIEWER_STDOUT_WORLD);
 
       MatDestroy(&tmp);
       schurPrimalAugmentedData.subSolver = NULL;

AMDiS/src/parallel/PetscSolverFeti.h

@@ -133,6 +133,8 @@ namespace AMDiS {
 
     void createMatAugmentedLagrange(vector<const FiniteElemSpace*> &feSpaces);
 
+    bool testWirebasketEdge(BoundaryObject &edge, const FiniteElemSpace *feSpace);
+
     ///
     void createPreconditionerMatrix(Matrix<DOFMatrix*> *mat);
 

AMDiS/src/parallel/PetscSolverFetiOperators.cc

@@ -66,6 +66,11 @@ namespace AMDiS {
       PetscInt muLocalSize = allLocalSize - primalLocalSize;
       PetscInt muSize = allSize - primalSize;
 
+      if (data->testMode) {
+	VecView(x, PETSC_VIEWER_STDOUT_WORLD);
+	MSG("LOCAL SIZE: %d %d %d\n", allLocalSize, primalLocalSize, muLocalSize);
+      }
+
       VecCreateMPI(PETSC_COMM_WORLD, muLocalSize, muSize, &x_mu);
       VecCreateMPI(PETSC_COMM_WORLD, muLocalSize, muSize, &y_mu);
 
@@ -76,10 +81,17 @@ namespace AMDiS {
       VecGetArray(x_primal, &primalValue);
       VecGetArray(x_mu, &muValue);
 
-      for (int i = 0; i < primalLocalSize; i++)
+      for (int i = 0; i < primalLocalSize; i++) {
+	if (data->testMode && allValue[i] != 1.0)
+	  MSG("ONE IS in %d ith local primal value!\n", i);
 	primalValue[i] = allValue[i];
-      for (int i = 0; i < muLocalSize; i++)
+      }
+      for (int i = 0; i < muLocalSize; i++) {
+	if (data->testMode && allValue[primalLocalSize + i] != 0.0)
+	  MSG("ONE IS in %d ith local mu value!\n", primalLocalSize + i);
+
 	muValue[i] = allValue[primalLocalSize + i];
+      }
 
       VecRestoreArray(x, &allValue);
       VecRestoreArray(x_primal, &primalValue);

AMDiS/src/parallel/PetscSolverFetiStructs.h

@@ -65,6 +65,8 @@ namespace AMDiS {
     PetscSolver* subSolver;
 
     bool nestedVec;
+
+    bool testMode;
   };