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

# ifndef _RHEO_SPACE_H
# define _RHEO_SPACE_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
/// 
/// =========================================================================

namespace rheolef { 
struct space_component;
struct const_space_component;
}// namespace rheolef

#include "rheolef/spacerep.h"

/*Class:space
NAME: @code{space} -- piecewise polynomial finite element space
@clindex space
DESCRIPTION:
@noindent
The @code{space} class contains some numbering 
for unknowns and blocked @emph{degrees of freedoms}
related to a given mesh and polynomial approximation.

DEGREE OF FREEDOM NUMBERING:
@cindex degree of freedom
@clindex field
@clindex form
@clindex geo
@fiindex @code{.field} field
@c
@noindent
Numbering of degrees of freedom (or shortly, @dfn{dof}) depends upon
the mesh and the piecewise polynomial approximation used.
This numbering is then used by the @code{field} and @code{form} classes.
See also @ref{field class} and @ref{form class}.

@noindent
The degree-of-freedom (or shortly, @dfn{dof}) follows some simple rules,
that depends of course of the approximation.
These rules are suitable for easy @code{.field} file retrieval 
(@pxref{field class}).
 
@apindex P0
@apindex bubble
@apindex P1
@apindex P2
@apindex P1d
@c
@table @code
   @itemx P0
   @itemx bubble
   	dof numbers follow element numbers in mesh.
   @itemx P1
   	dof numbers follow vertice numbers in mesh.
   @itemx P2
   	dof numbers related to vertices follow vertice numbers in mesh.
        Follow dof numbers related to edges, in the edge numbering order.
   @itemx P1d
   	dof numbers follow element numbers in mesh. In each element, we
	loop on vertices following the local vertice order provided in
	the mesh.
@end table
 
UNKNOWN AND BLOCKED DEGREE OF FREEDOM:
@cindex unknow degree of freedom
@cindex blocked degree of freedom
@c
@noindent
A second numbering is related to
the @dfn{unknown} and @dfn{blocked} degrees of freedom.
@example
      geo omega("square");
      space V(omega,"P1");
      V.block ("boundary");
@end example
@noindent
Here, all degrees of freedom located on the domain 
@code{"boundary"} are marked as blocked.

FILE FORMAT:
@fiindex @file{.space} space
@noindent
File output format for @code{space} is not of practical interest.
This file format is provided for completness.
@noindent
Here is a small example
@example
   space
   1 6
   square
   P1
      0 B
      1 B
      0
      1
      2 B
      2
@end example
@noindent
where the @samp{B} denotes blocked degres of freedom.
The method @code{set_dof(K, dof_array)} gives the numbering
in @code{dof_array}, supplying an element @code{K}.
The method @code{index(dof)} gives the unknown or blocked
numbering from the initial numbering, suppling a given @code{dof}
  
ADDITIONAL FEATURES:
@noindent
The method @code{xdof(K, i)} gives the geometrical location of
the i-th degree of freedom in the element @code{K}.

@noindent
Product spaces, for non-scalar fields, such as vector-valued 
or tensor-valued (symmetric) fields,
can also be defined:
@example
      space U (omega, "P1", "vector");
      space T (omega, "P1", "tensor");
@end example
Then, the number of component depends upon the geometrical
dimension of the mesh @code{omega}.

@noindent
Arbitrarily product spaces can also be constructed
@example
      space X = V0*V1;
@end example
Spaces for fields defined on a boundary domain can also be 
constructed
@noindent
@example
      space W (omega, omega["top"], "P1");
@end example
The correspondance between the degree-of-freedom numbering for
the space @code{W} and the global degree-of-freedom numbering
for the space @code{V} is supported.
The method @code{set_global_dof(S, dof_array)} gives the global numbering
in @code{dof_array}, supplying a boundary element @code{S}.
The method @code{domain_index(global_dof)} gives the unknown or blocked
numbering from the initial numbering, suppling a given @code{global_dof},
related to the space @code{V}.

AUTHORS: Pierre.Saramito@imag.fr
DATE:   14 may 1997   update: 28 jan 1998
METHODS: @space
End:
*/

