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

# ifndef _RHEO_SPACEREP_H
# define _RHEO_SPACEREP_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
/// 
/// =========================================================================
//
// representation for space object
//  apparent numbering = numbering of degrees of freedom (dof, d.d.l. in French)
//
// re-numbering = correspondance 
//    from global index "dof" to local one "num (dof)"
//    with blocked/not-bloked flag, i.e. imposed boundary
//    condition for this dof
//
// => we map in two vec<T> array
//
//     one for non-blocked dofs : u 	(stands for unknows)
//     the second for blocked:    b
//
// in addition, vector spaces can make use locally of tangential and normal
// dofs locally on a boundary rather than the native cartesian components,
// see has_locked function. --J.E., 25/08/2006
//
#include "rheolef/basis.h"
#include "rheolef/numbering.h"
#include "rheolef/geo.h"
#include "rheolef/tiny_matvec.h"
#include "rheolef/fem_helper.h"
#include "rheolef/basis_on_lattice.h"
namespace rheolef { 

class spacerep : public occurence
{
public:
// typedefs:

  typedef std::vector<basis>::size_type  size_type;

// data:
protected:
          geo                  _g;

  mutable bool                 _is_frozen  ; // TRUE if re-numbering is available
  mutable bool                 _has_locked ; // TRUE if local tangential-normal representation is used

          std::vector<basis>       _b;		      // basis per component
          std::vector<numbering>   _numbering;	      // global numbering
  
          // piola transformation basis
          // and precomputation of the the transformation basis polynoms
          // on nodes of the fem nodal basis 
	  basis		                    _tr_p1;
	  std::vector<basis_on_nodal_basis> _tr_p1_value;

          std::vector<std::string>  _name;	      // approx name per component
          std::vector<size_type>    _start;           // starting dof in components
  mutable std::vector<size_type>    _start_blocked;   // starting index for blocked in components
  mutable std::vector<size_type>    _start_unknown;   // starting index for unknown in components
  mutable std::vector<bool>         _is_blocked ;     // TRUE if dof (I_{global}) is blocked
  mutable std::vector<bool>         _is_periodic;     // TRUE if dof (I_{global}) is periodic
  mutable std::vector<size_type>    _period_assoc;    // tableau de correspondance pour les dof periodiques
  mutable std::vector<size_type>    _num        ;     // re-numbering from dof to index (u,b)
  mutable std::vector<size_type>    _dof_oriented_as_element; // dof-indexed array of corresponding element with same orientation 
  		// (these dofs have exactly 2 elements in common: 1D Hermite, Argyris) 
  mutable bool                 _dof_orientation_initialized ;

  // for spaces with normal-tangential local representation
  mutable std::vector<size_type>    _is_locked; // set if dof (I_{global}) has a locked direction
 	// set to numeric_limits if it has not
  	// e.g., the field must be 0 along the normal, but can take tangential values
	// value is the index to use for the _locked_dir vector.
  mutable std::vector<point>		_locked_dir;  // the direction which is locked
  mutable std::vector<size_type>	_locked_with; // dof of the other component with which a dof is locked


  // for spaces defined on a part of the boundary.
          bool                 _is_on_boundary_domain;
          const domain         _boundary_domain;
          const geo            _global_geo;
	  std::vector<size_type>    _global_dof;         // dof on domain-geo -> dof on full-geo
	  std::vector<size_type>    _boundary_dof;       // dof on full-geo   -> dof on domain-geo
	  std::vector<size_type>    _start_global_dof;   // starting global-geo dof in components

  // valued type: scalar, vector, tensor...

	  fem_helper::valued_field_type  _valued;

  // piola transformation F_K : hat_K --> K and it inverse
  meshpoint hatter   (const point& x, size_type K) const;
  point     dehatter (const meshpoint& S, bool is_a_vector=false) const;

// internals:
  size_type _n_dof       (size_type i_comp = 0);
  size_type _n_global_dof(size_type i_comp = 0);

public:
// allocators/deallocators:

  spacerep();
  spacerep(const spacerep& e);
  spacerep(const geo& g, const std::string& approx_name, const std::string& option);
  spacerep(const geo& g, const std::string& approx_name, 
	const domain& interface, const domain& subgeo, const std::string& option);
  spacerep(const geo& g, const domain& d, const std::string& approx_name);
  spacerep(const spacerep& X, const spacerep& Y);
  spacerep(const spacerep& X, size_type i_comp);
  ~spacerep();
  //
  // access
  //
  size_type  size()      const { return _num.size(); }
  size_type  dimension() const { return _g.dimension(); }
  std::string coordinate_system() const { return _g.coordinate_system(); }
  fem_helper::coordinate_type coordinate_system_type() const { return _g.coordinate_system_type(); }
  bool       is_frozen() const { return _is_frozen ; }
  
  size_type  n_component()                      const { return _start.size() - 1; }
  size_type  start (size_type i=0)              const { return _start[i]; }
  size_type  start_blocked (size_type i=0)      const { return _start_blocked[i]; }
  size_type  start_unknown (size_type i=0)      const { return _start_unknown[i]; }
  size_type  size_component(size_type i=0)      const { return start(i+1) - start(i); }
  size_type  n_blocked_component(size_type i=0) const { return start_blocked(i+1) - start_blocked(i); }
  size_type  n_unknown_component(size_type i=0) const { return start_unknown(i+1) - start_unknown(i); }

  size_type  n_blocked() const { return start_blocked(n_component()); }
  size_type  n_unknown() const { return start_unknown(n_component()); }

  size_type  start_global_dof       (size_type i=0) const { return _start_global_dof[i]; }
  size_type  n_global_dof_component (size_type i=0) const { return start_global_dof(i+1)-start_global_dof(i); }

  size_type  index (size_type dof)              const { return _num[dof]; }
  size_type  period_association (size_type dof) const { return _period_assoc[dof]; }
  bool       is_blocked(size_type dof)          const { return _is_blocked[dof]; }
  bool       is_periodic(size_type dof)         const { return _is_periodic[dof];}
  bool       is_locked(size_type dof)           const 
  	{ return _has_locked ? (_is_locked[dof]!=std::numeric_limits<size_type>::max()) : false; }
  size_type  locked_with  (size_type dof) const; 
  size_type  index_locked_with  (size_type dof) const; 
  	//! dof with the other component
  Float      locked_component   (size_type dof, size_type i_comp) const;  
  	//! i-th component of locked direction
  Float      unlocked_component (size_type dof, size_type i_comp) const;  
  	//! i-th component of unlocked direction
  bool       has_locked()                       const { return _has_locked; }

  const geo&    get_geo()                 const { return _g; }
  std::string   get_approx (size_type i=0) const { return _name[i] ; }  
  const basis&  get_basis  (size_type i=0) const { return _b[i] ; }  
  const numbering&  get_numbering  (size_type i=0) const { return _numbering[i] ; }  
  const basis&  get_p1_transformation  () const { return _tr_p1; }  
  std::string   get_valued () const;
  fem_helper::valued_field_type get_valued_type () const;

  static size_type inquire_size(const geo&, const numbering&);

  // domain-boundary spaces: extra infos
  bool          is_on_boundary_domain () const { return _is_on_boundary_domain; }
  const geo&    get_global_geo() const {
      				return _is_on_boundary_domain ? _global_geo : _g; }
  const domain& get_boundary_domain() const { return _boundary_domain; }
  size_type     global_size () const { return _boundary_dof.size(); }
  size_type     domain_dof (size_type global_dof) const { return _boundary_dof[global_dof]; }
  size_type     domain_index (size_type global_dof) const { 
				return _num[_boundary_dof[global_dof]]; } 

  //  
  // action
  //
  void block () const;
  void block (size_type i_comp) const;
  void block_Lagrange () const;
  void block_Lagrange (size_type i_comp) const;
  void block (const class domain& d) const;
  void block (const class domain& d, size_type i_comp) const;
  void block (const std::string& d_name)  const { block(_g[d_name]); }
  void block (const std::string& d_name, size_type i_comp) const { block(_g[d_name], i_comp); }
  void block_dof (size_type dof_idx) const;
  void unblock_dof (size_type dof_idx) const;
  void _build_global_dof_mapping();
  void freeze () const;
  void set_periodic (const class domain& d1, const class domain& d2) const;
  void set_periodic (const std::string& d1_name, const std::string& d2_name) const {
       set_periodic( _g[d1_name], _g[d2_name] );
  }
  void set_periodic (const class domain& d1, const class domain& d2, size_type i_comp) const;
  void lock_components (const class domain& d, point locked_dir) const;
  void lock_components (const std::string& d_name, point locked_dir) const 
  	{ lock_components(_g[d_name], locked_dir); }
  template <class Function>
  void lock_components (const std::string& d_name, Function locked_dir) const
  	{ lock_components(_g[d_name], locked_dir); }
  template <class Function>
  void lock_components (const class domain& d, Function locked_dir) const;
  
  void  set_dof       (const geo_element& K, tiny_vector<size_type>& idx) const;
  void  set_global_dof(const geo_element& K, tiny_vector<size_type>& idx) const;
  void  set_dof       (const geo_element& K, tiny_vector<size_type>& idx, size_type i_comp) const;
  void  set_dof       (const geo_element& K, tiny_vector<size_type>& idx, size_type i_comp,
  		reference_element::dof_family_type family) const;
  void  set_global_dof(const geo_element& K, tiny_vector<size_type>& idx, size_type i_comp) const;

  point x_dof      (const geo_element& K, size_type iloc) const ;
     //! Hermite/Argyris: returns 1 if the global dof contains the derivative along outward normal of K, -1 else  
  void  dof_orientation(const geo_element& K, tiny_vector<int>& orientation) const;
  void  initialize_dof_orientation() const;

  //
  // input/output
  //
  friend std::ostream& operator << (std::ostream&, const spacerep&);
  friend std::istream& operator >> (std::istream&, spacerep&);
  std::ostream& dump(std::ostream& s = std::cerr) const;
  void check() const;
  //
  // comparator
  //
  bool operator == (const spacerep&) const;
  friend class space;
}; // fin de la classe spacerep

inline
spacerep::size_type 
spacerep::locked_with  (size_type dof) const
 {
    check_macro(_is_locked[dof]!=std::numeric_limits<size_type>::max(), "Unexpected error, dof=" << dof);
    return _locked_with[dof];
 }
inline
spacerep::size_type 
spacerep::index_locked_with  (size_type dof) const
 {
    check_macro(_is_locked[dof]!=std::numeric_limits<size_type>::max(), "Unexpected error, dof=" << dof);
    return _num[_locked_with[dof]];
 }
inline
Float     
spacerep::locked_component   (size_type dof, size_type i_comp) const
 {
    check_macro(_is_locked[dof]!=std::numeric_limits<size_type>::max(), "Unexpected error");
    return _locked_dir[_is_locked[dof]][i_comp];
 }
inline
Float
spacerep::unlocked_component (size_type dof, size_type i_comp) const
 {
    check_macro(_is_locked[dof]!=std::numeric_limits<size_type>::max(), "Unexpected error");
    switch (i_comp) {
      case 0 : return -_locked_dir[_is_locked[dof]][1]; break;
      case 1 : return _locked_dir[_is_locked[dof]][0]; break;
      default: error_macro("Unexpected error");
    }
    return 0;
 }

template <class Function>
void
spacerep::lock_components (const domain& d, Function locked_dir) const
{
  _has_locked=true;
  check_macro(dimension()==2, "Components-locking only implemented in 2D");
  if (n_component()!=2) warning_macro("Components-locking only locks the first 2 components of vectors");
  size_type n_comp=2;
  check_macro(get_approx(0)==get_approx(1), "Components-locking only implemented when vector components are based on a the same approximation");
  if (_is_frozen) {
    error_macro("numbering already done -- cannot lock the components on any domain anymore");
  }
  resize(_is_locked, _is_blocked.size(), std::numeric_limits<size_type>::max());
  resize(_locked_with, _is_blocked.size(), std::numeric_limits<size_type>::max());
  size_type i_locked_domain=0;
  tiny_vector<size_type> idx0;
  tiny_vector<size_type> idx1;
  domain::const_iterator last_side = d.end(); 
  for (domain::const_iterator iter_side = d.begin(); iter_side != last_side; iter_side++) {
    const geo_element& S = *iter_side ;
    set_dof (S, idx0, 0);
    set_dof (S, idx1, 1);
    for (size_type i = 0 ; i < idx0.size(); i++) {
        _is_locked [idx0(i)] = i_locked_domain;
        _is_locked [idx1(i)] = i_locked_domain++;
	_is_blocked[idx1(i)] = true; // blocks only one of the two components
	_locked_with [idx0(i)] = idx1(i);
	_locked_with [idx1(i)] = idx0(i);
	_locked_dir.insert(_locked_dir.end(), locked_dir(x_dof(S, i)));
    }
  }
}

}// namespace rheolef
# endif // _RHEO_SPACEREP_H