Bugfix in element residual estimation, thanks to Katharina.

AMDiS/src/ProblemStat.cc

@@ -860,6 +860,7 @@ namespace AMDiS {
       solver->setMultipleRhs(true);
     }
 
+
 #ifdef HAVE_PARALLEL_DOMAIN_AMDIS
     MPI::COMM_WORLD.Barrier();
     INFO(info, 8)("buildAfterCoarsen needed %.5f seconds\n", 

AMDiS/src/est/ResidualEstimator.cc

@@ -304,16 +304,15 @@ namespace AMDiS {
 	for (it = dofMat->getOperatorsBegin(), itfac = dofMat->getOperatorEstFactorBegin();
 	     it != dofMat->getOperatorsEnd();  ++it, ++itfac) {
 	  if (*itfac == NULL || **itfac != 0.0) {
-	    if (num_rows(uhQP) == 0 && (*it)->zeroOrderTerms()) {
+	    if ((*it)->zeroOrderTerms()) {
 	      uhQP.change_dim(nPoints);
 	      uh[system]->getVecAtQPs(elInfo, NULL, quadFast[system], uhQP);
 	    }
-	    if (num_rows(grdUhQp) == 0 &&
-		((*it)->firstOrderTermsGrdPsi() || (*it)->firstOrderTermsGrdPhi())) {
+	    if ((*it)->firstOrderTermsGrdPsi() || (*it)->firstOrderTermsGrdPhi()) {
 	      grdUhQp.change_dim(nPoints);
 	      uh[system]->getGrdAtQPs(elInfo, NULL, quadFast[system], grdUhQp);
 	    }
-	    if (num_rows(D2UhQp) == 0 && degree > 2 && (*it)->secondOrderTerms()) {
+	    if (degree > 2 && (*it)->secondOrderTerms()) {
 	      D2UhQp.change_dim(nPoints);
 	      uh[system]->getD2AtQPs(elInfo, NULL, quadFast[system], D2UhQp);
 	    }

AMDiS/src/est/SimpleResidualEstimator.cc

@@ -41,6 +41,10 @@ namespace AMDiS {
   {
     FUNCNAME("SimpleResidualEstimator::init()");
 
+    double kappa = 0.0;
+    Parameters::get("kappa", kappa);
+    kappa_inv = 1.0 / kappa;
+
     // === Create data structures. ===
 
     basFcts = uh[0]->getFeSpace()->getBasisFcts();
@@ -188,12 +192,12 @@ namespace AMDiS {
     for (it = dofMat->getOperatorsBegin(), itfac = dofMat->getOperatorEstFactorBegin();
 	 it != dofMat->getOperatorsEnd();  ++it, ++itfac) {
       if (*itfac == NULL || **itfac != 0.0) {
-	if (num_rows(uhQP) == 0 && (*it)->zeroOrderTerms()) {
+	if ((*it)->zeroOrderTerms()) {
 	  uhQP.change_dim(nPoints);
 	  uh[0]->getVecAtQPs(elInfo, NULL, quadFast, uhQP);
 	}
 
-	if (num_rows(D2UhQp) == 0 && degree > 2 && (*it)->secondOrderTerms()) { 
+	if (degree > 2 && (*it)->secondOrderTerms()) { 
 	  D2UhQp.change_dim(nPoints);
 	  uh[0]->getD2AtQPs(elInfo, NULL, quadFast, D2UhQp);
 	}
@@ -208,12 +212,14 @@ namespace AMDiS {
     double result = 0.0;
     for (int iq = 0; iq < nPoints; iq++)
       result += quad->getWeight(iq) * riq[iq] * riq[iq];
+
+    double alpha0 = std::min(h2, kappa_inv);
    
     if (norm == NO_NORM || norm == L2_NORM)
-      result = C0 * h2 * h2 * det * result;
+      result = C0 * alpha0 * alpha0 * det * result;
     else
-      result = C0 * h2 * det * result;
-    
+      result = C0 * alpha0 * det * result;
+
     return result;
   }
 
@@ -228,7 +234,6 @@ namespace AMDiS {
     Element *el = elInfo->getElement();
     const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
     double det = elInfo->getDet();
-    double h2 = h2_from_det(det, dim);
 
     /// === Compute jump on all faces of the current element. ===
     for (int face = 0; face < nNeighbours; face++) {
@@ -320,23 +325,29 @@ namespace AMDiS {
 	val += surfaceQuad->getWeight(iq) * (jump[iq] * jump[iq]);
 
       double d = 0.5 * (det + detNeigh);
+      double h2 = h2_from_det(d, dim);
+
+      double alpha0 = std::min(h2, kappa_inv);
    
       if (norm == NO_NORM || norm == L2_NORM)
-	val *= C1 * h2_from_det(d, dim) * d;
+	val *= C1 * alpha0 * d;
       else
 	val *= C1 * d;
-      
+
       neigh->setEstimation(neigh->getEstimation(0) + val, 0);
       result += val;
     } // for face
-
     
-    double val = fh[0]->getBoundaryManager()->boundResidual(elInfo, matrix[0], uh[0]);    
+    /*
+    double val = 
+      fh[0]->getBoundaryManager()->boundResidual(elInfo, matrix[0], uh[0]);
     if (norm == NO_NORM || norm == L2_NORM)
       val *= C1 * h2;
     else
       val *= C1;   
     result += val;
+    */
+
 
     return result;
   }

AMDiS/src/est/SimpleResidualEstimator.h

@@ -193,6 +193,8 @@ namespace AMDiS {
 
     /// Maximal number of neighbours an element may have in the used dimension.
     int nNeighbours;
+
+    double kappa_inv;
   };
 }