namespace rheolef { 
//<space:  
class space : public smart_pointer<spacerep> {
public:
// typdefs:

    typedef spacerep::size_type size_type;
 
// allocator/deallocator:

    space ();
    space (const geo& g, const std::string& approx_name, const std::string& valued = "scalar");
    space (const geo& g, const std::string& approx_name, 
	const domain& interface, const domain& subgeo, const std::string& valued = "scalar");
	//! For spaces with a discontinuity through domain `interface', but otherwise continuous
    space (const geo& g, const domain& d, const std::string& approx_name);
    space(const const_space_component&);

    space operator = (const const_space_component&);
    friend space operator * (const space&, const space&);
    friend space pow (const space&, size_type);

// modifiers:

    void block () const;
    void block_Lagrange () const;
    void block (size_type i_comp) const;
    void block_Lagrange (size_type i_comp) const;
    void block (const std::string& d_name) const;
    void block (const domain& d) const;
    void block (const std::string& d_name,  size_type i_comp) const;
    void block (const domain& d, size_type i_comp) const;
    void block_dof (size_type dof_idx, size_type i_comp) const;
    void block_dof (size_type dof_idx) const;
    void unblock_dof (size_type dof_idx, size_type i_comp) const;
    void unblock_dof (size_type dof_idx) const;
    void set_periodic (const domain& d1, const domain& d2) const;
    void set_periodic (const std::string& d1_name, const std::string& d2_name) const;
    void lock_components (const domain& d, point locked_dir) const;
    void lock_components (const std::string& d_name, point locked_dir) const;
    template <class Function>
    void lock_components (const std::string& d_name, Function locked_dir) const;
    	//! Allows to enforce vector-boundary conditions.
	/*! Allows to enforce vector-boundary conditions.
	 *!  
	 *! Current limitation: only 2D vector fields with components having same FE approximation.
	 *! The space has a dof for the direction normal to locked_dir, while the value in the
	 *! direction locked_dir is blocked. 
	 *! The field::at function implements the reconstruction of the original cartesian components.
	 */
    template <class Function>
    void lock_components (const domain& d, Function locked_dir) const;

// accessors:

    size_type  size()      const;
    size_type  degree ()   const;
    size_type  n_blocked() const;
    size_type  n_unknown() const;
    size_type  dimension() const;
    std::string coordinate_system() const;
    fem_helper::coordinate_type coordinate_system_type() const;


    const geo&                get_geo   () const;
    const basis&              get_basis (size_type i_component = 0) const;
    const numbering&          get_numbering (size_type i_component = 0) const;
    std::string               get_approx(size_type i_component = 0) const;
    const basis&              get_p1_transformation () const;

    std::string               get_valued() const;
    fem_helper::valued_field_type get_valued_type() const;

    size_type  n_component() const;
    space_component operator [] (size_type i_comp);
    const_space_component operator [] (size_type i_comp) const;

    size_type  size_component      (size_type i_component = 0) const;
    size_type  n_unknown_component (size_type i_component = 0) const;
    size_type  n_blocked_component (size_type i_component = 0) const;
    size_type  start               (size_type i_component = 0) const;
    size_type  start_unknown       (size_type i_component = 0) const;
    size_type  start_blocked       (size_type i_component = 0) const;

    size_type  index              (size_type degree_of_freedom) const;
    size_type  period_association (size_type degree_of_freedom) const;
    bool       is_blocked         (size_type degree_of_freedom) const;
    bool       is_periodic        (size_type degree_of_freedom) const;
    
    bool       has_locked         () const; 
		 //! for (2D) vector spaces only, field is only blocked along locked_dir
    bool       is_locked          (size_type degree_of_freedom) const; 
    	//! for tangential-normal representation
    size_type  locked_with  (size_type degree_of_freedom) const; 
    size_type  index_locked_with  (size_type degree_of_freedom) const; 
        //! dof with the other component (thru space::index() as well)
    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

    void freeze () const;
    bool is_frozen() const;
	

