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

#ifndef _RHEO_FIELD_ELEMENT_REP_H
#define _RHEO_FIELD_ELEMENT_REP_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
/// 
/// =========================================================================

// used for assembly process of right-hand-sides

#include "rheolef/point.h"
#include "rheolef/basis.h"
#include "rheolef/numbering.h"
#include "rheolef/quadrature.h"
#include "rheolef/basis_on_lattice.h"
#include "rheolef/piola_v1.h"
#include "rheolef/field.h"
#include "rheolef/characteristic.h"
namespace rheolef { 

class field_element {
public:

// typedef:

    typedef size_t size_type;

// allocators:

    field_element();
    ~field_element();

// accessors

    template <class Function>
    void operator() (
	const geo_element& K,
	Function f,
	std::vector<Float>& uk) const;

    template <class T>
    void eval (const geo_element& K,
	const std::vector<T>& fq,
	std::vector<Float>& uk) const;

    template <class Function>
    void eval_tag (const geo_element& K, const Function& f,
	std::vector<Float>& uk, Float) const;
    template <class Function>
    void eval_tag (const geo_element& K, const Function& f,
	std::vector<Float>& uk, point) const;
    template <class Function>
    void eval_tag (const geo_element& K, const Function& f,
	std::vector<Float>& uk, tensor) const;

// modifiers (const, since data are mutable: avoid copies):

    void initialize (const space& X,
	const quadrature_option_type& qopt
	 = quadrature_option_type(quadrature_option_type::gauss_lobatto)) const;

// accessors (for building elementary fields):

    const space&     get_space() const;
    const basis&     get_basis() const;
    const numbering& get_numbering() const;
    const basis&     get_p1_transformation () const;

    const geo&       get_global_geo() const;
    size_type        coordinate_dimension() const;
    fem_helper::coordinate_type coordinate_system_type() const;
    geo::const_iterator_node    begin_vertex() const;

// data:
protected:
    mutable space                  _X;	// origin dof numbering
    mutable quadrature             _quad;
    mutable basis_on_quadrature	   _b_table;
    mutable piola_on_quadrature    _piola;
    mutable bool                   _initialized;
};
// -----------------------------------------------------------------------
// inlined
// -----------------------------------------------------------------------

inline
const space& 
field_element::get_space() const
{
  return _X;
}
inline
const basis& 
field_element::get_basis() const
{
  return get_space().get_basis();
}
inline
const numbering& 
field_element::get_numbering() const
{
  return get_space().get_numbering();
}
inline
fem_helper::coordinate_type
field_element::coordinate_system_type() const
{
    return get_global_geo().coordinate_system_type();
}
inline
field_element::size_type
field_element::coordinate_dimension() const
{
  return get_global_geo().dimension();
}
inline
geo::const_iterator_node
field_element::begin_vertex() const
{
  return get_global_geo().begin_node();
}
inline
field_element::field_element()
 : _X(),
   _quad(),
   _b_table(),
   _piola(),
   _initialized(false)
{
}
inline
field_element::~field_element()
{
}
template <class Function>
void 
field_element::eval_tag (const geo_element& K,
	const Function& f,
	std::vector<Float>& uk,
        Float) const
{
  std::vector<Float> fq(_quad.size(K));
  quadrature::const_iterator first = _quad.begin(K);
  quadrature::const_iterator last  = _quad.end  (K);
  for (size_type q = 0; first != last; first++, q++) {
    point xq = get_space().get_geo().dehatter ((*first).x, K.index());
    fq[q] = f(xq);
  }
  eval (K, fq, uk);
}
template <class Function>
void 
field_element::eval_tag (const geo_element& K,
	const Function& f,
	std::vector<Float>& uk,
        point) const
{
  // TODO: also vector and tensor valued fcts (see Function=field specialization)
  error_macro ("riez_representer not yet supported with vector-valued functions");
}
template <class Function>
void
field_element::eval_tag (const geo_element& K,
	const Function& f,
	std::vector<Float>& uk,
        tensor) const
{
  error_macro ("riez_representer not yet supported with tensor-valued functions");
}
template <class Function>
void 
field_element::operator() (const geo_element& K,
	Function f,
	std::vector<Float>& uk) const
{
  typedef typename Function::result_type result_type;
  eval_tag (K, f, uk, result_type());
}
// specializations for special Function classes:
//  * field 
//  * field_o_characteristic
//    -> more efficient: since we loop on element K
//       there is no need to localize x in K
//  * simple functions, for convenience
template<>
void 
field_element::operator() (const geo_element& K,
	const field_o_characteristic& f, std::vector<Float>& uk) const;

template <>
void 
field_element::operator() (const geo_element& K,
	const field& f,
	std::vector<Float>& uk) const;

template <>
inline
void 
field_element::operator() (const geo_element& K,
	Float (*f)(const point&),
	std::vector<Float>& uk) const
{
  eval_tag (K, std::ptr_fun(f), uk, Float());
}
template <>
inline
void 
field_element::operator() (const geo_element& K,
	point (*f)(const point&),
	std::vector<Float>& uk) const
{
  eval_tag (K, std::ptr_fun(f), uk, point());
}
template <>
inline
void 
field_element::operator() (const geo_element& K,
	tensor (*f)(const point&),
	std::vector<Float>& uk) const
{
  eval_tag (K, std::ptr_fun(f), uk, tensor());
}
}// namespace rheolef
#endif // _RHEO_FIELD_ELEMENT_REP_H