From 48e8963e949b8d84ade24d94101e2fcf36412eb0 Mon Sep 17 00:00:00 2001
From: Thomas Witkowski <thomas.witkowski@gmx.de>
Date: Mon, 17 Sep 2012 15:17:40 +0000
Subject: [PATCH] Fixed lumpred precondition for FETI-DP stokes mode.

---
 AMDiS/src/parallel/PetscSolverFeti.cc         | 59 +++++++++++++++++++
 AMDiS/src/parallel/PetscSolverFeti.h          |  4 ++
 .../src/parallel/PetscSolverFetiOperators.cc  |  5 +-
 AMDiS/src/parallel/PetscSolverFetiStructs.h   |  2 +
 4 files changed, 68 insertions(+), 2 deletions(-)

diff --git a/AMDiS/src/parallel/PetscSolverFeti.cc b/AMDiS/src/parallel/PetscSolverFeti.cc
index 34cb46b9..36899fe2 100644
--- a/AMDiS/src/parallel/PetscSolverFeti.cc
+++ b/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()");
diff --git a/AMDiS/src/parallel/PetscSolverFeti.h b/AMDiS/src/parallel/PetscSolverFeti.h
index 465f44dc..3313bab3 100644
--- a/AMDiS/src/parallel/PetscSolverFeti.h
+++ b/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);
 
diff --git a/AMDiS/src/parallel/PetscSolverFetiOperators.cc b/AMDiS/src/parallel/PetscSolverFetiOperators.cc
index 8a3dc324..d9c43ce3 100644
--- a/AMDiS/src/parallel/PetscSolverFetiOperators.cc
+++ b/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.
diff --git a/AMDiS/src/parallel/PetscSolverFetiStructs.h b/AMDiS/src/parallel/PetscSolverFetiStructs.h
index b8dd365f..820b2cf9 100644
--- a/AMDiS/src/parallel/PetscSolverFetiStructs.h
+++ b/AMDiS/src/parallel/PetscSolverFetiStructs.h
@@ -119,6 +119,8 @@ namespace AMDiS {
     Vec tmp_vec_b1;
 
     PetscSolver* subSolver;
+
+    Vec h2vec;
   };
 
 
-- 
GitLab