Merge to JUROPA.

AMDiS/src/parallel/PetscSolverFeti.cc

@@ -273,7 +273,14 @@ namespace AMDiS {
     for (DofContainerSet::iterator it = vertices.begin(); 
 	 it != vertices.end(); ++it) {
       if (meshLevel == 0) {
-	primals.insert(**it);
+	MSG("WARNING: Modified primal detection algorithm!\n");
+
+	double e = 1e-3;
+	WorldVector<double> c;
+	feSpace->getMesh()->getDofIndexCoords(*it, feSpace, c);
+
+	if (fabs(c[0]) > e && fabs(c[1]) > e && fabs(c[0] - 1.0) > e && fabs(c[1] - 1.0) > e)
+	  primals.insert(**it);
       } else {
 	double e = 1e-8;
 	WorldVector<double> c;
@@ -909,10 +916,11 @@ namespace AMDiS {
 	    lumpedData->localToDualMap[matIndexLocal] = matIndexDual;
 	  }
 	}
-	
+
 	if (stokesMode) {
 	  DOFVector<double> tmpvec(pressureFeSpace, "tmpvec");
 	  createH2vec(tmpvec);
+	  VtkWriter::writeFile(&tmpvec, "tmpvec.vtu");
 	  interfaceDofMap.createVec(fetiInterfaceLumpedPreconData.h2vec);
 
 	  DofMap& m = interfaceDofMap[pressureFeSpace].getMap();
@@ -926,14 +934,23 @@ namespace AMDiS {
 
 	  VecAssemblyBegin(fetiInterfaceLumpedPreconData.h2vec);
 	  VecAssemblyEnd(fetiInterfaceLumpedPreconData.h2vec);
-	  
+  
 	  localDofMap.createVec(fetiInterfaceLumpedPreconData.tmp_vec_b1);
+	  primalDofMap.createVec(fetiInterfaceLumpedPreconData.tmp_primal);
 	  fetiInterfaceLumpedPreconData.subSolver = subdomain;	
 	}
 	
 	KSPGetPC(ksp_feti, &precon_feti);
 	PCSetType(precon_feti, PCSHELL);
 	if (stokesMode) {
+	  KSPCreate(PETSC_COMM_WORLD, &fetiInterfaceLumpedPreconData.ksp_primal);
+	  KSPSetOperators(fetiInterfaceLumpedPreconData.ksp_primal,
+			  subdomain->getMatCoarse(0, 0),
+			  subdomain->getMatCoarse(0, 0), 
+			  SAME_NONZERO_PATTERN);
+	  KSPSetOptionsPrefix(fetiInterfaceLumpedPreconData.ksp_primal, "primal_");
+	  KSPSetFromOptions(fetiInterfaceLumpedPreconData.ksp_primal);
+
 	  PCShellSetContext(precon_feti, static_cast<void*>(&fetiInterfaceLumpedPreconData));
 	  PCShellSetApply(precon_feti, petscApplyFetiInterfaceLumpedPrecon);
 	} else {
@@ -1968,26 +1985,40 @@ namespace AMDiS {
 
     Mesh *mesh = vec.getFeSpace()->getMesh();    
     TEST_EXIT(mesh->getDim() == 2)("This works only in 2D!\n");
+    int n0 = vec.getFeSpace()->getAdmin()->getNumberOfPreDofs(VERTEX);
 
     TraverseStack stack;
     ElInfo *elInfo = 
-      stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS);
+      stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_DET);
     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);
+       	DegreeOfFreedom dof = el->getDof(i, n0);
+      	vec[dof] += elInfo->getDet() * 0.5;
+       	cvec[dof] += 1;
+      }
+
 
-	WorldVector<double> c0 = elInfo->getCoord(d0);
-	WorldVector<double> c1 = elInfo->getCoord(d1);
+      /*
+      
+      for (int i = 0; i < 3; i++) {
+	int idx0 = el->getVertexOfEdge(i, 0);
+	int idx1 = el->getVertexOfEdge(i, 1);
+	DegreeOfFreedom dof0 = el->getDof(idx0, n0);
+	DegreeOfFreedom dof1 = el->getDof(idx1, n0);
+
+	WorldVector<double> c0 = elInfo->getCoord(idx0);
+	WorldVector<double> c1 = elInfo->getCoord(idx1);
 	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;
+	vec[dof0] += length;
+	vec[dof1] += length;
+	cvec[dof0] += 1;
+	cvec[dof1] += 1;
       }
-
+      */
+      
+      
       elInfo = stack.traverseNext(elInfo);
     }
 

AMDiS/src/parallel/PetscSolverFetiOperators.cc

@@ -299,17 +299,22 @@ namespace AMDiS {
     VecNestGetSubVec(xvec, 1, &x_lagrange);
     VecNestGetSubVec(yvec, 0, &y_interface);
     VecNestGetSubVec(yvec, 1, &y_lagrange);
-    
-   
-    //    VecCopy(x_interface, y_interface);
-    //    VecScale(y_interface, 100.0);
-    VecPointwiseMult(x_interface, data->h2vec, y_interface);
 
+#if 0
+    MatMult(data->subSolver->getMatCoarse(0, 1), x_interface, data->tmp_primal);
+    KSPSolve(data->ksp_primal, data->tmp_primal, data->tmp_primal);
+    MatMult(data->subSolver->getMatCoarse(1, 0), data->tmp_primal, y_interface);
 
-    // Multiply with scaled Lagrange constraint matrix.
     MatMult(data->subSolver->getMatInteriorCoarse(1), x_interface, data->tmp_vec_b0);
     MatMultTranspose(*(data->mat_lagrange_scaled), x_lagrange, data->tmp_vec_b1);
     VecAXPY(data->tmp_vec_b0, 1.0, data->tmp_vec_b1);
+#else 
+    VecCopy(x_interface, y_interface);
+    VecScale(y_interface, 1024.0);
+    //    VecPointwiseMult(y_interface, x_interface, data->h2vec);
+
+    MatMultTranspose(*(data->mat_lagrange_scaled), x_lagrange, data->tmp_vec_b0);
+#endif
 
 
     // === Restriction of the B nodes to the boundary nodes. ===
@@ -352,7 +357,12 @@ namespace AMDiS {
     // Multiply with scaled Lagrange constraint matrix.
     MatMult(*(data->mat_lagrange_scaled), data->tmp_vec_b0, y_lagrange);
 
-    //  MatMult(data->subSolver->getMatCoarseInterior(1), data->tmp_vec_b0, y_interface);
+#if 0
+    Vec tmp_interface;
+    VecDuplicate(y_interface, &tmp_interface);
+    MatMult(data->subSolver->getMatCoarseInterior(1), data->tmp_vec_b0, tmp_interface);
+    VecAXPY(y_interface, 1.0, tmp_interface);
+#endif
 
     return 0;
   }

AMDiS/src/parallel/PetscSolverFetiStructs.h

@@ -121,6 +121,10 @@ namespace AMDiS {
     PetscSolver* subSolver;
 
     Vec h2vec;
+
+    KSP ksp_primal;
+
+    Vec tmp_primal;
   };