VecAtQP_SOT with NULL-Pointer fct

AMDiS/src/SecondOrderTerm.cc

@@ -209,8 +209,13 @@ namespace AMDiS {
     const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
     const int nPoints = static_cast<int>(LALt.size());
 
-    for (int iq = 0; iq < nPoints; iq++)
-      l1lt(grdLambda, LALt[iq], (*f)(vecAtQPs[iq]));    
+    if (f) {
+      for (int iq = 0; iq < nPoints; iq++)
+	l1lt(grdLambda, LALt[iq], (*f)(vecAtQPs[iq]));
+    } else {
+      for (int iq = 0; iq < nPoints; iq++)
+	l1lt(grdLambda, LALt[iq], vecAtQPs[iq]);
+    }
   }
 
   void VecAtQP_SOT::eval( int nPoints,
@@ -223,11 +228,20 @@ namespace AMDiS {
     int dow = Global::getGeo(WORLD);
 
     if (num_rows(D2UhAtQP) > 0) {
-      for (int iq = 0; iq < nPoints; iq++) {
-	double resultQP = 0.0;
-	for (int i = 0; i < dow; i++)
-	  resultQP += D2UhAtQP[iq][i][i];	
-	result[iq] += fac * (*f)(vecAtQPs[iq]) * resultQP;
+      if (f) {
+	for (int iq = 0; iq < nPoints; iq++) {
+	  double resultQP = 0.0;
+	  for (int i = 0; i < dow; i++)
+	    resultQP += D2UhAtQP[iq][i][i];
+	  result[iq] += fac * (*f)(vecAtQPs[iq]) * resultQP;
+	}
+      } else {
+	for (int iq = 0; iq < nPoints; iq++) {
+	  double resultQP = 0.0;
+	  for (int i = 0; i < dow; i++)
+	    resultQP += D2UhAtQP[iq][i][i];
+	  result[iq] += fac * vecAtQPs[iq] * resultQP;
+	}
       }
     }
   }
@@ -236,9 +250,16 @@ namespace AMDiS {
 			     std::vector<WorldVector<double> > &result) 
   {
     int nPoints = grdUhAtQP.size();
-    for (int iq = 0; iq < nPoints; iq++) {
-      double factor = (*f)(vecAtQPs[iq]);  
-      axpy(factor, grdUhAtQP[iq], result[iq]);      
+    if (f) {
+      for (int iq = 0; iq < nPoints; iq++) {
+	double factor = (*f)(vecAtQPs[iq]);
+	axpy(factor, grdUhAtQP[iq], result[iq]);
+      }
+    } else {
+      for (int iq = 0; iq < nPoints; iq++) {
+	double factor = vecAtQPs[iq];
+	axpy(factor, grdUhAtQP[iq], result[iq]);
+      }
     }
   }
 
@@ -1295,8 +1316,13 @@ namespace AMDiS {
     const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
     const int nPoints = static_cast<int>(LALt.size());
 
-    for (int iq = 0; iq < nPoints; iq++)
-      lalt_kl(grdLambda, xi, xj, LALt[iq], (*f)(vecAtQPs[iq]));
+    if (f) {
+      for (int iq = 0; iq < nPoints; iq++)
+	lalt_kl(grdLambda, xi, xj, LALt[iq], (*f)(vecAtQPs[iq]));
+    } else {
+      for (int iq = 0; iq < nPoints; iq++)
+	lalt_kl(grdLambda, xi, xj, LALt[iq], vecAtQPs[iq]);
+    }
   }  
 
   void VecAtQP_IJ_SOT::eval(int nPoints,
@@ -1307,9 +1333,16 @@ namespace AMDiS {
 			    double fac) 
   {
     if (num_rows(D2UhAtQP) > 0) {
-      for (int iq = 0; iq < nPoints; iq++) {
-	double factor = (*f)(vecAtQPs[iq]);
-	result[iq] += D2UhAtQP[iq][xi][xj] * factor * fac;
+      if (f) {
+	for (int iq = 0; iq < nPoints; iq++) {
+	  double factor = (*f)(vecAtQPs[iq]);
+	  result[iq] += D2UhAtQP[iq][xi][xj] * factor * fac;
+	}
+      } else {
+	for (int iq = 0; iq < nPoints; iq++) {
+	  double factor = vecAtQPs[iq];
+	  result[iq] += D2UhAtQP[iq][xi][xj] * factor * fac;
+	}
       }
     }
   }
@@ -1318,9 +1351,16 @@ namespace AMDiS {
 				std::vector<WorldVector<double> > &result) 
   {
     int nPoints = grdUhAtQP.size();
-    for (int iq = 0; iq < nPoints; iq++) {
-      double factor = (*f)(vecAtQPs[iq]);  
-      result[iq][xi] += grdUhAtQP[iq][xj] * factor;      
+    if (f) {
+      for (int iq = 0; iq < nPoints; iq++) {
+	double factor = (*f)(vecAtQPs[iq]);
+	result[iq][xi] += grdUhAtQP[iq][xj] * factor;
+      }
+    } else {
+      for (int iq = 0; iq < nPoints; iq++) {
+	double factor = vecAtQPs[iq];
+	result[iq][xi] += grdUhAtQP[iq][xj] * factor;
+      }
     }
   }