SimpleResidualEstimator ready and documented.

AMDiS/src/ResidualEstimator.cc → AMDiS/src/est/ResidualEstimator.cc

@@ -322,7 +322,7 @@ namespace AMDiS {
 	    }
 	    if (D2uhqp == NULL && degree > 2 && (*it)->secondOrderTerms()) { 
 	      D2uhqp = new WorldMatrix<double>[nPoints];
-	      uh[system]->getD2AtQPs(elInfo, NULL, quadFast[system], D2uhqp);	    
+	      uh[system]->getD2AtQPs(elInfo, NULL, quadFast[system], D2uhqp);  
 	    }
 	  }
 	}
@@ -533,11 +533,13 @@ namespace AMDiS {
 
       if (factor) {
 	if (D2UhIq)
-	  (*it)->evalSecondOrder(nPoints, uhIq, grdUhIq, D2UhIq, result, -factor);	
+	  (*it)->evalSecondOrder(nPoints, uhIq, grdUhIq, D2UhIq, result, -factor);
 
 	if (grdUhIq) {
-	  (*it)->evalFirstOrderGrdPsi(nPoints, uhIq, grdUhIq, D2UhIq, result, factor);
-	  (*it)->evalFirstOrderGrdPhi(nPoints, uhIq, grdUhIq, D2UhIq, result, factor);
+	  (*it)->evalFirstOrderGrdPsi(nPoints, uhIq, grdUhIq, D2UhIq, 
+				      result, factor);
+	  (*it)->evalFirstOrderGrdPhi(nPoints, uhIq, grdUhIq, D2UhIq, 
+				      result, factor);
 	}
 	
 	if (num_rows(uhIq) > 0)

AMDiS/src/ResidualEstimator.h → AMDiS/src/est/ResidualEstimator.h

@@ -59,7 +59,7 @@ namespace AMDiS {
    * \ingroup Estimator
    * 
    * \brief
-   * Estimator for scalar problems.
+   * Residual estimator.
    */
   class ResidualEstimator : public Estimator
   {
@@ -88,8 +88,8 @@ namespace AMDiS {
     virtual void init(double timestep);
 
     /** \brief
-     * Estimates the error on an element. For more information about the parameter,
-     * see the description \ref Estimator::estimateElement.
+     * Estimates the error on an element. For more information about the
+     * parameter, see the description \ref Estimator::estimateElement.
      */
     virtual void estimateElement(ElInfo *elInfo, DualElInfo *dualElInfo = NULL);
 

AMDiS/src/est/SimpleResidualEstimator.cc

+//
+// Software License for AMDiS
+//
+// Copyright (c) 2010 Dresden University of Technology 
+// All rights reserved.
+// Authors: Simon Vey, Thomas Witkowski et al.
+//
+// This file is part of AMDiS
+//
+// See also license.opensource.txt in the distribution.
+
+
+#include "est/SimpleResidualEstimator.h"
+#include "Operator.h"
+#include "DOFMatrix.h"
+#include "DOFVector.h"
+#include "Assembler.h"
+#include "Traverse.h"
+#include "Initfile.h"
+
+namespace AMDiS {
+
+  SimpleResidualEstimator::SimpleResidualEstimator(std::string name) 
+    : Estimator(name, 0),
+      C0(0.0), 
+      C1(0.0)
+  {
+    FUNCNAME("SimpleResidualEstimator::SimpleResidualEstimator()");
+
+    // === Read parameters C0 and C1 from init file. ===
+
+    Parameters::get(name + "->C0", C0);
+    Parameters::get(name + "->C1", C1);
+
+    C0 = C0 > 1.e-25 ? sqr(C0) : 0.0;
+    C1 = C1 > 1.e-25 ? sqr(C1) : 0.0;
+  }
+
+
+  void SimpleResidualEstimator::init(double)
+  {
+    FUNCNAME("SimpleResidualEstimator::init()");
+
+    // === Create data structures. ===
+
+    basFcts = uh[0]->getFeSpace()->getBasisFcts();
+    dim = mesh->getDim();
+    degree = basFcts->getDegree() * 2;
+    quad = Quadrature::provideQuadrature(dim, degree);
+    nPoints = quad->getNumPoints();
+
+    Flag flag = INIT_PHI;
+    if (degree > 2)
+      flag |= INIT_D2_PHI;    
+    quadFast = FastQuadrature::provideFastQuadrature(basFcts, *quad, flag);
+  
+    uhEl.change_dim(basFcts->getNumber()); 
+    uhNeigh.change_dim(basFcts->getNumber());
+    riq.change_dim(nPoints);
+    D2uhqp = NULL;
+
+    // === Clear error indicators and mark elements for jumpRes. ===
+
+    TraverseStack stack;
+    ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
+    while (elInfo) {
+      elInfo->getElement()->setEstimation(0.0, 0);
+      elInfo->getElement()->setMark(1);
+      elInfo = stack.traverseNext(elInfo);
+    }
+
+    est_sum = 0.0;
+    est_max = 0.0;
+
+    traverseFlag = 
+      Mesh::FILL_NEIGH      |
+      Mesh::FILL_COORDS     |
+      Mesh::FILL_OPP_COORDS |
+      Mesh::FILL_BOUND      |
+      Mesh::FILL_GRD_LAMBDA |
+      Mesh::FILL_DET        |
+      Mesh::CALL_LEAF_EL;
+    neighInfo = mesh->createNewElInfo();
+
+
+    // === Prepare date for computing jump residual. ===
+
+    if (C1 > 0.0 && dim > 1) {
+      surfaceQuad = Quadrature::provideQuadrature(dim - 1, degree);
+      nPointsSurface = surfaceQuad->getNumPoints();
+      grdUhEl.resize(nPointsSurface);
+      grdUhNeigh.resize(nPointsSurface);
+      jump.resize(nPointsSurface);
+      localJump.resize(nPointsSurface);
+      nNeighbours = Global::getGeo(NEIGH, dim);
+      lambdaNeigh = new DimVec<WorldVector<double> >(dim, NO_INIT);
+      lambda = new DimVec<double>(dim, NO_INIT);
+    }
+  }
+
+
+  void SimpleResidualEstimator::exit(bool output)
+  {
+    FUNCNAME("SimpleResidualEstimator::exit()");
+
+    // === Calculate the root of the estimations and make output. ===
+
+    est_sum = sqrt(est_sum);
+
+    if (output)
+      MSG("estimate = %.8e\n", est_sum);
+
+    /// === Delete data structures. ===
+
+    if (D2uhqp != NULL)
+      delete [] D2uhqp;
+
+    if (C1 && (dim > 1)) {
+      delete lambdaNeigh;
+      delete lambda;
+    }
+
+    delete neighInfo;
+  }
+
+
+  void SimpleResidualEstimator::estimateElement(ElInfo *elInfo, DualElInfo *)
+  {    
+    FUNCNAME("SimpleResidualEstimator::estimateElement()");
+
+    TEST_EXIT(matrix[0])("Should not happen!\n");
+
+    // Get pointer to element object.
+    Element *el = elInfo->getElement();
+    // Get error estimation of this element.
+    double est_el = el->getEstimation(0);
+    // Shortcut for the system matrix.
+    DOFMatrix *dofMat = const_cast<DOFMatrix*>(matrix[0]);
+    // Shortcut for the right hand side DOF vector
+    DOFVector<double> *dofVec = const_cast<DOFVector<double>*>(fh[0]);
+    
+    // === Init assembler ===
+
+    std::vector<Operator*>::iterator it;
+    std::vector<double*>::iterator itfac;
+    // Matrix assembler are only initialized with the corresponding multiply
+    // factors are != 0
+    for (it = dofMat->getOperatorsBegin(), itfac = dofMat->getOperatorEstFactorBegin();
+	 it != dofMat->getOperatorsEnd(); ++it, ++itfac)
+      if (*itfac == NULL || **itfac != 0.0)
+	(*it)->getAssembler()->initElement(elInfo, NULL, quad);	  
+
+
+    // Vector assembler are only initialized if C0 is set. Note that the jump 
+    // residual (thus C1) does not contain the right hand side.
+    if (C0 > 0.0)
+      for (it = dofVec->getOperatorsBegin(); it != dofVec->getOperatorsEnd(); ++it)
+	(*it)->getAssembler()->initElement(elInfo, NULL, quad);	  
+
+
+    // === Compute element residuals and time error estimation. ===
+    if (C0 > 0.0)
+      est_el += computeElementResidual(elInfo);
+
+    // === Compute jump residuals. ===
+    if (C1 > 0.0 && dim > 1)
+      est_el += computeJumpResidual(elInfo);
+
+    // === Update global residual variables. ===
+    el->setEstimation(est_el, 0);
+    el->setMark(0);
+    est_sum += est_el;
+    est_max = std::max(est_max, est_el);
+  }
+
+
+  double SimpleResidualEstimator::computeElementResidual(ElInfo *elInfo)
+  {
+    FUNCNAME("SimpleResidualEstimator::computeElementResidual()");
+
+    double det = elInfo->getDet();
+    double h2 = h2_from_det(det, dim);
+    riq = 0.0;
+    DOFMatrix *dofMat = const_cast<DOFMatrix*>(matrix[0]);
+    DOFVector<double> *dofVec = const_cast<DOFVector<double>*>(fh[0]);
+
+    // === If there is a valid left hand side operator get the solution  ===
+    // === vector or its derivations on the quadrature points of the     ===
+    // === element.                                                      ===
+    std::vector<Operator*>::iterator it;
+    std::vector<double*>::iterator itfac;      
+    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()) {
+	  uhQP.change_dim(nPoints);
+	  uh[0]->getVecAtQPs(elInfo, NULL, quadFast, uhQP);
+	}
+
+	if (D2uhqp == NULL && degree > 2 && (*it)->secondOrderTerms()) { 
+	  D2uhqp = new WorldMatrix<double>[nPoints];
+	  uh[0]->getD2AtQPs(elInfo, NULL, quadFast, D2uhqp);
+	}
+      }
+    }
+    
+    // === Compute the element residual and store it in irq. ===    
+    r(elInfo, nPoints, uhQP, D2uhqp, dofMat, dofVec, quad, riq);
+    
+    
+    // === Add integral over r square. ===
+    double result = 0.0;
+    for (int iq = 0; iq < nPoints; iq++)
+      result += quad->getWeight(iq) * riq[iq] * riq[iq];
+   
+    if (norm == NO_NORM || norm == L2_NORM)
+      result = C0 * h2 * h2 * det * result;
+    else
+      result = C0 * h2 * det * result;
+    
+    return result;
+  }
+
+
+  double SimpleResidualEstimator::computeJumpResidual(ElInfo *elInfo)
+  {
+    FUNCNAME("SimpleResidualEstimator::computeJumpResidual()");
+
+    // === Init temporary variables. ===
+    double result = 0.0;
+    int dow = Global::getGeo(WORLD);
+    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++) {
+      // Pointer to neighbouring element.
+      Element *neigh = const_cast<Element*>(elInfo->getNeighbour(face));
+
+      // If there is no neighbour, or we have already visited this pair, 
+      // then continue with next one.
+      if (!(neigh && neigh->getMark()))
+	continue;
+
+      // === The next lines are only to set all information about the ===
+      // === neighbouring element.                                    ===
+      int oppV = elInfo->getOppVertex(face);	
+      el->sortFaceIndices(face, &faceIndEl);
+      neigh->sortFaceIndices(oppV, &faceIndNeigh);	
+      neighInfo->setElement(const_cast<Element*>(neigh));
+      neighInfo->setFillFlag(Mesh::FILL_COORDS);
+      
+      for (int i = 0; i < dow; i++)
+	neighInfo->getCoord(oppV)[i] = elInfo->getOppCoord(face)[i];
+      
+      for (int i = 0; i < dim; i++) {
+	int i1 = faceIndEl[i];
+	int i2 = faceIndNeigh[i];
+	for (int j = 0; j < dow; j++)
+	  neighInfo->getCoord(i2)[j] = elInfo->getCoord(i1)[j];
+      }      
+	
+      // Compute determinant of the neighbouring element.
+      double detNeigh = abs(neighInfo->calcGrdLambda(*lambdaNeigh));
+           
+      // Reset data.
+      for (int iq = 0; iq < nPointsSurface; iq++)
+	jump[iq].set(0.0);     
+      
+      // Get solution vector on the nodes of the current element and 
+      // its neighbour.
+      uh[0]->getLocalVector(el, uhEl);
+      uh[0]->getLocalVector(neigh, uhNeigh);
+
+      // Compute the jump of the gradients of the solution over the face.
+      for (int iq = 0; iq < nPointsSurface; iq++) {
+	(*lambda)[face] = 0.0;
+	for (int i = 0; i < dim; i++)
+	  (*lambda)[faceIndEl[i]] = surfaceQuad->getLambda(iq, i);
+	
+	basFcts->evalGrdUh(*lambda, grdLambda, uhEl, &grdUhEl[iq]);
+	
+	(*lambda)[oppV] = 0.0;
+	for (int i = 0; i < dim; i++)
+	  (*lambda)[faceIndNeigh[i]] = surfaceQuad->getLambda(iq, i);
+	
+	basFcts->evalGrdUh(*lambda, *lambdaNeigh, uhNeigh, &grdUhNeigh[iq]);
+	
+	grdUhEl[iq] -= grdUhNeigh[iq];
+      } // for iq				
+      
+
+      // === Compute the action of the second order operators on the jump. ===
+      std::vector<double*>::iterator fac;
+      std::vector<Operator*>::iterator it;
+      DOFMatrix *mat = const_cast<DOFMatrix*>(matrix[0]);
+      for (it = mat->getOperatorsBegin(), fac = mat->getOperatorEstFactorBegin(); 
+	   it != mat->getOperatorsEnd(); ++it, ++fac) {
+	
+	if (*fac == NULL || **fac != 0.0) {
+	  for (int iq = 0; iq < nPointsSurface; iq++)
+	    localJump[iq].set(0.0);
+	  
+	  (*it)->weakEvalSecondOrder(grdUhEl, localJump);
+	  
+	  double factor = *fac ? **fac : 1.0;
+	  if (factor != 1.0)
+	    for (int i = 0; i < nPointsSurface; i++)
+	      localJump[i] *= factor;
+	  
+	  for (int i = 0; i < nPointsSurface; i++)
+	    jump[i] += localJump[i];
+	}		
+      } // for (it = ...
+      
+
+      // === Compute the squared integral of the jump and multiply with ===
+      // === the appropriate constants.                                 ===
+      
+      double val = 0.0;
+      for (int iq = 0; iq < nPointsSurface; iq++)
+	val += surfaceQuad->getWeight(iq) * (jump[iq] * jump[iq]);
+
+      double d = 0.5 * (det + detNeigh);
+   
+      if (norm == NO_NORM || norm == L2_NORM)
+	val *= C1 * h2_from_det(d, dim) * 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]);    
+    if (norm == NO_NORM || norm == L2_NORM)
+      val *= C1 * h2;
+    else
+      val *= C1;   
+    result += val;
+
+    return result;
+  }
+
+
+  void SimpleResidualEstimator::r(const ElInfo *elInfo,
+				  int nPoints,
+				  const ElementVector& uhIq,
+				  const WorldMatrix<double> *D2UhIq,
+				  DOFMatrix *A, 
+				  DOFVector<double> *fh,
+				  Quadrature *quad,
+				  ElementVector& result)
+  {
+    FUNCNAME("SimpleResidualEstimator::r()");
+
+    std::vector<Operator*>::iterator it;
+    std::vector<double*>::iterator fac;
+
+    // === Get left hand side on quadrature points. ===
+    for (it = A->getOperatorsBegin(), fac = A->getOperatorEstFactorBegin(); 
+	 it != A->getOperatorsEnd(); ++it, ++fac) {
+     
+      double factor = *fac ? **fac : 1.0;
+
+      if (factor) {
+	if (D2UhIq)
+	  (*it)->evalSecondOrder(nPoints, uhIq, NULL, D2UhIq, result, -factor);
+
+	if (num_rows(uhIq) > 0)
+	  (*it)->evalZeroOrder(nPoints, uhIq, NULL, D2UhIq, result, factor);
+      }
+    }
+    
+
+    // === Get right hand side on quadrature points. ===
+    for (it = fh->getOperatorsBegin(), fac = fh->getOperatorEstFactorBegin(); 
+	 it != fh->getOperatorsEnd(); ++it, ++fac) {
+
+      double factor = *fac ? **fac : 1.0;
+
+      if (factor) {
+	ElementVector fx(nPoints, 0.0);
+	(*it)->getC(elInfo, nPoints, fx);
+	
+	for (int iq = 0; iq < nPoints; iq++)
+	  result[iq] -= factor * fx[iq];
+      }
+    }    
+  }
+
+
+}

