Loading...
Searching...
No Matches
adapt.cc
Go to the documentation of this file.
#include "rheolef/adapt.h"
#include "rheolef/form.h"
#include "rheolef/field_wdof_sliced.h"
#include "rheolef/field_expr.h"
namespace rheolef {
// externs:
template<class T, class M>
geo_basic<T,M> adapt_gmsh (const field_basic<T,M>& uh, const adapt_option& opts);
template<class T, class M>
geo_basic<T,M> adapt_bamg (const field_basic<T,M>& uh, const adapt_option& opts);
// -----------------------------------------
// proj
// -----------------------------------------
template<class T, class M>
field_basic<T,M>
proj (const field_basic<T,M>& uh, const std::string& approx = "P1")
{
const space_basic<T,M>& Uh = uh.get_space();
if (Uh.get_approx() == approx) return uh;
space_basic<T,M> Vh (uh.get_geo(), approx, uh.valued());
form_basic<T,M> m (Vh, Vh, "lumped_mass");
form_basic<T,M> p (Uh, Vh, "mass");
solver_basic<T,M> sm (m.uu());
field_basic<T,M> vh (Vh, 0);
vh.set_u() = sm.solve((p*uh).u());
return vh;
}
// -----------------------------------------
// smooth
// -----------------------------------------
template<class T, class M>
field_basic<T,M>
smooth (const field_basic<T,M>& uh, size_t n = 1) {
const space_basic<T,M>& Xh = uh.get_space();
check_macro (Xh.get_approx() == "P1", "smooth: expect P1 approx, but have " << Xh.get_approx());
form_basic<T,M> m (Xh, Xh, "mass");
form_basic<T,M> md (Xh, Xh, "lumped_mass");
solver_basic<T,M> smd (md.uu());
field_basic<T,M> vh = uh;
for (size_t i = 0; i < n; i++) {
vh.set_u() = smd.solve ((m*vh).u());
}
return vh;
}
// -----------------------------------------
// hessian
// -----------------------------------------
template<class T, class M>
field_basic<T,M>
hessian (const field_basic<T,M>& uh)
{
// assume that uh is P1 and scalar
const geo_basic<T,M>& omega = uh.get_geo();
const space_basic<T,M>& Xh = uh.get_space();
check_macro (Xh.valued() == "scalar", "hessian: unexpected "<<Xh.valued()<<"-valued field");
check_macro (Xh.get_approx() == "P1", "hessian: unexpected "<<Xh.get_approx()<<" approximation");
space_basic<T,M> Vh (omega, "P1", "vector");
form_basic<T,M> bv (Xh, Vh, "grad");
form_basic<T,M> mv (Vh, Vh, "lumped_mass");
// TODO: inv_mass optimize: lumped and by components
solver_basic<T,M> smv (mv.uu());
field_basic<T,M> gh (Vh, 0);
gh.set_u() = smv.solve ((bv*uh).u());
space_basic<T,M> Th (omega, "P1", "tensor");
form_basic<T,M> bt (Vh, Th, "2D");
form_basic<T,M> mt (Th, Th, "lumped_mass");
solver_basic<T,M> smt (mt.uu());
field_basic<T,M> hh (Th, 0);
hh.set_u() = smt.solve ((bt*gh).u());
return hh;
}
// |hessian(uh)|
// M = ----------------------------
// err*hcoef^2*(sup(uh)-inf(uh))
// where
// |H| = same as H but with absolute value of eigenvalues
// = Q*diag(|lambda_i|)*Qt
//
template<class T, class M>
field_basic<T,M>
hessian_criterion (
const field_basic<T,M>& uh0,
const adapt_option& opts)
{
typedef typename geo_basic<T,M>::size_type size_type;
size_type d = uh0.get_geo().dimension();
size_type k = uh0.get_space().degree();
field_basic<T,M> uh = proj(uh0);
field_basic<T,M> hh = hessian(uh);
field_basic<T,M> mh (hh.get_space(), 0);
field_basic<T,M> sh (uh.get_space(), 0);
#ifdef TO_CLEAN
details::field_rdof_sliced_const<field_basic<T,M>> hh_comp [3][3];
for (size_type i_comp = 0; i_comp < d; i_comp++) {
for (size_type j_comp = 0; j_comp < d; j_comp++) {
hh_comp[i_comp][j_comp].proxy_assign (hh(i_comp,j_comp));
}}
#endif // TO_CLEAN
tensor_basic<T> h_value, m_value, Q, Qt;
point_basic<T> lambda, h_local;
T cut_off = 1e-5;
T uh_scale = std::max(cut_off, fabs(uh.max() - uh.min()));
T factor = opts.hcoef*sqrt(uh_scale)*pow(opts.err,1./(k+1));
T eps = std::numeric_limits<T>::epsilon();
for (size_type idof = 0, ndof = uh.ndof(); idof < ndof; idof++) {
for (size_type i_comp = 0; i_comp < d; i_comp++) {
for (size_type j_comp = 0; j_comp < d; j_comp++) {
h_value(i_comp,j_comp) = hh(i_comp,j_comp).dof (idof);
}}
const bool use_svd_when_2d = true;
// H = Q*diag(lambda)*Q^T
if (use_svd_when_2d && d == 2) {
// TODO: eig when d=2 is bad...
lambda = h_value.svd (Q, Qt, d);
} else {
lambda = h_value.eig (Q, d);
}
T h_min = opts.hmax;
for (size_type i_comp = 0; i_comp < d; i_comp++) {
if (fabs(lambda[i_comp]) < eps) continue;
// err = c*h^2
h_local[i_comp] = factor/sqrt(fabs(lambda[i_comp]));
// yield values
h_local[i_comp] = std::min (opts.hmax, h_local[i_comp]);
h_local[i_comp] = std::max (opts.hmin, h_local[i_comp]);
h_min = std::min (h_min, h_local[i_comp]);
}
// isotropic: takes h_local = h_min in all directions
sh.dof (idof) = h_min;
// anisotropic
m_value = Q*diag(h_local)*trans(Q);
for (size_type i_comp = 0; i_comp < d; i_comp++) {
for (size_type j_comp = 0; j_comp < d; j_comp++) {
mh(i_comp,j_comp).dof (idof) = m_value(i_comp,j_comp);
}}
}
#ifdef TO_CLEAN
trace_macro ("sh.min="<<sh.min()<<", sh.max="<<sh.max());
mh *= factor;
sh *= factor;
trace_macro ("uh_scale="<<uh_scale<<", factor="<<uh_scale<<" => sh.min="<<sh.min()<<", sh.max="<<sh.max());
#endif // TO_CLEAN
if (opts.isotropic) {
return smooth (sh, opts.n_smooth_metric);
} else {
return smooth (mh, opts.n_smooth_metric);
}
}
// -----------------------------------------
// adapt
// -----------------------------------------
template<class T, class M>
geo_basic<T,M>
adapt (
const field_basic<T,M>& uh,
const adapt_option& opts)
{
size_t d = uh.get_geo().dimension();
if (d == 2 && (opts.generator == "bamg" || opts.generator == "")) {
// default generator is bamg when d=2:
return adapt_bamg (uh, opts);
} else {
// use always gmsh when d != 2:
return adapt_gmsh (uh, opts);
}
}
// ----------------------------------------------------------------------------
// instanciation in library
// ----------------------------------------------------------------------------
#define _RHEOLEF_instanciation(T,M) \
template field_basic<T,M> proj (const field_basic<T,M>&, const std::string&); \
template field_basic<T,M> smooth (const field_basic<T,M>&, size_t); \
template field_basic<T,M> hessian (const field_basic<T,M>&); \
template field_basic<T,M> hessian_criterion (const field_basic<T,M>&, const adapt_option&); \
template geo_basic<T,M> adapt (const field_basic<T,M>&, const adapt_option&);
_RHEOLEF_instanciation(Float,sequential)
#ifdef _RHEOLEF_HAVE_MPI
_RHEOLEF_instanciation(Float,distributed)
#endif // _RHEOLEF_HAVE_MPI
} // namespace rheolef
_RHEOLEF_instanciation
#define _RHEOLEF_instanciation(T, M, A)
Definition asr.cc:223
size_type
field::size_type size_type
Definition branch.cc:430
gh
field gh(Float epsilon, Float t, const field &uh, const test &v)
Definition burgers_diffusion_operators.icc:37
Float
see the Float page for the full documentation
rheolef::details::field_rdof_sliced_const
Definition field_wdof_sliced.h:204
rheolef::field_basic
Definition field.h:219
rheolef::field_basic::dof
T & dof(size_type idof)
Definition field.h:738
rheolef::field_basic::set_u
vec< T, M > & set_u()
Definition field.h:284
rheolef::field_basic::min
T min() const
Definition field.h:786
rheolef::field_basic::max
T max() const
Definition field.h:802
rheolef::field_basic::get_geo
const geo_type & get_geo() const
Definition field.h:271
rheolef::field_basic::get_space
const space_type & get_space() const
Definition field.h:270
rheolef::field_basic::ndof
size_type ndof() const
Definition field.h:298
rheolef::field_basic::valued
const std::string & valued() const
Definition field.h:274
rheolef::form_basic
Definition form.h:146
rheolef::form_basic::uu
const csr< T, M > & uu() const
Definition form.h:204
rheolef::geo_basic
generic mesh with rerefence counting
Definition geo.h:1089
rheolef::point_basic
Definition point.h:86
rheolef::solver_basic
Definition solver.h:264
rheolef::solver_basic::solve
vec< T, M > solve(const vec< T, M > &b) const
Definition solver.h:345
rheolef::space_basic
the finite element space
Definition space.h:382
rheolef::tensor_basic
Definition tensor.h:90
rheolef::tensor_basic::eig
point_basic< T > eig(tensor_basic< T > &q, size_t dim=3) const
Definition tensor.cc:426
rheolef::tensor_basic::svd
point_basic< T > svd(tensor_basic< T > &u, tensor_basic< T > &v, size_t dim=3) const
Definition tensor.cc:470
trace_macro
#define trace_macro(message)
Definition dis_macros.h:111
T
Expr1::float_type T
Definition field_expr.h:230
check_macro
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
rheolef
This file is part of Rheolef.
Definition compiler_eigen.h:39
rheolef::hessian_criterion
field_basic< T, M > hessian_criterion(const field_basic< T, M > &uh0, const adapt_option &opts)
Definition adapt.cc:104
rheolef::adapt
geo_basic< T, M > adapt(const field_basic< T, M > &uh, const adapt_option &opts)
adapt(uh,opts): see the adapt page for the full documentation
Definition adapt.cc:179
rheolef::hessian
field_basic< T, M > hessian(const field_basic< T, M > &uh)
Definition adapt.cc:73
rheolef::pow
space_mult_list< T, M > pow(const space_basic< T, M > &X, size_t n)
Definition space_mult.h:120
rheolef::proj
field_basic< T, M > proj(const field_basic< T, M > &uh, const std::string &approx="P1")
Definition adapt.cc:39
rheolef::adapt_gmsh
geo_basic< T, M > adapt_gmsh(const field_basic< T, M > &uh, const adapt_option &opts)
Definition adapt_gmsh.cc:41
rheolef::smooth
field_basic< T, M > smooth(const field_basic< T, M > &uh, size_t n=1)
Definition adapt.cc:56
rheolef::h_local
details::field_expr_v2_nonlinear_terminal_function< details::h_local_pseudo_function< Float > > h_local()
h_local: see the expression page for the full documentation
Definition field_expr_terminal.h:520
rheolef::adapt_bamg
geo_basic< T, M > adapt_bamg(const field_basic< T, M > &uh, const adapt_option &opts)
Definition adapt_bamg.cc:63
rheolef::diag
csr< T, M > diag(const vec< T, M > &d)
Definition csr.cc:56
rheolef::trans
csr< T, sequential > trans(const csr< T, sequential > &a)
trans(a): see the form page for the full documentation
Definition csr.h:455
d
size_t d
Definition neumann-laplace-lambda.cc:28
lambda
Definition yield_slip_circle.h:34
p
Definition sphere.icc:25
rheolef::adapt_option
adapt_option: see the adapt page for the full documentation
Definition adapt.h:147
rheolef::adapt_option::err
Float err
Definition adapt.h:151
rheolef::adapt_option::generator
std::string generator
Definition adapt.h:149
rheolef::adapt_option::n_smooth_metric
size_type n_smooth_metric
Definition adapt.h:159
rheolef::adapt_option::hmax
Float hmax
Definition adapt.h:155
rheolef::adapt_option::hcoef
Float hcoef
Definition adapt.h:153
rheolef::adapt_option::hmin
Float hmin
Definition adapt.h:154
rheolef::adapt_option::isotropic
bool isotropic
Definition adapt.h:150
u
Definition leveque.h:25