Removed some static variables.

AMDiS/src/DOFMatrix.cc

@@ -34,8 +34,6 @@ namespace AMDiS {
 
   using namespace mtl;
 
-  DOFMatrix *DOFMatrix::traversePtr = NULL;
-
   DOFMatrix::DOFMatrix()
     : rowFeSpace(NULL),
       colFeSpace(NULL),

AMDiS/src/DOFMatrix.h

@@ -424,9 +424,6 @@ namespace AMDiS {
     /// default compressed2D<double>
     base_matrix_type matrix;
 
-    /// Used while mesh traversal
-    static DOFMatrix *traversePtr;
-  
     /// Pointers to all operators of the equation systems. Are used in the
     /// assembling process.
     vector<Operator*> operators;

AMDiS/src/DOFVector.cc

@@ -68,7 +68,8 @@ namespace AMDiS {
 
 
   template<>
-  const double& DOFVector<double>::evalAtPoint(WorldVector<double> &p, ElInfo *oldElInfo, double* values) const
+  double DOFVector<double>::evalAtPoint(WorldVector<double> &p,
+					ElInfo *oldElInfo) const
   {
     FUNCNAME("DOFVector<double>::evalAtCoords()");
   
@@ -82,7 +83,7 @@ namespace AMDiS {
     DimVec<double> lambda(dim, NO_INIT);
   
     ElInfo *elInfo = mesh->createNewElInfo();
-    static double value = 0.0;
+    double value = 0.0;
     bool inside = false;
   
     if (oldElInfo && oldElInfo->getMacroElement()) {
@@ -107,14 +108,13 @@ namespace AMDiS {
     if (oldElInfo == NULL)
       delete elInfo;
 
-    if (values != NULL)
-      *values = value;
     return value;
-  };
+  }
 
 
   template<>
-  const WorldVector<double>& DOFVector<WorldVector<double> >::evalAtPoint(WorldVector<double> &p, ElInfo *oldElInfo, WorldVector<double>* values) const
+  WorldVector<double> DOFVector<WorldVector<double> >::evalAtPoint(WorldVector<double> &p, 
+								   ElInfo *oldElInfo) const
   {
     FUNCNAME("DOFVector<double>::evalAtCoords()");
   
@@ -129,9 +129,7 @@ namespace AMDiS {
   
     ElInfo *elInfo = mesh->createNewElInfo();
 
-    static WorldVector<double> Values(DEFAULT_VALUE, 0.0);
-    WorldVector<double> *val = (NULL != values) ? values : &Values;
-
+    WorldVector<double> value(DEFAULT_VALUE, 0.0);
     bool inside = false;
   
     if (oldElInfo && oldElInfo->getMacroElement()) {
@@ -148,7 +146,7 @@ namespace AMDiS {
       mtl::dense_vector<WorldVector<double> > uh(nBasFcts);
       for (int i = 0; i < nBasFcts; i++)
         uh[i] = operator[](localIndices[i]);
-      *val = basFcts->evalUh(lambda, uh);
+      value = basFcts->evalUh(lambda, uh);
     } else
       throw(std::runtime_error("Can not eval DOFVector at point p, because point is outside geometry."));
 
@@ -156,8 +154,8 @@ namespace AMDiS {
     if (oldElInfo == NULL)
       delete elInfo;
 
-    return ((*val));
-  };
+    return value;
+  }
 
 
   template<>

AMDiS/src/DOFVector.h

@@ -613,13 +613,11 @@ namespace AMDiS {
     /// Eval DOFVector at given point p. If oldElInfo != NULL the search for 
     /// the element, where p is inside, starts from oldElInfo. implemented for:
     /// double, WorldVector< double >
-    inline const T& evalAtPoint(WorldVector<double> &p, 
-				ElInfo *oldElInfo = NULL, 
-				T* value = NULL) const 
+    T evalAtPoint(WorldVector<double> &p, 
+		  ElInfo *oldElInfo = NULL) const 
     {
       FUNCNAME("DOFVector::evalAtPoint())");
       TEST_EXIT(false)("Please implement your evaluation\n");
-      return *value;
     }
 
     /// Determine the DegreeOfFreedom that has coords with minimal euclidean 
@@ -674,12 +672,12 @@ namespace AMDiS {
     BoundaryType boundaryType, Quadrature* q) const;
 
   template<>
-  const double& DOFVector<double>::evalAtPoint(
-    WorldVector<double> &p, ElInfo *oldElInfo, double* value) const;
+  double DOFVector<double>::evalAtPoint(WorldVector<double> &p, 
+					 ElInfo *oldElInfo) const;
 
   template<>
-  const WorldVector<double>& DOFVector<WorldVector<double> >::evalAtPoint(
-    WorldVector<double> &p, ElInfo *oldElInfo, WorldVector<double>* value) const;
+  WorldVector<double> DOFVector<WorldVector<double> >::evalAtPoint(WorldVector<double> &p, 
+								   ElInfo *oldElInfo) const;
 
   template<>
   void DOFVector<double>::refineInterpol(RCNeighbourList&, int);

AMDiS/src/parallel/MeshDistributor.cc

@@ -1599,9 +1599,9 @@ namespace AMDiS {
       elObjDb.create(partitionMap, levelData);
     elObjDb.updateRankData();
 
-    unsigned long memsize = elObjDb.calculateMemoryUsage();
-    MSG("Memory usage of element object database = %5.f KByte\n", 
-	static_cast<double>(memsize / 1024));
+    //    unsigned long memsize = elObjDb.calculateMemoryUsage();
+    // MSG("Memory usage of element object database = %5.f KByte\n", 
+    // 	static_cast<double>(memsize / 1024));
 
     intBoundary.create(levelData, 0, elObjDb);
     ParallelDebug::printBoundaryInfo(intBoundary);

AMDiS/src/parallel/PetscSolverFeti.cc

@@ -680,7 +680,7 @@ namespace AMDiS {
 
 	VecAssemblyBegin(tmpVec);
 	VecAssemblyEnd(tmpVec);
-		
+	
 	subdomain->solve(tmpVec, tmpVec);
 
 	PetscScalar *tmpValues;
@@ -1018,9 +1018,8 @@ namespace AMDiS {
     MatMult(subdomain->getMatInterior(), ktest0, ktest1);
     
     PetscScalar *valarray;
-    VecGetArray(ktest1, &valarray);
-    
     Vec ktest2, ktest3;
+    VecGetArray(ktest1, &valarray);
     VecCreateMPIWithArray(PETSC_COMM_WORLD, 1, 
 			  localDofMap.getRankDofs(), nGlobalOverallInterior, 
 			  valarray, &ktest2);
@@ -1046,14 +1045,26 @@ namespace AMDiS {
     KSPSetNullSpace(ksp_feti, matNullSpace);
     MatNullSpaceDestroy(&matNullSpace);
 
-#if (DEBUG != 0)
+    //#if (DEBUG != 0)
     Vec vecSol;
-    VecDuplicate(nullSpaceBasis, &vecSol);
+    VecDuplicate(nullSpaceBasis, &vecSol);    
     MatMult(mat_feti, nullSpaceBasis, vecSol);
-    PetscReal norm;
+
+    Vec vecSol0, vecSol1;
+    VecNestGetSubVec(vecSol, 0, &vecSol0);
+    VecNestGetSubVec(vecSol, 1, &vecSol1);
+
+    PetscReal norm, norm0, norm1;
     VecNorm(vecSol, NORM_2, &norm);
-    MSG("NORM: %e\n", norm);  
-#endif
+    VecNorm(vecSol0, NORM_2, &norm0);
+    VecNorm(vecSol1, NORM_2, &norm1);
+
+    MSG("NORM: %e (%e %e)\n", norm, norm0, norm1);  
+    VecDestroy(&vecSol);    
+    //#endif
+
+    // AMDiS::finalize();
+    // exit(0);
 
     VecDestroy(&ktest0);
     VecDestroy(&ktest1);
@@ -1692,6 +1703,8 @@ namespace AMDiS {
   {
     FUNCNAME("PetscSolverFeti::solveReducedFetiMatrix()");
 
+    debugNullSpace(vec);
+
     // === Some temporary vectors. ===
 
     Vec tmp_b0, tmp_b1;
@@ -2017,4 +2030,96 @@ namespace AMDiS {
 		    &(nestMat[0]), &mat);
     }
   }
+
+
+  void PetscSolverFeti::debugNullSpace(SystemVector &vec)
+  {
+    FUNCNAME("PetscSolverFeti::debugNullSpace()");
+
+    TEST_EXIT(stokesMode)("This function works only for the stokes mode!\n");
+
+    Mat fetiMat;
+    createNestedFetiMat(fetiMat);
+
+    Vec vecArray[4];
+    
+
+    Vec ktest0, ktest1, ktest2;
+    localDofMap.createLocalVec(ktest0);
+    localDofMap.createLocalVec(ktest1);
+    localDofMap.createVec(ktest2, nGlobalOverallInterior);
+
+    DofMap& m = localDofMap[pressureFeSpace].getMap();
+    for (DofMap::iterator it = m.begin(); it != m.end(); ++it) {
+      if (meshDistributor->getDofMap()[pressureFeSpace].isRankDof(it->first)) {
+	int index = localDofMap.getLocalMatIndex(pressureComponent, it->first);
+	VecSetValue(ktest0, index, 1.0, INSERT_VALUES);
+      }
+    }    
+    VecAssemblyBegin(ktest0);
+    VecAssemblyEnd(ktest0);
+    MatMult(subdomain->getMatInterior(), ktest0, ktest1);
+    
+    PetscScalar *valarray;
+    VecGetArray(ktest0, &valarray);
+    VecCreateMPIWithArray(PETSC_COMM_WORLD, 1, 
+			  localDofMap.getRankDofs(), nGlobalOverallInterior, 
+			  valarray, &vecArray[0]);
+
+    Vec ktest3;
+    VecGetArray(ktest1, &valarray);
+    VecCreateMPIWithArray(PETSC_COMM_WORLD, 1, 
+			  localDofMap.getRankDofs(), nGlobalOverallInterior, 
+			  valarray, &ktest3);
+
+    primalDofMap.createVec(vecArray[1]);
+    VecSet(vecArray[1], 0.0);
+    
+    interfaceDofMap.createVec(vecArray[2]);
+    VecSet(vecArray[2], 1.0);
+
+    lagrangeMap.createVec(vecArray[3]);
+    MatMult(subdomain->getMatInteriorCoarse(1), vecArray[2], ktest2);
+    VecAXPY(ktest2, 1.0, ktest3);
+    MatMult(mat_lagrange_scaled, ktest2, vecArray[3]);
+    VecScale(vecArray[3], -1.0);
+
+    
+    Vec nullSpaceBasis;
+    VecCreateNest(mpiCommGlobal, 4, PETSC_NULL, vecArray, &nullSpaceBasis);
+
+    Vec vecSol;
+    VecDuplicate(nullSpaceBasis, &vecSol);
+    
+    MatMult(fetiMat, nullSpaceBasis, vecSol);
+    PetscReal norm;
+    VecNorm(vecSol, NORM_2, &norm);
+    MSG("Null space norm: %e\n", norm);  
+
+    Vec vec0;
+    Vec vec1;
+    Vec vec2;
+    Vec vec3;
+    VecNestGetSubVec(vecSol, 0, &vec0);
+    VecNestGetSubVec(vecSol, 1, &vec1);
+    VecNestGetSubVec(vecSol, 2, &vec2);
+    VecNestGetSubVec(vecSol, 3, &vec3);
+
+    PetscReal norm0, norm1, norm2, norm3;
+    VecNorm(vec0, NORM_2, &norm0);
+    VecNorm(vec1, NORM_2, &norm1);
+    VecNorm(vec2, NORM_2, &norm2);
+    VecNorm(vec3, NORM_2, &norm3);
+    MSG("Splitted null space norm: %e %e %e %e\n", norm0, norm1, norm2, norm3);  
+
+
+    recoverSolution(vec0, vec1, vec);
+    recoverInterfaceSolution(vec2, vec);
+
+    VtkWriter::writeFile(&vec, "nullspace.vtu");
+
+    MatDestroy(&fetiMat);
+    VecDestroy(&nullSpaceBasis);
+    VecDestroy(&vecSol);
+  }
 }

AMDiS/src/parallel/PetscSolverFeti.h

@@ -221,6 +221,9 @@ namespace AMDiS {
     /// For debugging only!
     void createNestedFetiMat(Mat &mat);
 
+    /// For debugging only!
+    void debugNullSpace(SystemVector &vec);
+
   protected:
     /// Mapping from primal DOF indices to a global index of primals.
     ParallelDofMapping primalDofMap;