geomap(9)

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,

std::string quadrature="default", size_t order=0) ;

// 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&
quadrature () const

{
// TO BE REMOVED: backward compatibility
if (!use_space)
return _quadrature; }

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
void init (const field& advection);
// 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;
field _advection;
bool advected;
std::string _quadrature;
size_t _order;
bool _allow_approximate_edges;

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

std::vector<point> quad_point;
std::vector<Float> quad_weight;

// 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) :

_Th_1 (Vh_1.get_geo()), _Th_2 (Th_2), _advection(advection), advected(true), _quadrature("default"), _order(0),
_allow_approximate_edges(true), use_space(true), _Vh_1(Vh_1),
_is_inside(), _barycentric_weight(0)

{ init (advection); };

inline
geomap::geomap (const space& Vh_1, const field& advection, const geomap_option_type& opt) :

_Th_1 (Vh_1.get_geo()), _Th_2 (Vh_1.get_geo()), _advection(advection), advected(true), _quadrature("default"), _order(0),
_allow_approximate_edges(true), use_space(true), _Vh_1(Vh_1),
_is_inside(), _barycentric_weight(0), _option(opt)

{ init (advection); };

inline
geomap::geomap (

const geo& Th_2,
const geo& Th_1,
const field& advection,
std::string quadrature,
size_t order)

:

_Th_1 (Th_1),
_Th_2 (Th_2),
_advection(advection),
advected(true),
_quadrature(quadrature),
_order(order),
_allow_approximate_edges(true),
use_space(false),
_Vh_1(),
_is_inside(),
_barycentric_weight(0),
_option()

{

if (_advection.get_geo() !=_Th_1)

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,
std::string quadrature,
size_t order,
bool allow_approximate_edges)

: _Th_1 (Th_1),

_Th_2 (Th_2),
_advection(),
advected(false),
_quadrature(quadrature),
_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>
quadrature_values (size_t iK) const;

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

const std::string&
quadrature() const
{ return g._quadrature; };

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

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

protected:

field f;
geomap g;

};

rheolef-5.89 geomap(9)