AMDiS/src/est/SimpleResidualEstimator.h

+// ============================================================================
+// ==                                                                        ==
+// == AMDiS - Adaptive multidimensional simulations                          ==
+// ==                                                                        ==
+// ==  http://www.amdis-fem.org                                              ==
+// ==                                                                        ==
+// ============================================================================
+//
+// Software License for AMDiS
+//
+// Copyright (c) 2010 Dresden University of Technology 
+// All rights reserved.
+// Authors: Simon Vey, Thomas Witkowski et al.
+//
+// This file is part of AMDiS
+//
+// See also license.opensource.txt in the distribution.
+
+
+
+/** \file SimpleResidualEstimator.h */
+
+/** \defgroup Estimator Estimator module
+ * @{ <img src="estimator.png"> @}
+ */
+
+#ifndef AMDIS_SIMPLE_RESIDUAL_ESTIMATOR_H
+#define AMDIS_SIMPLE_RESIDUAL_ESTIMATOR_H
+
+#include "Estimator.h"
+#include "FixVec.h"
+
+namespace AMDiS {
+
+  /**
+   * \ingroup Estimator
+   * 
+   * \brief
+   * Simple Residual estimator.
+   *
+   * The restrictions for using the simple residual estimator are the following:
+   *   - no system of PDEs, thus a PDE with only one unknown
+   *   - no first order terms
+   *   - no time dependent terms
+   */
+  class SimpleResidualEstimator : public Estimator
+  {
+  public:
+    /// Creator class used in the OEMSolverMap.
+    class Creator : public EstimatorCreator
+    {
+    public:
+      Creator() 
+	: EstimatorCreator() 
+      {}
+
+      virtual ~Creator() 
+      {}
+
+      /// Returns a new ODirSolver object.
+      Estimator* create() 
+      { 
+	return new SimpleResidualEstimator(name);
+      }
+    };
+  
+    /// Constructor.
+    SimpleResidualEstimator(std::string name);
+
+    /// Init the error estimator.
+    void init(double);
+
+    /// Run estimator on one element. 
+    /// \param[in]   elInfo      Info object for the element to be estimated.
+    /// \param[in]   dualElInfo  Not used here. In general, this may be used for
+    ///                          estimating with the multi-mesh technique.
+    void estimateElement(ElInfo *elInfo, DualElInfo *dualElInfo = NULL);
+
+    /// Finalize the error estimator, i.e., delete all temporary data structures.
+    void exit(bool output = true);
+
+    /// Returns pow(det,2.0/dim).
+    inline double h2_from_det(double det, int dim) 
+    {
+      return pow(det, 2.0 / dim);
+    }
+
+  protected:
+    /// Computes the element residual for a given element.
+    double computeElementResidual(ElInfo *elInfo);
+
+    /// Computes the jump residual for a given element.
+    double computeJumpResidual(ElInfo *elInfo);
+
+    /** \brief
+     * Returns residual square at quadrature point. Not Member of
+     * Estimator to avoid multiple instantiation.
+     *
+     * \param[in]  elInfo   Current element information object.
+     * \param[in]  nPoints  Number of quadrature points on this element.
+     * \param[in]  uhIq     Solution vector on element quadrature points.
+     * \param[in]  D2UhIq   Second derivations of the solution vector on
+     *                      element quadrature points.
+     * \param[in]  A        Matrix that contains the left hand side operators.
+     * \param[in]  fh       Vector that contains the right hand side operators.
+     * \param[in]  quad     Object for numerical quadrature.
+     * \param[out] result   Vector containing the residual on the quadrature 
+     *                      points.
+     */
+    void r(const ElInfo *elInfo,
+	   int nPoints,
+	   const ElementVector& uhIq,
+	   const WorldMatrix<double> *D2UhIq,
+	   DOFMatrix *A, 
+	   DOFVector<double> *fh,
+	   Quadrature *quad,
+	   ElementVector& result);
+ 
+  protected:
+    /// Constant in front of element residual.
+    double C0;
+ 
+    /// Constant in front of edge/face residual.
+    double C1;
+
+    /// Number of quadrature points.
+    int nPoints;
+
+    /// Dimension of the mesh.
+    int dim;
+
+    /// Polynomial degree.
+    int degree;
+
+    /// Object for numerical quadrature
+    Quadrature *quad;
+
+    /// Object for fast numerical quadrature
+    FastQuadrature *quadFast;
+
+    /// Pointer to the basis functions of the FE space.
+    const BasisFunction *basFcts;
+
+    /// Vector that stores all global DOFs of one element.
+    ElementVector uhEl;
+
+    /// Vector that stores all global DOFs of one element.
+    ElementVector uhNeigh;
+
+    /// Vector that stores values on all quadrature points (QP) on one element.
+    ElementVector uhQP;
+
+    /// Matrix that stores the second derivations on all quadrature points
+    /// on one element.
+    WorldMatrix<double> *D2uhqp; 
+
+    /// Stores the element residual computed at the quadrature points of 
+    /// the element.
+    ElementVector riq;
+
+    /// Pointer to the information object of some neighbouring element.