changed newest adaptation scheme to skeleton term proposed by diPietro and...

changed newest adaptation scheme to skeleton term proposed by diPietro and Ern. writing this and old error residuals into extra vtk files

dune/dorie/interface/adaptivity.hh

@@ -127,10 +127,10 @@ public:
 			U0 eta(p0gfs,0.0);
 
 			// compute new estimates
-			ESTLOP estlop_new(inifile,param,Impl::corrected);
-			ESTGO estgo_new(gfs,p0gfs,estlop_new,mbe);
-			U0 eta_new(p0gfs,0.0);
-			estgo_new.residual(uold,eta_new);
+			// ESTLOP estlop_new(inifile,param,Impl::corrected);
+			// ESTGO estgo_new(gfs,p0gfs,estlop_new,mbe);
+			// U0 eta_new(p0gfs,0.0);
+			// estgo_new.residual(uold,eta_new);
 
 			ESTLOP estlop_newer(inifile,param,Impl::renew);
 			ESTGO estgo_newer(gfs,p0gfs,estlop_newer,mbe);
@@ -146,12 +146,13 @@ public:
 			timer3.reset();
 
 			using Dune::PDELab::Backend::native;
-			maxeta = native(eta)[0];
+			maxeta = native(eta_newer)[0];
 			for (unsigned int i=0; i<eta.flatsize(); i++) {
 				native(eta)[i] = sqrt(native(eta)[i]); // eta contains squares
-				native(eta_new)[i] = sqrt(native(eta_new)[i]);
+				// native(eta_new)[i] = sqrt(native(eta_new)[i]);
 				native(eta_newer)[i] = sqrt(native(eta_newer)[i]);
-				if (native(eta)[i] > maxeta) maxeta = native(eta)[i];
+				if (native(eta_newer)[i] > maxeta)
+					maxeta = native(eta_newer)[i];
 			}
 			if (verbose>1)
 				std::cout << "Largest Local Error: " << maxeta << "   " << std::endl;
@@ -160,13 +161,13 @@ public:
 			Dune::VTKWriter<GV> vtkwriter(gv);
 			using ETADGF = typename Dune::PDELab::DiscreteGridFunction<AGFS,U0>;
 			ETADGF etadgf(p0gfs,eta);
-			ETADGF etadgf_new(p0gfs,eta_new);
+			// ETADGF etadgf_new(p0gfs,eta_new);
 			ETADGF etadgf_newer(p0gfs,eta_newer);
 			auto etadgf_vtk = std::make_shared<Dune::PDELab::VTKGridFunctionAdapter<ETADGF>>(etadgf,"err_est_old");
-			auto etadgf_new_vtk = std::make_shared<Dune::PDELab::VTKGridFunctionAdapter<ETADGF>>(etadgf_new,"err_est_corrected");
+			// auto etadgf_new_vtk = std::make_shared<Dune::PDELab::VTKGridFunctionAdapter<ETADGF>>(etadgf_new,"err_est_corrected");
 			auto etadgf_newer_vtk = std::make_shared<Dune::PDELab::VTKGridFunctionAdapter<ETADGF>>(etadgf_newer,"err_est_new");
 			vtkwriter.addCellData(etadgf_vtk);
-			vtkwriter.addCellData(etadgf_new_vtk);
+			// vtkwriter.addCellData(etadgf_new_vtk);
 			vtkwriter.addCellData(etadgf_newer_vtk);
 
 			vtkwriter.write("err_est" + std::to_string(count));

dune/dorie/solver/operator_error_indicator.hh

@@ -209,9 +209,9 @@ namespace Dune {
           if (impl == Impl::old) {
             sum += jump * jump * factor * factor;
           }
-          else if (impl == Impl::corrected) {
-            sum += jump * jump * factor;
-          }
+          // else if (impl == Impl::corrected) {
+          //   sum += jump * jump * factor;
+          // }
           else if (impl == Impl::renew) {
             const RF grad_jump = (satCond_s * condFactor_s * (n_F * gradu_s))
               - (satCond_n * condFactor_n * (n_F * gradu_n));
@@ -219,7 +219,7 @@ namespace Dune {
             const RF gamma = 2.0 * satCond_s * condFactor_s * satCond_n * condFactor_n / ( satCond_s * condFactor_s + satCond_n * condFactor_n );
             const RF head_jump = gamma * penalty * h_jump;
 
-            const RF total_jump = grad_jump + head_jump;
+            const RF total_jump = grad_jump/2.0 + head_jump;
 
             sum += total_jump * total_jump * factor;
           }
@@ -231,15 +231,27 @@ namespace Dune {
         }
 
         // what is this variable for?
-        DF h_T = 1.0;
+        DF h_T_s = 1.0;
+        DF h_T_n = 1.0;
         if (impl == Impl::corrected || impl == Impl::renew) {
-          h_T = std::max( diameter(ig.inside().geometry()),
-                         diameter(ig.outside().geometry()) );
+          h_T_s = diameter(ig.inside().geometry());
+          h_T_n = diameter(ig.outside().geometry());
         }
 
+        DF C_F_T_s = 1.0;
+        DF C_F_T_n = 1.0;
+        DF C_P_T = 1.0 / M_PI;
+        if (impl == Impl::corrected || impl == Impl::renew) {
+          C_F_T_s = h_T_s * ig.geometry().volume() / ig.inside().geometry().volume();
+          C_F_T_s *= C_P_T * (2.0 / dim + C_P_T);
+          C_F_T_n = h_T_n * ig.geometry().volume() / ig.outside().geometry().volume();
+          C_F_T_n *= C_P_T * (2.0 / dim + C_P_T);
+        }
+
+
         // accumulate indicator
-        r_s.accumulate(lfsv_s,0, h_T * sum );
-        r_n.accumulate(lfsv_n,0, h_T * sum );
+        r_s.accumulate(lfsv_s,0, h_T_s * C_F_T_s * sum );
+        r_n.accumulate(lfsv_n,0, h_T_n * C_F_T_n * sum );
       }
 
     private: