← Back to branch summary

~ubuntu-branches/ubuntu/saucy/rheolef/saucy

« back to all changes in this revision

Viewing changes to nfem/ptst/dirichlet.cc

  • Committer: Bazaar Package Importer
  • Author(s): Pierre Saramito
  • Date: 2011-03-23 11:14:26 UTC
  • mfrom: (1.1.4 upstream)
  • Revision ID: james.westby@ubuntu.com-20110323111426-cjvhey7lxt6077ty
Tags: 5.93-1
* New upstream release (minor changes):
  - some extra warning message deleted in heap_allocator
  - graphic output with mayavi2 fixed
  - add doc refman in .info and .pdf format

Show diffs side-by-side

added added

removed removed

Lines of Context:
nfem/ptst/dirichlet.cc
1
 
///
2
 
/// This file is part of Rheolef.
3
 
///
4
 
/// Copyright (C) 2000-2009 Pierre Saramito <Pierre.Saramito@imag.fr>
5
 
///
6
 
/// Rheolef is free software; you can redistribute it and/or modify
7
 
/// it under the terms of the GNU General Public License as published by
8
 
/// the Free Software Foundation; either version 2 of the License, or
9
 
/// (at your option) any later version.
10
 
///
11
 
/// Rheolef is distributed in the hope that it will be useful,
12
 
/// but WITHOUT ANY WARRANTY; without even the implied warranty of
13
 
/// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
14
 
/// GNU General Public License for more details.
15
 
///
16
 
/// You should have received a copy of the GNU General Public License
17
 
/// along with Rheolef; if not, write to the Free Software
18
 
/// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
19
 
///
20
 
/// =========================================================================
21
 
#include "rheolef.h"
22
 
#include "rounder.h"
23
 
using namespace rheolef;
24
 
using namespace std;
25
 
int main(int argc, char**argv) {
26
 
  environment parallel(argc, argv);
27
 
  geo omega (argv[1]);
28
 
  Float tol  = (argc > 2) ? atof(argv[2]) : 1e-15;
29
 
  Float prec = (argc > 3) ? atof(argv[3]) : 1e-10;
30
 
  space Xh (omega, "P1");
31
 
  if (omega.dimension() >= 3) {
32
 
    Xh.block("back");
33
 
    Xh.block("front");
34
 
  }
35
 
  if (omega.dimension() >= 2) {
36
 
    Xh.block("top");
37
 
    Xh.block("bottom");
38
 
  }
39
 
  Xh.block("left");
40
 
  Xh.block("right");
41
 
  form a (Xh, Xh, "grad_grad");
42
 
  form m (Xh, Xh, "mass");
43
 
  field fh (Xh, 1.0);
44
 
  field uh (Xh, 0.0);
45
 
#ifdef TODO
46
 
  uh ["boundary"] = 0;
47
 
#endif // TODO
48
 
#ifdef TODO
49
 
  uh.u() = 0;
50
 
#endif // TODO
51
 
  size_t max_iter = 10000;
52
 
  int status = pcg (a.uu(), uh.u(), m.uu()*fh.u() + m.ub()*fh.b() - a.ub()*uh.b(),
53
 
                  eye<Float>(), max_iter, tol, &pcerr);
54
 
  pcout << setprecision(15) << compose (rounder(prec),uh);
55
 
  pcerr << "uh.u.par_size = " << uh.u().par_size() << endl
56
 
        << uh.u();
57
 
  pcerr << "uh.b.par_size = " << uh.b().par_size() << endl
58
 
        << uh.b();
59
 
}