← Back to branch summary

~njansson/dolfin/hpc

« back to all changes in this revision

Viewing changes to bench/fem/convergence/main.cpp

  • Committer: Johannes Ring
  • Date: 2008-03-05 22:43:06 UTC
  • Revision ID: johannr@simula.no-20080305224306-2npsdyhfdpl2esji
The BIG commit!

Show diffs side-by-side

added added

removed removed

Lines of Context:
bench/fem/convergence/main.cpp
 
1
// Copyright (C) 2005-2007 Anders Logg.
 
2
// Licensed under the GNU LGPL Version 2.1.
 
3
//
 
4
// First added:  2005
 
5
// Last changed: 2007-08-20
 
6
 
 
7
#include <dolfin.h>
 
8
#include "Poisson2D_1.h"
 
9
#include "Poisson2D_2.h"
 
10
#include "Poisson2D_3.h"
 
11
#include "Poisson2D_4.h"
 
12
#include "Poisson2D_5.h"
 
13
#include "Poisson3D_1.h"
 
14
#include "Poisson3D_2.h"
 
15
#include "Poisson3D_3.h"
 
16
#include "Poisson3D_4.h"
 
17
#include "Poisson3D_5.h"
 
18
 
 
19
using namespace dolfin;
 
20
 
 
21
// Boundary condition
 
22
class DirichletBoundary : public SubDomain
 
23
{
 
24
  bool inside(const real* x, bool on_boundary) const
 
25
  {
 
26
    return on_boundary;
 
27
  }
 
28
};
 
29
 
 
30
// Right-hand side, 2D
 
31
class Source2D : public Function
 
32
{
 
33
public:
 
34
 
 
35
  Source2D(Mesh& mesh) : Function(mesh) {}
 
36
 
 
37
  real eval(const real* x) const
 
38
  {
 
39
    return 2.0*DOLFIN_PI*DOLFIN_PI*sin(DOLFIN_PI*x[0])*sin(DOLFIN_PI*x[1]);
 
40
  }
 
41
 
 
42
};
 
43
 
 
44
// Right-hand side, 3D
 
45
class Source3D : public Function
 
46
{
 
47
public:
 
48
  
 
49
  Source3D(Mesh& mesh) : Function(mesh) {}
 
50
 
 
51
  real eval(const real* x) const
 
52
  {
 
53
    return 3.0*DOLFIN_PI*DOLFIN_PI*sin(DOLFIN_PI*x[0])*sin(DOLFIN_PI*x[1])*sin(DOLFIN_PI*x[2]);
 
54
  }
 
55
 
 
56
};
 
57
 
 
58
// Solve equation and compute error, 2D
 
59
real solve2D(int q, int n)
 
60
{
 
61
  message("--------------------------------------------------");
 
62
  message("Solving Poisson's equation in 2D for q = %d, n = %d.", q, n);
 
63
  message("--------------------------------------------------");
 
64
 
 
65
  // Set up problem
 
66
  UnitSquare mesh(n, n);
 
67
  Source2D f(mesh);
 
68
  Function zero(mesh, 0.0);
 
69
  DirichletBoundary boundary;
 
70
  DirichletBC bc(zero, mesh, boundary);
 
71
 
 
72
  // Choose forms
 
73
  Form* a = 0;
 
74
  Form* L = 0;
 
75
  switch (q)
 
76
  {
 
77
  case 1:
 
78
    a = new Poisson2D_1BilinearForm();
 
79
    L = new Poisson2D_1LinearForm(f);
 
80
    break;
 
81
  case 2:
 
82
    a = new Poisson2D_2BilinearForm();
 
83
    L = new Poisson2D_2LinearForm(f);
 
84
    break;
 
85
  case 3:
 
86
    a = new Poisson2D_3BilinearForm();
 
87
    L = new Poisson2D_3LinearForm(f);
 
88
    break;
 
89
  case 4:
 
90
    a = new Poisson2D_4BilinearForm();
 
91
    L = new Poisson2D_4LinearForm(f);
 
92
    break;
 
93
  case 5:
 
94
    a = new Poisson2D_5BilinearForm();
 
95
    L = new Poisson2D_5LinearForm(f);
 
96
    break;
 
97
  default:
 
98
    error("Forms not compiled for q = %d.", q);
 
99
  }    
 
100
 
 
101
  // Discretize equation
 
102
  Matrix A;
 
103
  Vector x, b;
 
104
  assemble(A, *a, mesh);
 
105
  assemble(b, *L, mesh);
 
106
  bc.apply(A, b, *a);
 
107
 
 
108
  // Solve the linear system
 
109
  KrylovSolver solver(gmres);
 
110
  solver.set("Krylov relative tolerance", 1e-14); 
 
111
  solver.solve(A, x, b);
 
112
 
 
113
  // Compute maximum norm of error
 
114
  real emax = 0.0;
 
115
  real* U = new real[x.size()];
 
116
  x.get(U);
 
117
  for (VertexIterator v(mesh); !v.end(); ++v)
 
118
  {
 
119
    const Point p = v->point();
 
120
    const real u = sin(DOLFIN_PI*p.x())*sin(DOLFIN_PI*p.y());
 
121
    const real e = std::abs(U[v->index()] - u);
 
122
    emax = std::max(emax, e);
 
123
  }
 
124
  delete [] U;
 
125
 
 
126
  delete a;
 
127
  delete L;
 
128
 
 
129
  return emax;
 
130
}
 
