Bugfix in usage of FOT operators in multi-mesh method.

AMDiS/src/Assembler.cc

@@ -120,9 +120,11 @@ namespace AMDiS {
     }
 
     if (firstOrderAssemblerGrdPsi) {
+      //      std::cout << "PSI!" << std::endl;
+
       firstOrderAssemblerGrdPsi->calculateElementMatrix(smallElInfo, mat);
 
-      if (largeElInfo == colElInfo) {
+      if (largeElInfo == rowElInfo) {
 	ElementMatrix &m = 
 	  smallElInfo->getSubElemGradCoordsMat(rowFESpace->getBasisFcts()->getDegree());
 
@@ -140,7 +142,7 @@ namespace AMDiS {
     if (firstOrderAssemblerGrdPhi) {
       firstOrderAssemblerGrdPhi->calculateElementMatrix(smallElInfo, mat);
 
-      if (largeElInfo == rowElInfo) {
+      if (largeElInfo == colElInfo) {
 	ElementMatrix &m = 
 	  smallElInfo->getSubElemGradCoordsMat(rowFESpace->getBasisFcts()->getDegree());
 

AMDiS/src/ElInfo1d.cc

@@ -353,7 +353,7 @@ namespace AMDiS {
 
     using namespace mtl;
 
-    if (subElemMatrices[degree].count(refinementPath) == 0) {
+    if (subElemGradMatrices[degree].count(refinementPath) == 0) {
       dense2D<double> mat(mat_d1);
 
       for (int i = 0; i < refinementPathLength; i++)

AMDiS/src/ElInfo2d.cc

@@ -731,11 +731,29 @@ namespace AMDiS {
 
     using namespace mtl;
 
-    if (subElemMatrices[degree].count(refinementPath) == 0) {
+    if (subElemGradMatrices[degree].count(refinementPath) == 0) {
       dense2D<double> mat(mat_d1);
+      dense2D<double> tmpMat(num_rows(mat), num_rows(mat));
+
+      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++)
-	mat *= 0.5;
+	if (refinementPath & (1 << i)) {
+	  tmpMat = mat_right * mat;
+	  mat = tmpMat;
+	} else  {
+	  tmpMat = mat_left * mat;
+	  mat = tmpMat;
+	}
 
       subElemGradMatrices[1][refinementPath] = mat;
     }

AMDiS/src/FirstOrderAssembler.cc

@@ -85,15 +85,15 @@ namespace AMDiS {
 	dynamic_cast<FirstOrderAssembler*>(new Stand01(op, assembler, quad));    
     } else {
        if (pwConst) {
- 	newAssembler = 
- 	  (type == GRD_PSI) ?
- 	  dynamic_cast<FirstOrderAssembler*>(new Pre10(op, assembler, quad)) :
- 	  dynamic_cast<FirstOrderAssembler*>(new Pre01(op, assembler, quad));
+	 newAssembler = 
+	   (type == GRD_PSI) ?
+	   dynamic_cast<FirstOrderAssembler*>(new Pre10(op, assembler, quad)) :
+	   dynamic_cast<FirstOrderAssembler*>(new Pre01(op, assembler, quad));
        } else {
- 	newAssembler = 
- 	  (type == GRD_PSI) ?
- 	  dynamic_cast<FirstOrderAssembler*>(new Quad10(op, assembler, quad)) :
- 	  dynamic_cast<FirstOrderAssembler*>(new Quad01(op, assembler, quad));
+	 newAssembler = 
+	   (type == GRD_PSI) ?
+	   dynamic_cast<FirstOrderAssembler*>(new Quad10(op, assembler, quad)) :
+	   dynamic_cast<FirstOrderAssembler*>(new Quad01(op, assembler, quad));
        }
     }
 

AMDiS/src/FirstOrderTerm.cc

@@ -33,6 +33,14 @@ namespace AMDiS {
     vecAtQPs = getVectorAtQPs(vec, elInfo, subAssembler, quad);
   }
 
+  void VecAtQP_FOT::initElement(const ElInfo* smallElInfo,
+				const ElInfo* largeElInfo,
+				SubAssembler* subAssembler,
+				Quadrature *quad)
+  {
+    vecAtQPs = getVectorAtQPs(vec, smallElInfo, largeElInfo, subAssembler, quad);
+  }
+
   void VecAtQP_FOT::getLb(const ElInfo *elInfo, int nPoints, 
 			  VectorOfFixVecs<DimVec<double> >& Lb) const 
   {
@@ -64,7 +72,8 @@ namespace AMDiS {
 	double factor = (*f)(vecAtQPs[iq]);
 	double resultQP = 0.0;
 	for (int i = 0; i < dow; i++)
-	  resultQP += grdUhAtQP[iq][i];	
+	  resultQP += grdUhAtQP[iq][i];
+		
 	result[iq] += fac * factor * resultQP;
       } 
     }

AMDiS/src/FirstOrderTerm.h

@@ -266,6 +266,12 @@ namespace AMDiS {
     void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
 		     Quadrature *quad = NULL);
 
+    /// Implementation of \ref OperatorTerm::initElement() for multilpe meshes.
+    void initElement(const ElInfo* smallElInfo,
+		     const ElInfo* largeElInfo,
+		     SubAssembler* subAssembler,
+		     Quadrature *quad = NULL);
+
     /// Implements FirstOrderTerm::getLb().
     void getLb(const ElInfo *elInfo, int nPoints, 
 	       VectorOfFixVecs<DimVec<double> >& Lb) const;