add error estimator

- use L2 error norm of derivatives of the free energy
- the error does not dercease with final crack (mabye its better to user energynorm)
- working with brittle-fracture-half.py

dune/fracture-phasefields/l2projection.hh

@@ -34,7 +34,7 @@ namespace L2Projection {
     using type = std::vector<std::array<std::array<double, nProjectionData>, nQuadPoint> >;
   };
 
-  template<typename LQ, typename IDX>
+  template<typename LQ, typename IDX, typename VEC>
   class L2ProjectionLop
   : public PDELab::FullVolumePattern,
     public PDELab::LocalOperatorDefaultFlags
@@ -47,8 +47,9 @@ namespace L2Projection {
       // residual assembly flags
       enum { doAlphaVolume = true };
 
-      L2ProjectionLop (int intorder, LQ & localQuantity, IDX& idx)
-        : intorder_(intorder), localQuantity_(localQuantity), idx_(idx)
+      L2ProjectionLop (int intorder, LQ & localQuantity, IDX& idx, VEC& squaredElementErrorNorm, VEC& squaredElementPrStressNorm)
+        : intorder_(intorder), localQuantity_(localQuantity), idx_(idx),
+          squaredElementErrorNorm_(squaredElementErrorNorm),squaredElementPrStressNorm_(squaredElementPrStressNorm)
       {}
 
       template<typename EG, typename LFSU_HAT, typename X, typename LFSV, typename M>
@@ -116,6 +117,9 @@ namespace L2Projection {
         // loop over quadrature points
         int qp_index=0;  //index for gauss points (hack for plotting)
         const int element_index = idx_.index(eg.entity());
+        squaredElementErrorNorm_[element_index] = 0.0;
+        squaredElementPrStressNorm_[element_index] = 0.0;
+
         for (const auto& qp : quadratureRule(geo,intorder_))
         {
           // evaluate gradient of shape functions (we assume Galerkin method lfsu=lfsv)
@@ -127,12 +131,26 @@ namespace L2Projection {
             for (size_type i=0; i<lfsu.child(n).size(); i++)
               r.accumulate(lfsv.child(n), i,  localQuantity_[element_index][qp_index][n] * phi[i] * factor);
 
+          // error estimate
+          std::vector<RangeType> projectedStress(nstress,0.0);
+          for (size_t d=0; d<nstress; d++)
+            for (size_t i=0; i<lfsu.child(0).size(); i++)
+              projectedStress[d] += x(lfsu.child(d),i)*phi[i];
+
+          for (size_t d=0; d< nstress; d++){
+            squaredElementErrorNorm_[element_index] += (projectedStress[d]-localQuantity_[element_index][qp_index][d]) *
+                                                       (projectedStress[d]-localQuantity_[element_index][qp_index][d])*factor;
+            squaredElementPrStressNorm_[element_index] += projectedStress[d]*projectedStress[d]*factor;
+           }
+
           ++qp_index;
         }
       }
 
 
     protected:
+      VEC& squaredElementPrStressNorm_;
+      VEC& squaredElementErrorNorm_;
       IDX & idx_;
       int intorder_;
       LQ & localQuantity_;

problems/brittle-fracture-half.py

@@ -14,13 +14,16 @@ parameterSet = ParameterSet()
 parameterSet.structuredGrid = 'true'
 parameterSet.lower = '0 0'
 parameterSet.upper = '1 0.5'
-parameterSet.elements = '10 5'
+parameterSet.elements = '2 1'
 
 # parameterSet['structuredGrid'] = 'false'
 parameterSet.path = '../../doc/grids/'
 parameterSet.gridFile = 'square_with_hole.msh'
 
-parameterSet.numLevels = 4
+#tolerance for adaptivity
+parameterSet.tol = 0.25
+
+parameterSet.numLevels = 1 #not active with adaptive refinement
 
 #############################################
 #  Solver parameters
@@ -78,7 +81,7 @@ parameterSet.k_local = 0.0
 parameterSet.g_c = 2.7e-3
 
 # Fracture thickness
-parameterSet.l = 0.025
+parameterSet.l = 0.03125
 
 
 

src/brittle-fracture.cc

@@ -296,6 +296,8 @@ int main (int argc, char *argv[]) try
 
   int decompostion       = parameterSet.get("decompostion",0);  //decomposed(1)/non decomposed(0=default) energy
 
+  double tol             = parameterSet.get<double>("tol");
+
   ///////////////////////////////
   //   Generate the grid
   ///////////////////////////////
@@ -476,8 +478,10 @@ int main (int argc, char *argv[]) try
   L2ProjectionGFS l2ProjectionGFS(gridView,displacementFem);
   l2ProjectionGFS.name("stress");
 
-  typedef L2Projection::L2ProjectionLop<PD,IDX> L2ProjectionLop;
-  L2ProjectionLop l2ProjectionLop(intorder,projectionDataContainer,indexSet);
+  std::vector<double> squaredElementErrorNorm(gridView.size(0),0.0);
+  std::vector<double> squaredElementPrStressNorm(gridView.size(0),0.0);
+  typedef L2Projection::L2ProjectionLop<PD,IDX,std::vector<double>> L2ProjectionLop;
+  L2ProjectionLop l2ProjectionLop(intorder,projectionDataContainer,indexSet,squaredElementErrorNorm,squaredElementPrStressNorm);
 
   // Gridoperator for L2 projection
   typedef Dune::PDELab::GridOperator<L2ProjectionGFS,
@@ -668,8 +672,6 @@ int main (int argc, char *argv[]) try
   Vector xOld;
   Vector xDirichlet(0.0);
 
-  int maxRefinmentLevel = 0;
-
   for (int i=0; i<numTimeSteps; i++)
   {
     std::cout << "  ------  time step " << i << "  ------" << std::endl;
@@ -679,8 +681,6 @@ int main (int argc, char *argv[]) try
     bool adaptiveRefinement = false;
     int refinement=0;
     do{
-      if (adaptiveRefinement){
-      }
     /////////////////////////////
     //  Assemble load vector
     /////////////////////////////
@@ -948,22 +948,31 @@ int main (int argc, char *argv[]) try
     l2ProjectionGO.residual(projectedValues,l2Res);
     l2ProjectionGO.jacobian(projectedValues,lumpedMassMatrix);
     PDELab::ISTLBackend_SEQ_ExplicitDiagonal().apply(lumpedMassMatrix,projectedValues,l2Res,1e-10);
+    // error estimate
+    l2ProjectionGO.residual(projectedValues,l2Res);
+    double squaredErrorNorm = 0.0;
+    double squaredPrStressNorm = 0.0;
+    for (int i=0; i<gridView.size(0); ++i){
+      squaredErrorNorm += squaredElementErrorNorm[i];
+      squaredPrStressNorm += squaredElementPrStressNorm[i];
+    }
+
+    double error = std::sqrt(squaredErrorNorm / squaredPrStressNorm);
 
+    std::cout << "error " << error  << std::endl;
 
     adaptiveRefinement = false;
     std::map<GridType::LocalIdSet::IdType,Dune::FieldVector<double,dim+1>> persistentContainer;
     std::map<GridType::LocalIdSet::IdType,Dune::FieldVector<double,dim+1>> persistentContainerXOld;
 
-    for (int k=0; k<numLevels-1; ++k)
+     if (error > tol)
       {
-        // refinement criterion d>0.2
         for (const auto& e : Dune::elements(gridView))
           {
-            for(int i=0; i<e.subEntities(dim); i++)
-               if (x[indexSet.subIndex(e,i,dim)][dim] > 0.2 && e.level()!=numLevels-1){
+            if (std::sqrt(squaredElementErrorNorm[indexSet.index(e)]) > tol*std::sqrt(squaredPrStressNorm/gridView.size(0)))
+              {
                   grid->mark(1,e);
                   adaptiveRefinement = true;
-                  break;
               }
           }
 
@@ -982,8 +991,6 @@ int main (int argc, char *argv[]) try
 
         for(const auto& element : Dune::elements(gridView))
           {
-            if(element.level() > maxRefinmentLevel)
-              maxRefinmentLevel = element.level();
             if(element.isNew())
               {
                 for(int j=0; j < element.subEntities(dim); j++)
@@ -1073,6 +1080,8 @@ int main (int argc, char *argv[]) try
   lumpedMassMatrix = LumpedMassMatrix(l2ProjectionGO);
   lumpedMassMatrix = 0.0;
   projectionDataContainer.resize(gridView.size(0));
+  squaredElementErrorNorm.resize(gridView.size(0));
+  squaredElementPrStressNorm.resize(gridView.size(0));
 
   // reassemble quantity matrix
   quantityMatrix = quantityAssembler.assembleEvaluationMatrix<QuantityMatrix>();