RecoveryEstimator.cc

#include "RecoveryEstimator.h"
#include "Parameters.h"

RecoveryEstimator::RecoveryEstimator(std::string name, DOFVector<double> *uh, int r) 
  : Estimator(name, r), 
    uh_(uh), 
    relative_(0), 
    C(1.0), 
    method(0),
    feSpace(NULL), 
    f_vec(NULL), 
    f_scal(NULL), 
    aux_vec(NULL), 
    rec_struct(NULL)
{
  FUNCNAME("RecoveryEstimator::constructor()");

  GET_PARAMETER(0, name + "->rec method", "%d", &method);
  GET_PARAMETER(0, name + "->rel error", "%d", &relative_);

  GET_PARAMETER(0, name + "->C", "%f", &C);
  C = C > 1e-25 ? sqr(C) : 1.0;

  if (norm == H1_NORM) {
    feSpace = uh->getFESpace();
    degree  = feSpace->getBasisFcts()->getDegree();

    if ((degree <= 2) && (C != 1.0)) {
      WARNING("Recovery estimators in the H1_NORM usually very good for linear and quadratic finite element; normally you do not need to overwrite the default value of C\n");
      WAIT;
    }
  }
  else {
    degree  = uh->getFESpace()->getBasisFcts()->getDegree() + 1;
    
    feSpace = FiniteElemSpace::provideFESpace(NULL,
					      Lagrange::getLagrange(uh->getFESpace()->getMesh()->getDim(), degree),
					      uh->getFESpace()->getMesh(),
					      name + "->feSpace");


    if (method == 2) {
      ERROR("Simple averaging only for the H1_NORM; using SPR instead\n");
      WAIT;
      method = 0;
    }
  }

  if ((method == 2) && (degree !=1)) {
    ERROR("Simple averaging only for linear elements; using SPR instead\n");
    WAIT;
    method = 0;
  }

  rec_struct = new Recovery(norm, method);
}

