vecheat.cc 9.87 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 = x.getSize();
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
47
48
49
class Vecheat : public ProblemInstatVec
{
public:

50
  /// Constructor
51
52
53
54
55
56
57
58
59
60
61
62
63
  Vecheat(ProblemVec *vecheatSpace) 
    : ProblemInstatVec("vecheat", vecheatSpace)
  {
    // init theta scheme
    theta = -1.0;
    GET_PARAMETER(0, name + "->theta", "%f", &theta);
    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
71
72
73
  void setTime(AdaptInfo *adaptInfo) 
  {
    rhsTime = adaptInfo->getTime() - (1-theta) * adaptInfo->getTimestep();
    boundaryTime = adaptInfo->getTime();
    tau1 = 1.0 / adaptInfo->getTimestep();    
74
  }
75
76
77

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

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

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

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

  /// 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) {}
113

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

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

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

132
133
  /// Used by \ref problemInitial
  void beginIterationInitial(int) {}
134

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

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

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

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

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

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

160
  /// Returns pointer to \ref tau1
161
162
163
  double *getTau1Ptr() 
  { 
    return &tau1; 
164
  }
165

166
  /// Returns pointer to \ref rhsTime.
167
168
169
  double *getRHSTimePtr() 
  { 
    return &rhsTime; 
170
  }
171

172
  /// Returns pointer to \ref theta1.
173
174
175
  double *getBoundaryTimePtr() 
  { 
    return &boundaryTime; 
176
  }
177

178
  /// Returns \ref boundaryFct;
179
180
181
  AbstractFunction<double, WorldVector<double> > *getBoundaryFct() 
  {
    return boundaryFct;
182
  }
183
184

protected:
185
  /// Used for theta scheme.
186
187
  double theta;

188
  /// theta - 1
189
190
  double theta1;

191
  /// 1.0 / timestep
192
193
  double tau1;

194
  /// time for right hand side functions.
195
196
  double rhsTime;

197
  /// time for boundary functions.
198
199
  double boundaryTime;

200
  /// Pointer to boundary function. Needed for initial problem.
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
  AbstractFunction<double, WorldVector<double> > *boundaryFct;
};

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

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

  // ===== check for init file =====
  TEST_EXIT(argc > 1)("usage: vecheat initfile\n");

  // ===== init parameters =====
  Parameters::init(false, argv[1]);
  Parameters::readArgv(argc, argv);

  // ===== create and init stationary problem =====
220
  ProblemVec *vecheatSpace = new ProblemVec("vecheat->space");
221
222
223
  vecheatSpace->initialize(INIT_ALL);

  // ===== create instationary problem =====
224
  Vecheat *vecheat = new Vecheat(vecheatSpace);
225
226
227
228
  vecheat->initialize(INIT_ALL);

  // create adapt info
  AdaptInfo *adaptInfo = 
229
    new AdaptInfo("vecheat->adapt",
230
231
232
233
		  vecheatSpace->getNumComponents());

  // create initial adapt info
  AdaptInfo *adaptInfoInitial = 
234
    new AdaptInfo("vecheat->initial->adapt", 
235
236
237
238
		  vecheatSpace->getNumComponents());

  // create instationary adapt
  AdaptInstationary *adaptInstat =
239
    new AdaptInstationary("vecheat->adapt",
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
			  vecheatSpace,
			  adaptInfo,
			  vecheat,
			  adaptInfoInitial);

  double fac = (*(vecheat->getThetaPtr())) != 0 ? 1.0 : 1.0e-3;
  int nRefine = 0, dim = vecheatSpace->getMesh(0)->getDim();;
  GET_PARAMETER(0, vecheatSpace->getMesh(0)->getName() 
		+ "->global refinements", "%d", &nRefine);
  if((*(vecheat->getThetaPtr())) == 0.5) {
    (*(adaptInfo->getTimestepPtr())) *= 
      fac * pow(2.0, ((double) -nRefine)/dim);
  } else {
    (*(adaptInfo->getTimestepPtr())) *= 
      fac * pow(2.0, -nRefine);
  }


  // ===== create boundary functions =====
259
  G *boundaryFct = new G;
260
261
262
  boundaryFct->setTimePtr(vecheat->getBoundaryTimePtr());
  vecheat->setBoundaryFct(boundaryFct);

263
264
  vecheatSpace->addDirichletBC(DIRICHLET, 0, 0, boundaryFct);
  vecheatSpace->addDirichletBC(DIRICHLET, 1, 1, boundaryFct);
265
266

