librheolef-dev 6.7-6 / usr / include / rheolef / adapt.h

#ifndef _RHEO_ADAPT_OPTION_H
#define _RHEO_ADAPT_OPTION_H
///
/// This file is part of Rheolef.
///
/// Copyright (C) 2000-2009 Pierre Saramito <Pierre.Saramito@imag.fr>
///
/// Rheolef is free software; you can redistribute it and/or modify
/// it under the terms of the GNU General Public License as published by
/// the Free Software Foundation; either version 2 of the License, or
/// (at your option) any later version.
///
/// Rheolef is distributed in the hope that it will be useful,
/// but WITHOUT ANY WARRANTY; without even the implied warranty of
/// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
/// GNU General Public License for more details.
///
/// You should have received a copy of the GNU General Public License
/// along with Rheolef; if not, write to the Free Software
/// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
/// 
/// =========================================================================

#include "rheolef/field.h"
namespace rheolef { 

/*Class:
NAME: @code{adapt} - mesh adaptation
@findex  adapt
@clindex field
@clindex geo
@clindex adapt_option_type

SYNOPSYS:
 geo adapt (const field& phi);
 geo adapt (const field& phi, const adapt_option_type& opts);
DESCRIPTION:
@toindex  bamg
@toindex  gmsh
 The function @code{adapt} implements the mesh adaptation procedure,
 based on the @code{gmsh} (isotropic) or @code{bamg} (anisotropic) mesh generators.
 The @code{bamg} mesh generator is the default in two dimension.
 For dimension one or three, @code{gmsh} is the only generator supported yet.
 In the two dimensional case, the @code{gmsh} correspond to the @code{opts.generator="gmsh"}.

 The strategy based on a metric determined from the Hessian of
 a scalar governing field, denoted as @code{phi}, and that is supplied by the user.
 Let us denote by @code{H=Hessian(phi)} the Hessian tensor of the field @code{phi}.
 Then, @code{|H|} denote the tensor that has the same eigenvector as @code{H},
 but with absolute value of its eigenvalues:
 @example
  |H| = Q*diag(|lambda_i|)*Qt
 @end example
 The metric @code{M} is determined from @code{|H|}.
 Recall that an isotropic metric is such that @code{M(x)=hloc(x)^(-2)*Id}
 where @code{hloc(x)} is the element size field and @code{Id} is the
 identity @code{d*d} matrix, and @code{d=1,2,3} is the physical space dimension.

GMSH ISOTROPIC METRIC:
 @example
            max_(i=0..d-1)(|lambda_i(x)|)*Id
  M(x) = -----------------------------------------
         err*hcoef^2*(sup_y(phi(y))-inf_y(phi(y)))
 @end example
 Notice that the denominator involves a global (absolute) normalization
 @code{sup_y(phi(y))-inf_y(phi(y))} of the governing field @code{phi}
 and the two parameters @code{opts.err}, the target error,
 and @code{opts.hcoef}, a secondary normalization parameter (defaults to 1).

BAMG ANISOTROPIC METRIC:
 There are two approach for the normalization of the metric.
 The first one involves a global (absolute) normalization:
 @example
                         |H(x))|
  M(x) = -----------------------------------------
         err*hcoef^2*(sup_y(phi(y))-inf_y(phi(y)))
 @end example
 The first one involves a local (relative) normalization:
 @example
                         |H(x))|
  M(x) = -----------------------------------------
         err*hcoef^2*(|phi(x)|, cutoff*max_y|phi(y)|)
 @end example
 Notice that the denominator involves a local value @code{phi(x)}.
 The parameter is provided by the optional variable @code{opts.cutoff};
 its default value is @code{1e-7}.
 The default strategy is the local normalization.
 The global normalization can be enforced by setting
 @code{opts.additional="-AbsError"}.

 When choosing global or local normalization ?
 
 When the governing field @code{phi} is bounded, 
 i.e. when @code{err*hcoef^2*(sup_y(phi(y))-inf_y(phi(y)))}
 will converge versus mesh refinement to a bounded value,
 the global normalization defines a metric that is mesh-independent
 and thus the adaptation loop will converge.

@cindex corner singularity
 Otherwise, when @code{phi} presents singularities, with unbounded
 values (such as corner singularity, i.e. presents peacks when represented
 in elevation view), then the mesh adaptation procedure
 is more difficult. The global normalization
 divides by quantities that can be very large  and the mesh adaptation
 can diverges when focusing on the singularities.
 In that case, the local normalization is preferable.
 Moreover, the focus on singularities can also be controled 
 by setting @code{opts.hmin} not too small.

 The local normalization has been choosen as the default since it is 
 more robust. When your field @code{phi} does not present singularities,
 then you can swith to the global numbering that leads to a best
 equirepartition of the error over the domain.
End: */

//<adapt:
struct adapt_option_type {
    typedef std::vector<int>::size_type size_type;
    std::string generator;
    bool isotropic;
    Float err;
    Float errg;
    Float hcoef;
    Float hmin;
    Float hmax;
    Float ratio;
    Float cutoff;
    size_type n_vertices_max;
    size_type n_smooth_metric;
    bool splitpbedge;
    Float thetaquad;
    Float anisomax;
    bool clean;
    std::string additional;
    bool double_precision;
    Float anglecorner;	// angle below which bamg considers 2 consecutive edge to be part of
    			// the same spline
    adapt_option_type() :
        generator(""),
	isotropic(true), err(1e-2), errg(1e-1), hcoef(1), hmin(0.0001), hmax(0.3), ratio(0), cutoff(1e-7),
	n_vertices_max(50000), n_smooth_metric(1), 
	splitpbedge(true), thetaquad(std::numeric_limits<Float>::max()),
	anisomax(1e6), clean(false), additional("-RelError"), double_precision(false),
	anglecorner(0)
     {}
};
template <class T, class M>
geo_basic<T,M>
adapt (
  const field_basic<T,M>& phi,
  const adapt_option_type& options = adapt_option_type());
//>adapt:

}// namespace rheolef
#endif // _RHEO_ADAPT_OPTION_H