ResidualEstimator.h 4.78 KB
Newer Older
1 2 3 4
// ============================================================================
// ==                                                                        ==
// == AMDiS - Adaptive multidimensional simulations                          ==
// ==                                                                        ==
5
// ==  http://www.amdis-fem.org                                              ==
6 7
// ==                                                                        ==
// ============================================================================
8 9 10 11 12 13 14 15 16 17 18 19
//
// 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.


20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38

/** \file ResidualEstimator.h */

/** \defgroup Estimator Estimator module
 * @{ <img src="estimator.png"> @}
 */

#ifndef AMDIS_RESIDUALESTIMATOR_H
#define AMDIS_RESIDUALESTIMATOR_H

#include "Estimator.h"
#include "FixVec.h"

namespace AMDiS {

  /** \brief
   * Returns residual square at quadrature point. Not Member of
   * Estimator to avoid multiple instantiation.
   */
Thomas Witkowski's avatar
Thomas Witkowski committed
39 40
  void r(const ElInfo *elInfo,
	 int nPoints,
41
	 const ElementVector& uhIq,
42 43
	 const WorldVector<double> *grdUhIq,
	 const WorldMatrix<double> *D2UhIq,
44
	 const ElementVector& uhOldIq,
45 46 47 48 49
	 const WorldVector<double> *grdUhOldIq,
	 const WorldMatrix<double> *D2UhOldIq,
	 DOFMatrix *A, 
	 DOFVector<double> *fh,
	 Quadrature *quad,
50
	 ElementVector& result);
51
 
52 53 54
  /// Returns pow(det,2.0/dim). Not Member of Estimator to avoid multiple instantiation.
  inline double h2_from_det(double det, int dim) 
  {
55
    return pow(det, 2.0 / dim);
Thomas Witkowski's avatar
Thomas Witkowski committed
56
  }
57 58 59 60 61

  /**
   * \ingroup Estimator
   * 
   * \brief
62
   * Residual estimator.
63 64 65 66
   */
  class ResidualEstimator : public Estimator
  {
  public:
67
    /// Creator class used in the OEMSolverMap.
68 69 70
    class Creator : public EstimatorCreator
    {
    public:
71 72 73
      Creator() 
	: EstimatorCreator() 
      {}
74

75 76
      virtual ~Creator() 
      {}
77

78 79 80
      /// Returns a new ODirSolver object.
      Estimator* create() 
      { 
Thomas Witkowski's avatar
Thomas Witkowski committed
81
	return new ResidualEstimator(name, row);
Thomas Witkowski's avatar
Thomas Witkowski committed
82
      }
83 84
    };
  
85
    /// Constructor.
86
    ResidualEstimator(std::string name, int r);
87

88
    virtual void init(double timestep);
89

90
    /** \brief
91 92
     * Estimates the error on an element. For more information about the
     * parameter, see the description \ref Estimator::estimateElement.
93
     */
94
    virtual void estimateElement(ElInfo *elInfo, DualElInfo *dualElInfo = NULL);
95

96
    virtual void exit(bool output = true);
97 98 99 100 101 102 103

  protected:
    /// Computes the element residual for a given element.
    double computeElementResidual(ElInfo *elInfo, DualElInfo *dualElInfo);

    /// Computes the jump residual for a given element.
    double computeJumpResidual(ElInfo *elInfo, DualElInfo *dualElInfo);
104 105

  protected:
106
    /// Constant in front of element residual
107 108
    double C0;
 
109
    /// Constant in front of edge/face residual
110 111
    double C1;

112
    /// Not used! Was thought to be the constant in front of coarsening term.
113 114
    double C2;

115
    /// Constant in front of the time
116 117
    double C3;

118 119 120 121 122 123 124
    /** \brief
     * Is true, if C1 != 0, and C0 and C3 = 0, hence only the jump residual must be
     * calculated. In this case, some optimizations can be done, because the jump
     * residual is calculated only on second order terms.
     */
    bool jumpResidualOnly;

Thomas Witkowski's avatar
Thomas Witkowski committed
125 126
    /// Number of systems, e.g., number of variables in the equation.
    int nSystems;
127

Thomas Witkowski's avatar
Thomas Witkowski committed
128 129
    /// Number of quadrature points.
    int nPoints;
130 131 132 133 134 135 136 137 138 139 140

    int dim;

    int degree;

    Quadrature *quad;

    FastQuadrature **quadFast;

    const BasisFunction **basFcts;

141
    std::vector<ElementVector> uhEl;
142

143
    std::vector<ElementVector> uhOldEl;
144

145
    std::vector<ElementVector> uhNeigh;
146

147
    ElementVector uhQP;
148

149
    ElementVector uhOldQP;
150

151
    /// Stores the element residual computed at the quadrature points of the element.
152
    ElementVector riq;
153 154 155 156

    WorldVector<double> *grdUh_qp;

    WorldMatrix<double> *D2uhqp; 
157 158 159

    ElInfo *neighInfo;

160
    Quadrature *surfaceQuad;
161

162
    int nPointsSurface;
163

164
    std::vector<WorldVector<double> > grdUhEl;
165

166
    std::vector<WorldVector<double> > grdUhNeigh;
167

168
    std::vector<WorldVector<double> > jump;
169

170
    std::vector<WorldVector<double> > localJump;
171

172
    WorldVector<int> faceIndEl;
173
    
174
    WorldVector<int> faceIndNeigh;
175

176
    DimVec<WorldVector<double> > *lambdaNeigh;
177

178
    DimVec<double> *lambda;
179

180 181
    /// Maximal number of neighbours an element may have in the used dimension.
    int nNeighbours;
182 183 184 185 186 187 188

    /** \brief 
     * Defines for every system if there are second order terms. These values
     * are used to ommit computations of the jump residual that is defined
     * only on second order terms.
     */    
    std::vector<bool> secondOrderTerms;
189 190 191 192
  };
}

#endif // AMDIS_RESIDUALESTIMATOR_H