vecheat.cc 8.09 KB
Newer Older
1
2
3
4
5
6
7
8
9
#include "AMDiS.h"

using namespace std;
using namespace AMDiS;

// ===========================================================================
// ===== function definitions ================================================
// ===========================================================================

10
/// Dirichlet boundary function
11
12
13
14
15
class G : public AbstractFunction<double, WorldVector<double> >,
	  public TimedObject
{
public:

16
17
18
  /// Implementation of AbstractFunction::operator().
  double operator()(const WorldVector<double>& x) const 
  {
19
    return sin(M_PI * (*timePtr)) * exp(-10.0 * (x * x));
20
  }
21
22
};

23
/// RHS function
24
25
26
27
class F : public AbstractFunction<double, WorldVector<double> >,
	  public TimedObject
{
public:
28
  F(int degree) : AbstractFunction<double, WorldVector<double> >(degree) {}
29

30
31
32
  /// Implementation of AbstractFunction::operator().
  double operator()(const WorldVector<double>& x) const 
  {
33
    int dim = Global::getGeo(WORLD);
34
35
36
37
    double r2 = x * x;
    double ux = sin(M_PI * (*timePtr)) * exp(-10.0 * r2);
    double ut = M_PI * cos(M_PI * (*timePtr)) * exp(-10.0 * r2);
    return ut - (400.0 * r2 - 20.0 * dim) * ux;
38
39
40
41
42
43
44
  };
};

// ===========================================================================
// ===== instationary problem ================================================
// ===========================================================================

45
/// Instationary problem
46
class Vecheat : public ProblemInstat
47
48
49
{
public:

50
  /// Constructor
51
52
  Vecheat(ProblemStat &vecheatSpace) 
    : ProblemInstat("vecheat", vecheatSpace)
53
54
55
  {
    // init theta scheme
    theta = -1.0;
56
    Parameters::get(name + "->theta", theta);
57
58
59
60
61
62
63
    TEST_EXIT(theta >= 0)("theta not set!\n");
    if (theta == 0.0) {
      WARNING("You are using the explicit Euler scheme\n");
      WARNING("Use a sufficiently small time step size!!!\n");
    }
    MSG("theta = %f\n", theta);
    theta1 = theta - 1;
64
  }
65
66
67

  // ===== ProblemInstatBase methods ===================================

68
  /// set the time in all needed functions!
69
70
  void setTime(AdaptInfo *adaptInfo) 
  {
71
    ProblemInstat::setTime(adaptInfo);
72
    rhsTime = adaptInfo->getTime() - (1-theta) * adaptInfo->getTimestep();
73
  }
74
75
76

  // ===== initial problem methods =====================================

77
  /// Used by \ref problemInitial to solve the system of the initial problem
78
79
  void solveInitialProblem(AdaptInfo *adaptInfo) 
  {
80
    int size = problemStat->getNumComponents();
81
    rhsTime = 0.0;
82
    for (int i = 0; i < size; i++) {
83
84
85
      problemStat->getMesh(i)->dofCompress();
      problemStat->getSolution()->getDOFVector(i)->interpol(boundaryFct);
    }
86
  }
87

88
  /// Used by \ref problemInitial to do error estimation for the initial problem.
89
90
  void estimateInitial(AdaptInfo *adaptInfo) 
  {
91
    int size = problemStat->getNumComponents();
92
93
    double errMax, errSum;

94
    for (int i = 0; i < size; i++) {
95
96
97
98
99
100
      errSum = Error<double>::L2Err(*boundaryFct,
				    *(problemStat->getSolution()->
				      getDOFVector(i)), 
				    0, &errMax, false);
      adaptInfo->setEstSum(errSum, i);
    }
101
102
103
104
105
106
107
108
109
110
  }

  /// Used by \ref problemInitial to build before refinement.
  void buildBeforeRefineInitial(Flag) {}

  /// Used by \ref problemInitial to build before coarsening.
  void buildBeforeCoarsenInitial(Flag) {}

  /// Used by \ref problemInitial to build after coarsening.
  void buildAfterCoarsenInitial(Flag) {}
111

112
  /// Used by \ref problemInitial to mark elements.
113
114
115
  Flag markElementsInitial() 
  { 
    return 0; 
116
  }
117

118
  /// Used by \ref problemInitial
119
120
121
  Flag refineMeshInitial() 
  { 
    return 0; 
122
  }
123

124
  /// Used by \ref problemInitial
125
126
127
  Flag coarsenMeshInitial() 
  { 
    return 0; 
128
  }
129

130
131
  /// Used by \ref problemInitial
  void beginIterationInitial(int) {}
132

133
134
  /// Used by \ref problemInitial
  void endIterationInitial(int) {}
135
136
137

  // ===== setting methods ===============================================

138
  /// Sets \ref boundaryFct;
139
140
141
142
143
144
145
  void setBoundaryFct(AbstractFunction<double, WorldVector<double> > *fct) 
  {
    boundaryFct = fct;
  } 

  // ===== getting methods ===============================================

146
  /// Returns pointer to \ref theta.
147
148
149
  double *getThetaPtr() 
  { 
    return &theta; 
150
  }
151

152
  /// Returns pointer to \ref theta1.
153
154
155
  double *getTheta1Ptr() 
  { 
    return &theta1; 
156
  }
157

158
  /// Returns pointer to \ref rhsTime.
159
  double *getRhsTimePtr() 
160
161
  { 
    return &rhsTime; 
162
  }
163

164
  /// Returns \ref boundaryFct;
165
166
167
  AbstractFunction<double, WorldVector<double> > *getBoundaryFct() 
  {
    return boundaryFct;
168
  }
169
170

protected:
171
  /// Used for theta scheme.
172
173
  double theta;

174
  /// theta - 1
175
176
  double theta1;

177
  /// time for right hand side functions.
178
179
  double rhsTime;

180
  /// Pointer to boundary function. Needed for initial problem.
181
182
183
184
185
186
187
188
189
190
191
  AbstractFunction<double, WorldVector<double> > *boundaryFct;
};

