vecheat.cc 9.4 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
  Vecheat(ProblemVec &vecheatSpace) 
52
53
54
55
56
57
58
59
60
61
62
63
    : 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
221
  ProblemVec vecheatSpace("vecheat->space");
  vecheatSpace.initialize(INIT_ALL);
222
223

  // ===== create instationary problem =====
224
225
  Vecheat vecheat(vecheatSpace);
  vecheat.initialize(INIT_ALL);
226
227

  // create adapt info
228
  AdaptInfo adaptInfo("vecheat->adapt", vecheatSpace.getNumComponents());
229
230

  // create initial adapt info
231
  AdaptInfo adaptInfoInitial("vecheat->initial->adapt", vecheatSpace.getNumComponents());
232
233

  // create instationary adapt
234
235
236
237
238
239
240
241
242
243
  AdaptInstationary adaptInstat("vecheat->adapt",
				vecheatSpace,
				adaptInfo,
				vecheat,
				adaptInfoInitial);

  double fac = *(vecheat.getThetaPtr()) != 0.0 ? 1.0 : 1.0e-3;
  int nRefine = 0;
  int dim = vecheatSpace.getMesh(0)->getDim();
  GET_PARAMETER(0, vecheatSpace.getMesh(0)->getName() 
244
		+ "->global refinements", "%d", &nRefine);
245
246
247
248
249
  if (*(vecheat.getThetaPtr()) == 0.5)
    *(adaptInfo.getTimestepPtr()) *= fac * pow(2.0, static_cast<double>(-nRefine) / dim);
  else
    *(adaptInfo.getTimestepPtr()) *= fac * pow(2.0, -nRefine);
  
250
251

  // ===== create rhs functions =====
252
  F *rhsFct0 = new F(vecheatSpace.getFeSpace(0)->getBasisFcts()->getDegree());
253
  rhsFct0->setTimePtr(vecheat.getRHSTimePtr());
254
  F *rhsFct1 = new F(vecheatSpace.getFeSpace(1)->getBasisFcts()->getDegree());
255
  rhsFct1->setTimePtr(vecheat.getRHSTimePtr());
256
257
258
259
260
261

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

  // create laplace
262
  Operator A00 (Operator::MATRIX_OPERATOR | Operator::VECTOR_OPERATOR,
263
		vecheatSpace.getFeSpace(0), vecheatSpace.getFeSpace(0));
264
265
266
267
  A00.addSecondOrderTerm(new Laplace_SOT);
  A00.setUhOld(vecheat.getOldSolution()->getDOFVector(0));

  Operator A11(Operator::MATRIX_OPERATOR | Operator::VECTOR_OPERATOR,
268
	       vecheatSpace.getFeSpace(1), vecheatSpace.getFeSpace(1));
269
270
271
272
273
274
  A11.addSecondOrderTerm(new Laplace_SOT);
  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);
275
276
  }

277
278
279
  if (*(vecheat.getTheta1Ptr()) != 0.0) {
    vecheatSpace.addVectorOperator(A00, 0, vecheat.getTheta1Ptr(), &zero);
    vecheatSpace.addVectorOperator(A11, 1, vecheat.getTheta1Ptr(), &zero);
280
281
282
  }

  // create zero order operator
283
  Operator C00(Operator::MATRIX_OPERATOR | Operator::VECTOR_OPERATOR,
284
	       vecheatSpace.getFeSpace(0), vecheatSpace.getFeSpace(0));
285
286
287
288
289
290
291
292
  C00.addZeroOrderTerm(new Simple_ZOT);
  C00.setUhOld(vecheat.getOldSolution()->getDOFVector(0));
  vecheatSpace.addMatrixOperator(C00, 0, 0, 
				 vecheat.getTau1Ptr(), vecheat.getTau1Ptr());
  vecheatSpace.addVectorOperator(C00, 0, vecheat.getTau1Ptr(), vecheat.getTau1Ptr());
  
  
  Operator C11(Operator::MATRIX_OPERATOR | Operator::VECTOR_OPERATOR,
293
	       vecheatSpace.getFeSpace(1), vecheatSpace.getFeSpace(1));
294
295
296
297
298
299
  C11.addZeroOrderTerm(new Simple_ZOT);
  C11.setUhOld(vecheat.getOldSolution()->getDOFVector(1));
  vecheatSpace.addMatrixOperator(C11, 1, 1, 			    
				 vecheat.getTau1Ptr(), vecheat.getTau1Ptr());
  vecheatSpace.addVectorOperator(C11, 1,			    
				 vecheat.getTau1Ptr(), vecheat.getTau1Ptr());
300
301

  // create RHS operator
302
  Operator F0(Operator::VECTOR_OPERATOR, vecheatSpace.getFeSpace(0));
303
304
  F0.addZeroOrderTerm(new CoordsAtQP_ZOT(rhsFct0));
  vecheatSpace.addVectorOperator(F0, 0);
305

306
  Operator F1 = Operator(Operator::VECTOR_OPERATOR, vecheatSpace.getFeSpace(1));
307
308
  F1.addZeroOrderTerm(new CoordsAtQP_ZOT(rhsFct1));
  vecheatSpace.addVectorOperator(F1, 1);
309

310
311
312
313
314
315
316
317
  // ===== create boundary functions =====
  G *boundaryFct = new G;
  boundaryFct->setTimePtr(vecheat.getBoundaryTimePtr());
  vecheat.setBoundaryFct(boundaryFct);

  vecheatSpace.addDirichletBC(DIRICHLET, 0, 0, boundaryFct);
  vecheatSpace.addDirichletBC(DIRICHLET, 1, 1, boundaryFct);

318
  // ===== start adaption loop =====
319
  int errorCode = adaptInstat.adapt();
320
321
322

  return errorCode;
}