double RecoveryEstimator::estimate(double timestep)
{
  FUNCNAME("RecoveryEstimator::estimate()");

  const BasisFunction *basFcts = uh_->getFESpace()->getBasisFcts();
  Mesh                *mesh    = uh_->getFESpace()->getMesh();
  int                  dim     = mesh->getDim();
  int                  dow     = Global::getGeo(WORLD);
  double               h1Norm2 = 0.0;
  Flag                 flag    =
    Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_DET | Mesh::FILL_GRD_LAMBDA;
  TraverseStack        stack;
  ElInfo              *elInfo  = NULL;

  int    i, j;
  double det, estEl, errEl, normEl;
  Element *el = NULL;

  if (norm == H1_NORM)     // sets recovery gradient.
    {
      if (method == 2)
	rec_grd   = rec_struct->recovery(uh_, f_vec, f_scal, aux_vec);
      else
	rec_grd   = rec_struct->recovery(uh_, feSpace, f_vec, f_scal, aux_vec);
      rec_basFcts = rec_grd->getFESpace()->getBasisFcts();
    }
  else                     // sets higher-order recovery solution.
    {
      rec_uh      = rec_struct->recoveryUh(uh_, feSpace);
      rec_basFcts = rec_uh->getFESpace()->getBasisFcts();
    }

  int          deg = 2 * std::max(basFcts->getDegree(),
				    rec_basFcts->getDegree());
  Quadrature *quad = Quadrature::provideQuadrature(dim, deg);
  int    numPoints = quad->getNumPoints();

  WorldVector<double> quad_pt;
  WorldVector<double> *grdAtQP = new WorldVector<double>[numPoints];
  WorldVector<double> *recoveryGrdAtQP = new WorldVector<double>[numPoints];

  double *uhAtQP = new double[numPoints];
  double *recoveryUhAtQP = new double[numPoints];

  FastQuadrature *quadFast =
    FastQuadrature::provideFastQuadrature(basFcts, *quad,
					  INIT_PHI | INIT_GRD_PHI);
  FastQuadrature *rec_quadFast =
    FastQuadrature::provideFastQuadrature(rec_basFcts, *quad,
					  INIT_PHI | INIT_GRD_PHI);

  // clear error indicators
  elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
  while (elInfo)
    {
      elInfo->getElement()->setEstimation(0.0, row);
      elInfo = stack.traverseNext(elInfo);
    }

  // traverse mesh
  est_sum = est_max = est_t_sum = 0.0;

  elInfo = stack.traverseFirst(mesh, -1, flag);
  while (elInfo)
    {
      el    = elInfo->getElement();
      det   = elInfo->getDet();
      errEl = 0.0;

      if (norm == H1_NORM)
	{
	  // get gradient and recovery gradient at quadrature points
	  uh_->getGrdAtQPs(elInfo, NULL, quadFast, grdAtQP);
	  rec_grd->getVecAtQPs(elInfo, NULL, rec_quadFast, recoveryGrdAtQP);
	  if (f_scal)
	    {
	      if (aux_vec)
		aux_vec->getVecAtQPs(elInfo, NULL, quadFast, uhAtQP);
	      else
		uh_->getVecAtQPs(elInfo, NULL, quadFast, uhAtQP);
	    }

	  // calc h1 error
	  for (i=0; i<numPoints; i++)
	    {
	      double  err2 = 0.0;
	      double fAtQP = 1.0;
	      if (f_scal)
		fAtQP = (*f_scal)(uhAtQP[i]);
	      if (f_vec)
		{
		  elInfo->coordToWorld(quad->getLambda(i), quad_pt);
		  fAtQP = (*f_vec)(quad_pt);
		}

	      for (j=0; j<dow; j++)
		err2 += sqr(recoveryGrdAtQP[i][j] - fAtQP * grdAtQP[i][j]);
	      errEl += quad->getWeight(i) * err2;
	    }
	}
      else
	{
	  // get vector and recovery vector at quadrature points
	  uh_->getVecAtQPs(elInfo, NULL, quadFast, uhAtQP);
	  rec_uh->getVecAtQPs(elInfo, NULL, rec_quadFast, recoveryUhAtQP);

	  // calc l2 error
	  for (i=0; i<numPoints; i++)
	    errEl += quad->getWeight(i) * sqr(recoveryUhAtQP[i] - uhAtQP[i]);
	}

      estEl = C * det * errEl;
      el->setEstimation(estEl, row);
      est_sum += estEl;
      est_max = std::max(est_max, estEl);

      if (relative_)
	{
	  normEl = 0.0;

	  if (norm == H1_NORM)
	    for (i=0; i<numPoints; i++)
	      {
		double norm2 = 0.0;
		for (j = 0; j < dow; j++)
		  norm2 += sqr(recoveryGrdAtQP[i][j]);
		normEl += quad->getWeight(i) * norm2;
	      }
	  else
	    for (i=0; i<numPoints; i++)
	      normEl += quad->getWeight(i) * sqr(recoveryUhAtQP[i]);

	  h1Norm2 += det * normEl;
	}

      elInfo = stack.traverseNext(elInfo);
    }

  // Computing relative errors
  if (relative_)
    {
      elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
      while (elInfo)
	{
	  estEl  = elInfo->getElement()->getEstimation(row);
	  estEl /= h1Norm2;
	  elInfo->getElement()->setEstimation(estEl, row);
	  elInfo = stack.traverseNext(elInfo);
	}

      est_max /= h1Norm2;
      est_sum /= h1Norm2;
    }

  est_sum = sqrt(est_sum);

  MSG("estimate   = %.8e\n", est_sum);

  delete [] grdAtQP;
  delete [] recoveryGrdAtQP;
  delete [] uhAtQP;
  delete [] recoveryUhAtQP;

  return(est_sum);
}