heat.cc 6.41 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
class Heat : public ProblemInstatScal
{
public:
49
  Heat(ProblemScal &heatSpace)
50
51
52
53
54
55
56
57
58
59
60
61
    : 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
  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);

172

173
  // ===== create and init stationary problem =====
174
175
176
  ProblemScal heatSpace("heat->space");
  heatSpace.initialize(INIT_ALL);

177
178

  // ===== create instationary problem =====
179
180
  Heat heat(heatSpace);
  heat.initialize(INIT_ALL);
181
182

  // create adapt info
183
  AdaptInfo adaptInfo("heat->adapt");
184
185

  // create initial adapt info
186
  AdaptInfo adaptInfoInitial("heat->initial->adapt");
187
188

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

195
196

  // ===== create rhs functions =====
197
  int degree = heatSpace.getFeSpace()->getBasisFcts()->getDegree();
198
  F *rhsFct = new F(degree);
199
200
  rhsFct->setTimePtr(heat.getRHSTimePtr());

201
202
203
204
205
206

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

  // create laplace
207
  Operator A(Operator::MATRIX_OPERATOR | Operator::VECTOR_OPERATOR, 
208
	     heatSpace.getFeSpace());
209
210
211
212
213
214
  A.addSecondOrderTerm(new Laplace_SOT);
  A.setUhOld(heat.getOldSolution());
  if (*(heat.getThetaPtr()) != 0.0)
    heatSpace.addMatrixOperator(A, heat.getThetaPtr(), &one);
  if (*(heat.getTheta1Ptr()) != 0.0)
    heatSpace.addVectorOperator(A, heat.getTheta1Ptr(), &zero);
215
216

  // create zero order operator
217
  Operator C(Operator::MATRIX_OPERATOR | Operator::VECTOR_OPERATOR,
218
	     heatSpace.getFeSpace());
219
220
221
222
  C.addZeroOrderTerm(new Simple_ZOT);
  C.setUhOld(heat.getOldSolution());
  heatSpace.addMatrixOperator(C, heat.getTau1Ptr(), heat.getTau1Ptr());
  heatSpace.addVectorOperator(C, heat.getTau1Ptr(), heat.getTau1Ptr());
223
224

  // create RHS operator
225
  Operator F(Operator::VECTOR_OPERATOR, heatSpace.getFeSpace());
226
227
  F.addZeroOrderTerm(new CoordsAtQP_ZOT(rhsFct));
  heatSpace.addVectorOperator(F);
228
229


230
231
232
233
234
235
236
237
  // ===== create boundary functions =====
  G *boundaryFct = new G;
  boundaryFct->setTimePtr(heat.getBoundaryTimePtr());
  heat.setExactSolution(boundaryFct);

  heatSpace.addDirichletBC(DIRICHLET, boundaryFct);


238
  // ===== start adaption loop =====
239
  int errorCode = adaptInstat.adapt();
240
241
242

  return errorCode;
}