← Back to branch summary

~ubuntu-branches/ubuntu/trusty/rheolef/trusty-proposed

« back to all changes in this revision

Viewing changes to doc/usrman/embankment-adapt-2d.cc

  • Committer: Bazaar Package Importer
  • Author(s): Christophe Prud'homme
  • Date: 2010-06-12 09:08:59 UTC
  • Revision ID: james.westby@ubuntu.com-20100612090859-8gpm2gc7j3ab43et
Tags: upstream-5.89
Import upstream version 5.89

Show diffs side-by-side

added added

removed removed

Lines of Context:
doc/usrman/embankment-adapt-2d.cc
 
1
#include "rheolef.h"
 
2
#include "embankment.h"
 
3
using namespace std;
 
4
field criteria(Float lambda, const field& uh) {
 
5
  string approx = (uh.get_approx() == "P2") ? "P1d" : "P0";
 
6
  if (approx == "P0") return sqrt(sqr(uh[0])+sqr(uh[1]));
 
7
  space Th  (uh.get_geo(), approx, "tensor");
 
8
  space Xh  (uh.get_geo(), approx);
 
9
  form two_D  (uh.get_space(), Th, "2D");
 
10
  form div    (uh.get_space(), Xh, "div");
 
11
  form mt (Th, Th, "mass");
 
12
  form m  (Xh, Xh, "mass");
 
13
  form inv_mt (Th, Th, "inv_mass");
 
14
  form inv_m  (Xh, Xh, "inv_mass");
 
15
  field qh      = inv_m*(div*uh);
 
16
  field two_Duh = inv_mt*(two_D*uh);
 
17
  return sqrt(lambda*sqr(qh) +   sqr(field(two_Duh(0,0)))
 
18
         + sqr(field(two_Duh(1,1))) + 2*sqr(field(two_Duh(0,1))));
 
19
}
 
20
field solve(const geo& omega, const string& approx, Float lambda) {
 
21
  space Vh = embankment_space(omega, approx);
 
22
  field uh (Vh, 0.0), fh (Vh, 0.0);
 
23
  fh[1] = -1.0;
 
24
  form m  (Vh, Vh, "mass");
 
25
  form a1 (Vh, Vh, "div_div");
 
26
  form a2 (Vh, Vh, "2D_D");
 
27
  form a = lambda*a1 + a2;
 
28
  ssk<Float> fact = ldlt(a.uu);
 
29
  uh.u = fact.solve (m.uu*fh.u + m.ub*fh.b - a.ub*uh.b);
 
30
  return uh;
 
31
}
 
32
int main(int argc, char**argv) {
 
33
  const Float lambda = 1;
 
34
  geo  omega_h (argv[1]);
 
35
  string approx  = (argc > 2) ? argv[2] : "P1";
 
36
  size_t n_adapt = (argc > 3) ? atoi(argv[3]) : 5;
 
37
  adapt_option_type options;
 
38
  options.hcoef = 1.1;
 
39
  options.hmin  = 0.004;
 
40
  for (size_t i = 0; true; i++) {
 
41
      field uh = solve(omega_h, approx, lambda);
 
42
      orheostream o (omega_h.name(), "mfield");
 
43
      o << catchmark("lambda")  << lambda << endl
 
44
        << catchmark("u")       << uh;
 
45
      if (i == n_adapt) break;
 
46
      field ch = criteria(lambda,uh);
 
47
      omega_h  = geo_adapt(ch, options);
 
48
      omega_h.save();
 
49
  }
 
50
}