merge of branches/expr to trunk

AMDiS/AMDISConfig.cmake.in

@@ -59,12 +59,12 @@ unset(_AMDIS_LIB CACHE)
 # Boost libraries
 # ---------------
 set(AMDIS_HAS_PARALLEL_DOMAIN @ENABLE_PARALLEL_DOMAIN@)
-set(AMDIS_NEED_COMPRESSION @ENABLE_COMPRESSIN@)
+set(AMDIS_NEED_COMPRESSION @ENABLE_COMPRESSION@)
 
 set(BOOST_LIBS_REQUIRED system iostreams filesystem program_options date_time)
-if(AMDIS_NEED_COMPRESSION)
+if(WIN32 AND AMDIS_NEED_COMPRESSION)
   list(APPEND BOOST_LIBS_REQUIRED zlib bzip2)
-endif(AMDIS_NEED_COMPRESSION)
+endif(WIN32 AND AMDIS_NEED_COMPRESSION)
 
 set(BOOST_VERSION "1.42")
 if(AMDIS_HAS_PARALLEL_DOMAIN)

AMDiS/CMakeLists.txt

@@ -56,6 +56,8 @@ option(ENABLE_OUTPUT "AMDiS output printing, disable only for debugging!" true)
 
 
 mark_as_advanced(ENABLE_OUTPUT)
+mark_as_advanced(ENABLE_BDDCML)
+mark_as_advanced(USE_PETSC_DEV)
 
 
 find_package(Boost 1.44 REQUIRED)
@@ -81,6 +83,7 @@ SET(AMDIS_SRC ${SOURCE_DIR}/AdaptBase.cc
 	      ${SOURCE_DIR}/Boundary.cc
 	      ${SOURCE_DIR}/BoundaryManager.cc
 	      ${SOURCE_DIR}/BoundaryObject.cc
+	      ${SOURCE_DIR}/Bubble.cc
 	      ${SOURCE_DIR}/Cholesky.cc
 	      ${SOURCE_DIR}/CoarseningManager.cc
 	      ${SOURCE_DIR}/CoarseningManager1d.cc
@@ -125,7 +128,7 @@ SET(AMDIS_SRC ${SOURCE_DIR}/AdaptBase.cc
 	      ${SOURCE_DIR}/OperatorTerm.cc
 	      ${SOURCE_DIR}/Parametric.cc
 	      ${SOURCE_DIR}/PeriodicBC.cc
-	      ${SOURCE_DIR}/ProblemImplicit.cc
+	      #${SOURCE_DIR}/ProblemImplicit.cc
 	      ${SOURCE_DIR}/ProblemInstat.cc
 	      ${SOURCE_DIR}/ProblemInterpol.cc
 	      ${SOURCE_DIR}/ProblemStat.cc
@@ -342,7 +345,7 @@ endif(ENABLE_UMFPACK)
 
 
 if(ENABLE_HYPRE)
-  include(FindHYPRE.cmake)
+  include(HYPREConfig.cmake)
   message("have hypre: ${HAVE_HYPRE}")
   if (HAVE_HYPRE)
     if (NOT MPI_FOUND)
@@ -437,13 +440,15 @@ if(ENABLE_EXTENSIONS)
       if(ENABLE_SEQ_PETSC)
 	list(APPEND EXTENSIONS_SRC
 	      ${EXTENSIONS_DIR}/preconditioner/PetscPreconPfc.cc
+	      ${EXTENSIONS_DIR}/preconditioner/PetscPreconPfcDiag.cc
 	      ${EXTENSIONS_DIR}/preconditioner/PetscPreconCahnHilliard.cc)
       endif(ENABLE_SEQ_PETSC)
       
       if(ENABLE_PARALLEL_DOMAIN)
 	list(APPEND EXTENSIONS_SRC
 # 	      ${EXTENSIONS_DIR}/preconditioner/PetscSolverNavierStokes2.cc
-	      ${EXTENSIONS_DIR}/preconditioner/PetscSolverPfc.cc)
+	      ${EXTENSIONS_DIR}/preconditioner/PetscSolverPfc.cc
+	      ${EXTENSIONS_DIR}/preconditioner/PetscSolverPfc_diag.cc)
       endif(ENABLE_PARALLEL_DOMAIN)
       
       list(APPEND COMPILEFLAGS "-DHAVE_EXTENSIONS=1")
@@ -502,7 +507,7 @@ if(ENABLE_EXTENSIONS)
       
       if(ENABLE_BASE_PROBLEMS)
 	  SET(BASE_PROBLEMS_SRC
-# 		  ${EXTENSIONS_DIR}/base_problems/CahnHilliard.cc
+		  ${EXTENSIONS_DIR}/base_problems/CahnHilliard.cc
 		  ${EXTENSIONS_DIR}/base_problems/CahnHilliard_RB.cc
 		  ${EXTENSIONS_DIR}/base_problems/CahnHilliardNavierStokes.cc
     # 	      ${EXTENSIONS_DIR}/base_problems/DiffuseDomainFsi.cc
@@ -511,13 +516,15 @@ if(ENABLE_EXTENSIONS)
     # 	      ${EXTENSIONS_DIR}/base_problems/NavierStokes_Chorin.cc
 		  ${EXTENSIONS_DIR}/base_problems/NavierStokesCahnHilliard.cc
 		  ${EXTENSIONS_DIR}/base_problems/NavierStokesPhase_TaylorHood.cc
-# 		  ${EXTENSIONS_DIR}/base_problems/NavierStokes_TaylorHood.cc
+		  ${EXTENSIONS_DIR}/base_problems/NavierStokes_TaylorHood.cc
 		  ${EXTENSIONS_DIR}/base_problems/NavierStokes_TaylorHood_RB.cc
-# 		  ${EXTENSIONS_DIR}/base_problems/NavierStokes_TH_MultiPhase.cc
+		  ${EXTENSIONS_DIR}/base_problems/NavierStokes_TH_MultiPhase.cc
 		  ${EXTENSIONS_DIR}/base_problems/NavierStokes_TH_MultiPhase_RB.cc
 		  ${EXTENSIONS_DIR}/base_problems/PhaseFieldCrystal.cc
 		  ${EXTENSIONS_DIR}/base_problems/PhaseFieldCrystal_Phase.cc
-		  ${EXTENSIONS_DIR}/base_problems/PhaseFieldCrystal_RB.cc)
+		  ${EXTENSIONS_DIR}/base_problems/PhaseFieldCrystal_RB.cc
+		  ${EXTENSIONS_DIR}/base_problems/PolarizationField.cc
+		  ${EXTENSIONS_DIR}/base_problems/QuasiCrystal.cc)
 	list(APPEND COMPILEFLAGS "-DHAVE_BASE_PROBLEMS=1")
 	list(APPEND AMDIS_INCLUDE_DIRS ${EXTENSIONS_DIR}/base_problems)
 	if(WIN32)