  // ===== create rhs functions =====
267
  F *rhsFct0 = new F(vecheatSpace->getFESpace(0)->getBasisFcts()->getDegree());
268
  rhsFct0->setTimePtr(vecheat->getRHSTimePtr());
269
  F *rhsFct1 = new F(vecheatSpace->getFESpace(1)->getBasisFcts()->getDegree());
270
271
272
273
274
275
276
  rhsFct1->setTimePtr(vecheat->getRHSTimePtr());

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

  // create laplace
277
  Operator *A00 = new Operator(Operator::MATRIX_OPERATOR | 
278
279
280
281
282
283
284
			       Operator::VECTOR_OPERATOR,
			       vecheatSpace->getFESpace(0),
			       vecheatSpace->getFESpace(0));

  A00->addSecondOrderTerm(new Laplace_SOT);
  A00->setUhOld(vecheat->getOldSolution()->getDOFVector(0));

285
  Operator *A11 = new Operator(Operator::MATRIX_OPERATOR |
286
287
288
289
290
291
292
			       Operator::VECTOR_OPERATOR,
			       vecheatSpace->getFESpace(1),
			       vecheatSpace->getFESpace(1));

  A11->addSecondOrderTerm(new Laplace_SOT);
  A11->setUhOld(vecheat->getOldSolution()->getDOFVector(1));

293
  if ((*(vecheat->getThetaPtr())) != 0.0) {
294
295
296
297
298
299
300
301
    vecheatSpace->addMatrixOperator(A00, 0, 0, 
				      vecheat->getThetaPtr(), 
				      &one);
    vecheatSpace->addMatrixOperator(A11, 1, 1, 
				      vecheat->getThetaPtr(), 
				      &one);
  }

302
303
304
  if ((*(vecheat->getTheta1Ptr())) != 0.0) {
    vecheatSpace->addVectorOperator(A00, 0, vecheat->getTheta1Ptr(), &zero);
    vecheatSpace->addVectorOperator(A11, 1, vecheat->getTheta1Ptr(), &zero);
305
306
307
  }

  // create zero order operator
308
  Operator *C00 = new Operator(Operator::MATRIX_OPERATOR | 
309
310
311
312
			       Operator::VECTOR_OPERATOR,
			       vecheatSpace->getFESpace(0),
			       vecheatSpace->getFESpace(0));

313
  C00->addZeroOrderTerm(new Simple_ZOT);
314
315
316
317
318
319
320
321
322
323
324

  C00->setUhOld(vecheat->getOldSolution()->getDOFVector(0));

  vecheatSpace->addMatrixOperator(C00, 0, 0, 
				  vecheat->getTau1Ptr(), 
				  vecheat->getTau1Ptr());

  vecheatSpace->addVectorOperator(C00, 0,
				  vecheat->getTau1Ptr(), 
				  vecheat->getTau1Ptr());

325
  Operator *C11 = new Operator(Operator::MATRIX_OPERATOR | 
326
327
328
329
			       Operator::VECTOR_OPERATOR,
			       vecheatSpace->getFESpace(1),
			       vecheatSpace->getFESpace(1));

330
  C11->addZeroOrderTerm(new Simple_ZOT);
331
332
333
334
335
336
337
338
339
340
341
342

  C11->setUhOld(vecheat->getOldSolution()->getDOFVector(1));

  vecheatSpace->addMatrixOperator(C11, 1, 1, 
				    vecheat->getTau1Ptr(), 
				    vecheat->getTau1Ptr());

  vecheatSpace->addVectorOperator(C11, 1,
				    vecheat->getTau1Ptr(), 
				    vecheat->getTau1Ptr());

  // create RHS operator
343
  Operator *F0 = new Operator(Operator::VECTOR_OPERATOR,
344
345
			      vecheatSpace->getFESpace(0));

346
  F0->addZeroOrderTerm(new CoordsAtQP_ZOT(rhsFct0));
347
348
349

  vecheatSpace->addVectorOperator(F0, 0);

350
  Operator *F1 = new Operator(Operator::VECTOR_OPERATOR,
351
352
			      vecheatSpace->getFESpace(1));

353
  F1->addZeroOrderTerm(new CoordsAtQP_ZOT(rhsFct1));
354
355
356
357
358
359
360
361

  vecheatSpace->addVectorOperator(F1, 1);

  // ===== start adaption loop =====
  int errorCode = adaptInstat->adapt();

  return errorCode;
}