Fixed lagrange scaling for general 3D problems.

AMDiS/src/ProblemStat.cc

@@ -65,7 +65,8 @@ namespace AMDiS {
       writeAsmInfo(false)
   {
     Parameters::get(name + "->components", nComponents);
-    TEST_EXIT(nComponents > 0)("components not set!\n");    
+    TEST_EXIT(nComponents > 0)("No value set for parameter \"%s->components\"!\n",
+			       name.c_str());    
     estimator.resize(nComponents, NULL);
     marker.resize(nComponents, NULL);
     

AMDiS/src/parallel/PetscSolverFeti.cc

@@ -254,8 +254,9 @@ namespace AMDiS {
 	    interfaceDofMap[feSpace].nRankDofs, 
 	    interfaceDofMap[feSpace].nOverallDofs);
       } else {
-	MSG("  nRankPrimals = %d   nOverallPrimals = %d\n", 
+	MSG("  nRankPrimals = %d   nLocalPrimals = %d  nOverallPrimals = %d\n", 
 	    primalDofMap[feSpace].nRankDofs, 
+	    primalDofMap[feSpace].nLocalDofs,
 	    primalDofMap[feSpace].nOverallDofs);
 	
 	MSG("  nRankDuals = %d  nOverallDuals = %d\n",
@@ -606,6 +607,9 @@ namespace AMDiS {
 		 &mat_lagrange);
     MatSetOption(mat_lagrange, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
 
+    Vec vec_scale_lagrange;
+    lagrangeMap.createVec(vec_scale_lagrange);
+
     // === Create for all duals the corresponding Lagrange constraints. On ===
     // === each rank we traverse all pairs (n, m) of ranks, with n < m,    ===
     // === that contain this node. If the current rank number is r, and    ===
@@ -631,16 +635,30 @@ namespace AMDiS {
 
 	TEST_EXIT_DBG(degree > 1)("Should not happen!\n");
 	
+	int counter = 0;
 	for (int i = 0; i < degree; i++) {
 	  for (int j = i + 1; j < degree; j++) {
 	    if (W[i] == mpiRank || W[j] == mpiRank) {
-	      int colIndex = 
-		localDofMap.getMatIndex(k, it->first) + rStartInterior;
-
-	      double value = (W[i] == mpiRank ? 1.0 : -1.0);
-	      MatSetValue(mat_lagrange, index, colIndex, value, INSERT_VALUES);
+	      MatSetValue(mat_lagrange, 
+			  index + counter, 
+			  localDofMap.getMatIndex(k, it->first) + rStartInterior,
+			  (W[i] == mpiRank ? 1.0 : -1.0),
+			  INSERT_VALUES);
 	    }
-	    index++;	      
+	    counter++;
+	  }
+	}
+
+	// === Create scaling factors for scaling the lagrange matrix, which ===
+	// === is required for FETI-DP preconditioners.                      ===
+	
+	if (meshDistributor->getDofMap()[feSpaces[k]].isRankDof(it->first)) {
+	  int nConstraints = (degree * (degree - 1)) / 2;
+	  for (int i = 0; i < nConstraints; i++) {
+	    VecSetValue(vec_scale_lagrange,
+			index + i,
+			1.0 / static_cast<double>(degree),
+			INSERT_VALUES);
 	  }
 	}
       }
@@ -649,6 +667,24 @@ namespace AMDiS {
     MatAssemblyBegin(mat_lagrange, MAT_FINAL_ASSEMBLY);
     MatAssemblyEnd(mat_lagrange, MAT_FINAL_ASSEMBLY);
 
+
+    // === If required, create \ref mat_lagrange_scaled ===
+
+    VecAssemblyBegin(vec_scale_lagrange);
+    VecAssemblyEnd(vec_scale_lagrange);
+
+    VecView(vec_scale_lagrange, PETSC_VIEWER_STDOUT_WORLD);
+
+    if (fetiPreconditioner != FETI_NONE || stokesMode) {
+      MatDuplicate(mat_lagrange, MAT_COPY_VALUES, &mat_lagrange_scaled);
+      MatDiagonalScale(mat_lagrange_scaled, vec_scale_lagrange, PETSC_NULL);
+    }
+
+    VecDestroy(&vec_scale_lagrange);
+
+
+    // === Print final timings. ===
+
     if (printTimings) {
       MPI::COMM_WORLD.Barrier();
       MSG("FETI-DP timing 05: %.5f seconds (creation of lagrange constraint matrix)\n", 
@@ -971,16 +1007,6 @@ namespace AMDiS {
     if (stokesMode)
       KSPMonitorSet(ksp_feti, KSPMonitorFetiStokes, &fetiKspData, PETSC_NULL);
 
-
-
-    // === Create FETI-DP preconditioner object. ===
-
-    if (fetiPreconditioner != FETI_NONE || stokesMode) {
-      MatDuplicate(mat_lagrange, MAT_COPY_VALUES, &mat_lagrange_scaled);
-      MatScale(mat_lagrange_scaled, 0.5);
-    }
-
-
     // === Create null space objects. ===
     
     createNullSpace();