@@ -647,11 +654,21 @@ INSTALL(FILES ${HEADERS}
 	DESTINATION include/amdis/est/)
 list(APPEND deb_add_dirs "include/amdis/est")
 
-FILE(GLOB HEADERS "${SOURCE_DIR}/expressions/*.h")
+FILE(GLOB HEADERS "${SOURCE_DIR}/expressions/*.h*")
 INSTALL(FILES ${HEADERS} 
 	DESTINATION include/amdis/expressions/)
 list(APPEND deb_add_dirs "include/amdis/expressions")
 
+FILE(GLOB HEADERS "${SOURCE_DIR}/operations/*.h*")
+INSTALL(FILES ${HEADERS} 
+	DESTINATION include/amdis/operations/)
+list(APPEND deb_add_dirs "include/amdis/operations")
+
+FILE(GLOB HEADERS "${SOURCE_DIR}/traits/*.h*")
+INSTALL(FILES ${HEADERS} 
+	DESTINATION include/amdis/traits/)
+list(APPEND deb_add_dirs "include/amdis/traits")
+
 FILE(GLOB HEADERS "${SOURCE_DIR}/time/*.h")
 INSTALL(FILES ${HEADERS} 
 	DESTINATION include/amdis/time/)

AMDiS/src/AMDiS.h

@@ -37,6 +37,7 @@
 #include "Boundary.h"
 #include "BoundaryCondition.h"
 #include "BoundaryManager.h"
+#include "Bubble.h"
 #include "CoarseningManager.h"
 #include "CoarseningManager1d.h"
 #include "CoarseningManager2d.h"
@@ -73,7 +74,7 @@
 #include "Line.h"
 #include "MacroElement.h"
 #include "Marker.h"
-#include "MathFunctions.h"
+// #include "MathFunctions.h"
 #include "MatrixVector.h"
 #include "Mesh.h"
 #include "MeshStructure.h"
@@ -85,7 +86,7 @@
 #include "ProblemStat.h"
 #include "ProblemInstat.h"
 #include "ProblemTimeInterface.h"
-#include "ProblemImplicit.h"
+// #include "ProblemImplicit.h"
 #include "ProblemInterpol.h"
 #include "ProblemStatBase.h"
 #include "SecondOrderTerm.h"
