ResidualEstimator.h 5.27 KB
Newer Older
1
2
3
4
5
6
// ============================================================================
// ==                                                                        ==
// == AMDiS - Adaptive multidimensional simulations                          ==
// ==                                                                        ==
// ============================================================================
// ==                                                                        ==
7
// ==  TU Dresden                                                            ==
8
// ==                                                                        ==
9
10
11
// ==  Institut fr Wissenschaftliches Rechnen                               ==
// ==  Zellescher Weg 12-14                                                  ==
// ==  01069 Dresden                                                         ==
12
13
14
15
// ==  germany                                                               ==
// ==                                                                        ==
// ============================================================================
// ==                                                                        ==
16
// ==  https://gforge.zih.tu-dresden.de/projects/amdis/                      ==
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
// ==                                                                        ==
// ============================================================================

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

  /**
   * \ingroup Estimator
   * 
   * \brief
   * Estimator for scalar problems.
   */
  class ResidualEstimator : public Estimator
  {
  public:
66
    /// Creator class used in the OEMSolverMap.
67
68
69
    class Creator : public EstimatorCreator
    {
    public:
Thomas Witkowski's avatar
Thomas Witkowski committed
70
      Creator() : EstimatorCreator() {}
71

Thomas Witkowski's avatar
Thomas Witkowski committed
72
      virtual ~Creator() {}
73

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

84
    void init(double timestep);
85

86
87
88
89
90
    /** \brief
     * Estimates the error on an element. For more information about the parameter,
     * see the description \ref Estimator::estimateElement.
     */
    void estimateElement(ElInfo *elInfo, DualElInfo *dualElInfo = NULL);
91

92
93
94
95
96
97
98
99
    void exit(bool output = true);

  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);
100
101

  protected:
102
    /// Constant in front of element residual
103
104
    double C0;
 
105
    /// Constant in front of edge/face residual
106
107
    double C1;

108
    /// Not used! Was thought to be the constant in front of coarsening term.
109
110
    double C2;

111
    /// Constant in front of the time
112
113
    double C3;

114
115
116
117
118
119
120
    /** \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
121
122
    /// Number of systems, e.g., number of variables in the equation.
    int nSystems;
123

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

    int dim;

    int degree;

    Quadrature *quad;

    FastQuadrature **quadFast;

    const BasisFunction **basFcts;

    double **uhEl;

    double **uhOldEl;

141
142
    double **uhNeigh;

143
144
145
146
    double *uhQP;

    double *uhOldQP;

147
    /// Stores the element residual computed at the quadrature points of the element.
148
    double *riq;
149
150
151
152

    WorldVector<double> *grdUh_qp;

    WorldMatrix<double> *D2uhqp; 
153
154
155

    ElInfo *neighInfo;

156
    Quadrature *surfaceQuad;
157

158
    int nPointsSurface;
159

160
    std::vector<WorldVector<double> > grdUhEl;
161

162
    std::vector<WorldVector<double> > grdUhNeigh;
163

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

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

168
    WorldVector<int> faceIndEl;
169
    
170
    WorldVector<int> faceIndNeigh;
171

172
    DimVec<WorldVector<double> > *lambdaNeigh;
173

174
    DimVec<double> *lambda;
175

176
177
    /// Maximal number of neighbours an element may have in the used dimension.
    int nNeighbours;
178
179
180
181
182
183
184

    /** \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;
185
186
187
188
  };
}

#endif // AMDIS_RESIDUALESTIMATOR_H