← Back to branch summary

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

« back to all changes in this revision

Viewing changes to doc/pexamples/p_laplacian1.icc

  • Committer: Package Import Robot
  • Author(s): Pierre Saramito
  • Date: 2012-04-06 09:12:21 UTC
  • mfrom: (1.1.5)
  • Revision ID: package-import@ubuntu.com-20120406091221-m58me99p1nxqui49
Tags: 6.0-1
http://bugs.debian.org/667356
http://bugs.debian.org/661689
* New upstream release 6.0 (major changes):
  - massively distributed and parallel support
  - full FEM characteristic method (Lagrange-Gakerkin method) support
  - enhanced users documentation 
  - source code supports g++-4.7 (closes: #667356)
* debian/control: dependencies for MPI distributed solvers added
* debian/rules: build commands simplified
* debian/librheolef-dev.install: man1/* to man9/* added
* debian/changelog: package description rewritted (closes: #661689)

Show diffs side-by-side

added added

removed removed

Lines of Context:
doc/pexamples/p_laplacian1.icc
 
1
#include "dirichlet.icc"
 
2
#include "not_too_small.h"
 
3
p_laplacian::p_laplacian(Float p1, const geo& omega_h, string approx1)
 
4
 : p(p1), Xh(), Kh(), fh(),
 
5
  m(), inv_mt(), grad(), sm(), a1(), sa1() {
 
6
  reset(omega_h, approx1);
 
7
}
 
8
void p_laplacian::reset(const geo& omega_h1, string approx1) {
 
9
  if (approx1 == "previous") approx1 = Xh.get_approx();
 
10
  Xh = space(omega_h1, approx1);
 
11
  Xh.block ("boundary");
 
12
  fh = field(Xh, 1);
 
13
  m = form (Xh, Xh, "mass");
 
14
  sm = solver(m.uu());
 
15
  string grad_approx = "P" + itos(Xh.degree()-1) + "d";
 
16
  space Th (fh.get_geo(), grad_approx, "vector");
 
17
  inv_mt = form (Th, Th, "inv_mass");
 
18
  grad = form (Xh, Th, "grad");
 
19
  Kh = space(fh.get_geo(), grad_approx, "tensor");
 
20
}
 
21
field p_laplacian::initial () const {
 
22
  field uh(Xh);
 
23
  uh [Xh.get_geo()["boundary"]] = 0;
 
24
  dirichlet (fh, uh);
 
25
  return uh;
 
26
}
 
27
void p_laplacian::update_derivative (const field& uh) const {
 
28
  field grad_uh = inv_mt*(grad*uh); 
 
29
  field norm2_grad_uh = norm2(grad_uh);
 
30
  if (p/2 - 4 <= 0) norm2_grad_uh = compose (not_too_small(1e-10), norm2_grad_uh);
 
31
  field w0h = pow(norm2_grad_uh, p/2)/norm2_grad_uh;
 
32
  field w1h = pow(norm2_grad_uh, p/2)/sqr(norm2_grad_uh);
 
33
  field eta_h (Kh);
 
34
  eta_h(0,0) = w0h + (p-2)*w1h*sqr(grad_uh[0]);
 
35
  size_t d = uh.get_geo().dimension();
 
36
  if (d >= 2) {
 
37
    eta_h(1,1) = w0h + (p-2)*w1h*sqr(grad_uh[1]);
 
38
    eta_h(0,1) =       (p-2)*w1h*grad_uh[0]*grad_uh[1];
 
39
  }
 
40
  if (d == 3) {
 
41
    eta_h(2,2) = w0h + (p-2)*w1h*sqr(grad_uh[2]);
 
42
    eta_h(1,2) =       (p-2)*w1h*grad_uh[1]*grad_uh[2];
 
43
    eta_h(0,2) =       (p-2)*w1h*grad_uh[0]*grad_uh[2];
 
44
  }
 
45
  a1 = form (Xh, Xh, "grad_grad", eta_h);
 
46
  sa1 = solver(a1.uu());
 
47
}
 
48
field p_laplacian::residue (const field& uh) const {
 
49
  field grad_uh = inv_mt*(grad*uh); 
 
50
  field norm2_grad_uh = norm2(grad_uh);
 
51
  field w0h = pow(norm2_grad_uh, p/2-1);
 
52
  form a (Xh, Xh, "grad_grad", w0h);
 
53
  field mrh = a*uh - m*fh;
 
54
  mrh.set_b() = 0;
 
55
  return mrh;
 
56
}
 
57
field p_laplacian::derivative_solve (const field& mrh) const {
 
58
  field delta_uh (Xh,0);
 
59
  delta_uh.set_b() = 0;
 
60
  delta_uh.set_u() = sa1.solve(mrh.u());
 
61
  return delta_uh;
 
62
}