heat.cc 6.62 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
class G : public AbstractFunction<double, WorldVector<double> >,
	  public TimedObject
{
public:
15
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));
Thomas Witkowski's avatar
Thomas Witkowski committed
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;
Thomas Witkowski's avatar
Thomas Witkowski committed
38
  }
39
40
41
42
43
44
};

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

45
/// Instationary problem
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
class Heat : public ProblemInstatScal
{
public:
  Heat(ProblemScal *heatSpace)
    : ProblemInstatScal("heat", heatSpace)
  {
    // 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;
Thomas Witkowski's avatar
Thomas Witkowski committed
62
  }
63
64
65

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

66
67
68
  /// set the time in all needed functions!
  void setTime(AdaptInfo *adaptInfo) 
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
69
    rhsTime = adaptInfo->getTime() - (1 - theta) * adaptInfo->getTimestep();
70
71
    boundaryTime = adaptInfo->getTime();
    tau1 = 1.0 / adaptInfo->getTimestep();    
Thomas Witkowski's avatar
Thomas Witkowski committed
72
  }
73

74
75
  void closeTimestep(AdaptInfo *adaptInfo) 
  {
76
77
    ProblemInstatScal::closeTimestep(adaptInfo);
    WAIT;
Thomas Witkowski's avatar
Thomas Witkowski committed
78
  }
79
80
81

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

82
83
84
  /// Used by \ref problemInitial to solve the system of the initial problem
  void solve(AdaptInfo *adaptInfo) 
  {
85
    problemStat->getSolution()->interpol(exactSolution);
Thomas Witkowski's avatar
Thomas Witkowski committed
86
  }
87

88
89
90
  /// Used by \ref problemInitial to do error estimation for the initial problem.
  void estimate(AdaptInfo *adaptInfo) 
  {
91
92
93
94
95
96
    double errMax, errSum;

    errSum = Error<double>::L2Err(*exactSolution,
				  *(problemStat->getSolution()), 0, &errMax, false);
    adaptInfo->setEstSum(errSum, 0);
    adaptInfo->setEstMax(errMax, 0);
Thomas Witkowski's avatar
Thomas Witkowski committed
97
  }
98
99
100

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

101
102
103
  /// Sets \ref exactSolution;
  void setExactSolution(AbstractFunction<double, WorldVector<double> > *fct) 
  {
104
105
106
107
108
    exactSolution = fct;
  } 

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

109
110
111
  /// Returns pointer to \ref theta.
  double *getThetaPtr() 
  { 
112
    return &theta; 
Thomas Witkowski's avatar
Thomas Witkowski committed
113
  }
114

115
116
117
  /// Returns pointer to \ref theta1.
  double *getTheta1Ptr() 
  { 
118
    return &theta1; 
Thomas Witkowski's avatar
Thomas Witkowski committed
119
  }
120

121
122
123
  /// Returns pointer to \ref tau1
  double *getTau1Ptr() 
  { 
124
    return &tau1; 
Thomas Witkowski's avatar
Thomas Witkowski committed
125
  }
126

127
128
129
  /// Returns pointer to \ref rhsTime.
  double *getRHSTimePtr() 
  { 
130
    return &rhsTime; 
Thomas Witkowski's avatar
Thomas Witkowski committed
131
  }
132

133
134
135
  /// Returns pointer to \ref theta1.
  double *getBoundaryTimePtr() 
  {
136
    return &boundaryTime; 
Thomas Witkowski's avatar
Thomas Witkowski committed
137
  }
138
139

private:
140
  /// Used for theta scheme.
141
142
  double theta;

143
  /// theta - 1
144
145
  double theta1;

146
  /// 1.0 / timestep
147
148
  double tau1;

149
  /// time for right hand side functions.
150
151
  double rhsTime;

152
  /// time for boundary functions.
153
154
  double boundaryTime;

155
  /// Pointer to boundary function. Needed for initial problem.
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
  AbstractFunction<double, WorldVector<double> > *exactSolution;
};

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

int main(int argc, char** argv)
{
  // ===== check for init file =====
  TEST_EXIT(argc > 1)("usage: heat initfile\n");

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

  // ===== create and init stationary problem =====
173
  ProblemScal *heatSpace = new ProblemScal("heat->space");
174
175
176
177
178
179
180
  heatSpace->initialize(INIT_ALL);

  // ===== create instationary problem =====
  Heat *heat = new Heat(heatSpace);
  heat->initialize(INIT_ALL);

  // create adapt info
181
  AdaptInfo *adaptInfo = new AdaptInfo("heat->adapt");
182
183

  // create initial adapt info
184
  AdaptInfo *adaptInfoInitial = new AdaptInfo("heat->initial->adapt");
185
186
187

  // create instationary adapt
  AdaptInstationary *adaptInstat =
188
    new AdaptInstationary("heat->adapt",
189
190
191
192
193
194
			  heatSpace,
			  adaptInfo,
			  heat,
			  adaptInfoInitial);

  // ===== create boundary functions =====
195
  G *boundaryFct = new G;
196
197
198
199
200
201
202
  boundaryFct->setTimePtr(heat->getBoundaryTimePtr());
  heat->setExactSolution(boundaryFct);

  heatSpace->addDirichletBC(DIRICHLET, boundaryFct);

  // ===== create rhs functions =====
  int degree = heatSpace->getFESpace()->getBasisFcts()->getDegree();
203
  F *rhsFct = new F(degree);
204
205
206
207
208
209
210
  rhsFct->setTimePtr(heat->getRHSTimePtr());

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

  // create laplace
211
  Operator *A = new Operator(Operator::MATRIX_OPERATOR | 
212
213
214
215
216
217
218
			     Operator::VECTOR_OPERATOR,
			     heatSpace->getFESpace());

  A->addSecondOrderTerm(new Laplace_SOT);

  A->setUhOld(heat->getOldSolution());

219
  if ((*(heat->getThetaPtr())) != 0.0)
220
221
    heatSpace->addMatrixOperator(A, heat->getThetaPtr(), &one);

222
  if ((*(heat->getTheta1Ptr())) != 0.0)
223
224
225
    heatSpace->addVectorOperator(A, heat->getTheta1Ptr(), &zero);

  // create zero order operator
226
  Operator *C = new Operator(Operator::MATRIX_OPERATOR | 
227
228
229
			     Operator::VECTOR_OPERATOR,
			     heatSpace->getFESpace());

230
  C->addZeroOrderTerm(new Simple_ZOT);
231
232
233
234
235
236
237

  C->setUhOld(heat->getOldSolution());

  heatSpace->addMatrixOperator(C, heat->getTau1Ptr(), heat->getTau1Ptr());
  heatSpace->addVectorOperator(C, heat->getTau1Ptr(), heat->getTau1Ptr());

  // create RHS operator
238
  Operator *F = new Operator(Operator::VECTOR_OPERATOR,
239
240
			     heatSpace->getFESpace());

241
  F->addZeroOrderTerm(new CoordsAtQP_ZOT(rhsFct));
242
243
244
245
246
247
248
249

  heatSpace->addVectorOperator(F);

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

  return errorCode;
}