Cleanup of last big merge!

AMDiS/bin/Makefile.am

@@ -65,25 +65,17 @@ endif
 
 libamdis_la_SOURCES =  \
 $(PARALLEL_AMDIS_SOURCES) \
-$(SOURCE_DIR)/MultiGridPreconWrapper.h $(SOURCE_DIR)/MultiGridPreconWrapper.cc \
-$(SOURCE_DIR)/LagrangeInterpolRestrict.h \
-$(SOURCE_DIR)/InterpolRestrictMatrix.h $(SOURCE_DIR)/InterpolRestrictMatrix.cc \
 $(SOURCE_DIR)/DOFIndexed.h $(SOURCE_DIR)/DOFIndexed.cc \
 $(SOURCE_DIR)/GNUPlotWriter.h $(SOURCE_DIR)/GNUPlotWriter.cc \
-$(SOURCE_DIR)/SmootherBase.h \
-$(SOURCE_DIR)/GSSMoother.h $(SOURCE_DIR)/GSSmoother.cc \
-$(SOURCE_DIR)/JacobiSMoother.h $(SOURCE_DIR)/JacobiSmoother.cc \
-$(SOURCE_DIR)/SparseVector.h $(SOURCE_DIR)/SparseVector.hh $(SOURCE_DIR)/SparseVector.cc \
 $(SOURCE_DIR)/VertexVector.h $(SOURCE_DIR)/VertexVector.cc \
 $(SOURCE_DIR)/PeriodicBC.h $(SOURCE_DIR)/PeriodicBC.cc \
-$(SOURCE_DIR)/MultiGridSolverBase.h $(SOURCE_DIR)/MultiGridSolverBase.hh \
-$(SOURCE_DIR)/MultiGridSolver.h $(SOURCE_DIR)/MultiGridSolver.cc \
 $(SOURCE_DIR)/Recovery.h $(SOURCE_DIR)/Recovery.cc \
 $(SOURCE_DIR)/RecoveryEstimator.h $(SOURCE_DIR)/RecoveryEstimator.cc \
 $(SOURCE_DIR)/ResidualEstimator.h $(SOURCE_DIR)/ResidualEstimator.cc \
 $(SOURCE_DIR)/Cholesky.h $(SOURCE_DIR)/Cholesky.cc \
 $(SOURCE_DIR)/AdaptBase.h $(SOURCE_DIR)/AdaptBase.cc \
-$(SOURCE_DIR)/ProblemIterationInterface.h $(SOURCE_DIR)/StandardProblemIteration.h $(SOURCE_DIR)/StandardProblemIteration.cc \
+$(SOURCE_DIR)/ProblemIterationInterface.h \
+$(SOURCE_DIR)/StandardProblemIteration.h $(SOURCE_DIR)/StandardProblemIteration.cc \
 $(SOURCE_DIR)/ProblemScal.h $(SOURCE_DIR)/ProblemScal.cc \
 $(SOURCE_DIR)/ProblemVec.h $(SOURCE_DIR)/ProblemVec.cc \
 $(SOURCE_DIR)/DualTraverse.h $(SOURCE_DIR)/DualTraverse.cc \
@@ -95,7 +87,9 @@ $(SOURCE_DIR)/CreatorInterface.h \
 $(SOURCE_DIR)/ElementFunction.h \
 $(SOURCE_DIR)/ProblemInterpolScal.h $(SOURCE_DIR)/ProblemInterpolScal.cc \
 $(SOURCE_DIR)/ProblemInterpolVec.h $(SOURCE_DIR)/ProblemInterpolVec.cc \