@@ -106,6 +107,7 @@
 #include "TimedObject.h"
 #include "TransformDOF.h"
 #include "Traverse.h"
+#include "Traits.h"
 #include "Triangle.h"
 #include "VertexVector.h"
 #include "ZeroOrderTerm.h"
@@ -153,7 +155,7 @@
 #include "parallel/ParallelProblemStat.h"
 
 #if HAVE_PARALLEL_MTL4
-#include "parallel/PMTL4Solver.h"
+// #include "parallel/PMTL4Solver.h"
 #else
 #include "parallel/PetscSolver.h"
 #include "parallel/PetscSolverNavierStokes.h"

AMDiS/src/AMDiS_fwd.h

@@ -104,7 +104,7 @@ namespace AMDiS {
     class ElementFileWriter;
     class GNUPlotWriter;
     class MacroReader;
-    class MacroWriter;
+    struct MacroWriter;
     class PngWriter;
     class PovrayWriter;
     class Spreadsheet;

AMDiS/src/AdaptInfo.h

@@ -57,6 +57,7 @@ namespace AMDiS {
 	  fac_sum(1.0),
 	  spaceTolerance(0.0),
 	  timeTolerance(0.0),
+	  timeRelativeTolerance(0.0),
 	  timeErrLow(0.0),
 	  coarsenAllowed(0),
 	  refinementAllowed(1),
@@ -65,6 +66,7 @@ namespace AMDiS {
       {
 	Parameters::get(prefix + "->tolerance", spaceTolerance);
 	Parameters::get(prefix + "->time tolerance", timeTolerance);
+	Parameters::get(prefix + "->time relative tolerance", timeRelativeTolerance);
 	Parameters::get(prefix + "->coarsen allowed", coarsenAllowed);
 	Parameters::get(prefix + "->refinement allowed", refinementAllowed);
 	Parameters::get(prefix + "->refine bisections", refineBisections);
@@ -72,6 +74,8 @@ namespace AMDiS {
 	Parameters::get(prefix + "->sum factor", fac_sum);
 	Parameters::get(prefix + "->max factor", fac_max);
 
+	if (timeTolerance == 0.0 && timeRelativeTolerance == 0.0)
+	  timeTolerance = 1.0;
 	timeErrLow = timeTolerance * 0.3;
       }
 
@@ -95,6 +99,9 @@ namespace AMDiS {
 
       /// Time tolerance.
       double timeTolerance;
+      
+      /// Relative time tolerance
+      double timeRelativeTolerance;
 
       /// Lower bound for the time error.
       double timeErrLow;
@@ -143,6 +150,7 @@ namespace AMDiS {
 	maxSolverIterations(0),
 	solverTolerance(1e-8),
 	solverResidual(0.0),
+	globalTimeTolerance(1.0),
         scalContents(size),
 	deserialized(false),
 	rosenbrockMode(false)
@@ -175,6 +183,7 @@ namespace AMDiS {
       Parameters::get(name + "->min timestep", minTimestep);
       Parameters::get(name + "->max timestep", maxTimestep);
       Parameters::get(name + "->number of timesteps", nTimesteps);
+      Parameters::get(name + "->time tolerance", globalTimeTolerance);
     }
 
     /// Resets all variables to zero (or something equivalent)
@@ -481,6 +490,12 @@ namespace AMDiS {
       return scalContents[index]->timeTolerance; 
     }  
 
+    /// Returns \ref timeRelativeTolerance.
+    inline double getTimeRelativeTolerance(int index) 
+    { 
+      return scalContents[index]->timeRelativeTolerance; 
+    }  
+
     /// Sets \ref time
     inline double setTime(double t) 
     { 
@@ -690,6 +705,18 @@ namespace AMDiS {
     {
       return solverResidual;
     }
+  
+    inline void setGlobalTimeTolerance(double tol) 
+    {
+      globalTimeTolerance = tol;
+    }
+  
+    inline double getGlobalTimeTolerance() 
+    {
+      return globalTimeTolerance;
+    }
+    
+    
 
     /// Returns true, if the adaptive procedure was deserialized from a file.
     bool isDeserialized() const 
@@ -796,6 +823,9 @@ namespace AMDiS {
 
     ///
     double solverResidual;
+    
+    /// tolerance for the overall time error
+    double globalTimeTolerance;
 
     /// Scalar adapt infos.
     std::vector<ScalContent*> scalContents;

AMDiS/src/AdaptInstationary.cc

@@ -162,12 +162,6 @@ namespace AMDiS {
     problemIteration->oneIteration(adaptInfo, FULL_ITERATION);
     problemIteration->endIteration(adaptInfo);
     adaptInfo->setLastProcessedTimestep(adaptInfo->getTimestep()); 
-    if (dbgMode) {
-	std::cout << "=== ADAPT INFO DEBUG MODE           ===\n";
-	std::cout << "=== explicitTimeStrategy() ===\n";
-	adaptInfo->printTimeErrorLowInfo();
-    }
-
   }
 
 
@@ -207,15 +201,24 @@ namespace AMDiS {
 	// === Space iterations. ===
 	do {
 	  problemIteration->beginIteration(adaptInfo);
-	      if (dbgMode) {
-		std::cout << "=== ADAPT INFO DEBUG MODE           ===\n";
-		std::cout << "=== in implicitTimeStrategy() ===\n";
-		std::cout << "=== space/time iteration "<< adaptInfo->getSpaceIteration()
-			  <<" : "<< adaptInfo->getTimestepIteration() <<" ===\n";
-		adaptInfo->printTimeErrorLowInfo();
-	      }
-
-	  if (problemIteration->oneIteration(adaptInfo, FULL_ITERATION)) {
+	      
+	  if (dbgMode) {
+	    std::cout << "=== ADAPT INFO DEBUG MODE           ===\n";
+	    std::cout << "=== in implicitTimeStrategy() ===\n";
+	    std::cout << "=== space/time iteration "<< adaptInfo->getSpaceIteration()
+	      <<" : "<< adaptInfo->getTimestepIteration() <<" ===\n";
+	    adaptInfo->printTimeErrorLowInfo();
+	  }
+	  
+	  Flag adapted = problemIteration->oneIteration(adaptInfo, FULL_ITERATION);
+#if HAVE_PARALLEL_DOMAIN_AMDIS
+          int recvAllValues = 0;
+          int isAdapted = static_cast<bool>(adapted);
+          MPI::COMM_WORLD.Allreduce(&isAdapted, &recvAllValues, 1, MPI_INT, MPI_SUM);
+          if (recvAllValues) {
+#else
+          if (adapted) {
+#endif
 	    if (!fixedTimestep && 
 		!adaptInfo->timeToleranceReached() &&
 		!(adaptInfo->getTimestep() <= adaptInfo->getMinTimestep())) {
@@ -244,7 +247,6 @@ namespace AMDiS {
 
     adaptInfo->setLastProcessedTimestep(adaptInfo->getTimestep()); 
 
-
     // After successful iteration/timestep the timestep will be changed according 
     // adaption rules for next timestep. 
     // First, check for increase of timestep

AMDiS/src/AdaptStationary.cc

@@ -99,7 +99,7 @@ namespace AMDiS {
 
 
   void AdaptStationary::initialize()
-  {  
+  {
     Parameters::get(name + "->info", info);
   }
 

AMDiS/src/BasisFunction.cc

@@ -25,6 +25,7 @@
 #include "DOFVector.h"
 #include "BasisFunction.h"
 #include "Lagrange.h"
+#include "Bubble.h"
 
 namespace AMDiS {
 

AMDiS/src/BasisFunction.h

@@ -27,6 +27,7 @@
 
 #include <string>
 #include "AMDiS_fwd.h"
+#include "CreatorInterface.h"
 #include "Global.h"
 #include "Boundary.h"
 #include "MatrixVector.h"
@@ -311,6 +312,12 @@ namespace AMDiS {
 					const ElementVector& uh,
 					WorldMatrix<double>* val) const;
 
+    /**
+    * override this method, if the base of your finite element space is not
+    * nodal
+    */
+    virtual bool isNodal() const = 0;
+
   protected:
     /// Textual description
     std::string name;     
@@ -339,6 +346,26 @@ namespace AMDiS {
     /// Vector of second derivatives
     std::vector<D2BasFctType*> *d2Phi;
   };
+  
+  /** 
+   * \brief
+   * Interface for creators of concrete BasisFunctions. 
+   */
+  class BasisFunctionCreator : public CreatorInterface<BasisFunction>
+  { 
+  public:
+    virtual ~BasisFunctionCreator() {}
+
+    /// Sets \ref problem
+    void setDim(int dim_) 
+    { 
+      dim = dim_; 
+    }
+
+  protected:
+    /// dimension of the mesh
+    int dim;
+  };
 
 }
 #include "BasisFunction.hh"

AMDiS/src/Cholesky.cc

@@ -21,6 +21,11 @@
 
 #include "Cholesky.h"
 
+namespace AMDiS {
+
+using namespace std;
+using namespace AMDiS;
+
 bool Cholesky::factorization(Matrix<double> *A, Vector<double> *p)
 {
   FUNCNAME("Cholesky::factorization()");
@@ -183,3 +188,5 @@ bool Cholesky::solve(Matrix<double> *A, Vector<WorldVector<double> > *b,
 
   return success;
 }
+
+}

AMDiS/src/Cholesky.h

@@ -27,6 +27,8 @@
 
 #include "FixVec.h"
 
+namespace AMDiS {
+
 using namespace std;
 using namespace AMDiS;
 
@@ -63,4 +65,6 @@ class Cholesky
 		    Vector<double> *p = nullptr);
 };
 
+}
+
 #endif

AMDiS/src/ComponentTraverseInfo.cc

@@ -94,9 +94,9 @@ namespace AMDiS {
     TEST_EXIT_DBG(rowFeSpace != nullptr)("No row FE space!\n");
     
     for (int i = 0; i < nComponents; i++) {
-      if (matrixComponents[row][i].getColFeSpace() != rowFeSpace)
+      if (matrixComponents[row][i].getColFeSpace() && matrixComponents[row][i].getColFeSpace() != rowFeSpace)
 	return matrixComponents[row][i].getColFeSpace();
-      if (matrixComponents[row][i].getAuxFeSpace() != rowFeSpace)
+      if (matrixComponents[row][i].getAuxFeSpace() && matrixComponents[row][i].getAuxFeSpace() != rowFeSpace)
 	return matrixComponents[row][i].getAuxFeSpace();
     }
     

AMDiS/src/CreatorMap.cc

@@ -19,6 +19,8 @@
  ******************************************************************************/
 
 
+#include "BasisFunction.h"
+#include "Bubble.h"
 #include "CreatorMap.h"
 #include "MTL4Types.h"
 #include "solver/LinearSolver.h"
@@ -28,6 +30,7 @@
 #include "solver/KrylovPreconditioner.h"
 #include "MatrixVector.h"
 #include "SystemVector.h"
+#include "Lagrange.h"
 #include "LeafData.h"
 #include "SurfaceRegion_ED.h"
 #include "ElementRegion_ED.h"
@@ -172,6 +175,25 @@ namespace AMDiS {
   }
 
 
+  template<>
+  void CreatorMap<BasisFunction>::addDefaultCreators()
+  {
+    BasisFunctionCreator *creator;
+
+    // Lagrange-functions up to order 4
+    for (int i = 1; i <= MAX_DEGREE; i++) {
+      creator = new Lagrange::Creator(i);
+      addCreator("P" + boost::lexical_cast<std::string>(i), creator);
+      addCreator("Lagrange" + boost::lexical_cast<std::string>(i), creator);
+      addCreator("CG" + boost::lexical_cast<std::string>(i), creator);
+    }
+
+    // linear Lagrange functions plus element bubble function
+    creator = new Bubble::Creator;
+    addCreator("P1+bubble", creator);
+  }
+
+
   template<>
   void CreatorMap<ElementData>::addDefaultCreators()
   {

AMDiS/src/DOFAdmin.h

@@ -251,6 +251,11 @@ namespace AMDiS {
     }
 
     int calcMemoryUsage();
+    
+    void reset()
+    {
+      init();
+    }
 
     /** \} */
 

AMDiS/src/DOFIndexed.h

@@ -91,6 +91,8 @@ namespace AMDiS {
     typedef DegreeOfFreedom	size_type;
     typedef value_type&		reference;
     typedef value_type const&	const_reference;
+    typedef value_type*		pointer;
+    typedef value_type const*	const_pointer;
 
     typedef typename std::vector<value_type>::iterator		iterator;
     typedef typename std::vector<value_type>::const_iterator	const_iterator;
@@ -116,6 +118,42 @@ namespace AMDiS {
     /// Returns container element at index i
     virtual const_reference operator[](size_type i) const = 0;
   };
+  
+  template<>
+  class DOFIndexed<bool> : public DOFIndexedBase
+  {
+  public: // typedefs
+    typedef bool 				value_type;
+    typedef DegreeOfFreedom			size_type;
+    typedef std::vector<bool>::reference	reference;
+    typedef std::vector<bool>::const_reference	const_reference;
+    typedef std::vector<bool>::pointer		pointer;
+    typedef std::vector<bool>::const_pointer	const_pointer;
+
+    typedef std::vector<bool>::iterator		iterator;
+    typedef std::vector<bool>::const_iterator	const_iterator;
+    
+  public:
+    virtual ~DOFIndexed() {}
+
+    /// Returns iterator to the begin of container
+    virtual iterator begin() = 0;
+
+    /// Returns iterator to the end of container
+    virtual iterator end() = 0;
+
+    /// Returns iterator to the begin of container
+    virtual