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

#ifndef _RHEOLEF_INV_PIOLA_H
#define _RHEOLEF_INV_PIOLA_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
/// 
/// =========================================================================
// 
// invert piola tranformation on nonlinear elements :
//  quadrangles, high-order, etc
// by using the newton generic algorithm
//
#include "rheolef/geo.h"
#include "rheolef/piola.h"
namespace rheolef {

template<class T>
class inv_piola {
public:
  typedef T                              float_type;
  typedef point_basic<T>                 value_type;
  typedef typename value_type::size_type size_type;
  inv_piola();
#ifdef TO_CLEAN
  template<class M>
  void reset (const geo_basic<T,M>& omega, const geo_element& K);
#endif // TO_CLEAN
  template<class M>
  void reset (const geo_basic<T,M>& omega, const reference_element& hat_K, const std::vector<size_t>& dis_inod);
  void set_x (const value_type& x1) { x0 = x1; }
  value_type initial() const;
  value_type residue     (const value_type& hat_xh) const;
  void update_derivative (const value_type& hat_xh) const;
  value_type derivative_solve      (const value_type& r) const;
  value_type derivative_trans_mult (const value_type& r) const;
  float_type space_norm       (const value_type& hat_xh) const;
  float_type dual_space_norm  (const value_type& r) const;
  float_type duality_product  (const value_type& r, const value_type& s) const;
protected:
  size_t                  dim, map_dim;
  basis_basic<T>          b;
  reference_element       hat_K;
  std::vector<value_type> node;
  value_type              x0;
  mutable std::vector<float_type> value;
  mutable std::vector<value_type> grad_value;
  mutable tensor_basic<T>         DF, inv_DF;
};
template<class T>
inv_piola<T>::inv_piola()
{
}
#ifdef TO_CLEAN
template<class T>
template<class M>
void
inv_piola<T>::reset (const geo_basic<T,M>& omega, const geo_element& K) {
  dim = omega.dimension(); 
  b = omega.get_piola_basis();
  hat_K   = K.variant();
  map_dim = K.dimension(); 
  node.resize (K.n_node());
  for (size_t loc_inod = 0, loc_nnod = K.n_node(); loc_inod < loc_nnod; ++loc_inod) {
    node[loc_inod] = omega.dis_node (K[loc_inod]);
  }
}
#endif // TO_CLEAN
template<class T>
template<class M>
void
inv_piola<T>::reset (const geo_basic<T,M>& omega, const reference_element& hat_K1, const std::vector<size_t>& dis_inod) {
  dim = omega.dimension(); 
  b = omega.get_piola_basis();
  hat_K   = hat_K1;
  map_dim = hat_K1.dimension(); 
  node.resize (dis_inod.size());
  for (size_t loc_inod = 0, loc_nnod = node.size(); loc_inod < loc_nnod; ++loc_inod) {
    node[loc_inod] = omega.dis_node (dis_inod[loc_inod]);
  }
}
template<class T>
typename inv_piola<T>::value_type
inv_piola<T>::initial() const {
  switch (hat_K.variant()) {
    case reference_element::e : return value_type(0.5);
    case reference_element::t : return value_type(1/float_type(3),1/float_type(3));
    case reference_element::q : return value_type(0,0);
    case reference_element::T : return value_type(1/float_type(3),1/float_type(3),1/float_type(3));
    case reference_element::P : return value_type(1/float_type(3),1/float_type(3),0);
    case reference_element::H : return value_type(0,0,0);
  }
  return value_type(0);
}
template<class T>
typename inv_piola<T>::value_type
inv_piola<T>::residue (const value_type& hat_x) const {
  b.eval (hat_K, hat_x, value);
  value_type r;
  for (size_t loc_inod = 0, loc_nnod = node.size(); loc_inod < loc_nnod; ++loc_inod) {
    r = r + value[loc_inod]*node[loc_inod];
  }
  return r - x0;
}
template<class T>
void
inv_piola<T>::update_derivative (const value_type& hat_x) const {
  b.grad_eval (hat_K, hat_x, grad_value);
  DF.reset();
  for (size_t loc_inod = 0, loc_nnod = node.size(); loc_inod < loc_nnod; ++loc_inod) {
    cumul_otimes(DF, node[loc_inod], grad_value[loc_inod], dim, map_dim);
  }
  inv_DF = pseudo_inverse_jacobian_piola_transformation (DF, dim, map_dim);
}
template<class T>
typename inv_piola<T>::value_type
inv_piola<T>::derivative_solve (const value_type& r) const {
  return inv_DF*r;
}
template<class T>
typename inv_piola<T>::float_type
inv_piola<T>::dual_space_norm (const value_type& r) const {
  return norm(r);
}
template<class T>
typename inv_piola<T>::float_type
inv_piola<T>::space_norm (const value_type& hat_xh) const {
  return norm(hat_xh);
}
template<class T>
typename inv_piola<T>::float_type
inv_piola<T>::duality_product (const value_type& r, const value_type& s) const {
  return dot (r, s);
}
template<class T>
typename inv_piola<T>::value_type
inv_piola<T>::derivative_trans_mult (const value_type& r) const {
  return DF.trans_mult(r);
}

}// namespace rheolef
#endif // _RHEOLEF_PIOLA_H