    // informations related to boundary spaces
    bool  is_on_boundary_domain() const;
    const domain& get_boundary_domain() const;
    const geo&    get_global_geo     () const;
    size_type  global_size () const;
    size_type  domain_dof   (size_type global_dof) const;
    size_type  domain_index (size_type global_dof) const;

    // get indices of degrees of freedom associated to an element
	//  the tiny_vector is a "local index" -> "global index" table
	//	"global" variants differ only for boundary spaces (non-global uses domain-indices, global use
	//  mesh-wide indices)
    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_global_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;
    point x_dof      (const geo_element& K, size_type i_local) 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;

    // piola transformation
    meshpoint hatter (const point& x, size_type K) const;
    point dehatter (const meshpoint& S) const;
    point dehatter (const point& x_hat, size_type e) const;

// comparator

    bool operator == (const space&) const;
    bool operator != (const space&) const;

// input/output

   friend std::ostream& operator << (std::ostream&, const space&);
   friend std::istream& operator >> (std::istream&, space&);
   std::ostream& dump(std::ostream& s = std::cerr) const;
   void check() const;

// inquires:

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

// implementation:
protected:
   space(smart_pointer<spacerep>);
   space (const space&, const space&); // X*Y
   space (const space&, size_type i_comp); // X[i_comp]
   friend class field;
};
//>space:

// ---------------------------------------------------------
// V[1].block("boundary");        multicomponent spaces...
// ---------------------------------------------------------
struct space_component {
    typedef space::size_type size_type;
    space_component();
    space_component(space& V, size_type i);
    void block (const domain& d);
    void block (const std::string& d_name);
    space*    _pV;
    size_type _i_comp;
};
struct const_space_component {
    typedef space::size_type size_type;
    const_space_component();
    const_space_component(const space& V, size_type i);
    const_space_component(const space_component&);
    const space*    _pV;
    size_type _i_comp;
};
inline
space_component::space_component()
 : _pV(), _i_comp(0) {}

inline
const_space_component::const_space_component()
 : _pV(), _i_comp(0) {}

inline
const_space_component::const_space_component(const space_component& Vi)
 : _pV(Vi._pV), _i_comp(Vi._i_comp) {}

inline
space_component::space_component(space& V, size_type i)
 : _pV(&V), _i_comp(i) {}
inline
const_space_component::const_space_component(const space& V, size_type i)
 : _pV(&V), _i_comp(i) {}

inline
space_component
space::operator [] (size_type i_comp)
{
    return space_component (*this, i_comp);
}
inline
const_space_component
space::operator [] (size_type i_comp) const
{
    return const_space_component (*this, i_comp);
}
inline
void
space_component::block (const domain& d)
{
    (*_pV).block (d, _i_comp);
}
inline
void
space_component::block (const std::string& d_name)
{
    (*_pV).block (d_name, _i_comp);
}
// ---------------------------------------------------------

inline
space::space () 
  : smart_pointer<spacerep> (new_macro(spacerep)) {}

inline
space::space (smart_pointer<spacerep> x) 
  : smart_pointer<spacerep> (x) {}

inline
space::space(const geo& g, const std::string& approx_name, const std::string& option)
  : smart_pointer<spacerep> (new_macro(spacerep(g, approx_name, option))) {}

inline
space::space (const geo& g, const std::string& approx_name, 
	const domain& interface, const domain& subgeo, const std::string& option)
  : smart_pointer<spacerep> (new_macro(spacerep(g, approx_name, interface, subgeo, option))) {}

inline
space::space (const geo& g, const domain& d, const std::string& approx_name) 
  : smart_pointer<spacerep> (new_macro(spacerep(g, d, approx_name))) {}

inline
space::space (const space& X, const space& Y) 
  : smart_pointer<spacerep> (new_macro(spacerep(X.data(),Y.data()))) {}

inline
space::space (const space& X, size_type i_comp) 
  : smart_pointer<spacerep> (new_macro(spacerep(X.data(),i_comp))) {}

inline
space::size_type
space::size() const
{
  return data().size();
}
inline
space::size_type
space::degree() const
{
  return get_basis().degree();
}
inline
space::size_type
space::dimension() const
{
  return data().dimension();
}
inline
std::string
space::coordinate_system() const
{
  return data().coordinate_system();
}
inline
fem_helper::coordinate_type
space::coordinate_system_type() const
{
  return data().coordinate_system_type();
}
inline
space::size_type
space::n_blocked() const
{
  return data().n_blocked();
}

inline
space::size_type
space::n_unknown() const
{
  return data().n_unknown();
}
inline
space::size_type
space::size_component(space::size_type i) const 
{
  return data().size_component(i) ;
}
inline
space::size_type
space::start(space::size_type i) const 
{
  return data().start(i) ;
}
inline
space::size_type
space::start_unknown(space::size_type i) const 
{
  return data().start_unknown(i);
}
inline
space::size_type
space::start_blocked(space::size_type i) const 
{
  return data().start_blocked(i);
}
inline 
space::size_type
space::n_blocked_component(space::size_type i) const 
{
  return data().n_blocked_component(i) ;
}
inline 
space::size_type
space::n_unknown_component(space::size_type i) const 
{
  return data().n_unknown_component(i) ;
}

inline 
space::size_type
space::n_component() const 
{
  return data().n_component() ;
}

inline
space::size_type
space::index (space::size_type degree_of_freedom) const
{
  return data().index (degree_of_freedom);
}

inline
space::size_type
space::period_association (space::size_type degree_of_freedom) const
{
  return data().period_association (degree_of_freedom);
}

inline
bool
space::is_blocked (space::size_type degree_of_freedom) const
{
  return data().is_blocked (degree_of_freedom);
}

inline
bool
space::is_periodic (space::size_type degree_of_freedom) const
{
  return data().is_periodic (degree_of_freedom);
}

inline
bool
space::is_locked (space::size_type degree_of_freedom) const
{
  return data().is_locked (degree_of_freedom);
}
inline
space::size_type  
space::index_locked_with  (size_type degree_of_freedom) const
{
  return data().index_locked_with  (degree_of_freedom);
}
inline
space::size_type  
space::locked_with  (size_type degree_of_freedom) const
{
  return data().locked_with  (degree_of_freedom);
}
inline
Float
space::unlocked_component (size_type dof, size_type i_comp) const
{
  return data().unlocked_component (dof, i_comp);
}
inline
Float
space::locked_component (size_type dof, size_type i_comp) const
{
  return data().locked_component (dof, i_comp);
}
inline
bool
space::has_locked () const
{
  return data().has_locked ();
}

inline
std::string
space::get_approx (space::size_type i) const
{
  return data().get_approx(i);
}
inline
const basis&
space::get_p1_transformation () const
{
  return data().get_p1_transformation();
}
inline
fem_helper::valued_field_type 
space::get_valued_type() const
{
  return data().get_valued_type();
}
inline
std::string 
space::get_valued() const
{
  return data().get_valued();
}
inline
const geo&
space::get_geo () const
{
  return data().get_geo();
}

inline
const basis&
space::get_basis (space::size_type i) const
{
  return data().get_basis(i);
}
inline
const numbering&
space::get_numbering (space::size_type i) const
{
  return data().get_numbering(i);
}

inline
space::size_type
space::global_size() const
{
  return data().global_size();
}

inline
void 
space::block () const
{
  data().block();
}
inline
void 
space::block (size_type i_comp) const
{
  data().block(i_comp);
}
inline
void 
space::block_Lagrange () const
{
  data().block_Lagrange();
}
inline
void 
space::block_Lagrange (size_type i_comp) const
{
  data().block_Lagrange(i_comp);
}
inline
void 
space::block (const domain& d) const
{
  data().block(d);
}
inline
void 
space::block (const domain& d, size_type i_comp) const
{
  data().block(d, i_comp);
}
inline
void 
space::block (const std::string& d_name) const
{
  data().block(d_name);
}
inline
void 
space::block (const std::string& d_name, size_type i_comp) const
{
  data().block(d_name, i_comp);
}  
inline
void 
space::block_dof (size_type dof_idx) const
{
  data().block_dof(dof_idx);
}
inline
void 
space::unblock_dof (size_type dof_idx) const
{
  data().unblock_dof(dof_idx);
}

inline
void 
space::set_periodic (const domain& d1, const domain& d2) const
{
  data().set_periodic(d1,d2);
}
inline
void 
space::set_periodic (const std::string& d1_name, const std::string& d2_name) const
{
  data().set_periodic(d1_name, d2_name);
}
inline
void 
space::lock_components (const domain& d, point locked_dir) const
{
  data().lock_components(d, locked_dir);
}
inline
void 
space::lock_components (const std::string& d_name, point locked_dir) const
{
  data().lock_components(d_name, locked_dir);
}
template <class Function>
inline
void 
space::lock_components (const domain& d, Function locked_dir) const
{
  data().lock_components(d, locked_dir);
}
template <class Function>
inline
void 
space::lock_components (const std::string& d_name, Function locked_dir) const
{
  data().lock_components(d_name, locked_dir);
}

inline
space::size_type
space::domain_dof (size_type global_degree_of_freedom) const
{
  return data().domain_dof (global_degree_of_freedom);
}

inline
bool
space::is_on_boundary_domain() const
{
  return data().is_on_boundary_domain();
}
inline
const domain&
space::get_boundary_domain() const
{
  return data().get_boundary_domain();
}
inline
const geo&
space::get_global_geo () const
{
  return data().get_global_geo();
}

inline
space::size_type
space::domain_index (size_type global_degree_of_freedom) const
{
  return data().domain_index (global_degree_of_freedom);
}


inline
void
space::freeze () const
{
  data().freeze();
}
inline
bool
space::is_frozen () const
{
  return data().is_frozen();
}
inline
void
space::set_dof( const geo_element& K, tiny_vector<size_type>& idx) const
{
  data().set_dof(K, idx);
}
inline
void
space::set_dof( const geo_element& K, tiny_vector<size_type>& idx,
		size_type i_comp) const
{
  data().set_dof(K, idx, i_comp);
}
inline
void
space::set_dof( const geo_element& K, tiny_vector<size_type>& idx,
		size_type i_comp, reference_element::dof_family_type family) const
{
  data().set_dof(K, idx, i_comp, family);
}
inline
void
space::set_global_dof(
		const geo_element& K, 
		tiny_vector<size_type>& idx) const
{
  data().set_global_dof(K, idx);
}
inline
void
space::set_global_dof(
		const geo_element& K, 
		tiny_vector<size_type>& idx,
		size_type i_comp) const
{
  data().set_global_dof(K, idx, i_comp);
}
inline
point
space::x_dof(const geo_element& K, size_type iloc) const
{
  return data().x_dof(K, iloc) ;
}
inline
void   
//! Hermite/Argyris: returns 1 if the global dof contains the derivative along outward normal of K, -1 else  
space::dof_orientation(const geo_element& K, tiny_vector<int>& orientation) const
{
  return data().dof_orientation(K, orientation);
}

inline
space 
operator * (const space& X, const space& Y)
{
  return space(X,Y);
}

inline
std::ostream&
operator << (std::ostream& s, const space& x)
{
  return s << x.data();
}
inline
std::istream&
operator >> (std::istream& s, space& x)
{
  return s >> x.data();
}
inline
std::ostream& 
space::dump(std::ostream& s) const
{
    return data().dump(s);
}
inline
void
space::check() const
{
    data().check();
}
inline
bool
space::operator == (const space& y) const
{
    return data().operator== (y.data());
}
inline
bool
space::operator != (const space& y) const
{
    return !operator== (y);
}
/* static */
inline
space::size_type
space::inquire_size(const geo& g, const numbering& numb)
{
    return spacerep::inquire_size(g,numb);
}
inline
meshpoint 
space::hatter (const point& x, size_type K) const
{
    return data().hatter(x,K);
}
inline
point 
space::dehatter (const meshpoint& S) const
{
    return data().dehatter(S);
}
inline
point
space::dehatter (const point& x_hat, size_type e) const
{
    return data().dehatter(meshpoint(e,x_hat));
}
}// namespace rheolef
# endif /* _RHEO_SPACE_H */