Fixed 0-diagonal problem in parallel computations.

AMDiS/libtool

@@ -44,7 +44,7 @@ available_tags=" CXX F77"
 
 # ### BEGIN LIBTOOL CONFIG
 
-# Libtool was configured on host deimos103:
+# Libtool was configured on host deimos102:
 
 # Shell to use when invoking shell scripts.
 SHELL="/bin/sh"
@@ -6760,7 +6760,7 @@ build_old_libs=`case $build_libtool_libs in yes) $echo no;; *) $echo yes;; esac`
 # End:
 # ### BEGIN LIBTOOL TAG CONFIG: CXX
 
-# Libtool was configured on host deimos103:
+# Libtool was configured on host deimos102:
 
 # Shell to use when invoking shell scripts.
 SHELL="/bin/sh"
@@ -7065,7 +7065,7 @@ include_expsyms=""
 
 # ### BEGIN LIBTOOL TAG CONFIG: F77
 
-# Libtool was configured on host deimos103:
+# Libtool was configured on host deimos102:
 
 # Shell to use when invoking shell scripts.
 SHELL="/bin/sh"

AMDiS/src/parallel/GlobalMatrixSolver.cc

@@ -4,7 +4,7 @@
 #include "DOFVector.h"
 #include "Debug.h"
 #include "SystemVector.h"
-
+#include "VtkWriter.h"
 #include "petscksp.h"
 
 namespace AMDiS {
@@ -94,6 +94,8 @@ namespace AMDiS {
       int globalRowDof = meshDistributor->mapLocalToGlobal(*cursor);
       // Test if the current row dof is a periodic dof.
       bool periodicRow = meshDistributor->isPeriodicDof(globalRowDof);
+      // Calculate petsc row index.
+      int rowIndex = globalRowDof * dispMult + dispAddRow;
 
 
       // === Traverse all non zero entries of the row and produce vector cols ===
@@ -103,13 +105,14 @@ namespace AMDiS {
       for (icursor_type icursor = begin<nz>(cursor), icend = end<nz>(cursor); 
 	   icursor != icend; ++icursor) {
 	
-	// Set only non null values.
-	if (value(*icursor) != 0.0) {
 	// Global index of the current column index.
 	int globalColDof = meshDistributor->mapLocalToGlobal(col(*icursor));
 	// Calculate the exact position of the column index in the petsc matrix.
 	int colIndex = globalColDof * dispMult + dispAddCol;
 
+	// Set only non null values.
+	if (value(*icursor) != 0.0 || rowIndex == colIndex) {
+
 	  // If the current row is not periodic, but the current dof index is periodic,
 	  // we have to duplicate the value to the other corresponding periodic columns.
  	  if (!periodicRow && meshDistributor->isPeriodicDof(globalColDof)) {
@@ -141,9 +144,6 @@ namespace AMDiS {
       // === Up to now we have assembled on row. Now, the row must be send to the ===
       // === corresponding rows to the petsc matrix.                              ===
      
-      // Calculate petsc row index.
-      int rowIndex = globalRowDof * dispMult + dispAddRow;
-      
       if (periodicRow) {
 	// The row dof is periodic, so send dof to all the corresponding rows.
 	std::set<int>& perAsc = meshDistributor->getPerDofAssociations(globalRowDof);
@@ -174,6 +174,7 @@ namespace AMDiS {
 	  }
 
 	  int perRowIndex = (meshDistributor->getPeriodicMapping(*perIt, globalRowDof) * dispMult) + dispAddRow;
+
 	  MatSetValues(petscMatrix, 1, &perRowIndex, perCols.size(), 
 		       &(perCols[0]), &(perValues[0]), ADD_VALUES);
 	}

AMDiS/src/parallel/MeshDistributor.cc

@@ -1815,10 +1815,13 @@ namespace AMDiS {
     debug::getAllDofs(feSpace, testDofs);
     for (std::set<const DegreeOfFreedom*>::iterator it = testDofs.begin();
 	 it != testDofs.end(); ++it) {
+//       WorldVector<double> coords;
+//       mesh->getDofIndexCoords(**it, feSpace, coords);
       MSG("DOF %d:   mapLocalGlobalDofs = %d   vertexDof = %d    isRankDof = %d\n", **it, 
 	  mapLocalGlobalDofs[**it], 
 	  vertexDof[*it],
 	  isRankDof[**it]);
+      //      MSG("coord = %f %f\n", coords[0], coords[1]);
     }
 #endif
 #endif