librheolef-dev 6.6-1build2 / usr / include / rheolef / field_component.h

#ifndef _RHEOLEF_FIELD_COMPONENT_H
#define _RHEOLEF_FIELD_COMPONENT_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
///
/// =========================================================================

/* TODO: when slip B.C. or Raviart-Thomas RT H(div) approx :
   uh[0] values are computed on the fly in a temporary
 
struct field_component {	
// data:
   field_basic<T,M>::iterator _start, _last;
   const space_constit&       _constit;         // a ce niveau
   field*                     _tmp_comp_ptr;    // si necesite un temporaire
// actions:
   field_component (field& uh, i_comp) {
	constit = uh.get_space().get_constitution();
	_constit = constit[i_comp]);
	if (! constit.has_slip_bc() && ! constit.has_hdiv()) {
          _tmp_comp_ptr = 0;
	  _star & _last : pointe dans le field uh
	} else {
	  space Vi (_constit);
          _tmp_comp_ptr = new field (Vi);
	  uh.compute_slip_comp (i_comp, *_tmp_comp_ptr);
	  _star & _last : pointe dans *_tmp_comp_ptr
	}
   }
   ~field_component () {
      if (_tmp_comp_ptr) delete _tmp_comp_ptr;
   }
};
*/

#include "rheolef/field.h"
#include "rheolef/field_indirect.h"

namespace rheolef {

template<class Expr> class field_expr; // forward declaration

// =========================================================================
// field_component
// =========================================================================
template<class T, class M>
class field_component {
public:

// typename:

  typedef typename field_basic<T,M>::size_type      size_type;
  typedef T                                         scalar_type;
  typedef T                                         value_type;
  typedef typename field_basic<T,M>::iterator       iterator;
  typedef typename field_basic<T,M>::const_iterator const_iterator;

// allocators:

  field_component ();
  field_component (field_basic<T,M>&     uh, size_type i_comp);
  field_component (field_component<T,M>& uh, size_type i_comp);
  field_component<T,M>& operator= (const field_component<T,M>& uh_comp);
  field_component<T,M>& operator= (const field_basic<T,M>& uh);
  template<class Expr>
  field_component<T,M>& operator= (const field_expr<Expr>& expr);
  const T& operator= (const T& alpha);

// high-level accessors & modifiers:

  field_indirect<T,M> operator[] (const geo_basic<T,M>& dom);
  field_indirect<T,M> operator[] (const std::string&    dom_name);

// recursive call:

  field_component<T,M> operator[] (size_t i_comp);

// low-level accessors & modifiers:

  const space_constitution<T,M>& get_constitution() const { return _constit; }
  std::string stamp() const { return _constit.stamp(); }
  geo_basic<T,M>   get_geo() const { return _constit.get_geo(); }
  space_basic<T,M> get_space() const { return space_basic<T,M>(_constit); }
  const distributor& ownership() const { return _constit.ownership(); }	
  const communicator& comm() const { return ownership().comm(); }
  size_type     ndof() const { return ownership().size(); }
  size_type dis_ndof() const { return ownership().dis_size(); }
        T& dof (size_type idof)       { return _start [idof]; }
  const T& dof (size_type idof) const { return _start [idof]; }
  iterator begin_dof() { return _start; }
  iterator end_dof()   { return _last; }
  const_iterator begin_dof() const { return const_iterator(_start); }
  const_iterator end_dof() const   { return const_iterator(_last); }

protected:
// data:
  space_constitution<T,M>        _constit;     // at this level
  iterator                       _start;
  iterator                       _last;
};
// =========================================================================
// field_component_const
// =========================================================================
template<class T, class M>
class field_component_const {
public:

// typename:

  typedef typename field_basic<T,M>::size_type      size_type;
  typedef T                                         scalar_type;
  typedef T                                         value_type;
  typedef typename field_basic<T,M>::const_iterator const_iterator;

// allocators:

  field_component_const ();
  field_component_const (const field_basic<T,M>& uh, size_type i_comp);
  field_component_const (const field_component<T,M>& uh_comp);
  field_component_const<T,M>& operator= (const field_component_const<T,M>& uh_comp);
  const space_constitution<T,M>& get_constitution() const { return _constit; }
  std::string stamp() const { return _constit.stamp(); }
  const distributor& ownership() const { return _constit.ownership(); }	
  const communicator& comm() const { return ownership().comm(); }
  geo_basic<T,M>   get_geo() const { return _constit.get_geo(); }
  space_basic<T,M> get_space() const { return space_basic<T,M>(_constit); }
  size_type     ndof() const { return ownership().size(); }
  size_type dis_ndof() const { return ownership().dis_size(); }
  const T& dof (size_type idof) const { return _start [idof]; }
  const_iterator begin_dof() const { return const_iterator(_start); }
  const_iterator end_dof() const   { return const_iterator(_last); }
protected:
// data:
  space_constitution<T,M>        _constit;     // at this level
  const_iterator                 _start;
  const_iterator                 _last;
};
// =========================================================================
// inlined
// =========================================================================
template<class T, class M>
inline
field_component<T,M>::field_component()
 : _constit(),
   _start(),
   _last()
{
}
template<class T, class M>
inline
field_component_const<T,M>::field_component_const()
 : _constit(),
   _start(),
   _last()
{
}
template<class T, class M>
field_component<T,M>::field_component (field_basic<T,M>& uh, size_type i_comp)
 : _constit (uh.get_space().get_constitution() [i_comp]),
   _start(uh.begin_dof()),
   _last(uh.begin_dof())
{
  check_macro (i_comp < uh.size(),
      "field component index "<<i_comp<<" is out of range [0:"<<uh.size()<<"[");
  const space_constitution<T,M>& sup_constit = uh.get_space().get_constitution();
  size_type shift = 0;
  for (size_type j_comp = 0; j_comp < i_comp; j_comp++) {
    shift += sup_constit [j_comp].ndof();
  }
  size_type sz = sup_constit [i_comp].ndof();
  _start += shift;	
  _last  += shift + sz;	
}
template<class T, class M>
field_component_const<T,M>::field_component_const (const field_basic<T,M>& uh, size_type i_comp)
 : _constit (uh.get_space().get_constitution() [i_comp]),
   _start(uh.begin_dof()),
   _last(uh.begin_dof())
{
  check_macro (i_comp < uh.size(),
      "field component index "<<i_comp<<" is out of range [0:"<<uh.size()<<"[");
  const space_constitution<T,M>& sup_constit = uh.get_space().get_constitution();
  size_type shift = 0;
  for (size_type j_comp = 0; j_comp < i_comp; j_comp++) {
    shift += sup_constit [j_comp].ndof();
  }
  size_type sz = sup_constit [i_comp].ndof();
  _start += shift;	
  _last  += shift + sz;	
}
template<class T, class M>
field_component<T,M>::field_component (field_component<T,M>& uh, size_type i_comp)
 : _constit (),
   _start(uh._start),
   _last(uh._last)
{
  const space_constitution<T,M>& sup_constit = uh._constit;
  check_macro (i_comp < sup_constit.size(),
      "field component index "<<i_comp<<" is out of range [0:"<<sup_constit.size()<<"[");
  size_type shift = 0;
  for (size_type j_comp = 0; j_comp < i_comp; j_comp++) {
    shift += sup_constit [j_comp].ndof();
  }
  _constit = sup_constit [i_comp];
  size_type sz = _constit.ndof();
  _start += shift;	
  _last  += shift + sz;	
}
template<class T, class M>
inline
field_component_const<T,M>::field_component_const (const field_component<T,M>& uh_comp)
 : _constit (uh_comp.get_constitution()),
   _start(uh_comp.begin_dof()),
   _last(uh_comp.end_dof())
{
}
template<class T, class M>
inline
field_component<T,M>&
field_component<T,M>::operator= (const field_component<T,M>& uh_comp)
{
  _constit = uh_comp._constit;
  _start   = uh_comp._start;
  _last    = uh_comp._last;
  return *this;
}
template<class T, class M>
inline
field_component_const<T,M>&
field_component_const<T,M>::operator= (const field_component_const<T,M>& uh_comp)
{
  _constit = uh_comp._constit;
  _start   = uh_comp._start;
  _last    = uh_comp._last;
  return *this;
}
template<class T, class M>
inline
field_component<T,M>
field_basic<T,M>::operator[] (size_type i_comp) 
{
  dis_dof_update_needed();
  return field_component<T,M> (*this, i_comp);
}
template<class T, class M>
inline
field_component<T,M>
field_component<T,M>::operator[] (size_t i_comp)
{
  return field_component<T,M> (*this, i_comp);
}
template<class T, class M>
inline
field_component_const<T,M>
field_basic<T,M>::operator[] (size_type i_comp) const
{
  return field_component_const<T,M> (*this, i_comp);
}
template<class T, class M>
inline
field_basic<T,M>&
field_basic<T,M>::operator=  (const field_component_const<T,M>& uh_comp)
{
  if (stamp() == "") {
    resize (uh_comp.get_space());
  } else {
    check_macro (stamp() == uh_comp.stamp(), "incompatible spaces "
        << stamp() << " and " << uh_comp.stamp()
        << " in field = field[i_comp]");
  }
  dis_dof_update_needed();
  std::copy (uh_comp.begin_dof(), uh_comp.end_dof(), begin_dof());
  return *this;
}
template<class T, class M>
inline
field_component<T,M>&
field_component<T,M>::operator= (const field_basic<T,M>& uh)
{
  check_macro (stamp() == uh.stamp(), "incompatible spaces "
        << stamp() << " and " << uh.stamp()
        << " in field[i_comp] = field");
  std::copy (uh.begin_dof(), uh.end_dof(), begin_dof());
  return *this;
}
template<class T, class M>
inline
const T&
field_component<T,M>::operator= (const T& alpha)
{
  std::fill (begin_dof(), end_dof(), alpha);
  return alpha;
}
template<class T, class M>
inline
field_basic<T,M>&
field_basic<T,M>::operator=  (const field_component<T,M>& uh_comp)
{
  return operator= (field_component_const<T,M>(uh_comp));
}
template<class T, class M>
inline
field_basic<T,M>::field_basic (const field_component<T,M>& uh_comp)
 : _V(),
   _u(),
   _b(),
   _dis_dof_update_needed(true)
{
  operator= (uh_comp);
}
template<class T, class M>
inline
field_basic<T,M>::field_basic (const field_component_const<T,M>& uh_comp)
 : _V(),
   _u(),
   _b(),
   _dis_dof_update_needed(true)
{
  operator= (uh_comp);
}
#ifdef TODO
// TODO: as: field vh; din >> vh ; uh[i_comp] = vh;
template <class T, class M>
inline
idiststream& operator >> (odiststream& ids, field_basic<T,M>& u);
#endif // TODO

template <class T, class M>
odiststream& operator << (odiststream& ods, const field_component<T,M>& uh_comp)
{
  return ods << field_basic<T,M>(uh_comp);
}
template <class T, class M>
inline
odiststream& operator << (odiststream& ods, const field_component_const<T,M>& uh_comp)
{
  return ods << field_basic<T,M>(uh_comp);
}
// =========================================================================
// uh[0]["boundary"] = alpha;
// => field_component & field_indirect together
// =========================================================================
template <class T, class M>
inline
field_indirect<T,M>::field_indirect (field_component<T,M>& uh_comp, const geo_basic<T,M>& dom)
  : _V(uh_comp.get_space()),
    _W(dom, _V.get_numbering().name()),
    _dom(dom),
    _indirect(_V.build_indirect_array (_W, _dom.name())),
    _val(uh_comp.begin_dof())
{
}
template <class T, class M>
inline
field_indirect<T,M>
field_component<T,M>::operator[] (const geo_basic<T,M>& dom)
{
  return field_indirect<T,M> (*this, dom);
}
template <class T, class M>
inline
field_indirect<T,M>
field_component<T,M>::operator[] (const std::string& dom_name)
{
  return field_indirect<T,M> (*this, get_geo().operator[](dom_name));
}

} // namespace rheolef
#endif // _RHEOLEF_FIELD_COMPONENT_H