131
 
 
132
// Solve equation and compute error, 3D
 
133
real solve3D(int q, int n)
 
134
{
 
135
  message("--------------------------------------------------");
 
136
  message("Solving Poisson's equation in 3D for q = %d, n = %d.", q, n);
 
137
  message("--------------------------------------------------");
 
138
 
 
139
  // Set up problem
 
140
  UnitCube mesh(n, n, n);
 
141
  Source3D f(mesh);
 
142
  Function zero(mesh, 0.0);
 
143
  DirichletBoundary boundary;
 
144
  DirichletBC bc(zero, mesh, boundary);
 
145
 
 
146
  // Choose forms
 
147
  Form* a = 0;
 
148
  Form* L = 0;
 
149
  switch (q)
 
150
  {
 
151
  case 1:
 
152
    a = new Poisson3D_1BilinearForm();
 
153
    L = new Poisson3D_1LinearForm(f);
 
154
    break;
 
155
  case 2:
 
156
    a = new Poisson3D_2BilinearForm();
 
157
    L = new Poisson3D_2LinearForm(f);
 
158
    break;
 
159
  case 3:
 
160
    a = new Poisson3D_3BilinearForm();
 
161
    L = new Poisson3D_3LinearForm(f);
 
162
    break;
 
163
  case 4:
 
164
    a = new Poisson3D_4BilinearForm();
 
165
    L = new Poisson3D_4LinearForm(f);
 
166
    break;
 
167
  case 5:
 
168
    a = new Poisson3D_5BilinearForm();
 
169
    L = new Poisson3D_5LinearForm(f);
 
170
    break;
 
171
  default:
 
172
    error("Forms not compiled for q = %d.", q);
 
173
  }    
 
174
 
 
175
  // Discretize equation
 
176
  Matrix A;
 
177
  Vector x, b;
 
178
  assemble(A, *a, mesh);
 
179
  assemble(b, *L, mesh);
 
180
  bc.apply(A, b, *a);
 
181
 
 
182
  // Solve the linear system
 
183
  KrylovSolver solver(gmres);
 
184
  solver.set("Krylov relative tolerance", 1e-14); 
 
185
  solver.solve(A, x, b);
 
186
 
 
187
  // Compute maximum norm of error
 
188
  real emax = 0.0;
 
189
  real* U = new real[x.size()];
 
190
  x.get(U);
 
191
  for (VertexIterator v(mesh); !v.end(); ++v)
 
192
  {
 
193
    const Point p = v->point();
 
194
    const real u = sin(DOLFIN_PI*p.x())*sin(DOLFIN_PI*p.y())*sin(DOLFIN_PI*p.z());
 
195
    const real e = std::abs(U[v->index()] - u);
 
196
    emax = std::max(emax, e);
 
197
  }
 
198
  delete [] U;
 
199
 
 
200
  delete a;
 
201
  delete L;
 
202
 
 
203
  return emax;
 
204
}
 
205
 
 
206
int main()
 
207
{
 
208
  const int qmax = 5;
 
209
  const int num_meshes = 3;
 
210
  real e2D[qmax][num_meshes];
 
211
  real e3D[qmax][num_meshes];
 
212
 
 
213
  // Compute errors in 2D
 
214
  for (int q = 1; q <= qmax; q++)
 
215
  {
 
216
    int n = 2;
 
217
    for (int i = 0; i < num_meshes; i++)
 
218
    {
 
219
      e2D[q - 1][i] = solve2D(q, n);
 
220
      n *= 2;
 
221
    }
 
222
  }
 
223
 
 
224
  // Compute errors in 3D
 
225
  for (int q = 1; q <= qmax; q++)
 
226
  {
 
227
    int n = 2;
 
228
    for (int i = 0; i < num_meshes; i++)
 
229
    {
 
230
      e3D[q - 1][i] = solve3D(q, n);
 
231
      n *= 2;
 
232
    }
 
233
  }
 
234
 
 
235
  // Write errors in 2D
 
236
  printf("\nMaximum norm error in 2D:\n");
 
237
  printf("-------------------------\n");
 
238
  for (int q = 1; q <= qmax; q++)
 
239
  {
 
240
    printf("q = %d:", q);
 
241
    for (int i = 0; i < num_meshes; i++)
 
242
      printf(" %.3e", e2D[q - 1][i]);
 
243
    printf("\n");
 
244
  }
 
245
 
 
246
  // Write errors in 3D
 
247
  printf("\nMaximum norm error in 3D:\n");
 
248
  printf("-------------------------\n");
 
249
  for (int q = 1; q <= qmax; q++)
 
250
  {
 
251
    printf("q = %d:", q);
 
252
    for (int i = 0; i < num_meshes; i++)
 
253
      printf(" %.3e", e3D[q - 1][i]);
 
254
    printf("\n");
 
255
  }
 
256
}