-$(SOURCE_DIR)/Serializable.h $(SOURCE_DIR)/BallProject.h $(SOURCE_DIR)/CylinderProject.h \
+$(SOURCE_DIR)/Serializable.h \
+$(SOURCE_DIR)/BallProject.h \
+$(SOURCE_DIR)/CylinderProject.h \
 $(SOURCE_DIR)/MacroReader.h $(SOURCE_DIR)/MacroReader.cc \
 $(SOURCE_DIR)/ValueReader.h $(SOURCE_DIR)/ValueReader.cc \
 $(SOURCE_DIR)/Projection.h $(SOURCE_DIR)/Projection.cc \
@@ -110,18 +104,17 @@ $(SOURCE_DIR)/SystemVector.h \
 $(SOURCE_DIR)/MatrixVector.h $(SOURCE_DIR)/MatrixVector.cc \
 $(SOURCE_DIR)/SurfaceQuadrature.h $(SOURCE_DIR)/SurfaceQuadrature.cc \
 $(SOURCE_DIR)/LeafData.h $(SOURCE_DIR)/LeafData.cc \
-$(SOURCE_DIR)/BoundaryManager.h $(SOURCE_DIR)/BoundaryManager.cc $(SOURCE_DIR)/BoundaryCondition.h \
+$(SOURCE_DIR)/BoundaryManager.h $(SOURCE_DIR)/BoundaryManager.cc \
+$(SOURCE_DIR)/BoundaryCondition.h \
 $(SOURCE_DIR)/DirichletBC.h $(SOURCE_DIR)/DirichletBC.cc \
 $(SOURCE_DIR)/RobinBC.h $(SOURCE_DIR)/RobinBC.cc \
 $(SOURCE_DIR)/AbstractFunction.h \
 $(SOURCE_DIR)/ProblemStatBase.h \
-$(SOURCE_DIR)/MatVecMultiplier.cc \
 $(SOURCE_DIR)/DOFContainer.h \
 $(SOURCE_DIR)/FileWriter.h $(SOURCE_DIR)/FileWriter.hh $(SOURCE_DIR)/FileWriter.cc \
 $(SOURCE_DIR)/ElementFileWriter.h $(SOURCE_DIR)/ElementFileWriter.cc \
 $(SOURCE_DIR)/ElInfo.cc \
 $(SOURCE_DIR)/ElInfoStack.h $(SOURCE_DIR)/ElInfoStack.cc \
-$(SOURCE_DIR)/MatVecMultiplier.h \
 $(SOURCE_DIR)/Operator.h $(SOURCE_DIR)/Operator.cc \
 $(SOURCE_DIR)/Mesh.cc \
 $(SOURCE_DIR)/AMDiS.h \

AMDiS/src/AMDiS.h

@@ -48,7 +48,6 @@
 #include "MacroElement.h"
 #include "MacroWriter.h"
 #include "Marker.h"
-#include "MatVecMultiplier.h"
 #include "MatrixVector.h"
 #include "MemoryManager.h"
 #include "MemoryPool.h"
@@ -83,9 +82,7 @@
 #include "RefinementManager2d.h"
 #include "RefinementManager3d.h"
 #include "RobinBC.h"
-#include "SmootherBase.h"
 #include "SolutionDataStorage.h"
-#include "SparseVector.h"
 #include "SurfaceOperator.h"
 #include "SurfaceQuadrature.h"
 #include "SystemVector.h"

AMDiS/src/AMDiS_fwd.h

