Added more nullspace support.

AMDiS/src/Mesh.h

@@ -402,11 +402,9 @@ namespace AMDiS {
     /// Adds a macro element to the mesh
     void addMacroElement(MacroElement* me);
 
-    /* \brief
-     * Removes a set of macro elements from the mesh. This works only for the 
-     * case, that there are no global or local refinements, i.e., all macro 
-     * elements have no children.
-     */
+    /// Removes a set of macro elements from the mesh. This works only for the 
+    /// case, that there are no global or local refinements, i.e., all macro 
+    /// elements have no children.
     void removeMacroElements(std::set<MacroElement*>& macros,
 			     vector<const FiniteElemSpace*>& feSpaces);
 
@@ -484,10 +482,8 @@ namespace AMDiS {
     }
 
 
-    /** \brief
-     * This function is equal to \ref getDofIndexCoords as defined above, but takes
-     * a DOF index instead of a DOF pointer.
-     */
+    /// This function is equal to \ref getDofIndexCoords as defined above, but
+    /// takes a DOF index instead of a DOF pointer.
     bool getDofIndexCoords(DegreeOfFreedom dof, 
 			   const FiniteElemSpace* feSpace,
 			   WorldVector<double>& coords);

AMDiS/src/ProblemStat.h

@@ -49,11 +49,9 @@ namespace AMDiS {
   };
 
 
-  /** \brief
-   * This class defines the stationary problem definition in sequential
-   * computations. For parallel computations, see 
-   * \ref ParallelProblemStatBase.
-   */
+  /// This class defines the stationary problem definition in sequential
+  /// computations. For parallel computations, see 
+  /// \ref ParallelProblemStatBase.
   class ProblemStatSeq : public ProblemStatBase,
 			 public StandardProblemIteration
   {
@@ -149,11 +147,9 @@ namespace AMDiS {
     void dualAssemble(AdaptInfo *adaptInfo, Flag flag, 
 		      bool asmMatrix = true, bool asmVector = true);
 
-    /** \brief
-     * Determines the execution order of the single adaption steps. If adapt is
-     * true, mesh adaption will be performed. This allows to avoid mesh adaption,
-     * e.g. in timestep adaption loops of timestep adaptive strategies.
-     */
+    /// Determines the execution order of the single adaption steps. If adapt is
+    /// true, mesh adaption will be performed. This allows to avoid mesh adaption,
+    /// e.g. in timestep adaption loops of timestep adaptive strategies.
     virtual Flag oneIteration(AdaptInfo *adaptInfo, Flag toDo = FULL_ITERATION);
 
     /// Returns number of managed problems
@@ -168,10 +164,8 @@ namespace AMDiS {
       return nComponents; 
     }
 
-    /** \brief
-     * Returns the problem with the given number. If only one problem
-     * is managed by this master problem, the number hasn't to be given.
-     */
+    /// Returns the problem with the given number. If only one problem
+    /// is managed by this master problem, the number hasn't to be given.
     virtual ProblemStatBase *getProblem(int number = 0) 
     { 
       return this; 
@@ -260,10 +254,8 @@ namespace AMDiS {
       addBoundaryVectorOperator(type, &op, row);
     }
 
-    /** \brief
-     * This function assembles a DOFMatrix and a DOFVector for the case,
-     * the meshes from row and col FE-space are equal.
-     */
+    /// This function assembles a DOFMatrix and a DOFVector for the case,
+    /// the meshes from row and col FE-space are equal.
     void assembleOnOneMesh(FiniteElemSpace *feSpace, 
 			   Flag assembleFlag,
 			   DOFMatrix *matrix, DOFVector<double> *vector);
@@ -533,10 +525,8 @@ namespace AMDiS {
     }
     /** \} */
 
-    /** \brief
-     * Outputs the mesh of the given component, but the values are taken from the 
-     * residual error estimator. 
-     */
+    /// Outputs the mesh of the given component, but the values are taken from
+    /// the residual error estimator. 
     void writeResidualMesh(int comp, AdaptInfo *adaptInfo, string name);
 
     /// Function that implements the serialization procedure.
@@ -553,10 +543,8 @@ namespace AMDiS {
     }
 
   protected:
-    /** \brief
-     * If the exact solution is known, the problem can compute the exact
-     * error instead of the error estimation. This is done in this function.
-     */
+    /// If the exact solution is known, the problem can compute the exact
+    /// error instead of the error estimation. This is done in this function.
     void computeError(AdaptInfo *adaptInfo);
 
   protected:
@@ -571,11 +559,9 @@ namespace AMDiS {
     /// vectors, \ref solution.
     vector<string> componentNames;
 
-    /** \brief
-     * Number of problem meshes. If all components are defined on the same mesh,
-     * this number is 1. Otherwise, this variable is the number of different meshes
-     * within the problem.
-     */
+    /// Number of problem meshes. If all components are defined on the same mesh,
+    /// this number is 1. Otherwise, this variable is the number of different meshes
+    /// within the problem.
     int nMeshes;
 
     /// FE spaces of this problem.
@@ -590,11 +576,9 @@ namespace AMDiS {
     /// Pointer to the meshes for the different problem components
     vector<Mesh*> componentMeshes;
 
-    /** \brief
-     * Stores information about which meshes must be traversed to assemble the
-     * specific components. I.e., it was implemented to make use of different
-     * meshes for different components.
-     */
+    /// Stores information about which meshes must be traversed to assemble the
+    /// specific components. I.e., it was implemented to make use of different
+    /// meshes for different components.
     ComponentTraverseInfo traverseInfo;
     
     /// Responsible for element marking.
@@ -618,21 +602,17 @@ namespace AMDiS {
     /// Composed system matrix
     SolverMatrix<Matrix<DOFMatrix*> > solverMatrix;
 
-    /** \brief
-     * Some DOFMatrices of the systemMatrix may be assembled only once (for
-     * example if they are independent of the time or older solutions). If
-     * [i][j] of this field is set to true, the corresponding DOFMatrix will
-     * be assembled only once. All other matrices will be assembled at every
-     * time step.
-     */
+    /// Some DOFMatrices of the systemMatrix may be assembled only once (for
+    /// example if they are independent of the time or older solutions). If
+    /// [i][j] of this field is set to true, the corresponding DOFMatrix will
+    /// be assembled only once. All other matrices will be assembled at every
+    /// time step.
     vector<vector<bool> > assembleMatrixOnlyOnce;
 
-    /** \brief
-     * If [i][j] of this field is set to true, the corresponding DOFMatrix of
-     * the systemMatrix has been assembled at least once. This field is used
-     * to determine, if assembling of a matrix can be ommitted, if it is set
-     * to be assembled only once.
-     */
+    /// If [i][j] of this field is set to true, the corresponding DOFMatrix of
+    /// the systemMatrix has been assembled at least once. This field is used
+    /// to determine, if assembling of a matrix can be ommitted, if it is set
+    /// to be assembled only once.
     vector<vector<bool> > assembledMatrix;
 
     /// Determines whether domain boundaries should be considered at assembling.
@@ -641,47 +621,40 @@ namespace AMDiS {
     /// Writes the meshes and solution after the adaption loop.
     vector<FileWriterInterface*> fileWriters;
 
-    /** \brief
-     * All actions of mesh refinement are performed by refinementManager.
-     * If new refinement algorithms should be realized, one has to override
-     * RefinementManager and give one instance of it to AdaptStationary.
-     */
+    /// All actions of mesh refinement are performed by refinementManager.
+    /// If new refinement algorithms should be realized, one has to override
+    /// RefinementManager and give one instance of it to AdaptStationary.
     RefinementManager *refinementManager;
 
-    /** \brief
-     * All actions of mesh coarsening are performed by coarseningManager.
-     * If new coarsening algorithms should be realized, one has to override
-     * CoarseningManager and give one instance of it to AdaptStationary.
-     */
+    /// All actions of mesh coarsening are performed by coarseningManager.
+    /// If new coarsening algorithms should be realized, one has to override
+    /// CoarseningManager and give one instance of it to AdaptStationary.
     CoarseningManager *coarseningManager;
   
     /// Info level.
     int info;
 
-    /// If true, the stationary problem was deserialized from some serialization file.
+    /// If true, the stationary problem was deserialized from some serialization
+    /// file.
     bool deserialized;
 
-    /** \brief
-     * This vectors stores pointers to functions defining the exact solution of
-     * the problem. This may be used to compute the real error of the computed
-     * solution.
-     */
+    /// This vectors stores pointers to functions defining the exact solution of
+    /// the problem. This may be used to compute the real error of the computed
+    /// solution.
     vector<AbstractFunction<double, WorldVector<double> >*> exactSolutionFcts;
 
-    /** \brief
-     * If true, the error is not estimated but computed from the exact solution
-     * defined by \ref exactSolutionFcts.
-     */
+    /// If true, the error is not estimated but computed from the exact solution
+    /// defined by \ref exactSolutionFcts.
     bool computeExactError;
     
-    /** \brief
-     * If at least on boundary condition is set, this variable is true. It is used
-     * to ensure that no operators are added after boundary condition were set. If
-     * this would happen, boundary conditions could set wrong on off diagonal matrices.
-     */
+    /// If at least on boundary condition is set, this variable is true. It is 
+    /// used to ensure that no operators are added after boundary condition were
+    /// set. If this would happen, boundary conditions could set wrong on off 
+    /// diagonal matrices.
     bool boundaryConditionSet;
 
-    /// If true, AMDiS prints information about the assembling procedure to the screen.
+    /// If true, AMDiS prints information about the assembling procedure to 
+    /// the screen.
     bool writeAsmInfo;
 
     map<Operator*, vector<OperatorPos> > operators;

AMDiS/src/parallel/ElementObjectDatabase.cc

@@ -243,6 +243,7 @@ namespace AMDiS {
 	("Should not happen (%d)!\n", multPeriodicDof2.size());
       TEST_EXIT_DBG(multPeriodicDof3.size() == 0)("Should not happen!\n");
     }
+
     if (mesh->getDim() == 3) {
       TEST_EXIT_DBG(multPeriodicDof3.size() == 0 || 
 		    multPeriodicDof3.size() == 8)

AMDiS/src/parallel/MeshDistributor.h

@@ -295,6 +295,14 @@ namespace AMDiS {
       return levelData;
     }
 
+    void updateLocalGlobalNumbering();
+
+    /// Updates the local and global DOF numbering after the mesh has been 
+    /// changed.
+    void updateLocalGlobalNumbering(ParallelDofMapping &dmap,
+				    DofComm &dcom,
+				    const FiniteElemSpace *feSpace);
+
   protected:
     void addProblemStat(ProblemStatSeq *probStat);
 
@@ -308,14 +316,6 @@ namespace AMDiS {
     /// partition.
     void removeMacroElements();
 
-    void updateLocalGlobalNumbering();
-
-    /// Updates the local and global DOF numbering after the mesh has been 
-    /// changed.
-    void updateLocalGlobalNumbering(ParallelDofMapping &dmap,
-				    DofComm &dcom,
-				    const FiniteElemSpace *feSpace);
-
     /// Calls \ref createPeriodicMap(feSpace) for all FE spaces that are
     /// handled by the mesh distributor.
     void createPeriodicMap();

AMDiS/src/parallel/PetscSolver.cc

@@ -37,6 +37,8 @@ namespace AMDiS {
 
     Parameters::get("parallel->remove rhs null space", removeRhsNullspace);
     Parameters::get("parallel->has constant null space", hasConstantNullspace);
+    Parameters::get("parallel->nullspace->const in comp", 
+		    constNullspaceComponent);
   }
 
 

AMDiS/src/parallel/PetscSolver.h

@@ -256,6 +256,8 @@ namespace AMDiS {
     bool removeRhsNullspace;
 
     bool hasConstantNullspace;
+
+    vector<int> constNullspaceComponent;
   };
 
 

AMDiS/src/parallel/PetscSolverGlobalMatrix.cc

@@ -119,9 +119,8 @@ namespace AMDiS {
     KSPSetType(kspInterior, KSPBCGS);
     KSPSetOptionsPrefix(kspInterior, kspPrefix.c_str());
     KSPSetFromOptions(kspInterior);
-    PCSetFromOptions(pcInterior);
 
-    createFieldSplit(pcInterior);
+    initPreconditioner(pcInterior);
 
     // Do not delete the solution vector, use it for the initial guess.
     if (!zeroStartVector)
@@ -131,6 +130,14 @@ namespace AMDiS {
   }
 
 
+  void PetscSolverGlobalMatrix::fillPetscMatrix(DOFMatrix *mat)
+  {
+    Matrix<DOFMatrix*> m(1, 1);
+    m[0][0] = mat;
+    fillPetscMatrix(&m);
+  }
+
+
   void PetscSolverGlobalMatrix::fillPetscMatrixWithCoarseSpace(Matrix<DOFMatrix*> *mat)
   {
     FUNCNAME("PetscSolverGlobalMatrix::fillPetscMatrixWithCoarseSpace()");
@@ -317,7 +324,6 @@ namespace AMDiS {
     KSPSetOperators(kspInterior, matIntInt, matIntInt, SAME_NONZERO_PATTERN);
     KSPSetOptionsPrefix(kspInterior, "interior_");
     KSPSetType(kspInterior, KSPPREONLY);
-    PC pcInterior;
     KSPGetPC(kspInterior, &pcInterior);
     PCSetType(pcInterior, PCLU);
     if (subdomainLevel == 0)
@@ -393,9 +399,21 @@ namespace AMDiS {
     MatNullSpace matNullspace;
     Vec nullspaceBasis;
 
-    if (nullspace.size() > 0 || hasConstantNullspace) {
+    if (nullspace.size() > 0 || 
+	hasConstantNullspace ||
+	constNullspaceComponent.size() > 0) {
       TEST_EXIT_DBG(nullspace.size() <= 1)("Not yet implemented!\n");
 
+      if (constNullspaceComponent.size() > 0) {
+	nullspace.clear();
+	SystemVector *basisVec = new SystemVector(vec);
+	basisVec->set(0.0);
+	for (unsigned int i = 0; i < constNullspaceComponent.size(); i++)
+	  basisVec->getDOFVector(constNullspaceComponent[i])->set(1.0);
+
+	nullspace.push_back(basisVec);
+      } 
+
       if (nullspace.size() > 0) {
 	VecDuplicate(petscSolVec, &nullspaceBasis);
 	for (int i = 0; i < nComponents; i++)
@@ -403,6 +421,8 @@ namespace AMDiS {
 	
 	VecAssemblyBegin(nullspaceBasis);
 	VecAssemblyEnd(nullspaceBasis);
+
+	VecNormalize(nullspaceBasis, PETSC_NULL);
 	
 	MatNullSpaceCreate(mpiCommGlobal, (hasConstantNullspace ? PETSC_TRUE : PETSC_FALSE), 
 			   1, &nullspaceBasis, &matNullspace);
@@ -415,6 +435,7 @@ namespace AMDiS {
 	MatNullSpaceCreate(mpiCommGlobal, PETSC_TRUE, 0, PETSC_NULL, &matNullspace);
       }
 
+      MatSetNullSpace(matIntInt, matNullspace);
       KSPSetNullSpace(kspInterior, matNullspace);
 
       // === Remove null space, if requested. ===
@@ -455,11 +476,6 @@ namespace AMDiS {
     }
 
     VecRestoreArray(petscSolVec, &vecPointer);
-
-
-    // === Synchronize DOFs at common DOFs, i.e., DOFs that correspond to ===
-    // === more than one partition.                                       ===
-    meshDistributor->synchVector(vec);
   }
 
 
@@ -517,6 +533,8 @@ namespace AMDiS {
   {
     FUNCNAME("PetscSolverGlobalMatrix::destroyMatrixData()");
 
+    exitPreconditioner(pcInterior);
+
     MatDestroy(&matIntInt);
     KSPDestroy(&kspInterior);
 
@@ -580,6 +598,21 @@ namespace AMDiS {
   }
 
 
+  void PetscSolverGlobalMatrix::initPreconditioner(PC pc)
+  {
+    FUNCNAME("PetscSolverGlobalMatrix::initPreconditioner()");
+
+    PCSetFromOptions(pc);
+    createFieldSplit(pc);
+  }
+
+ 
+  void PetscSolverGlobalMatrix::exitPreconditioner(PC pc)
+  {
+    FUNCNAME("PetscSolverGlobalMatrix::exitPreconditioner()");
+  }
+
+
   void PetscSolverGlobalMatrix::setDofMatrix(DOFMatrix* mat,
 					     int nRowMat, int nColMat)
   {
@@ -766,7 +799,6 @@ namespace AMDiS {
 	      else
 		perColDof = entry[i].second;	      	      
 	      
-
 	      entry.push_back(make_pair(perRowDof, perColDof));
 	    }
 	  }
@@ -815,6 +847,7 @@ namespace AMDiS {
     // Traverse all used DOFs in the dof vector.
     DOFVector<double>::Iterator dofIt(vec, USED_DOFS);
     for (dofIt.reset(); !dofIt.end(); ++dofIt) {
+
       if (rankOnly && !(*interiorMap)[feSpace].isRankDof(dofIt.getDOFIndex()))
 	continue;
 
@@ -836,7 +869,7 @@ namespace AMDiS {
 	else
 	  index =
 	    interiorMap->getMatIndex(nRowVec, dofIt.getDOFIndex()) + rStartInterior;
-	
+
 	if (perMap.isPeriodic(feSpace, globalRowDof)) {
 	  std::set<int>& perAsc = perMap.getAssociations(feSpace, globalRowDof);
 	  double value = *dofIt / (perAsc.size() + 1.0);
@@ -846,6 +879,7 @@ namespace AMDiS {
 	       perIt != perAsc.end(); ++perIt) {
 	    int mappedDof = perMap.map(feSpace, *perIt, globalRowDof);
 	    int mappedIndex = interiorMap->getMatIndex(nRowVec, mappedDof);
+
 	    VecSetValue(vecInterior, mappedIndex, value, ADD_VALUES);
 	  }	  
 	} else {	  

AMDiS/src/parallel/PetscSolverGlobalMatrix.h

@@ -50,6 +50,10 @@ namespace AMDiS {
 
     void fillPetscMatrix(Matrix<DOFMatrix*> *mat);
 
+    /// This function is just a small wrapper that creates a 1x1 matrix that
+    /// contains exactly one DOFMatrix and than calls \ref fillPetscMatrix
+    void fillPetscMatrix(DOFMatrix *mat);
+
     void fillPetscMatrixWithCoarseSpace(Matrix<DOFMatrix*> *mat);
 
     void fillPetscRhs(SystemVector *vec);
@@ -62,9 +66,13 @@ namespace AMDiS {
 
     void destroyVectorData();
 
+  protected:
     void createFieldSplit(PC pc);
 
-  protected:
+    virtual void initPreconditioner(PC pc);
+
+    virtual void exitPreconditioner(PC pc);
+
     /// Creates a new non zero pattern structure for the PETSc matrix. 
     void createPetscNnzStructure(Matrix<DOFMatrix*> *mat);