Provided by: librheolef-dev_5.93-2_i386 bug

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;
        };