some baseproblems extended

extensions/POperators.cc

@@ -52,7 +52,7 @@ void Phase_SOT::getLALt(const ElInfo *elInfo,
   const int nPoints = static_cast<int>(LALt.size());
   
   for (int iq = 0; iq < nPoints; iq++)
-    l1lt(Lambda, LALt[iq], f(iq) * phase[iq] * fac);
+    l1lt(Lambda, LALt[iq], f(iq) * fac);
 }
 
 void Phase_SOT::eval(int nPoints,
@@ -64,7 +64,7 @@ void Phase_SOT::eval(int nPoints,
 {	
   if (num_rows(D2UhAtQP) > 0) {
     for (int iq = 0; iq < nPoints; iq++) {
-      double feval = f(iq) * phase[iq] * opFactor * fac;
+      double feval = f(iq) * opFactor * fac;
       double resultQP = 0.0;
       for (int i = 0; i < dimOfWorld; i++)
 	resultQP += D2UhAtQP[iq][i][i];	
@@ -78,14 +78,14 @@ void Phase_SOT::weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
 {
   int nPoints = grdUhAtQP.size();
   for (int iq = 0; iq < nPoints; iq++) {
-    double factor = f(iq) * phase[iq] * fac;
+    double factor = f(iq) * fac;
     axpy(factor, grdUhAtQP[iq], result[iq]);
   }
 }
 
 double Phase_SOT::f(const int iq) const
 {
-  return 1.0;
+  return std::max(0.0, std::min( 1.0, phase[iq] ) );
 }
 
 /* ----------------------------------------------------------- */
@@ -343,7 +343,6 @@ WorldVec_FOT::WorldVec_FOT(WorldVector<DOFVector<double>*> vecs_, double fac_)
   fac(fac_)
 {
   numVecs=vecs_.size();
-  TEST_EXIT(numVecs==2 || numVecs==3)("Only Dim=2 or Dim=3 possible\n");
   
   for (int i = 0; i < numVecs; i++) {
     TEST_EXIT(vecs_[i])("One vector is NULL!\n");	
@@ -351,7 +350,7 @@ WorldVec_FOT::WorldVec_FOT(WorldVector<DOFVector<double>*> vecs_, double fac_)
   }
   
   vec0DV=vecs_[0];
-  vec1DV=vecs_[1];
+  if(numVecs>=2) vec1DV=vecs_[1];
   if(numVecs>=3) vec2DV=vecs_[2];
 }
 
@@ -360,7 +359,7 @@ void WorldVec_FOT::initElement(const ElInfo* elInfo,
 					Quadrature *quad) 
 {
   getVectorAtQPs(vec0DV, elInfo, subAssembler, quad, vec0);
-  getVectorAtQPs(vec1DV, elInfo, subAssembler, quad, vec1);
+  if(numVecs>=2) getVectorAtQPs(vec1DV, elInfo, subAssembler, quad, vec1);
   if(numVecs>=3) getVectorAtQPs(vec2DV, elInfo, subAssembler, quad, vec2);
 }
 
@@ -369,7 +368,7 @@ void WorldVec_FOT::initElement(const ElInfo* largeElInfo, const ElInfo* smallElI
 					Quadrature *quad) 
 {
   getVectorAtQPs(vec0DV, smallElInfo, largeElInfo, subAssembler, quad, vec0);
-  getVectorAtQPs(vec1DV, smallElInfo, largeElInfo, subAssembler, quad, vec1);
+  if(numVecs>=2) getVectorAtQPs(vec1DV, smallElInfo, largeElInfo, subAssembler, quad, vec1);
   if(numVecs>=3) getVectorAtQPs(vec2DV, smallElInfo, largeElInfo, subAssembler, quad, vec2);
 }
 
@@ -382,7 +381,7 @@ void WorldVec_FOT::getLb(const ElInfo *elInfo,
   for (int iq = 0; iq < nPoints; iq++) {
     WorldVector<double> vec;
     vec[0]=vec0[iq];
-    vec[1]=vec1[iq];
+    if(numVecs>=2) vec[1]=vec1[iq];
     if(numVecs>=3) vec[2]=vec2[iq];
     lb(Lambda, vec, result[iq], fac);
   }
@@ -400,7 +399,7 @@ void WorldVec_FOT::eval(int nPoints,
     double resultQP = 0.0;
 		    
     resultQP += grdUhAtQP[iq][0] * vec0[iq];
-    resultQP += grdUhAtQP[iq][1] * vec1[iq];
+    if(numVecs>=2) resultQP += grdUhAtQP[iq][1] * vec1[iq];
     if(numVecs>=3) resultQP += grdUhAtQP[iq][2] * vec2[iq];
     
     result[iq] += factor * resultQP;