extensions extended by unit-vector/-matrix; FixVec,Vector and Matrix corrected...

extensions extended by unit-vector/-matrix; FixVec,Vector and Matrix corrected and extended; lots of other changes

AMDiS/CMakeLists.txt

@@ -122,6 +122,7 @@ SET(AMDIS_SRC ${SOURCE_DIR}/AdaptBase.cc
 	      ${SOURCE_DIR}/MacroElement.cc
 	      ${SOURCE_DIR}/Marker.cc
 	      ${SOURCE_DIR}/MatrixVector.cc
+# 	      ${SOURCE_DIR}/Test_MatrixVectorOperations.cc
 	      ${SOURCE_DIR}/Mesh.cc
 	      ${SOURCE_DIR}/MeshStructure.cc
 	      ${SOURCE_DIR}/Operator.cc
@@ -175,6 +176,7 @@ SET(AMDIS_SRC ${SOURCE_DIR}/AdaptBase.cc
               ${SOURCE_DIR}/io/detail/Arh2Writer.cc
 	      ${SOURCE_DIR}/io/DofWriter.cc
 	      ${SOURCE_DIR}/io/ElementFileWriter.cc
+	      ${SOURCE_DIR}/io/FileWriterInterface.cc
 	      ${SOURCE_DIR}/io/FileWriter.cc
 	      ${SOURCE_DIR}/io/GNUPlotWriter.cc
 	      ${SOURCE_DIR}/io/MacroInfo.cc
@@ -433,9 +435,7 @@ if(ENABLE_EXTENSIONS)
 	      ${EXTENSIONS_DIR}/POperators.cc
 	      ${EXTENSIONS_DIR}/SingularDirichletBC2.cc
 	      ${EXTENSIONS_DIR}/time/ExtendedRosenbrockStationary.cc
-	      ${EXTENSIONS_DIR}/pugixml/src/pugixml.cpp
-	      ${EXTENSIONS_DIR}/preconditioner/PhaseFieldCrystal_.cc
-	      ${EXTENSIONS_DIR}/preconditioner/CahnHilliard_.cc)
+	      ${EXTENSIONS_DIR}/pugixml/src/pugixml.cpp)
 	      
       if(ENABLE_SEQ_PETSC)
 	list(APPEND EXTENSIONS_SRC
@@ -510,11 +510,15 @@ if(ENABLE_EXTENSIONS)
 		  ${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
+		  ${EXTENSIONS_DIR}/base_problems/CahnHilliardNavierStokes_RB.cc
+		  ${EXTENSIONS_DIR}/base_problems/CahnHilliardNavierStokes_TwoPhase.cc
+		  ${EXTENSIONS_DIR}/base_problems/CahnHilliardNavierStokes_TwoPhase_RB.cc
+		  ${EXTENSIONS_DIR}/base_problems/DiffuseDomainFsi.cc
 		  ${EXTENSIONS_DIR}/base_problems/LinearElasticity.cc
 		  ${EXTENSIONS_DIR}/base_problems/LinearElasticityPhase.cc
-    # 	      ${EXTENSIONS_DIR}/base_problems/NavierStokes_Chorin.cc
+# 		  ${EXTENSIONS_DIR}/base_problems/NavierStokes_Chorin.cc
 		  ${EXTENSIONS_DIR}/base_problems/NavierStokesCahnHilliard.cc
+# 		  ${EXTENSIONS_DIR}/base_problems/NavierStokesPhase_Chorin.cc
 		  ${EXTENSIONS_DIR}/base_problems/NavierStokesPhase_TaylorHood.cc
 		  ${EXTENSIONS_DIR}/base_problems/NavierStokes_TaylorHood.cc
 		  ${EXTENSIONS_DIR}/base_problems/NavierStokes_TaylorHood_RB.cc
@@ -524,7 +528,10 @@ if(ENABLE_EXTENSIONS)
 		  ${EXTENSIONS_DIR}/base_problems/PhaseFieldCrystal_Phase.cc
 		  ${EXTENSIONS_DIR}/base_problems/PhaseFieldCrystal_RB.cc
 		  ${EXTENSIONS_DIR}/base_problems/PolarizationField.cc
-		  ${EXTENSIONS_DIR}/base_problems/QuasiCrystal.cc)
+		  ${EXTENSIONS_DIR}/base_problems/QuasiCrystal.cc
+		  ${EXTENSIONS_DIR}/base_problems/QuasiCrystal_RB.cc
+# 		  ${EXTENSIONS_DIR}/base_problems/VacancyPhaseFieldCrystal.cc
+		  )
 	list(APPEND COMPILEFLAGS "-DHAVE_BASE_PROBLEMS=1")
 	list(APPEND AMDIS_INCLUDE_DIRS ${EXTENSIONS_DIR}/base_problems)
 	if(WIN32)
