Fixed some bugs with dual mesh. Some more code refactoring.

AMDiS/src/Assembler.cc

@@ -80,6 +80,7 @@ namespace AMDiS {
 					 const ElInfo *colElInfo,
 					 const ElInfo *smallElInfo,
 					 const ElInfo *largeElInfo,
+					 bool rowColFeSpaceEqual,
 					 ElementMatrix& userMat,
 					 double factor)
   {
@@ -93,9 +94,9 @@ namespace AMDiS {
    
     if ((el != lastMatEl && el != lastVecEl) || !operat->isOptimized()) {
       if (smallElInfo == colElInfo)
-	initElement(smallElInfo);
+	initElement(smallElInfo);	
       else
-	initElement(smallElInfo, largeElInfo);
+	initElement(smallElInfo, largeElInfo);      
     }      
 
     if (el != lastMatEl || !operat->isOptimized()) {
@@ -118,65 +119,70 @@ namespace AMDiS {
       ElementMatrix &m = 
 	smallElInfo->getSubElemGradCoordsMat(rowFeSpace->getBasisFcts()->getDegree());
       
-      if (smallElInfo == colElInfo)
-	tmpMat = m * mat;	
-      else
-  	tmpMat = mat * trans(m);
-      
-      mat = tmpMat;
+      if (!rowColFeSpaceEqual) {
+	if (smallElInfo == colElInfo)
+	  tmpMat = m * mat;	
+	else
+	  tmpMat = mat * trans(m);
+	
+	mat = tmpMat;
+      }
     }
 
     if (firstOrderAssemblerGrdPsi) {
       firstOrderAssemblerGrdPsi->calculateElementMatrix(smallElInfo, mat);
 
-      if (largeElInfo == rowElInfo) {
-	ElementMatrix &m = 
-	  smallElInfo->getSubElemGradCoordsMat(rowFeSpace->getBasisFcts()->getDegree());
-
-	tmpMat = m * mat;
-      } else {
-	ElementMatrix &m = 
-	  smallElInfo->getSubElemCoordsMat(rowFeSpace->getBasisFcts()->getDegree());
+      if (!rowColFeSpaceEqual) {
+	if (largeElInfo == rowElInfo) {
+	  ElementMatrix &m = 
+	    smallElInfo->getSubElemGradCoordsMat(rowFeSpace->getBasisFcts()->getDegree());
+	  
+	  tmpMat = m * mat;
+	} else {
+	  ElementMatrix &m = 
+	    smallElInfo->getSubElemCoordsMat(rowFeSpace->getBasisFcts()->getDegree());
+	  
+	  tmpMat = mat * trans(m);
+	}
 	
-	tmpMat = mat * trans(m);
+	mat = tmpMat;
       }
-	
-      mat = tmpMat;
     }
 
     if (firstOrderAssemblerGrdPhi) {
       firstOrderAssemblerGrdPhi->calculateElementMatrix(smallElInfo, mat);
 
-      if (largeElInfo == colElInfo) {
-	ElementMatrix &m = 
-	  smallElInfo->getSubElemGradCoordsMat(rowFeSpace->getBasisFcts()->getDegree());
-
-	tmpMat = mat * trans(m);
-      } else {
-	ElementMatrix &m = 
-	  smallElInfo->getSubElemCoordsMat(rowFeSpace->getBasisFcts()->getDegree());
+      if (!rowColFeSpaceEqual) {
+	if (largeElInfo == colElInfo) {
+	  ElementMatrix &m = 
+	    smallElInfo->getSubElemGradCoordsMat(rowFeSpace->getBasisFcts()->getDegree());
+	  
+	  tmpMat = mat * trans(m);
+	} else {
+	  ElementMatrix &m = 
+	    smallElInfo->getSubElemCoordsMat(rowFeSpace->getBasisFcts()->getDegree());
+	  
+	  tmpMat = m * mat;	
+	}
 	
-	tmpMat = m * mat;	
+	mat = tmpMat;
       }
-
-      mat = tmpMat;
     }
 
     if (zeroOrderAssembler) {
       zeroOrderAssembler->calculateElementMatrix(smallElInfo, mat);
-      
-      ElementMatrix &m = 
-	smallElInfo->getSubElemCoordsMat(rowFeSpace->getBasisFcts()->getDegree());
-      
-      if (smallElInfo == colElInfo)
-	tmpMat = m * mat;
-      else 
-	tmpMat = mat;
 
-//       else
-//   	tmpMat = mat * trans(m);
-      
-      mat = tmpMat;
+      if (!rowColFeSpaceEqual) {
+	ElementMatrix &m = 
+	  smallElInfo->getSubElemCoordsMat(rowFeSpace->getBasisFcts()->getDegree());
+	
+	if (smallElInfo == colElInfo)
+	  tmpMat = m * mat;
+	else 
+	  tmpMat = mat * trans(m);
+	
+	mat = tmpMat;
+      }
     }
 
     if (rememberElMat && &userMat != &elementMatrix)
@@ -342,7 +348,7 @@ namespace AMDiS {
     if (mainEl != lastMatEl) {
       set_to_zero(elementMatrix);
       calculateElementMatrix(mainElInfo, auxElInfo, smallElInfo, largeElInfo, 
- 			     elementMatrix);    
+			     false, elementMatrix);    
     }
 
     for (int i = 0; i < nBasFcts; i++) {

AMDiS/src/Assembler.h

@@ -69,6 +69,7 @@ namespace AMDiS {
 				const ElInfo *colElInfo,
 				const ElInfo *smallElInfo,
 				const ElInfo *largeElInfo,
+				bool rowColFeSpaceEqual,
 				ElementMatrix& userMat,
 				double factor = 1.0);
 

AMDiS/src/DOFMatrix.cc

@@ -286,13 +286,16 @@ namespace AMDiS {
 
     if (op) {
       op->getElementMatrix(rowElInfo, colElInfo, smallElInfo, largeElInfo, 
-			   elementMatrix);
+			   false, elementMatrix);
     } else {
       std::vector<Operator*>::iterator it = operators.begin();
       std::vector<double*>::iterator factorIt = operatorFactor.begin();
       for (; it != operators.end(); ++it, ++factorIt)
-	(*it)->getElementMatrix(rowElInfo, colElInfo,
-				smallElInfo, largeElInfo,
+	(*it)->getElementMatrix(rowElInfo, 
+				colElInfo,
+				smallElInfo, 
+				largeElInfo,
+				false,
 				elementMatrix, 
 				*factorIt ? **factorIt : 1.0);	     
     }
@@ -328,8 +331,11 @@ namespace AMDiS {
 	   it != operators.end(); 
 	   ++it, ++factorIt) {
 	if ((*it)->getNeedDualTraverse()) {
-	  (*it)->getElementMatrix(mainElInfo, auxElInfo,
-				  smallElInfo, largeElInfo,
+	  (*it)->getElementMatrix(mainElInfo, 
+				  auxElInfo,
+				  smallElInfo, 
+				  largeElInfo,
+				  rowFeSpace == colFeSpace,
 				  elementMatrix, 
 				  *factorIt ? **factorIt : 1.0);
 	}

AMDiS/src/DOFVector.cc

@@ -776,7 +776,8 @@ namespace AMDiS {
 	double val = 0.0;
 	for (int k = 0; k < nBasFcts; k++)
 	  val += m[j][k] * (*(basFcts->getPhi(k)))(quad->getLambda(i));
-	vecAtQPs[i] += val;
+
+	vecAtQPs[i] += localVec[j] * val;
       }
     }
   }

AMDiS/src/ElInfo2d.cc

@@ -108,19 +108,11 @@ namespace AMDiS {
 
   ElInfo2d::ElInfo2d(Mesh *aMesh) 
     : ElInfo(aMesh) 
-  {
-    e1 = new WorldVector<double>;
-    e2 = new WorldVector<double>;
-    normal = new WorldVector<double>;
-  }
+  {}
 
 
   ElInfo2d::~ElInfo2d()
-  {
-    delete e1;
-    delete e2;
-    delete normal;
-  }
+  {}
 
 
   void ElInfo2d::fillMacroInfo(const MacroElement * mel)
@@ -585,7 +577,7 @@ namespace AMDiS {
     }
   }
 
-  double ElInfo2d::calcGrdLambda(DimVec<WorldVector<double> >& grd_lam)
+  double ElInfo2d::calcGrdLambda(DimVec<WorldVector<double> >& grd)
   {
     FUNCNAME("ElInfo2d::calcGrdLambda()");
 
@@ -595,31 +587,30 @@ namespace AMDiS {
     int dim = mesh->getDim();
 
     for (int i = 0; i < dimOfWorld; i++) {
-      (*e1)[i] = coord[1][i] - coord[0][i];
-      (*e2)[i] = coord[2][i] - coord[0][i];
+      e1[i] = coord[1][i] - coord[0][i];
+      e2[i] = coord[2][i] - coord[0][i];
     }
 
     if (dimOfWorld == 2) {
-      double sdet = (*e1)[0] * (*e2)[1] - (*e1)[1] * (*e2)[0];
+      double sdet = e1[0] * e2[1] - e1[1] * e2[0];
       adet = abs(sdet);
 
       if (adet < 1.0E-25) {
 	MSG("abs(det) = %f\n", adet);
 	for (int i = 0; i <= dim; i++)
-	  for (int j = 0; j < dimOfWorld; j++)
-	    grd_lam[i][j] = 0.0;
+	  grd[i].set(0.0);
       } else {
 	double det1 = 1.0 / sdet;
 
-	grd_lam[1][0] = (*e2)[1] * det1;  // a11: (a_ij) = A^{-T}
-	grd_lam[1][1] = -(*e2)[0] * det1; // a21
-	grd_lam[2][0] = -(*e1)[1] * det1; // a12
-	grd_lam[2][1] = (*e1)[0] * det1;  // a22
-	grd_lam[0][0] = - grd_lam[1][0] - grd_lam[2][0];
-	grd_lam[0][1] = - grd_lam[1][1] - grd_lam[2][1];
+	grd[1][0] = e2[1] * det1;  // a11: (a_ij) = A^{-T}
+	grd[1][1] = -e2[0] * det1; // a21
+	grd[2][0] = -e1[1] * det1; // a12
+	grd[2][1] = e1[0] * det1;  // a22
+	grd[0][0] = -grd[1][0] - grd[2][0];
+	grd[0][1] = -grd[1][1] - grd[2][1];
       }
     } else {  
-      vectorProduct(*e1, *e2, *normal);
+      vectorProduct(e1, e2, normal);
 
       adet = norm(normal);
 
@@ -627,21 +618,21 @@ namespace AMDiS {
 	MSG("abs(det) = %lf\n", adet);
 	for (int i = 0; i <= dim; i++)
 	  for (int j = 0; j < dimOfWorld; j++)
-	    grd_lam[i][j] = 0.0;
+	    grd[i][j] = 0.0;
       } else {
-	vectorProduct(*e2, *normal, grd_lam[1]);
-	vectorProduct(*normal, *e1, grd_lam[2]);
+	vectorProduct(e2, normal, grd[1]);
+	vectorProduct(normal, e1, grd[2]);
       
 	double adet2 = 1.0 / (adet * adet);
 
 	for (int i = 0; i < dimOfWorld; i++) {
-	  grd_lam[1][i] *= adet2;
-	  grd_lam[2][i] *= adet2;
+	  grd[1][i] *= adet2;
+	  grd[2][i] *= adet2;
 	}
 
-	grd_lam[0][0] = - grd_lam[1][0] - grd_lam[2][0];
-	grd_lam[0][1] = - grd_lam[1][1] - grd_lam[2][1];
-	grd_lam[0][2] = - grd_lam[1][2] - grd_lam[2][2];
+	grd[0][0] = -grd[1][0] - grd[2][0];
+	grd[0][1] = -grd[1][1] - grd[2][1];
+	grd[0][2] = -grd[1][2] - grd[2][2];
       }
     }
 
@@ -710,7 +701,7 @@ namespace AMDiS {
       normal[0] = coord[i1][1] - coord[i0][1];
       normal[1] = coord[i0][0] - coord[i1][0];
     } else { // dow == 3
-      WorldVector<double> e0, e1,e2, elementNormal;
+      WorldVector<double> e0, elementNormal;
 
       e0 = coord[i1]; 
       e0 -= coord[i0];

AMDiS/src/ElInfo2d.h

@@ -64,7 +64,7 @@ namespace AMDiS {
 
   protected:
     /// Temp vectors for function \ref calcGrdLambda.
-    WorldVector<double> *e1, *e2, *normal;
+    WorldVector<double> e1, e2, normal;
 
     static double mat_d1_left_val[3][3];
     static mtl::dense2D<double> mat_d1_left;

AMDiS/src/FirstOrderTerm.cc

@@ -44,17 +44,17 @@ namespace AMDiS {
   void VecAtQP_FOT::getLb(const ElInfo *elInfo, int nPoints, 
 			  VectorOfFixVecs<DimVec<double> >& Lb) const 
   {
-    const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
+    const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
 
     if (bOne > -1) {
       for (int iq = 0; iq < nPoints; iq++)
-	lb_one(Lambda, Lb[iq], (*f)(vecAtQPs[iq]));
+	lb_one(grdLambda, Lb[iq], (*f)(vecAtQPs[iq]));
     } else if (b) {
       for (int iq = 0; iq < nPoints; iq++)
-	lb(Lambda, *b, Lb[iq], (*f)(vecAtQPs[iq]));
+	lb(grdLambda, *b, Lb[iq], (*f)(vecAtQPs[iq]));
     } else {
       for (int iq = 0; iq < nPoints; iq++)
-	l1(Lambda, Lb[iq], (*f)(vecAtQPs[iq]));      
+	l1(grdLambda, Lb[iq], (*f)(vecAtQPs[iq]));      
     }
   }
 
@@ -92,9 +92,9 @@ namespace AMDiS {
   void CoordsAtQP_FOT::getLb(const ElInfo *elInfo, int nPoints, 
 			     VectorOfFixVecs<DimVec<double> >& Lb) const 
   {
-    const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
+    const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
     for (int iq = 0; iq < nPoints; iq++)
-      l1(Lambda, Lb[iq], (*g)(coordsAtQPs[iq]));
+      l1(grdLambda, Lb[iq], (*g)(coordsAtQPs[iq]));
   }
 
   void CoordsAtQP_FOT::eval(int nPoints,
@@ -130,9 +130,9 @@ namespace AMDiS {
   void VecCoordsAtQP_FOT::getLb(const ElInfo *elInfo, int nPoints,
 				VectorOfFixVecs<DimVec<double> >& Lb) const
   {
-    const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
+    const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
     for (int iq = 0; iq < nPoints; iq++)
-      lb(Lambda, *b, Lb[iq], (*g)(coordsAtQPs[iq]));
+      lb(grdLambda, *b, Lb[iq], (*g)(coordsAtQPs[iq]));
   }
 
   void VecCoordsAtQP_FOT::eval(int nPoints,
@@ -178,13 +178,13 @@ namespace AMDiS {
   void VectorGradient_FOT::getLb(const ElInfo *elInfo, int nPoints, 
 				 VectorOfFixVecs<DimVec<double> >& Lb) const 
   {
-    const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
+    const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
     if (f) {
       for (int iq = 0; iq < nPoints; iq++)
-	lb(Lambda, (*f)(gradAtQPs[iq]), Lb[iq], 1.0);      
+	lb(grdLambda, (*f)(gradAtQPs[iq]), Lb[iq], 1.0);      
     } else {
       for (int iq = 0; iq < nPoints; iq++)
-	lb(Lambda, gradAtQPs[iq], Lb[iq], 1.0);      
+	lb(grdLambda, gradAtQPs[iq], Lb[iq], 1.0);      
     }
   }
 
@@ -230,9 +230,9 @@ namespace AMDiS {
   void VectorFct_FOT::getLb(const ElInfo *elInfo, int nPoints, 
 			    VectorOfFixVecs<DimVec<double> >& Lb) const 
   {
-    const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
+    const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
     for (int iq = 0; iq < nPoints; iq++)
-      lb(Lambda, (*vecFct)(vecAtQPs[iq]), Lb[iq], 1.0);    
+      lb(grdLambda, (*vecFct)(vecAtQPs[iq]), Lb[iq], 1.0);    
   }
 
   void VectorFct_FOT::eval(int nPoints,
@@ -264,9 +264,9 @@ namespace AMDiS {
   void VecFctAtQP_FOT::getLb(const ElInfo *elInfo, int nPoints, 
 			     VectorOfFixVecs<DimVec<double> >& Lb) const 
   {
-    const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
+    const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
     for (int iq = 0; iq < nPoints; iq++)
-      lb(Lambda, (*g)(coordsAtQPs[iq]), Lb[iq], 1.0);    
+      lb(grdLambda, (*g)(coordsAtQPs[iq]), Lb[iq], 1.0);    
   }
 
   void VecFctAtQP_FOT::eval(int nPoints,
@@ -322,9 +322,9 @@ namespace AMDiS {
 
   void VecGrad_FOT::getLb(const ElInfo *elInfo, int nPoints, 
 			  VectorOfFixVecs<DimVec<double> >& Lb) const {
-    const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
+    const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
     for (int iq = 0; iq < nPoints; iq++)
-      lb(Lambda, (*vecFct)(vecAtQPs[iq], gradAtQPs[iq]), Lb[iq], 1.0);
+      lb(grdLambda, (*vecFct)(vecAtQPs[iq], gradAtQPs[iq]), Lb[iq], 1.0);
   }
   
   void VecGrad_FOT::eval(int nPoints,
@@ -356,9 +356,9 @@ namespace AMDiS {
   void FctGrad2_FOT::getLb(const ElInfo *elInfo, int nPoints, 
 			   VectorOfFixVecs<DimVec<double> >& Lb) const 
   {
-    const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
+    const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
     for (int iq = 0; iq < nPoints; iq++)
-      lb(Lambda, (*vecFct)(grad1AtQPs[iq], grad2AtQPs[iq]), Lb[iq], 1.0);
+      lb(grdLambda, (*vecFct)(grad1AtQPs[iq], grad2AtQPs[iq]), Lb[iq], 1.0);
   }
   
   void FctGrad2_FOT::eval(int nPoints,
@@ -402,12 +402,12 @@ namespace AMDiS {
   }
     
   void Vec2AndGradAtQP_FOT::getLb(const ElInfo *elInfo, int nPoints, VectorOfFixVecs<DimVec<double> >& Lb) const {
-    const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
+    const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
     for (int iq = 0; iq < nPoints; iq++) {
       if (b)
-	lb(Lambda, *b, Lb[iq], (*f)(vec1AtQPs[iq], vec2AtQPs[iq],gradAtQPs1[iq]));
+	lb(grdLambda, *b, Lb[iq], (*f)(vec1AtQPs[iq], vec2AtQPs[iq],gradAtQPs1[iq]));
       else
-	l1(Lambda, Lb[iq], (*f)(vec1AtQPs[iq], vec2AtQPs[iq], gradAtQPs1[iq]));
+	l1(grdLambda, Lb[iq], (*f)(vec1AtQPs[iq], vec2AtQPs[iq], gradAtQPs1[iq]));
     }
   }
   
@@ -448,12 +448,12 @@ namespace AMDiS {
   void FctVecAtQP_FOT::getLb(const ElInfo *elInfo, int nPoints, 
 			     VectorOfFixVecs<DimVec<double> >& Lb) const