← Back to branch summary

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

« back to all changes in this revision

Viewing changes to nfem/plib/adapt.cc

  • 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:
nfem/plib/adapt.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/adapt.h"
 
22
#include "rheolef/form.h"
 
23
#include "rheolef/field_component.h"
 
24
#include "rheolef/field_expr_ops.h"
 
25
 
 
26
namespace rheolef {
 
27
 
 
28
// externs:
 
29
template<class T, class M>
 
30
geo_basic<T,M> adapt_gmsh (const field_basic<T,M>&  uh, const adapt_option_type& opts);
 
31
template<class T, class M>
 
32
geo_basic<T,M> adapt_bamg (const field_basic<T,M>&  uh, const adapt_option_type& opts);
 
33
 
 
34
// -----------------------------------------
 
35
// proj
 
36
// -----------------------------------------
 
37
template<class T, class M>
 
38
field_basic<T,M>
 
39
proj (const field_basic<T,M>& uh, const std::string& approx = "P1")
 
40
{
 
41
  const space_basic<T,M>& Uh = uh.get_space();
 
42
  if (Uh.get_approx() == approx) return uh;
 
43
  space_basic<T,M> Vh (uh.get_geo(), approx, uh.valued());
 
44
  form_basic<T,M>  m  (Vh, Vh, "lumped_mass");
 
45
  form_basic<T,M>  p  (Uh, Vh, "mass");
 
46
  solver_basic<T,M> sm (m.uu());
 
47
  field_basic<T,M>  vh (Vh, 0);
 
48
  vh.set_u() = sm.solve((p*uh).u());
 
49
  return vh;
 
50
}
 
51
// -----------------------------------------
 
52
// smooth
 
53
// -----------------------------------------
 
54
template<class T, class M>
 
55
field_basic<T,M>
 
56
smooth (const field_basic<T,M>& uh, size_t n = 1) {
 
57
  const space_basic<T,M>& Xh = uh.get_space();
 
58
  check_macro (Xh.get_approx() == "P1", "smooth: expect P1 approx, but have " << Xh.get_approx());
 
59
  form_basic<T,M> m  (Xh, Xh, "mass");
 
60
  form_basic<T,M> md (Xh, Xh, "lumped_mass");
 
61
  solver_basic<T,M> smd (md.uu());
 
62
  field_basic<T,M> vh = uh;
 
63
  for (size_t i = 0; i < n; i++) {
 
64
    vh.set_u() = smd.solve ((m*vh).u());
 
65
  }
 
66
  return vh;
 
67
}
 
68
// -----------------------------------------
 
69
// hessian
 
70
// -----------------------------------------
 
71
template<class T, class M>
 
72
field_basic<T,M>
 
73
hessian (const field_basic<T,M>&  uh)
 
74
{
 
75
  // assume that uh is P1 and scalar
 
76
  const geo_basic<T,M>& omega = uh.get_geo();
 
77
  const space_basic<T,M>& Xh = uh.get_space();
 
78
  check_macro (Xh.valued() == "scalar", "hessian: unexpected "<<Xh.valued()<<"-valued field");
 
79
  check_macro (Xh.get_approx() == "P1", "hessian: unexpected "<<Xh.get_approx()<<" approximation");
 
80
  space_basic<T,M> Vh (omega, "P1", "vector");
 
81
  form_basic<T,M>  bv (Xh, Vh, "grad");
 
82
  form_basic<T,M>  mv (Vh, Vh, "lumped_mass");
 
83
  // TODO: inv_mass optimize: lumped and by components
 
84
  solver_basic<T,M> smv (mv.uu());
 
85
  field_basic<T,M> gh (Vh, 0);
 
86
  gh.set_u() = smv.solve ((bv*uh).u());
 
87
  space_basic<T,M> Th (omega, "P1", "tensor");
 
88
  form_basic<T,M>  bt (Vh, Th, "2D");
 
89
  form_basic<T,M>  mt (Th, Th, "lumped_mass");
 
90
  solver_basic<T,M> smt (mt.uu());
 
91
  field_basic<T,M> hh (Th, 0);
 
92
  hh.set_u() = smt.solve ((bt*gh).u());
 
93
  return hh;
 
94
}
 
95
//            |hessian(uh)|
 
96
//  M = ----------------------------
 
97
//      err*hcoef^2*(sup(uh)-inf(uh))
 
98
// where
 
99
//  |H| = same as H but with absolute value of eigenvalues
 
100
//      = Q*diag(|lambda_i|)*Qt
 
101
//
 
102
template<class T, class M>
 
103
field_basic<T,M>
 
104
hessian_criterion (
 
105
    const field_basic<T,M>&  uh0,
 
106
    const adapt_option_type& opts)
 
107
{
 
108
  typedef typename geo_basic<T,M>::size_type size_type;
 
109
  size_type d = uh0.get_geo().dimension();
 
110
  size_type k = uh0.get_space().degree();
 
111
  field_basic<T,M> uh = proj(uh0);
 
112
  field_basic<T,M> hh = hessian(uh);
 
113
  field_basic<T,M> mh (hh.get_space(), 0);
 
114
  field_basic<T,M> sh (uh.get_space(), 0);
 
115
  field_component_const<T,M> hh_comp [3][3];
 
116
  field_component<T,M>       mh_comp [3][3];
 
117
  for (size_type i_comp = 0; i_comp < d; i_comp++) {
 
118
  for (size_type j_comp = 0; j_comp < d; j_comp++) {
 
119
    hh_comp[i_comp][j_comp] = hh(i_comp,j_comp);
 
120
    mh_comp[i_comp][j_comp] = mh(i_comp,j_comp);
 
121
  }}
 
122
  tensor h_value;
 
123
  tensor m_value;
 
124
  tensor Q, Qt;
 
125
  point_basic<T> lambda;
 
126
  point_basic<T> h_local;
 
127
  for (size_type idof = 0, ndof = uh.ndof(); idof < ndof; idof++) {
 
128
    for (size_type i_comp = 0; i_comp < d; i_comp++) {
 
129
    for (size_type j_comp = 0; j_comp < d; j_comp++) {
 
130
      h_value(i_comp,j_comp) = hh_comp[i_comp][j_comp].dof (idof);
 
131
    }}
 
132
    const bool use_svd_when_2d = true;
 
133
    // H = Q*diag(lambda)*Q^T
 
134
    if (use_svd_when_2d && d == 2) {
 
135
      // TODO: eig when d=2 is bad...
 
136
      lambda = h_value.svd (Q, Qt, d);
 
137
    } else {
 
138
      lambda = h_value.eig (Q, d);
 
139
    }
 
140
    T h_min = std::numeric_limits<T>::max();
 
141
    for (size_type i_comp = 0; i_comp < d; i_comp++) {
 
142
      // err = c*h^2
 
143
      h_local[i_comp] = 1/sqrt(fabs(lambda[i_comp]));
 
144
      // yield value:
 
145
      h_local[i_comp] = std::max (opts.hmin, h_local[i_comp]);
 
146
      h_min = std::min (h_min, h_local[i_comp]);
 
147
    }
 
148
    // isotropic: takes h_local = h_min in all directions
 
149
    sh.dof (idof) = h_min;
 
150
 
 
151
    // anisotropic
 
152
    m_value = Q*diag(h_local)*trans(Q);
 
153
    for (size_type i_comp = 0; i_comp < d; i_comp++) {
 
154
    for (size_type j_comp = 0; j_comp < d; j_comp++) {
 
155
      mh_comp[i_comp][j_comp].dof (idof) = m_value(i_comp,j_comp);
 
156
    }}
 
157
  }
 
158
  T cut_off = 1e-5;
 
159
  T uh_scale  = std::max(cut_off, fabs(uh.max() - uh.min()));
 
160
  T factor = opts.hcoef*sqrt(uh_scale)*pow(opts.err,1./(k+1));
 
161
  mh = factor*mh;
 
162
  sh = factor*sh;
 
163
  if (opts.isotropic) {
 
164
    return smooth (sh, opts.n_smooth_metric);
 
165
  } else {
 
166
    return smooth (mh, opts.n_smooth_metric);
 
167
  }
 
168
}
 
169
// -----------------------------------------
 
170
// adapt
 
171
// -----------------------------------------
 
172
template<class T, class M>
 
173
geo_basic<T,M>
 
174
adapt (
 
175
    const field_basic<T,M>&  uh,
 
176
    const adapt_option_type& opts)
 
177
{
 
178
  size_t d = uh.get_geo().dimension();
 
179
  if (d == 2 && (opts.generator == "bamg" || opts.generator == "")) {
 
180
    // default generator is bamg when d=2:
 
181
    return adapt_bamg (uh, opts);
 
182
  } else {
 
183
    // use always gmsh when d != 2:
 
184
    return adapt_gmsh (uh, opts);
 
185
  }
 
186
}
 
187
// ----------------------------------------------------------------------------
 
188
// instanciation in library
 
189
// ----------------------------------------------------------------------------
 
190
#define _RHEOLEF_instanciation(T,M)                                                     \
 
191
template field_basic<T,M> proj (const field_basic<T,M>&, const std::string&);           \
 
192
template field_basic<T,M> smooth (const field_basic<T,M>&, size_t);                     \
 
193
template field_basic<T,M> hessian (const field_basic<T,M>&);                            \
 
194
template field_basic<T,M> hessian_criterion (const field_basic<T,M>&, const adapt_option_type&); \
 
195
template geo_basic<T,M> adapt (const field_basic<T,M>&, const adapt_option_type&);
 
196
 
 
197
_RHEOLEF_instanciation(Float,sequential)
 
198
#ifdef _RHEOLEF_HAVE_MPI
 
199
_RHEOLEF_instanciation(Float,distributed)
 
200
#endif // _RHEOLEF_HAVE_MPI
 
201
 
 
202
} // namespace rheolef