// ===========================================================================
// ===== main program ========================================================
// ===========================================================================

int main(int argc, char** argv)
{
  FUNCNAME("main");

192
  AMDiS::init(argc, argv);
193
194

  // ===== create and init stationary problem =====
195
  ProblemStat vecheatSpace("vecheat->space");
196
  vecheatSpace.initialize(INIT_ALL);
197
198

  // ===== create instationary problem =====
199
200
  Vecheat vecheat(vecheatSpace);
  vecheat.initialize(INIT_ALL);
201
202

  // create adapt info
203
  AdaptInfo adaptInfo("vecheat->adapt", vecheatSpace.getNumComponents());
204
205

  // create initial adapt info
206
  AdaptInfo adaptInfoInitial("vecheat->initial->adapt", vecheatSpace.getNumComponents());
207
208

  // create instationary adapt
209
210
211
212
213
214
215
  AdaptInstationary adaptInstat("vecheat->adapt",
				vecheatSpace,
				adaptInfo,
				vecheat,
				adaptInfoInitial);

  
216
  // ===== create rhs functions =====
217
  F *rhsFct0 = new F(vecheatSpace.getFeSpace(0)->getBasisFcts()->getDegree());
218
219
  rhsFct0->setTimePtr(vecheat.getRhsTimePtr());

220
  F *rhsFct1 = new F(vecheatSpace.getFeSpace(1)->getBasisFcts()->getDegree());
221
222
  rhsFct1->setTimePtr(vecheat.getRhsTimePtr());

223
224
225
226
227
228

  // ===== create operators =====
  double one = 1.0;
  double zero = 0.0;

  // create laplace
229
  Operator A00 (vecheatSpace.getFeSpace(0), vecheatSpace.getFeSpace(0));
230
  A00.addTerm(new Simple_SOT);
231
232
  A00.setUhOld(vecheat.getOldSolution()->getDOFVector(0));

233
  Operator A11(vecheatSpace.getFeSpace(1), vecheatSpace.getFeSpace(1));
234
  A11.addTerm(new Simple_SOT);
235
236
237
238
239
  A11.setUhOld(vecheat.getOldSolution()->getDOFVector(1));

  if (*(vecheat.getThetaPtr()) != 0.0) {
    vecheatSpace.addMatrixOperator(A00, 0, 0, vecheat.getThetaPtr(), &one);
    vecheatSpace.addMatrixOperator(A11, 1, 1, vecheat.getThetaPtr(), &one);
240
241
  }

242
243
244
  if (*(vecheat.getTheta1Ptr()) != 0.0) {
    vecheatSpace.addVectorOperator(A00, 0, vecheat.getTheta1Ptr(), &zero);
    vecheatSpace.addVectorOperator(A11, 1, vecheat.getTheta1Ptr(), &zero);
245
246
247
  }

  // create zero order operator
248
  Operator C00(vecheatSpace.getFeSpace(0), vecheatSpace.getFeSpace(0));
249
  C00.addTerm(new Simple_ZOT);
250
251
  C00.setUhOld(vecheat.getOldSolution()->getDOFVector(0));
  vecheatSpace.addMatrixOperator(C00, 0, 0, 
252
253
				 vecheat.getInvTau(), vecheat.getInvTau());
  vecheatSpace.addVectorOperator(C00, 0, vecheat.getInvTau(), vecheat.getInvTau());
254
255
  
  
256
  Operator C11(vecheatSpace.getFeSpace(1), vecheatSpace.getFeSpace(1));
257
  C11.addTerm(new Simple_ZOT);
258
259
  C11.setUhOld(vecheat.getOldSolution()->getDOFVector(1));
  vecheatSpace.addMatrixOperator(C11, 1, 1, 			    
260
				 vecheat.getInvTau(), vecheat.getInvTau());
261
  vecheatSpace.addVectorOperator(C11, 1,			    
262
				 vecheat.getInvTau(), vecheat.getInvTau());
263
264

  // create RHS operator
265
  Operator F0(vecheatSpace.getFeSpace(0));
266
  F0.addTerm(new CoordsAtQP_ZOT(rhsFct0));
267
  vecheatSpace.addVectorOperator(F0, 0);
268

269
  Operator F1 = Operator(vecheatSpace.getFeSpace(1));
270
  F1.addTerm(new CoordsAtQP_ZOT(rhsFct1));
271
  vecheatSpace.addVectorOperator(F1, 1);
272

273

274
275
  // ===== create boundary functions =====
  G *boundaryFct = new G;
276
  boundaryFct->setTimePtr(vecheat.getTime());
277
278
  vecheat.setBoundaryFct(boundaryFct);

279
280
281
  vecheatSpace.addDirichletBC(1, 0, 0, boundaryFct);
  vecheatSpace.addDirichletBC(1, 1, 1, boundaryFct);

282

283
  // ===== start adaption loop =====
284
  int errorCode = adaptInstat.adapt();
285

286
287
  AMDiS::finalize();

288
289
  return errorCode;
}