~ubuntu-branches/ubuntu/maverick/dolfin/maverick

« back to all changes in this revision

Viewing changes to sandbox/mtl4/poisson_solve/main.cpp

  • Committer: Bazaar Package Importer
  • Author(s): Johannes Ring
  • Date: 2008-09-16 08:41:20 UTC
  • Revision ID: james.westby@ubuntu.com-20080916084120-i8k3u6lhx3mw3py3
Tags: upstream-0.9.2
ImportĀ upstreamĀ versionĀ 0.9.2

Show diffs side-by-side

added added

removed removed

Lines of Context:
 
1
#include <dolfin.h>
 
2
#include "Poisson.h"
 
3
 
 
4
using namespace dolfin;
 
5
 
 
6
int main()
 
7
{
 
8
  // Create mesh and forms
 
9
  UnitSquare mesh(128, 128);
 
10
  DomainBoundary boundary;
 
11
  Function g(mesh, 0.0);
 
12
  DirichletBC bc(g, mesh, boundary); 
 
13
 
 
14
  Function f(mesh, 2.0);
 
15
  PoissonBilinearForm a;
 
16
  PoissonLinearForm L(f);
 
17
 
 
18
  Table table;
 
19
 
 
20
  // Linear solver
 
21
  KrylovSolver solver(bicgstab, ilu);
 
22
  //LUSolver solver;
 
23
  
 
24
  // Assemble and solve using uBLAS
 
25
  cout << "uBLAS ------------------------------------------------------" << endl;
 
26
  Assembler ass_ublas(mesh);
 
27
  uBLASMatrix<ublas_sparse_matrix> A_ublas;
 
28
  uBLASVector b_ublas, x_ublas;
 
29
 
 
30
  tic(); 
 
31
  ass_ublas.assemble(A_ublas, a); 
 
32
  ass_ublas.assemble(b_ublas, L); 
 
33
  bc.apply(A_ublas, b_ublas, a);
 
34
  solver.solve(A_ublas, x_ublas, b_ublas);
 
35
  table("uBLAS", "solve") =  toc();
 
36
 
 
37
  //x_ublas.disp();
 
38
 
 
39
 
 
40
  // Assemble and solve using MTL4
 
41
  printf("MTL4 with specified number of nonzeros----------------------\n");
 
42
  Assembler ass_mtl4(mesh); 
 
43
  MTL4Matrix A_mtl4(A_ublas.size(0), A_ublas.size(1), 10);
 
44
  MTL4Vector b_mtl4(A_ublas.size(0)), x_mtl4;
 
45
 
 
46
  tic();
 
47
  ass_mtl4.assemble(A_mtl4, a, false); 
 
48
  ass_mtl4.assemble(b_mtl4, L, false);
 
49
  bc.apply(A_mtl4, b_mtl4, a);
 
50
  solver.solve(A_mtl4, x_mtl4, b_mtl4);
 
51
  table("MTL4", "solve") =  toc();
 
52
 
 
53
 
 
54
  //Function u;
 
55
  //u.init(mesh, x_mtl4, a, 1);
 
56
  //File file("test.pvd");
 
57
  //file << u;
 
58
  //x_mtl4.disp();
 
59
 
 
60
  table.disp();
 
61
 
 
62
  return 0;
 
63
}