ResidualEstimator.h 4.79 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
62
63
64
65
66

  /**
   * \ingroup Estimator
   * 
   * \brief
   * Estimator for scalar problems.
   */
  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
91
92
93
    /** \brief
     * Estimates the error on an element. For more information about the parameter,
     * see the description \ref Estimator::estimateElement.
     */
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