@@ -44,13 +44,11 @@ namespace AMDiS {
   class FileWriterInterface;
   class FiniteElemSpace; 
   class IdentityPreconditioner; 
-  class InterpolRestrictMatrix;
   class ITL_BasePreconditioner;
   class LeafDataPeriodic;
   class LevelAdmin;
   class Marker;
   class Mesh; 
-  class MultiGridSortSolver; 
   class OEMSolver;
   class Operator;
   class ProblemInstat;
@@ -65,14 +63,12 @@ namespace AMDiS {
   class VertexVector;
 
   template<typename ReturnType, typename ArgumentType> class AbstractFunction;
-  template<typename VectorType>                        class MatVecMultiplier;
   template<typename T>                                 class DOFVectorBase;
   template<typename T>                                 class DOFVector;
   template<typename T>                                 class ITL_Preconditioner;
   template<typename T>                                 class Matrix;
   template<typename T>                                 class MatVecMultiplier;
   template<typename T>                                 class SolverMatrix;
-  template<typename T>                                 class SparseVector;
   template<typename T>                                 class Vector;
 
   template <typename ITLSolver>    class ITL_OEMSolver;

AMDiS/src/BoxSmoother.h

-#ifndef AMDIS_BOXSMOOTHER_H
-#define AMDIS_BOXSMOOTHER_H
-
-#include "SmootherBase.h"
-#include "MemoryManager.h"
-#include "DOFMatrix.h"
-#include "SparseVector.h"
-#include "MatrixVector.h"
-#include "Global.h"
-#include "GaussElimination.h"
-#include "FiniteElemSpace.h"
-#include "GaussElimination.h"
-#include <set>
-
-namespace AMDiS {
-
-  class BoxSmoother 
-    : public SmootherBase<Matrix<DOFMatrix*>, 
-			  Vector<SparseVector<double>*>,
-			  Vector< std::set<DegreeOfFreedom>*> >
-  {
-  public:
-    MEMORY_MANAGED(BoxSmoother);
-
-    class Creator 
-      : public SmootherCreator<Matrix<DOFMatrix*>, 
-			       Vector<SparseVector<double>*>,
-			       Vector< std::set<DegreeOfFreedom>*> >
-    {
-    public:
-      MEMORY_MANAGED(Creator);
-
-      SmootherBase<Matrix<DOFMatrix*>, 
-		   Vector<SparseVector<double>*>,
-		   Vector< std::set<DegreeOfFreedom>*> >*
-      create() {
-	return NEW BoxSmoother(this->name);
-      };
-    };
-
-    BoxSmoother(std::string name)
-      : SmootherBase<Matrix<DOFMatrix*>, 
-		     Vector<SparseVector<double>*>,
-		     Vector< std::set<DegreeOfFreedom>*> >(name),
-	boxMatrix_(NULL),
-	boxSolution_(NULL),
-	boxRHS_(NULL),
-	numComponents_(0),
-	sumDOFs_(0)
-    {
-      directSolver_ = 
-	NEW GaussElimination<Matrix<double>, Vector<double> >();
-    };
- 
-    virtual ~BoxSmoother()
-    {
-      if(boxMatrix_) DELETE boxMatrix_;
-      if(boxSolution_) DELETE boxSolution_;
-      if(boxRHS_) DELETE boxRHS_;
-      DELETE directSolver_;
-    };
-
-    void smooth(Matrix<DOFMatrix*>                       *matrix,
-		Vector<SparseVector<double>*>            *solution,
-		Vector<SparseVector<double>*>            *rhs,
-		int                                       iterations,
-		const Vector< std::set<DegreeOfFreedom>*> &dofSet)
-    {
-      FUNCNAME("BoxSmoother::smooth()");
-      ERROR_EXIT("Will be rewritten - Peter\n");
-
-#if 0
-      int i, j, k;
-
-      // allocate memory
-      if(!boxMatrix_) {
-	numComponents_ = solution->getSize();
-	numDOFs_.resize(numComponents_);
-	sumDOFs_ = 0;
-	for(i = 0; i < numComponents_; i++) {
-	  numDOFs_[i] = (*solution)[i]->getFESpace()->getBasisFcts()->getNumber();
-	  sumDOFs_ += numDOFs_[i];
-	}
-
-	boxMatrix_ = NEW Matrix<double>(sumDOFs_, sumDOFs_);
-	boxSolution_ = NEW Vector<double>(sumDOFs_);
-	boxRHS_ = NEW Vector<double>(sumDOFs_);
-      }
-
-      // create local equation system
-      int boxSize = 0;
-      std::vector<int> offset(numComponents_, 0);
-      std::vector< std::map<DegreeOfFreedom, int> > indexOfDOF(numComponents_);
-      std::set<DegreeOfFreedom>::const_iterator dofIt, dofBegin, dofEnd;
-      std::vector<MatEntry>::iterator rowIt, rowBegin, rowEnd;
-
-      int pos;
-      for(i = 0; i < numComponents_; i++) {
-	if(dofSet[i]) {
-	  offset[i] = boxSize;
-	  boxSize += numDOFs_[i];
-	  dofBegin = dofSet[i]->begin();
-	  dofEnd = dofSet[i]->end();
-	  pos = 0;
-	  for(dofIt = dofBegin; dofIt != dofEnd; ++dofIt) {
-	    indexOfDOF[i][*dofIt] = pos + 1;
-// 	    MSG("%d : %d, %d \n", *dofIt, i, pos);
-	    ++pos;
-	  }
-	}
-      }
-
-      for(i = 0; i < boxSize; i++) {
-	for(j = 0; j < boxSize; j++) {
-	  (*boxMatrix_)[i][j] = 0.0;
-	}
-	(*boxRHS_)[i] = 0.0;
-      }
-
-      pos = 0;
-      for(i = 0; i < numComponents_; i++) {
-	if(dofSet[i]) {
-	  dofBegin = dofSet[i]->begin();
-	  dofEnd = dofSet[i]->end();
-	  for(dofIt = dofBegin; dofIt != dofEnd; ++dofIt) {
-
-	    // fill box rhs
-	    (*boxRHS_)[pos] = (*((*rhs)[i]))[*dofIt];
-
-	    // fill box matrix
-	    for(j = 0; j < numComponents_; j++) {
-	      if((*matrix)[i][j] != NULL) {
-		rowBegin = (*((*matrix)[i][j]))[*dofIt].begin();
-		rowEnd = (*((*matrix)[i][j]))[*dofIt].end();
-		for(rowIt = rowBegin; rowIt != rowEnd; ++rowIt) {
-		  if(rowIt->col == DOFMatrix::UNUSED_ENTRY) continue;
-		  if(rowIt->col == DOFMatrix::NO_MORE_ENTRIES) break;
-		  int ind = indexOfDOF[j][rowIt->col] - 1;
-		  if(ind != -1) {
-		    (*boxMatrix_)[pos][offset[j] + ind] = rowIt->entry;
-		  } else {
-		    (*boxRHS_)[pos] -= 
-		      rowIt->entry * (*((*solution)[j]))[rowIt->col];
-		  }
-		}
-	      } else {
-		for(k = 0; k < numDOFs_[j]; k++) {
-		  (*boxMatrix_)[pos][offset[j] + k] = 0.0;
-		}
-	      }
-	    }
-
-	    ++pos;
-	  }
-	}
-      }
-
-//       int a,b;
-
-//       for(i = 0; i < numComponents_; i++) {
-// 	MSG("DOFs %d: ", i);
-// 	dofBegin = dofSet[i]->begin();
-// 	dofEnd = dofSet[i]->end();
-// 	for(dofIt = dofBegin; dofIt != dofEnd; ++dofIt) {
-// 	  MSG("%d ", *dofIt);
-// 	}
-// 	MSG("\n");
-//       }
-
-//       for(a = 0; a < numComponents_; a++) {
-// 	for(b = 0; b < numComponents_; b++) {
-// 	  MSG("matrix[%d][%d]:\n", a, b);
-// 	  if((*matrix)[a][b]) {
-// 	    (*matrix)[a][b]->print();
-// 	    MSG("solution[%d]\n", b);
-// 	    (*solution)[b]->print();
-// 	    MSG("rhs[%d]\n", a);
-// 	    (*rhs)[a]->print();
-// 	  }
-// 	}
-//       }
-
-//       MSG("matrix:\n");
-//       for(a = 0; a < boxSize; a++) {
-// 	for(b = 0; b < boxSize; b++) {
-// 	  MSG("%e ", (*boxMatrix_)[a][b]);
-// 	}
-// 	MSG("\n");
-//       }
-
-//       MSG("rhs:\n");
-//       for(a = 0; a < boxSize; a++) {
-// 	MSG("%e\n", (*boxRHS_)[a]);
-//       }
-
-      // apply direct solver
-      directSolver_->solve(*boxMatrix_, *boxRHS_, *boxSolution_, boxSize);
-
-//       MSG("solution:\n");
-//       for(a = 0; a < boxSize; a++) {
-// 	MSG("%e\n", (*boxSolution_)[a]);
-//       }
-
-//       WAIT_REALLY;
-
-      // copy solution
-      pos = 0;
-      for(i = 0; i < numComponents_; i++) {
-	dofBegin = dofSet[i]->begin();
-	dofEnd = dofSet[i]->end();
-	for(dofIt = dofBegin; dofIt != dofEnd; ++dofIt) {
-	  (*((*solution)[i]))[*dofIt] = (*boxSolution_)[pos];	    
-	  ++pos;
-	}
-      }
-#endif
-    };
-
-  protected:
-    Matrix<double> *boxMatrix_;
-    Vector<double> *boxSolution_;
-    Vector<double> *boxRHS_;
-    int numComponents_;
-    std::vector<int> numDOFs_;
-    int sumDOFs_;
-    DirectSolverInterface<Matrix<double>, Vector<double> > *directSolver_;
-  };
-}
-
-#endif

AMDiS/src/CreatorMap.cc

@@ -10,14 +10,8 @@
 #include "ResidualEstimator.h"
 #include "LeafData.h"
 #include "SurfaceRegion_ED.h"
-#include "MultiGridWrapper.h"
 #include "DOFMatrix.h"
-#include "SparseVector.h"
-#include "GSSmoother.h"
-#include "JacobiSmoother.h"
 #include "ElementRegion_ED.h"
-#include "BoxSmoother.h"
-#include "MultiGridPreconWrapper.h"
 #include "UmfPackSolver.h"
 
 namespace AMDiS {
@@ -64,12 +58,6 @@ namespace AMDiS {
 
     creator = NEW ICPreconditioner::Creator;
     addCreator("ic", creator);
-
-
-#if 0
-    creator =  NEW MGPreconWrapperScal::Creator;
-    addCreator("mg", creator);
-#endif
   }
 
 

AMDiS/src/DOFMatrix.cc

@@ -43,9 +43,8 @@ namespace AMDiS {
   {
     TEST_EXIT(rowFESpace)("no rowFESpace\n");
   
-    if (!colFESpace) {
+    if (!colFESpace)
       colFESpace = rowFESpace;
-    }
 
     if (rowFESpace && rowFESpace->getAdmin())
       (const_cast<DOFAdmin*>(rowFESpace->getAdmin()))->addDOFIndexed(this);
@@ -149,19 +148,22 @@ namespace AMDiS {
 
   DOFMatrix& DOFMatrix::operator=(const DOFMatrix& rhs)
   {
-    TEST_EXIT(rhs.inserter == 0 && inserter == 0)("Cannot copy during insertion");
-
     rowFESpace = rhs.rowFESpace;
     colFESpace = rhs.colFESpace;
     operators = rhs.operators;
     operatorFactor = rhs.operatorFactor;
-    matrix = rhs.matrix;
     coupleMatrix = rhs.coupleMatrix;
+
+    /// The matrix values may only be copyed, if there is no active insertion.
+    if (rhs.inserter == 0 && inserter == 0)
+      matrix = rhs.matrix;
+
     if (rhs.boundaryManager) {
       boundaryManager = NEW BoundaryManager(*rhs.boundaryManager);
     } else {
       boundaryManager = NULL;
     }
+
     if (rhs.elementMatrix) {
       elementMatrix = NEW ElementMatrix(rowFESpace->getBasisFcts()->getNumber(),
 					colFESpace->getBasisFcts()->getNumber());

AMDiS/src/DOFMatrix.h

@@ -22,10 +22,6 @@
 #ifndef AMDIS_DOFMATRIX_H
 #define AMDIS_DOFMATRIX_H
 
-// ============================================================================
-// ===== includes =============================================================
-// ============================================================================
-
 #include <vector>
 #include <set>
 #include <memory>
@@ -49,110 +45,6 @@
 
 namespace AMDiS {
 
-
-  // ===========================================================================
-  // ===== struct MatEntry =====================================================
-  // ===========================================================================
-
-  /** \brief
-   * Represents one entry of a DOFMatrix
-   */
-  struct MatEntry
-  {
-    /** \brief
-     * column index of entry (if >= 0; else unused)
-     */
-    DegreeOfFreedom col;
-
-    /** \brief
-     * matrix entry
-     */
-    double entry;
-  };
-
-
-  /** \brief
-   * Is used to search for all entries of a matrix which column index is 
-   * smaller than a given value.
-   */
-  class MatEntryValueLess : public std::unary_function<MatEntry, bool> {
-  private:
-    const double value_;
-  public:
-    MatEntryValueLess(const double& value) 
-      : value_(value) 
-    {}
-
-    bool operator()(const MatEntry& m) {
-      return (fabs(m.entry) < value_);
-    }      
-  };
-
-
-  /** \brief
-   * Is used to search for all entries of a matrix which value is smaller 
-   * than a given value.
-   */
-  class MatEntryColLess : public std::unary_function<MatEntry, bool> {
-  private:
-    const int col_;
-  public:
-    MatEntryColLess(const int& col) 
-      : col_(col) 
-    {}
-
-    bool operator()(const MatEntry& m) {
-      return (fabs(m.col) < col_);
-    }      
-  };
-
-
-  /** \brief
-   * This function is required if two MatEntries are compared by their col
-   * entry (for example when sorting a vector of MatEntries).
-   */
-  struct CmpMatEntryCol : public std::binary_function<MatEntry, MatEntry, bool> {
-    bool operator()(const MatEntry& m1, const MatEntry m2) const {
-      return m1.col < m2.col;
-    }
-  };
-  
-
-  /** \brief
-   * This function compares two matrix entries by their values. It returns true,
-   * if the value of m2 is greater than the value of m1.
-   */
-  struct CmpMatEntryValue : public std::binary_function<MatEntry, MatEntry, bool> {
-    bool operator()(const MatEntry& m1, const MatEntry m2) const {
-      return m1.entry < m2.entry;
-    }
-  };
-
-
-  /** \brief
-   * This function compares two matrix entries by their values. It returns true,
-   * if the value of m2 is greater than the value of m1.
-   */
-  struct CmpMatEntryAbsValueLess : public std::binary_function<MatEntry, MatEntry, bool> {
-    bool operator()(const MatEntry& m1, const MatEntry m2) const {
-      return fabs(m1.entry) < fabs(m2.entry);
-    }
-  };
-  
-  /** \brief
-   * This function compares two matrix entries by their values. It returns true,
-   * if the value of m1 is greater than the value of m2.
-   */
-  struct CmpMatEntryAbsValueGreater : public std::binary_function<MatEntry, MatEntry, bool> {
-    bool operator()(const MatEntry& m1, const MatEntry m2) const {
-      return fabs(m1.entry) > fabs(m2.entry);
-    }
-  };
-
-  // ============================================================================
-  // ===== class DOFMatrix ======================================================
-  // ============================================================================
-
   /** \ingroup DOFAdministration
    * \brief
    * A DOFMatrix is a sparse matrix representation for matrices that work
@@ -175,28 +67,7 @@ namespace AMDiS {
     typedef mtl::matrix::inserter<base_matrix_type>  inserter_type;
 
   private: