librheolef-dev 5.93-2 / usr / include / rheolef / geomap.h

#ifndef _RHEO_GEOMAP_H
#define _RHEO_GEOMAP_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
/// 
/// =========================================================================
// 
// Composition of a field by an interpolation
// or an advection operator defined on a space
// or on a lattice over a mesh.
//
// author: jocelyn.etienne@imag.fr
//
// date: 4 Apr 2002
#include "rheolef/field.h"
#include "rheolef/meshpoint.h"
#include "rheolef/hazel.h"
namespace rheolef { 

/*Class:geomap
NAME: @code{geomap} - discrete mesh advection by a field: gh(x)=fh(x-dt*uh(x))
@cindex  mesh
@cindex  advection
@cindex  method of characteristic
@clindex geomap
@clindex space
@clindex geo

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:
@noindent
 The following code compute gh(x)=fh(x-dt*uh(x)):
 @example
  field uh = interpolate (Vh, u);
  field fh = interpolate (Fh, f);
  geomap X (Fh, -dt*uh);
  field gh = compose (fh, X);
 @end example
 For a complete example, see @file{convect.cc} in the
 example directory.
End: */

//! Maps the set of Vh space dof's or a lattice of quadrature points onto the meshpoints
//! obtained through the advection of the mesh.

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

//! Field seen through a composition with a lattice-based geomap.
/*! The composition remains a dynamic one, upon demand. */
//<fieldog:
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;
 };
//>fieldog:

template<class Function>
field compose (const field& f, const geomap& g, Function f_outside)
{
	check_macro(g.no_barycentric_weight(), "Unsupported barycentric weight in geomap");
	// multi-component: assume same approx for each component !
	typedef field::size_type size_type;
	const space& _Vh_1 =g.get_space();
	check_macro (g.get_target_triangulation() == f.get_space().get_geo(),
		"compose : Meshes do not match : " << g.get_target_triangulation().name() 
		<< " and " << f.get_space().get_geo().name() << ".");
	field fog(_Vh_1);
	size_type offset = 0;
	for (size_type i_comp = 0; i_comp < _Vh_1.n_component(); i_comp++) {
		geomap::const_iterator ig = g.begin();
		size_type sz_comp = _Vh_1.size_component(i_comp);
        check_macro (sz_comp == g.size(), "incompatible " << i_comp << "-th component size "
		<< sz_comp << " and geomap size " << g.size());
		for (size_type i_dof = 0; i_dof < sz_comp; i_dof++, ig++) {
			Float value;
			if (g.is_inside (i_dof)) {
			  value = f.evaluate (*ig, i_comp);
			} else {
			/* HAVE TO USE TARGET TRIANGULATION, NOT ORIGIN ! */
			  point x = g.get_target_triangulation().dehatter (*ig);
			  value = f_outside (x);
			}
			fog.at(offset+i_dof) = value;
		}
		offset += sz_comp;
	}
	return fog;
}
}// namespace rheolef
#endif // _RHEO_GEOMAP_H