Introduction of optimized dual mesh traverse procedure.

AMDiS/src/Assembler.cc

@@ -23,15 +23,17 @@ namespace AMDiS {
       rememberElVec(false),
       elementMatrix(nRow, nCol),
       elementVector(nRow),
+      tmpMat(nRow, nCol),
       lastMatEl(NULL),
       lastVecEl(NULL),
       lastTraverseId(-1)
-  
   {}
 
+
   Assembler::~Assembler()
   {}
 
+
   void Assembler::calculateElementMatrix(const ElInfo *elInfo, 
 					 ElementMatrix& userMat,
 					 double factor)
@@ -110,10 +112,9 @@ namespace AMDiS {
       ERROR_EXIT("Da muss i noch ma testen!\n");
       secondOrderAssembler->calculateElementMatrix(smallElInfo, mat);
 
-      ElementMatrix m(nRow, nRow);
-      smallElInfo->getSubElemCoordsMat_so(m, rowFESpace->getBasisFcts()->getDegree());
+      ElementMatrix &m =       
+	smallElInfo->getSubElemCoordsMat_so(rowFESpace->getBasisFcts()->getDegree());
 
-      ElementMatrix tmpMat(nRow, nRow);
       tmpMat = m * mat;
       mat = tmpMat;      
     }
@@ -131,10 +132,8 @@ namespace AMDiS {
     if (zeroOrderAssembler) {
       zeroOrderAssembler->calculateElementMatrix(smallElInfo, mat);
 
-      ElementMatrix m(nRow, nRow);
-      smallElInfo->getSubElemCoordsMat(m, rowFESpace->getBasisFcts()->getDegree());
-
-      ElementMatrix tmpMat(nRow, nRow);
+      ElementMatrix &m = 
+	smallElInfo->getSubElemCoordsMat(rowFESpace->getBasisFcts()->getDegree());
 
       if (smallElInfo == colElInfo) 
 	tmpMat = m * mat;	
@@ -366,22 +365,24 @@ namespace AMDiS {
     }
   }
 
+
   void Assembler::finishAssembling()
   {
     lastVecEl = NULL;
     lastMatEl = NULL;
   }
 
+
   OptimizedAssembler::OptimizedAssembler(Operator  *op,
 					 Quadrature *quad2,
 					 Quadrature *quad1GrdPsi,
 					 Quadrature *quad1GrdPhi,
 					 Quadrature *quad0,
-					 const FiniteElemSpace *row,
-					 const FiniteElemSpace *col) 
-    : Assembler(op, row, col)
+					 const FiniteElemSpace *rowFeSpace,
+					 const FiniteElemSpace *colFeSpace) 
+    : Assembler(op, rowFeSpace, colFeSpace)
   {
-    bool opt = (row == col);
+    bool opt = (rowFeSpace->getBasisFcts() == colFeSpace->getBasisFcts());
 
     // create sub assemblers
     secondOrderAssembler = 
@@ -396,14 +397,15 @@ namespace AMDiS {
     checkQuadratures();
   }
 
+
   StandardAssembler::StandardAssembler(Operator *op,
 				       Quadrature *quad2,
 				       Quadrature *quad1GrdPsi,
 				       Quadrature *quad1GrdPhi,
 				       Quadrature *quad0,
-				       const FiniteElemSpace *row,
-				       const FiniteElemSpace *col) 
-    : Assembler(op, row, col)
+				       const FiniteElemSpace *rowFeSpace,
+				       const FiniteElemSpace *colFeSpace) 
+    : Assembler(op, rowFeSpace, colFeSpace)
   {
     remember = false;
 

AMDiS/src/Assembler.h

@@ -247,6 +247,9 @@ namespace AMDiS {
 
     /// Locally stored element vector
     ElementVector elementVector;
+
+    ///
+    ElementMatrix tmpMat;
   
     /** \brief
      * Used to check whether \ref initElement() must be called, because

AMDiS/src/ComponentTraverseInfo.h

@@ -87,7 +87,7 @@ namespace AMDiS {
       return NULL;
     }
 
-    int getStatus() 
+    inline int getStatus() const
     {
       return status;
     }
@@ -166,16 +166,30 @@ namespace AMDiS {
       return vectorComponents[row];
     }
 
-    int getStatus(int row, int col) 
+    inline int getStatus(int row, int col) const
     {
       return matrixComponents[row][col].getStatus();
     }
 
-    int getStatus(int row) 
+    inline int getStatus(int row) const
     {
       return vectorComponents[row].getStatus();
     }
 
+    /** \brief
+     * Returns true, if for the given matrix component the row and the col FE spaces
+     * are equal. Note that this is also the case, if there is another aux FE space,
+     * that may be different from both, the row and the col FE spaces.
+     */
+    inline bool eqSpaces(int row, int col) const
+    {
+      int status = matrixComponents[row][col].getStatus();
+
+      return (status == SingleComponentInfo::EQ_SPACES_NO_AUX ||
+	      status == SingleComponentInfo::EQ_SPACES_WITH_AUX ||
+	      status == SingleComponentInfo::EQ_SPACES_WITH_DIF_AUX);
+    }
+
     const FiniteElemSpace* getAuxFeSpace(int row, int col) 
     {
       return matrixComponents[row][col].getAuxFeSpace();
@@ -185,6 +199,24 @@ namespace AMDiS {
     {
       return vectorComponents[row].getAuxFeSpace();
     }
+
+    /// Returns true if there is an aux FE that is different from row and col FE spaces.
+    inline bool difAuxSpace(int row, int col) const
+    {
+      int status = matrixComponents[row][col].getStatus();
+
+      return (status == SingleComponentInfo::EQ_SPACES_WITH_DIF_AUX ||
+	      status == SingleComponentInfo::DIF_SPACES_WITH_DIF_AUX);
+    }
+
+    /// Returns true if there is an aux FE that is different from row and col FE spaces.
+    inline bool difAuxSpace(int row) const
+    {
+      int status = vectorComponents[row].getStatus();
+
+      return (status == SingleComponentInfo::EQ_SPACES_WITH_DIF_AUX ||
+	      status == SingleComponentInfo::DIF_SPACES_WITH_DIF_AUX);
+    }
     
     /** \brief
      * Returns the row FE space for a given row number, i.e., the FE space

AMDiS/src/DOFMatrix.cc

@@ -221,6 +221,7 @@ namespace AMDiS {
 	  DegreeOfFreedom col = colIndices[j];
 	  double entry = elMat[i][j];
   
+	  //	  std::cout << "ADD at " << row << " " << col << std::endl;
 	  ins[row][col] += entry;
 	}
       }

AMDiS/src/DOFVector.cc

@@ -337,8 +337,7 @@ namespace AMDiS {
     getLocalVector(largeElInfo->getElement(), localVec);
 
     const BasisFunction *basFcts = feSpace->getBasisFcts();    
-    mtl::dense2D<double> m(nBasFcts, nBasFcts);
-    smallElInfo->getSubElemCoordsMat(m, basFcts->getDegree());
+    mtl::dense2D<double> &m = smallElInfo->getSubElemCoordsMat(basFcts->getDegree());
 
     DimVec<double> &grd1 = *grdTmp[myRank];
     int parts = Global::getGeo(PARTS, dim);
@@ -779,8 +778,7 @@ namespace AMDiS {
     double *localVec = localVectors[omp_get_thread_num()];
     getLocalVector(largeElInfo->getElement(), localVec);
 
-    mtl::dense2D<double> m(nBasFcts, nBasFcts);    
-    smallElInfo->getSubElemCoordsMat(m, basFcts->getDegree());
+    mtl::dense2D<double> &m = smallElInfo->getSubElemCoordsMat(basFcts->getDegree());
 
     for (int i = 0; i < nPoints; i++) {
       result[i] = 0.0;

AMDiS/src/ElInfo.cc

@@ -14,7 +14,7 @@
 
 namespace AMDiS {
 
-  std::map<unsigned long, mtl::dense2D<double> > ElInfo::subElemMatrices;
+  std::vector<std::map<unsigned long, mtl::dense2D<double> > > ElInfo::subElemMatrices(4);
 
   ElInfo::ElInfo(Mesh *aMesh) 
     : mesh(aMesh),

AMDiS/src/ElInfo.h

@@ -239,12 +239,15 @@ namespace AMDiS {
       return refinementPathLength;
     }
 
+    virtual mtl::dense2D<double>& getSubElemCoordsMat(int degree) const 
+    {
+      return subElemMatrices[degree][refinementPath];
+    }
 
-    virtual void getSubElemCoordsMat(mtl::dense2D<double>& mat,
-				     int basisFctsDegree) const {}
-
-    virtual void getSubElemCoordsMat_so(mtl::dense2D<double>& mat,
-					int basisFctsDegree) const {}
+    virtual mtl::dense2D<double>& getSubElemCoordsMat_so(int degree) const 
+    {
+      return subElemMatrices[degree][refinementPath];
+    }
 
     /** \} */ 
 
@@ -544,7 +547,7 @@ namespace AMDiS {
     mtl::dense2D<double> subElemCoordsMat_so;
 
   public:
-    static std::map<unsigned long, mtl::dense2D<double> > subElemMatrices;
+    static std::vector<std::map<unsigned long, mtl::dense2D<double> > > subElemMatrices;
 
     /** \brief 
      * child_vertex[el_type][child][i] = father's local vertex index of new 

AMDiS/src/ElInfo1d.cc

@@ -295,45 +295,69 @@ namespace AMDiS {
   }
 
 
-  void ElInfo1d::getSubElemCoordsMat(mtl::dense2D<double>& mat, int degree) const
+  mtl::dense2D<double>& ElInfo1d::getSubElemCoordsMat(int degree) const
   {
     FUNCNAME("ElInfo1d::getSubElemCoordsMat()");
 
-    switch (degree) {
-    case 1:
-      mat = mat_d1;
+    using namespace mtl;
+
+    if (subElemMatrices[degree].count(refinementPath) == 0) {
+      switch (degree) {
+      case 1:
+	{
+	  dense2D<double> mat(mat_d1);
+	  dense2D<double> tmpMat(num_rows(mat), num_rows(mat));
+	  
+	  for (int i = 0; i < refinementPathLength; i++) {
+	    if (refinementPath & (1 << i)) {
+	      tmpMat = mat_d1_right * mat;
+	      mat = tmpMat;
+	    } else  {
+	      tmpMat = mat_d1_left * mat;
+	      mat = tmpMat;
+	    }
+	  }
 
-      for (int i = 0; i < refinementPathLength; i++) {
-	if (refinementPath & (1 << i)) 
-	  mat *= mat_d1_left;
-	else 
-	  mat *= mat_d1_right;
-      }
+	  subElemMatrices[1][refinementPath] = mat;  
+	}
+	
+	break;	
+      case 2:
+	{
+	  dense2D<double> mat(mat_d2);
+	  dense2D<double> tmpMat(num_rows(mat), num_rows(mat));
+	  
+	  for (int i = 0; i < refinementPathLength; i++) {
+	    if (refinementPath & (1 << i)) {
+	      tmpMat = mat_d2_right * mat;
+	      mat = tmpMat;
+	    } else  {
+	      tmpMat = mat_d2_left * mat;
+	      mat = tmpMat;
+	    }
+	  }
 
-      break;
-      
-    case 2:
-      mat = mat_d2;
-
-      for (int i = 0; i < refinementPathLength; i++) {
-	if (refinementPath & (1 << i)) 
-	  mat *= mat_d2_left;
-	else 
-	  mat *= mat_d2_right;
+	  subElemMatrices[2][refinementPath] = mat;  
+	}
+	break;
+      default:
+	ERROR_EXIT("Not supported for basis function degree: %d\n", degree);
       }
-
-      break;
-
-    default:
-      ERROR_EXIT("Not supported for basis function degree: %d\n", degree);
     }
+
+    return subElemMatrices[degree][refinementPath];
   }
 
 
-  void ElInfo1d::getSubElemCoordsMat_so(mtl::dense2D<double>& mat, int degree) const
+  mtl::dense2D<double>& ElInfo1d::getSubElemCoordsMat_so(int degree) const
   {
     FUNCNAME("ElInfo1d::getSubElemCoordsMat_so()");
 
+    ERROR_EXIT("Reimplement!\n");
+
+    return subElemMatrices[degree][refinementPath];
+
+#if 0
     switch (degree) {
     case 1:
       mat = test2_val;
@@ -345,6 +369,7 @@ namespace AMDiS {
     default:
       ERROR_EXIT("Not supported for basis function degree: %d\n", degree);
     }
+#endif
   }
 
 

AMDiS/src/ElInfo1d.h

@@ -62,9 +62,9 @@ namespace AMDiS {
       return (i + 1) % 2; 
     }
 
-    void getSubElemCoordsMat(mtl::dense2D<double>& mat, int degree) const;
+    mtl::dense2D<double>& getSubElemCoordsMat(int degree) const;
 
-    void getSubElemCoordsMat_so(mtl::dense2D<double>& mat, int degree) const;
+    mtl::dense2D<double>& getSubElemCoordsMat_so(int degree) const;
 
   protected:
     static double mat_d1_val[2][2];

AMDiS/src/ElInfo2d.cc

@@ -688,46 +688,47 @@ namespace AMDiS {
   }
 
 
-  void ElInfo2d::getSubElemCoordsMat(mtl::dense2D<double>& mat, int degree) const
+  mtl::dense2D<double>& ElInfo2d::getSubElemCoordsMat(int degree) const
   {
     FUNCNAME("ElInfo2d::getSubElemCoordsMat()");
 
     using namespace mtl;
 
-    switch (degree) {
-    case 1:
-      {
-	if (subElemMatrices.count(refinementPath) > 0) {
-	  mat = subElemMatrices[refinementPath];
-	  return;
-	}
-
-	mat = mat_d1;
-	dense2D<double> tmpMat(num_rows(mat), num_rows(mat));
-	
-	for (int i = 0; i < refinementPathLength; i++) {
-	  if (refinementPath & (1 << i)) {
-	    tmpMat = mat_d1_right * mat;
-	    mat = tmpMat;
-	  } else  {
-	    tmpMat = mat_d1_left * mat;
-	    mat = tmpMat;
+    if (subElemMatrices[degree].count(refinementPath) == 0) {
+      switch (degree) {
+      case 1:
+	{
+	  dense2D<double> mat(mat_d1);
+	  dense2D<double> tmpMat(num_rows(mat), num_rows(mat));
+	  
+	  for (int i = 0; i < refinementPathLength; i++) {
+	    if (refinementPath & (1 << i)) {
+	      tmpMat = mat_d1_right * mat;
+	      mat = tmpMat;
+	    } else  {
+	      tmpMat = mat_d1_left * mat;
+	      mat = tmpMat;
+	    }
 	  }
-	}
 
-	subElemMatrices[refinementPath] = mat;
+	  subElemMatrices[1][refinementPath] = mat;  
+	}
+	break;
+      default:
+	ERROR_EXIT("Not supported for basis function degree: %d\n", degree);
       }
-      break;
-    default:
-      ERROR_EXIT("Not supported for basis function degree: %d\n", degree);
     }
+
+    return subElemMatrices[degree][refinementPath];
   }
 
 
-  void ElInfo2d::getSubElemCoordsMat_so(mtl::dense2D<double>& mat, int degree) const
+  mtl::dense2D<double>& ElInfo2d::getSubElemCoordsMat_so(int degree) const
   {
     FUNCNAME("ElInfo2d::getSubElemCoordsMat_so()");
 
     ERROR_EXIT("Not yet implemented!\n");
+
+    return subElemMatrices[degree][refinementPath];
   }
 }

AMDiS/src/ElInfo2d.h

@@ -58,9 +58,9 @@ namespace AMDiS {
     /// 2-dimensional realisation of ElInfo's getElementNormal method.
     double getElementNormal(WorldVector<double> &normal) const;
 
-    void getSubElemCoordsMat(mtl::dense2D<double>& mat, int degree) const;
+    mtl::dense2D<double>& getSubElemCoordsMat(int degree) const;
 
-    void getSubElemCoordsMat_so(mtl::dense2D<double>& mat, int degree) const;
+    mtl::dense2D<double>& getSubElemCoordsMat_so(int degree) const;
 
   protected:
     /// Temp vectors for function \ref calcGrdLambda.

AMDiS/src/ElInfo3d.cc

@@ -591,19 +591,23 @@ namespace AMDiS {
   }
 
 
-  void ElInfo3d::getSubElemCoordsMat(mtl::dense2D<double>& mat, int degree) const
+  mtl::dense2D<double>& ElInfo3d::getSubElemCoordsMat(int degree) const
   {
     FUNCNAME("ElInfo3d::getSubElemCoordsMat()");
 
     ERROR_EXIT("Not yet implemented!\n");
+
+    return subElemMatrices[degree][refinementPath];
   }
 
 
-  void ElInfo3d::getSubElemCoordsMat_so(mtl::dense2D<double>& mat, int degree) const
+  mtl::dense2D<double>& ElInfo3d::getSubElemCoordsMat_so(int degree) const
   {
     FUNCNAME("ElInfo3d::getSubElemCoordsMat_so()");
 
     ERROR_EXIT("Not yet implemented!\n");
+