Provided by: librheolef-dev_5.93-2_i386

NAME

```       geomap - discrete mesh advection by a field: gh(x)=fh(x-dt*uh(x))

```

SYNOPSYS

```        The class geomap is a fundamental class used for the correspondance
between fields defined on different meshes or for advection problems.
This class is used for the method of characteristic.

```

EXAMPLE

```        The following code compute gh(x)=fh(x-dt*uh(x)):

field uh = interpolate (Vh, u);
field fh = interpolate (Fh, f);
geomap X (Fh, -dt*uh);
field gh = compose (fh, X);

For a complete example, see `convect.cc' in the
example directory.

```

IMPLEMENTATION

```       struct geomap_option_type {
size_t n_track_step; // loop while tracking: y = X_u(x)
geomap_option_type() : n_track_step(1) {}
};
class geomap : public Vector<meshpoint> {
public:
geomap () : advected(false), use_space(true) {}
//  Maps a quadrature-dependent lattice over Th_1 onto Th_2 triangulation.
geomap (const geo& Th_2, const geo& Th_1,
std::string quadrature="default", size_t order=0, bool allow_approximate_edges=true);

//  Maps a quadrature-dependent lattice over Th_1 onto Th_2 triangulation
//  thro' an offset vector field u to be defined on Th_1.
geomap (const geo& Th_2, const geo& Th_1, const field& advection_h_1,

//  TO BE REMOVED: backward compatibility
//  Maps space Vh_1 dof's onto the meshpoints of Th_2 triangulation
/*  geomap : ( Vh_1.dof's ) |--> ( Th_2 )
*/
geomap (const geo& Th_2, const space& Vh_1,
bool allow_approximate_edges=true );
//  Same but for P1d, using points inside the triangle to preserve discontinuity
geomap (const geo& Th_2, const space& Vh_1,
Float barycentric_weight );

//  TO BE REMOVED: backward compatibility
//  Maps space Vh_1 dof's onto the meshpoints of Th_2 triangulation
//  thro' an offset u to be defined on Vh_1 dof's.
/*  geomap : ( Vh_1.dof's ) |--> ( Th_2 )
*/
geomap (const geo& Th_2, const space& Vh_1, const field& advection_h_1);

//  Does the same, but assumes that Th_2 is the triangulation for Vh_1
geomap (const space& Vh_2, const field& advection_h_1, const geomap_option_type& opts = geomap_option_type());

~geomap () {};

//accessors
//  TO BE REMOVED: backward compatibility
const space&
get_space () const
{ if (!use_space) error_macro("Lattice-based geomaps use no space");
return _Vh_1;
};

const geo&
get_origin_triangulation () const
{ return _Th_1; };

const geo&
get_target_triangulation () const
{ return _Th_2; };

size_t
order () const
{
//  TO BE REMOVED: backward compatibility
if (!use_space)
return _order;
}

const std::string&
{
//  TO BE REMOVED: backward compatibility
if (!use_space)

bool is_inside (size_t dof) const { return _is_inside [dof]; }
bool no_barycentric_weight() const { return _barycentric_weight == Float(0); }
Float barycentric_weight() const { return _barycentric_weight; }

protected:
friend class fieldog;
//  use_space mode
//  non use_space mode
void init ();
meshpoint advect (const point& q, size_t iK_Th_1);
meshpoint robust_advect_1 (const point& x0, const point& va, bool& is_inside) const;
meshpoint robust_advect_N (const point& x0, const point& v0, const field& vh, size_t n, bool& is_inside) const;

// data:
geo _Th_1;
geo _Th_2;
size_t _order;
bool _allow_approximate_edges;

// TO BE REMOVED:
bool use_space;
space _Vh_1;

// flag when a dof go outside of the domain
std::vector<bool> _is_inside;

Float _barycentric_weight;
geomap_option_type _option;
};

inline
geomap::geomap (const geo& Th_2, const space& Vh_1, const field& advection) :
_allow_approximate_edges(true), use_space(true), _Vh_1(Vh_1),
_is_inside(), _barycentric_weight(0)

inline
geomap::geomap (const space& Vh_1, const field& advection, const geomap_option_type& opt) :
_allow_approximate_edges(true), use_space(true), _Vh_1(Vh_1),
_is_inside(), _barycentric_weight(0), _option(opt)

inline
geomap::geomap (
const geo& Th_2,
const geo& Th_1,
size_t order)
:
_Th_1 (Th_1),
_Th_2 (Th_2),
_order(order),
_allow_approximate_edges(true),
use_space(false),
_Vh_1(),
_is_inside(),
_barycentric_weight(0),
_option()
{
error_macro("The advection field should be defined on the original mesh Th_1="
<< Th_1.name() << " and not " << _advection.get_geo().name());
init();
}

inline
geomap::geomap (
const geo& Th_2,
const geo& Th_1,
size_t order,
bool allow_approximate_edges)
: _Th_1 (Th_1),
_Th_2 (Th_2),
_order(order),
_allow_approximate_edges(allow_approximate_edges),
use_space(false),
_Vh_1(),
_is_inside(),
_barycentric_weight(0),
_option()
{
init();
}

//  Composition with a field
/*  f being a field of space V, X a geomap on this space, foX=compose(f,X) is their
*  composition.
*/
field compose (const field& f, const geomap& g);
/*  send also a callback function when advection goes outside a domain
*  this is usefull when using upstream boundary conditions
*/
field compose (const field& f, const geomap& g, Float (*f_outside)(const point&));
template <class Func>
field compose (const field& f, const geomap& g, Func f_outside);

```

IMPLEMENTATION

```       class fieldog: public field // and friend of geomap.
{
public:
// Constructor
fieldog (const field& _f, const geomap& _g) : f(_f), g(_g) {};

// Accessors
Float
operator() (const point& x) const;

Float
operator() (const meshpoint& x) const;

Float
operator() (const point& x_hat, size_t e) const
{ return operator() (meshpoint(x_hat,e)); };

std::vector<Float>

std::vector<Float>
quadrature_d_dxi_values (size_t i, size_t iK) const;

const std::string&

size_t
order() const
{ return g._order; };

const space&
get_space() const
{ return f.get_space(); };

protected:
field f;
geomap g;
};
```