Fixed lumpred precondition for FETI-DP stokes mode.

AMDiS/src/parallel/PetscSolverFeti.cc

@@ -911,6 +911,22 @@ namespace AMDiS {
 	}
 	
 	if (stokesMode) {
+	  DOFVector<double> tmpvec(pressureFeSpace, "tmpvec");
+	  createH2vec(tmpvec);
+	  interfaceDofMap.createVec(fetiInterfaceLumpedPreconData.h2vec);
+
+	  DofMap& m = interfaceDofMap[pressureFeSpace].getMap();
+	  for (DofMap::iterator it = m.begin(); it != m.end(); ++it) {
+	    if (meshDistributor->getDofMap()[pressureFeSpace].isRankDof(it->first)) {
+	      int index = interfaceDofMap.getMatIndex(pressureComponent, it->first);
+	      VecSetValue(fetiInterfaceLumpedPreconData.h2vec, 
+			  index, tmpvec[it->first], INSERT_VALUES);
+	    }
+	  }  
+
+	  VecAssemblyBegin(fetiInterfaceLumpedPreconData.h2vec);
+	  VecAssemblyEnd(fetiInterfaceLumpedPreconData.h2vec);
+	  
 	  localDofMap.createVec(fetiInterfaceLumpedPreconData.tmp_vec_b1);
 	  fetiInterfaceLumpedPreconData.subSolver = subdomain;	
 	}
@@ -1941,6 +1957,49 @@ namespace AMDiS {
   }
 
 
+  void PetscSolverFeti::createH2vec(DOFVector<double> &vec)
+  {
+    FUNCNAME("PetscSolverFeti::createH2vec()");
+    
+    vec.set(0.0);
+    
+    DOFVector<int> cvec(vec.getFeSpace(), "cvec");
+    cvec.set(0.0);
+
+    Mesh *mesh = vec.getFeSpace()->getMesh();    
+    TEST_EXIT(mesh->getDim() == 2)("This works only in 2D!\n");
+
+    TraverseStack stack;
+    ElInfo *elInfo = 
+      stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS);
+    while (elInfo) {
+      Element *el = elInfo->getElement();
+      for (int i = 0; i < 3; i++) {
+	int d0 = el->getVertexOfEdge(i, 0);
+	int d1 = el->getVertexOfEdge(i, 1);
+
+	WorldVector<double> c0 = elInfo->getCoord(d0);
+	WorldVector<double> c1 = elInfo->getCoord(d1);
+	c0 -= c1;
+	double length = norm(c0);
+	vec[el->getDof(d0, 0)] += length;
+	vec[el->getDof(d1, 0)] += length;
+	cvec[el->getDof(d0, 0)] += 1;
+	cvec[el->getDof(d1, 0)] += 1;
+      }
+
+      elInfo = stack.traverseNext(elInfo);
+    }
+
+    DOFIterator<double> dofIt(&vec, USED_DOFS);
+    for (dofIt.reset(); !dofIt.end(); ++dofIt) {
+      DegreeOfFreedom d = dofIt.getDOFIndex();
+      vec[d] /= static_cast<double>(cvec[d]);
+      vec[d] = 1.0 / pow(vec[d], 2.0);
+    }
+  }
+
+
   void PetscSolverFeti::createInteriorMat(Mat &mat)
   {
     FUNCNAME("PetscSolverFeti::createInteriorMat()");

AMDiS/src/parallel/PetscSolverFeti.h

@@ -215,6 +215,10 @@ namespace AMDiS {
       return false;
     }
 
+    /// Traverse the mesh and create a DOF vector where each DOF is in the
+    /// order of: 1/h^2
+    void createH2vec(DOFVector<double> &vec);
+
     /// For debugging only!
     void createInteriorMat(Mat &mat);
 

AMDiS/src/parallel/PetscSolverFetiOperators.cc

@@ -301,8 +301,9 @@ namespace AMDiS {
     VecNestGetSubVec(yvec, 1, &y_lagrange);
     
    
-    VecCopy(x_interface, y_interface);
-    //    VecScale(y_interface, 0.16);
+    //    VecCopy(x_interface, y_interface);
+    //    VecScale(y_interface, 100.0);
+    VecPointwiseMult(x_interface, data->h2vec, y_interface);
 
 
     // Multiply with scaled Lagrange constraint matrix.

AMDiS/src/parallel/PetscSolverFetiStructs.h

@@ -119,6 +119,8 @@ namespace AMDiS {
     Vec tmp_vec_b1;
 
     PetscSolver* subSolver;
+
+    Vec h2vec;
   };