@@ -635,6 +642,11 @@ INSTALL(FILES ${HEADERS}
 	DESTINATION include/amdis/)
 list(APPEND deb_add_dirs "include/amdis")
 
+FILE(GLOB HEADERS "${SOURCE_DIR}/config/*.h*")
+INSTALL(FILES ${HEADERS} 
+	DESTINATION include/amdis/config/)
+list(APPEND deb_add_dirs "include/amdis/config")
+
 FILE(GLOB HEADERS "${SOURCE_DIR}/*.hh")
 INSTALL(FILES ${HEADERS} 
 	DESTINATION include/amdis/)

AMDiS/src/AMDiS.h

@@ -76,6 +76,7 @@
 #include "Marker.h"
 // #include "MathFunctions.h"
 #include "MatrixVector.h"
+#include "MatrixVectorOperations.h"
 #include "Mesh.h"
 #include "MeshStructure.h"
 #include "ComponentTraverseInfo.h"

AMDiS/src/Assembler.cc

@@ -64,7 +64,7 @@ namespace AMDiS {
     Element *el = elInfo->getElement();
 
     if (el != lastMatEl || !operat->isOptimized()) {
-      initElement(elInfo);
+      initElement(elInfo, elInfo);
 
       if (rememberElMat)
 	set_to_zero(elementMatrix);
@@ -125,10 +125,12 @@ namespace AMDiS {
     ElementMatrix& mat = rememberElMat ? elementMatrix : userMat;
     
     if (secondOrderAssembler) {
+      // calculate element matrices always on smallest element
       secondOrderAssembler->calculateElementMatrix(smallElInfo, mat);
 
+      // smallElInfo stores refinement-relation to largeElInfo
       ElementMatrix &m = 
-	smallElInfo->getSubElemGradCoordsMat(rowFeSpace->getBasisFcts()->getDegree());
+	smallElInfo->getSubElemGradCoordsMat(rowFeSpace->getBasisFcts()->getDegree());  // muste be moved to next if-else block when generalized for multiple polynomial degrees
       
       if (!rowColFeSpaceEqual) {
 	if (smallElInfo == colElInfo)
@@ -321,14 +323,7 @@ namespace AMDiS {
       calculateElementMatrix(elInfo, elementMatrix);
     }
 
-    // vec += elementMatrix*uhOldLoc;
-    for (int i = 0; i < nRow; i++) {
-      double val = 0.0;
-      for (int j = 0; j < nCol; j++)
-	val += elementMatrix[i][j] * uhOldLoc[j];
-      
-      vec[i] += val;
-    }   
+    vec += elementMatrix*uhOldLoc;
   }
 
 
@@ -358,15 +353,10 @@ namespace AMDiS {
     if (mainElInfo->getElement() != lastMatEl) {
       set_to_zero(elementMatrix);
       calculateElementMatrix(mainElInfo, auxElInfo, smallElInfo, largeElInfo, 
-			     false, elementMatrix);    
+			     rowFeSpace == operat->uhOld->getFeSpace(), elementMatrix);  
     }
 
-    for (int i = 0; i < nBasFcts; i++) {
-      double val = 0.0;
-      for (int j = 0; j < nBasFcts; j++)
- 	val += elementMatrix[i][j] * uhOldLoc[j];
-      vec[i] += val;
-    }   
+    vec += elementMatrix * uhOldLoc;
   }
 
 

AMDiS/src/BallProject.h

@@ -25,6 +25,8 @@
 #ifndef AMDIS_BALLPROJECT_H
 #define AMDIS_BALLPROJECT_H
 
+#include "MatrixVectorOperations.h"
+
 namespace AMDiS {
 
   /** \brief

AMDiS/src/Cholesky.cc

@@ -20,6 +20,7 @@
 
 
 #include "Cholesky.h"
+#include "MatrixVectorOperations.h"
 
 namespace AMDiS {
 

AMDiS/src/Config.h

+#pragma once
+
+/** \brief current AMDiS version */
+#ifndef AMDIS_VERSION
+#define AMDIS_VERSION  "AMDiS: Version 0.9.1"
+#endif
+
+#include <boost/config.hpp>
+
+#define CACHE_LINE 16
+
+#if defined(__clang__)					// Clang/LLVM.
+  #include "config/Config_clang.h"
+  
+#elif defined(__ICC) || defined(__INTEL_COMPILER)	// Intel ICC/ICPC. 
+  #include "config/Config_intel.h"
+  
+#elif defined(__GNUC__) || defined(__GNUG__)		// GNU GCC/G++.
+  #include "config/Config_gcc.h"
+  
+#elif defined(__HP_cc) || defined(__HP_aCC)
+  error: not supported compiler
+  
+#elif defined(__IBMC__) || defined(__IBMCPP__)
+  error: not supported compiler
+  
+#elif defined(_MSC_VER)					// Microsoft Visual Studio. 
+  #include "config/Config_msc.h"
+  
+#elif defined(__PGI)					// Portland Group PGCC/PGCPP.
+  error: not supported compiler
+//   #include "Config_pgi.h"  
+
+#elif defined(__SUNPRO_C) || defined(__SUNPRO_CC)
+  error: not supported compiler
+#endif
+
+#include "config/Config_defaults.h"
\ No newline at end of file

AMDiS/src/CylinderProject.h

@@ -25,6 +25,8 @@
 #ifndef AMDIS_CYLINDERPROJECT_H
 #define AMDIS_CYLINDERPROJECT_H
 
+#include "MatrixVectorOperations.h"
+
 namespace AMDiS {
 
   /** \brief

AMDiS/src/DOFVector.h

@@ -350,10 +350,10 @@ namespace AMDiS {
     {}
 
     /// Constructs a DOFVector with name n belonging to FiniteElemSpace f
-    DOFVector(const FiniteElemSpace* f, std::string n, bool addToSynch = true); 
+    DOFVector(const FiniteElemSpace* f, std::string n, bool addToSynch = false); 
 
     /// Initialization.
-    void init(const FiniteElemSpace* f, std::string n, bool addToSynch = true);
+    void init(const FiniteElemSpace* f, std::string n, bool addToSynch = false);
 
     /// Copy Constructor
     DOFVector(const DOFVector& rhs) : DOFVectorBase<T>()

AMDiS/src/DOFVector.hh

@@ -130,7 +130,7 @@ namespace AMDiS {
       (this->feSpace->getAdmin())->addDOFIndexed(this);
     this->boundaryManager = new BoundaryManager(f);
 #ifdef HAVE_PARALLEL_DOMAIN_AMDIS
-    if ( Parallel::MeshDistributor::globalMeshDistributor != NULL) 
+    if (addToSynch && Parallel::MeshDistributor::globalMeshDistributor != NULL) 
       Parallel::MeshDistributor::globalMeshDistributor->addInterchangeVector(this);
 #endif
   }

AMDiS/src/DualTraverse.h

@@ -37,10 +37,10 @@ namespace AMDiS {
    */
   struct DualElInfo 
   {
-    ElInfo *rowElInfo;
-    ElInfo *colElInfo;
-    ElInfo *smallElInfo;
-    ElInfo *largeElInfo;
+    ElInfo *rowElInfo;    ///< elInfo related to testfunction
+    ElInfo *colElInfo;    ///< elInfo related to trialfunction
+    ElInfo *smallElInfo;  ///< the smaller element of (rowElInfo, colElInfo) with refinementPath relative to largeElInfo
+    ElInfo *largeElInfo;  ///< the larger element of (rowElInfo, colElInfo)
   };
 
   /// Parallel traversal of two meshes. 

AMDiS/src/ElInfo.h

@@ -252,7 +252,7 @@ namespace AMDiS {
 
     virtual mtl::dense2D<double>& getSubElemGradCoordsMat(int degree) const 
     {
-      return subElemGradMatrices[degree][std::make_pair(refinementPathLength, refinementPath)];
+      return getSubElemCoordsMat(degree);
     }
 
     /** \} */ 

AMDiS/src/ElInfo1d.cc

@@ -367,22 +367,7 @@ namespace AMDiS {
 
   mtl::dense2D<double>& ElInfo1d::getSubElemGradCoordsMat(int degree) const
   {
-    FUNCNAME("ElInfo1d::getSubElemGradCoordsMat()");
-
-    TEST_EXIT(degree == 1)("Not supported for basis functions with degree > 1!\n");
-
-    using namespace mtl;
-
-    if (subElemGradMatrices[degree].count(std::make_pair(refinementPathLength, refinementPath)) == 0) {
-      dense2D<double> mat(mat_d1);
-
-      for (int i = 0; i < refinementPathLength; i++)
-	mat *= 0.5;
-
-      subElemGradMatrices[1][std::make_pair(refinementPathLength, refinementPath)] = mat;
-    }
-
-    return subElemGradMatrices[degree][std::make_pair(refinementPathLength, refinementPath)];
+    return getSubElemCoordsMat(degree);
   }
 
 

AMDiS/src/ElInfo2d.cc

@@ -66,63 +66,63 @@ namespace AMDiS {
 
 
 
-  double ElInfo2d::mat_d3_left_val[10][10] = {{0.0, 1.0, -0.0625, 0.3125, 0.0, 0.0, 0.0625, 0.0, 0.0, -0.0625},
-					      {0.0, 0.0, -0.0625, 0.0625, 0.0, 0.0, 0.0625, 0.0, 0.0, 0.0625},
+  double ElInfo2d::mat_d3_left_val[10][10] = {{0.0,  1.0, -6.25e-02,  3.125e-01,  0.0,  0.0,  6.25e-02,  0.0,  0.0, -6.25e-02},
+					      {0.0,  0.0, -6.25e-02,  6.25e-02,  0.0,  0.0,  6.25e-02,  0.0,  0.0, 6.25e-02},
 					      {1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
-					      {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -0.25, 0.0, 0.0, -0.125},
+					      {0.0,  0.0,  0.0,  0.0,  0.0,  0.0, -2.5e-01,  0.0,  0.0, -0.125},
 					      {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.0},
 					      {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 1.0, 0.0, 0.0},
-					      {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -0.25, 0.0, 1.0, 0.375},
-					      {0.0, 0.0, 0.5625, 0.9375, 1.0, 0.0, -0.0625, 0.0, 0.0, 0.1875},
-					      {0.0, 0.0, 0.5625, -0.3125, 0.0, 0.0, -0.0625, 0.0, 0.0, 0.0},
-					      {0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.5, 0.0, 0.0, 0.75}};
+					      {0.0,  0.0,  0.0,  0.0,  0.0,  0.0, -2.5e-01,  0.0,  1.0, 0.375},
+					      {0.0,  0.0,  5.625e-01,  9.375e-01,  1.0,  0.0, -6.25e-02,  0.0,  0.0, 1.875e-01},
+					      {0.0,  0.0,  5.625e-01, -3.125e-01,  0.0,  0.0, -6.25e-02,  0.0,  0.0, -1.875e-01},
+					      {0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.5, 0.0, 0.0, 7.5e-01}};
   mtl::dense2D<double> ElInfo2d::mat_d3_left(mat_d3_left_val);
 
-  double ElInfo2d::mat_d3_right_val[10][10] = {{0.0, 0.0, -0.0625, 0.0625, 0.0, 0.0, 0.0625, 0.0, 0.0, 0.0625},
-					       {1.0, 0.0, -0.0625, 0.0625, 0.0, 0.0, 0.3125, 0.0, 0.0, -0.0625},
+  double ElInfo2d::mat_d3_right_val[10][10] = {{0.0,  0.0, -6.25e-02,  6.25e-02,  0.0,  0.0,  6.25e-02,  0.0,  0.0, 6.25e-02},
+					      {1.0,  0.0, -6.25e-02,  6.25e-02,  0.0,  0.0,  3.125e-01,  0.0,  0.0, -6.25e-02},
 					      {0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
-					       {0.0, 0.0, 0.0, -0.25, 0.0, 0.0, 0.0, 1.0, 0.0, 0.375},
+					      {0.0,  0.0,  0.0, -2.5e-01,  0.0,  0.0,  0.0,  1.0,  0.0, 0.375},
 					      {0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0},
 					      {0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
-					       {0.0, 0.0, 0.0, -0.25, 0.0, 0.0, 0.0, 0.0, 0.0, -0.125},
-					       {0.0, 0.0, 0.5625, -0.0625, 0.0, 0.0, -0.3125, 0.0, 0.0, -0.1875},
-					       {0.0, 0.0, 0.5625, -0.0625, 0.0, 1.0, 0.9375, 0.0, 0.0, 0.1875},
-					       {0.0, 0.0, 0.0, 0.5, 1.0, 0.0, 0.0, 0.0, 0.0, 0.75}};
+					      {0.0,  0.0,  0.0, -2.5e-01,  0.0,  0.0,  0.0,  0.0,  0.0, -0.125},
+					      {0.0,  0.0,  5.625e-01, -6.25e-02,  0.0,  0.0, -3.125e-01,  0.0,  0.0, -1.875e-01},
+					      {0.0,  0.0,  5.625e-01, -6.25e-02,  0.0,  1.0,  9.375e-01,  0.0,  0.0, 1.875e-01},
+					      {0.0, 0.0, 0.0, 0.5, 1.0, 0.0, 0.0, 0.0, 0.0, 7.5e-01}};
   mtl::dense2D<double> ElInfo2d::mat_d3_right(mat_d3_right_val);
 
 
 
-  double ElInfo2d::mat_d4_left_val[15][15] = {{0.0, 1.0, 0.0, 2.734375e-01, 0.0, -3.90625e-02, 2.34375e-02, 0.0, -3.90625e-02, 0.0, 0.0, 0.0, 2.34375e-02, -3.90625e-02, 0.0},
-					      {0.0, 0.0, 0.0, -3.90625e-02, 0.0, 2.34375e-02, 2.34375e-02, 0.0, -3.90625e-02, 0.0, 0.0, 0.0, -3.90625e-02, -3.90625e-02, 0.0},
+  double ElInfo2d::mat_d4_left_val[15][15] = {{0.0,  1.0,  0.0,  2.734375e-01,  0.0, -3.906250e-02,  2.343750e-02,  0.0, -3.906250e-02,  0.0,  0.0,  0.0,  2.343750e-02, -3.906250e-02, 0.0},
+					      {0.0,  0.0,  0.0, -3.906250e-02,  0.0,  2.343750e-02,  2.343750e-02,  0.0, -3.906250e-02,  0.0,  0.0,  0.0, -3.906250e-02, -3.906250e-02, 0.0},
 					      {1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
-					      {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -6.25e-02, 0.0, 1.875e-01, 0.0, 0.0, 0.0, 1.25e-01, 6.25e-02, 0.0},
-					      {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -3.75e-01, 0.0, 0.0, 0.0, -1.25e-01, 0.0, 0.0},
-					      {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 5.e-01, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
-					      {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 5.e-01, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0},
-					      {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -3.75e-01, 0.0, 1.0, 0.0, 3.75e-01, 0.0, 0.0},
-					      {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -6.25e-02, 0.0, 1.875e-01, 0.0, 0.0, 1.0, -1.25e-01, 3.125e-01, 0.0},
-					      {0.0, 0.0, 0.0, 1.09375e+00, 1.0, 4.6875e-01, -9.375e-02, 0.0, 3.125e-02, 0.0, 0.0, 0.0, -3.125e-02, 1.5625e-01, 0.0},
-					      {0.0, 0.0, 1.0, -5.46875e-01, 0.0, 7.03125e-01, 1.40625e-01, 0.0, 1.5625e-02, 0.0, 0.0, 0.0, -4.6875e-02, -2.34375e-01, 0.0},
-					      {0.0, 0.0, 0.0, 2.1875e-01, 0.0, -1.5625e-01, -9.375e-02, 0.0, 3.125e-02, 0.0, 0.0, 0.0, 9.375e-02, 1.5625e-01, 0.0},
-					      {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 5.625e-01, 0.0, -1.875e-01, 0.0, 0.0, 0.0, 3.75e-01, 9.375e-01, 1.0},
-					      {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 5.625e-01, 0.0, -1.875e-01, 0.0, 0.0, 0.0, -3.75e-01, -3.125e-01, 0.0},
+					      {0.0,  0.0,  0.0,  0.0,  0.0,  0.0, -6.25e-02,  0.0,  1.875e-01,  0.0,  0.0,  0.0,  0.125,  6.25e-02, 0.0},
+					      {0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0, -0.375,  0.0,  0.0,  0.0, -0.125,  0.0, 0.0},
+					      {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
+					      {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0},
+					      {0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0, -0.375,  0.0,  1.0,  0.0,  0.375,  0.0, 0.0},
+					      {0.0,  0.0,  0.0,  0.0,  0.0,  0.0, -6.25e-02,  0.0,  1.875e-01,  0.0,  0.0,  1.0, -0.125,  3.125e-01, 0.0},
+					      {0.0,  0.0,  0.0,  1.093750e+00,  1.0,  4.687500e-01, -9.375e-02,  0.0,  3.125e-02,  0.0,  0.0,  0.0, -3.125e-02,  1.562500e-01, 0.0},
+					      {0.0,  0.0,  1.0, -5.468750e-01,  0.0,  7.031250e-01,  1.406250e-01,  0.0,  1.562500e-02,  0.0,  0.0,  0.0, -4.687500e-02, -2.343750e-01, 0.0},
+					      {0.0,  0.0,  0.0,  2.187500e-01,  0.0, -1.562500e-01, -9.375e-02,  0.0,  3.125e-02,  0.0,  0.0,  0.0,  9.375e-02,  1.562500e-01, 0.0},
+					      {0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  5.625e-01,  0.0, -1.875e-01,  0.0,  0.0,  0.0,  0.375,  9.375e-01, 1.0},
+					      {0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  5.625e-01,  0.0, -1.875e-01,  0.0,  0.0,  0.0, -0.375, -3.125e-01, 0.0},
 					      {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 7.5e-01, 0.0, 0.0, 0.0, 7.5e-01, 0.0, 0.0}};
   mtl::dense2D<double> ElInfo2d::mat_d4_left(mat_d4_left_val);
 
-  double ElInfo2d::mat_d4_right_val[15][15] = {{0.0, 0.0, 0.0, -3.90625e-02, 0.0, 2.34375e-02, 2.34375e-02, 0.0, -3.90625e-02, 0.0, 0.0, 0.0, -3.90625e-02, -3.90625e-02, 0.0},
-					       {1.0, 0.0, 0.0, -3.90625e-02, 0.0, 2.34375e-02, -3.90625e-02, 0.0, 2.734375e-01, 0.0, 0.0, 0.0, -3.90625e-02, 2.34375e-02, 0.0},
+  double ElInfo2d::mat_d4_right_val[15][15] = {{0.0,  0.0,  0.0, -3.906250e-02,  0.0,  2.343750e-02,  2.343750e-02,  0.0, -3.906250e-02,  0.0,  0.0,  0.0, -3.906250e-02, -3.906250e-02, 0.0},
+						{1.0,  0.0,  0.0, -3.906250e-02,  0.0,  2.343750e-02, -3.906250e-02,  0.0,  2.734375e-01,  0.0,  0.0,  0.0, -3.906250e-02,  2.343750e-02, 0.0},
 						{0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
-					       {0.0, 0.0, 0.0, 1.875e-01, 0.0, -6.25e-02, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 3.125e-01, -1.25e-01, 0.0},
-					       {0.0, 0.0, 0.0, -3.75e-01, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 3.75e-01, 0.0},
-					       {0.0, 0.0, 0.0, 5.0e-01, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0},
-					       {0.0, 0.0, 0.0, 5.0e-01, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
-					       {0.0, 0.0, 0.0, -3.75e-01, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -1.25e-01, 0.0},
-					       {0.0, 0.0, 0.0, 1.875e-01, 0.0, -6.25e-02, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 6.25e-02, 1.25e-01, 0.0},
-					       {0.0, 0.0, 0.0, 3.125e-02, 0.0, -9.375e-02, -1.5625e-01, 0.0, 2.1875e-01, 0.0, 0.0, 0.0, 1.5625e-01, 9.375e-02, 0.0},
-					       {0.0, 0.0, 1.0, 1.5625e-02, 0.0, 1.40625e-01, 7.03125e-01, 0.0, -5.46875e-01, 0.0, 0.0, 0.0, -2.34375e-01, -4.6875e-02, 0.0},
-					       {0.0, 0.0, 0.0, 3.125e-02, 0.0, -9.375e-02, 4.6875e-01, 1.0, 1.09375e+00, 0.0, 0.0, 0.0, 1.5625e-01, -3.125e-02, 0.0},
-					       {0.0, 0.0, 0.0, -1.875e-01, 0.0, 5.625e-01, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -3.125e-01, -3.75e-01, 0.0},
-					       {0.0, 0.0, 0.0, -1.875e-01, 0.0, 5.625e-01, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 9.375e-01, 3.75e-01, 1.0},
+						{0.0,  0.0,  0.0,  1.875e-01,  0.0, -6.25e-02,  0.0,  0.0,  0.0,  1.0,  0.0,  0.0,  3.125e-01, -0.125, 0.0},
+						{0.0,  0.0,  0.0, -0.375,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  1.0,  0.0,  0.0,  0.375, 0.0},
+						{0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0},
+						{0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
+						{0.0,  0.0,  0.0, -0.375,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0, -0.125, 0.0},
+						{0.0,  0.0,  0.0,  1.875e-01,  0.0, -6.25e-02,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  6.25e-02,  0.125, 0.0},
+						{0.0,  0.0,  0.0,  3.125e-02,  0.0, -9.375e-02, -1.562500e-01,  0.0,  2.187500e-01,  0.0,  0.0,  0.0,  1.562500e-01,  9.375e-02, 0.0},
+						{0.0,  0.0,  1.0,  1.562500e-02,  0.0,  1.406250e-01,  7.031250e-01,  0.0, -5.468750e-01,  0.0,  0.0,  0.0, -2.343750e-01, -4.687500e-02, 0.0},
+						{0.0,  0.0,  0.0,  3.125e-02,  0.0, -9.375e-02,  4.687500e-01,  1.0,  1.093750e+00,  0.0,  0.0,  0.0,  1.562500e-01, -3.125e-02, 0.0},
+						{0.0,  0.0,  0.0, -1.875e-01,  0.0,  5.625e-01,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0, -3.125e-01, -0.375, 0.0},
+						{0.0,  0.0,  0.0, -1.875e-01,  0.0,  5.625e-01,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  9.375e-01,  0.375, 1.0},
 						{0.0, 0.0, 0.0, 7.5e-01, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 7.5e-01, 0.0}};
   mtl::dense2D<double> ElInfo2d::mat_d4_right(mat_d4_right_val);
 
@@ -873,39 +873,6 @@ namespace AMDiS {
 
   mtl::dense2D<double>& ElInfo2d::getSubElemGradCoordsMat(int degree) const
   {
-    FUNCNAME("ElInfo2d::getSubElemGradCoordsMat()");
-
-    TEST_EXIT(degree == 1)("Not supported for basis functions with degree > 1!\n");
-
-    using namespace mtl;
-
-    if (subElemGradMatrices[degree].count(std::make_pair(refinementPathLength, refinementPath)) == 0) {
-      dense2D<double> mat(3, 3), tmpMat(3, 3);
-      mat = 1;
-
-      double test_left[3][3] = {{0.0, 0.0, 0.5},
-				{-0.5, -0.5, 0.0},
-				{1.0, 0.0, 0.0}};
-      double test_right[3][3] = {{0.0, 0.0, 0.5},
-				 {0.5, -0.5, 0.0},
-				 {0.0, 1.0, 0.0}};
-      
-      mtl::dense2D<double> mat_left(test_left);
-      mtl::dense2D<double> mat_right(test_right);
-
-
-      for (int i = 0; i < refinementPathLength; i++)
-	if (refinementPath & (1 << i)) {
-	  tmpMat = mat_right * mat;
-	  mat = tmpMat;
-	} else  {
-	  tmpMat = mat_left * mat;
-	  mat = tmpMat;
-	}
-
-      subElemGradMatrices[1][std::make_pair(refinementPathLength, refinementPath)] = mat;
-    }
-
-    return subElemGradMatrices[degree][std::make_pair(refinementPathLength, refinementPath)];
+    return getSubElemCoordsMat(degree);
   }
 }

AMDiS/src/ElInfo3d.cc

@@ -824,11 +824,7 @@ namespace AMDiS {
 
   mtl::dense2D<double>& ElInfo3d::getSubElemGradCoordsMat(int degree) const
   {
-    FUNCNAME("ElInfo3d::getSubElemGradCoordsMat()");
-
-    ERROR_EXIT("Not yet implemented!\n");
-
-    return subElemGradMatrices[degree][std::make_pair(refinementPathLength, refinementPath)];
+    return getSubElemCoordsMat(degree);
   }
 
 }

AMDiS/src/Expressions.h

@@ -28,6 +28,7 @@
 #include "AMDiS_fwd.h"
 #include "OperatorTerm.h"
 #include "Functors.h"
+#include "MatrixVectorOperations.h"
 
 #include <boost/static_assert.hpp>
 #include <boost/type_traits.hpp>
@@ -94,6 +95,7 @@
 
 namespace AMDiS {
 
+/// helper class to adopt the correct OperatorTerm based on the term order
 template<int Order>
 struct GetTerm {
   typedef typename boost::mpl::if_c<Order == 0, ZeroOrderTerm, 
@@ -103,19 +105,27 @@ struct GetTerm {
 	  >::type >::type >::type type;
 };
 
-
+/// basic interface for OperatorTerms based on expressions
 template<typename Term, int Order = -1>
 struct GenericOperatorTerm : public GetTerm<Order>::type
 {
   typedef typename GetTerm<Order>::type super;
   
+  /// Expression term stored as copy
   Term term;
+  
+  /// constructor
+  /// adds all feSpaces provided by the expression term to auxFeSpaces liste
   GenericOperatorTerm(const Term& term_)
     : super(term_.getDegree()), term(term_) 
   {
     term.insertFeSpaces(this->auxFeSpaces);
+#ifndef NDEBUG
+    test_auxFeSpaces(this->auxFeSpaces);
+#endif
   }
 
+  /// calls initElement() for \ref term
   void initElement(const ElInfo* elInfo,
 		   SubAssembler* subAssembler,
 		   Quadrature *quad)
@@ -123,6 +133,7 @@ struct GenericOperatorTerm : public GetTerm<Order>::type
     term.initElement(this, elInfo, subAssembler, quad, NULL);
   }
 
+  /// calls initElement() for \ref term
   void initElement(const ElInfo* smallElInfo,
 		   const ElInfo* largeElInfo,
 		   SubAssembler* subAssembler,
@@ -130,9 +141,22 @@ struct GenericOperatorTerm : public GetTerm<Order>::type
   {
     term.initElement(this, smallElInfo, largeElInfo, subAssembler, quad, NULL);
   }
-};
-
 
+  /// test for only one mesh allowed in expressions
+  template<typename FeSpaceList>
+  void test_auxFeSpaces(FeSpaceList const& auxFeSpaces)
+  {
+    typedef typename FeSpaceList::const_iterator fe_iter;
+    if (auxFeSpaces.size() > 0) {
+      Mesh* mesh0 = (*auxFeSpaces.begin())->getMesh();
+      for (fe_iter it = auxFeSpaces.begin(); it != auxFeSpaces.end(); it++) {
+	if ((*it)->getMesh() != mesh0) {
+	  ERROR_EXIT("Only one mesh allowed in expression.\n");
+	}
+      }
+    }
+  }
+};
 
 template<typename Term>
 struct GenericOperatorTerm<Term, -1> : public GenericOperatorTerm<Term, -2>
@@ -687,7 +711,7 @@ template<typename Term>
 inline typename boost::enable_if<typename traits::is_expr<Term>::type, typename Term::value_type>::type
 max(Term term)
 {
-  typename Term::value_type value0 = -1.e25;
+  typename Term::value_type value0 = std::numeric_limits<typename Term::value_type>::min();
   value0 = accumulate(term, functors::max<typename Term::value_type>(), value0);
   
 #ifdef HAVE_PARALLEL_DOMAIN_AMDIS
@@ -701,7 +725,7 @@ template<typename Term>
 inline typename boost::enable_if<typename traits::is_expr<Term>::type, typename Term::value_type>::type
 min(Term term)
 {
-  typename Term::value_type value0 = 1.e25;
+  typename Term::value_type value0 = std::numeric_limits<typename Term::value_type>::max();
   value0 = accumulate(term, functors::min<typename Term::value_type>(), value0);
   
 #ifdef HAVE_PARALLEL_DOMAIN_AMDIS
@@ -715,7 +739,7 @@ template<typename Term>
 inline typename boost::enable_if<typename traits::is_expr<Term>::type, typename Term::value_type>::type
 abs_max(Term term)
 {
-  typename Term::value_type value0 = 0.0;
+  typename Term::value_type value0 = 0;
   value0 = accumulate(term, functors::abs_max<typename Term::value_type>(), value0);
   
 #ifdef HAVE_PARALLEL_DOMAIN_AMDIS
@@ -729,7 +753,7 @@ template<typename Term>
 inline typename boost::enable_if<typename traits::is_expr<Term>::type, typename Term::value_type>::type
 abs_min(Term term)
 {
-  typename Term::value_type value0 = 1.e25;
+  typename Term::value_type value0 = std::numeric_limits<typename Term::value_type>::max();
   value0 = accumulate(term, functors::abs_min<typename Term::value_type>(), value0);
   
 #ifdef HAVE_PARALLEL_DOMAIN_AMDIS
@@ -744,19 +768,61 @@ abs_min(Term term)
 template<typename T, typename Term>
 inline typename boost::enable_if<
   typename boost::mpl::and_<typename traits::is_expr<Term>::type, 
-			    typename boost::is_convertible<typename Term::value_type, T>::type
+			    typename traits::is_convertible<typename Term::value_type, T>::type
 			    >::type
   >::type
 transformDOF(Term term, DOFVector<T>* result);
 
+/// Assign an expression to a DOFVector (using multi-mesh if term and result vector are on different meshes)
+template<typename T, typename Term>
+inline typename boost::enable_if<
+  typename boost::mpl::and_<typename traits::is_expr<Term>::type, 
+			    typename traits::is_convertible<typename Term::value_type, T>::type
+			    >::type
+  >::type
+transformDOF_mm(Term term, DOFVector<T>* result);
+
 /// Assign an expression to a DOFVector
 template<typename T, typename Term>
 typename boost::enable_if<