/usr/include/rheolef/piola_algo

#ifndef _RHEO_PIOLA_ALGO_H
#define _RHEO_PIOLA_ALGO_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
/// 
#include "rheolef/tensor.h"
#include "rheolef/geo.h"
namespace rheolef { 
// ---------------------------------------------------------------------------
// piola transform and its inverse on simplex
// ---------------------------------------------------------------------------
// piola transformation
// F_K : hat_K --> K
//       hat_x --> x = F_K(hat_x)
// ------------------------------------------
// TODO: should be merged with
//    piola.h  = general piola on all elements (non-simplex)
// this file is actually used only by ZZ posteriori error estimators
//
inline
point
piola_e (
  const point& hat_x, 
  const point& a, 
  const point& b,
  bool is_a_vector = false) 
{
  return point((!is_a_vector?a[0]:0)
               + (b[0]-a[0])*hat_x[0]);
}
inline
point
piola_t (
  const point& hat_x, 
  const point& a, 
  const point& b,
  const point& c,
  bool is_a_vector = false)
{
  return ((!is_a_vector)?a:point(0,0))
          + (b-a)*hat_x[0]
	  + (c-a)*hat_x[1];
}
inline
point
piola_T (
  const point& hat_x, 
  const point& a, 
  const point& b,
  const point& c,
  const point& d,
  bool is_a_vector = false)
{
  return ((!is_a_vector)?a:point(0,0))
          + (b-a)*hat_x[0] 
          + (c-a)*hat_x[1]
          + (d-a)*hat_x[2];
}
// ------------------------------------------
// piola transformation matrix: M_K where
//  F_K(hat_x) = M_K*hat_x + t_K
// ------------------------------------------
inline
void
set_piola_matrix_e (
  tensor& m,
  const point& a, 
  const point& b)
{
  m(0,0) = b[0]-a[0];
}
inline
void
set_piola_matrix_t (
  tensor& m,
  const point& a, 
  const point& b,
  const point& c)
{
   m.set_column (b-a, 0, 2);
   m.set_column (c-a, 1, 2);
}
inline
void
set_piola_matrix_T (
  tensor& m,
  const point& a, 
  const point& b,
  const point& c,
  const point& d)
{
   m.set_column (b-a, 0, 3);
   m.set_column (c-a, 1, 3);
   m.set_column (d-a, 2, 3);
}
inline
void
set_piola_matrix (
  tensor& M_K,
  const geo& omega,
  const geo_element& K)
{
    switch (K.type()) {
      case reference_element::e:
	set_piola_matrix_e (M_K,
	   omega.vertex(K[0]),
	   omega.vertex(K[1]));
	break;
      case reference_element::t:
	set_piola_matrix_t (M_K,
	   omega.vertex(K[0]),
	   omega.vertex(K[1]),
	   omega.vertex(K[2]));
	break;
      case reference_element::T:
	set_piola_matrix_T (M_K,
	   omega.vertex(K[0]),
	   omega.vertex(K[1]),
	   omega.vertex(K[2]),
	   omega.vertex(K[3]));
	break;
      default:
	error_macro ("not yet supported element");
    }
}
// ------------------------------------------
// inverse piola transformation
// F_K^{-1} : K --> hat(K)
//            x --> hat(x)
// ------------------------------------------
inline
point 
inv_piola_e (
  const point& x, 
  const point& a, 
  const point& b) 
{
  return point((x[0]-a[0])/(b[0]-a[0]));
}
inline
point 
inv_piola_t (
  const point& x, 
  const point& a, 
  const point& b, 
  const point& c) 
{
  Float t9 = 1/(-b[0]*c[1]+b[0]*a[1]+a[0]*c[1]+c[0]*b[1]-c[0]*a[1]-a[0]*b[1]);
  Float t11 = -a[0]+x[0];
  Float t15 = -a[1]+x[1];
  return point((-c[1]+a[1])*t9*t11-(-c[0]+a[0])*t9*t15,
               (b[1]-a[1])*t9*t11-(b[0]-a[0])*t9*t15);
}
inline
point 
inv_piola_T (
  const point& x, 
  const point& a, 
  const point& b, 
  const point& c, 
  const point& d) 
{
  tensor A;
  point ax;
  for (size_t i = 0; i < 3; i++) {
    ax[i]  = x[i]-a[i];
    A(i,0) = b[i]-a[i];
    A(i,1) = c[i]-a[i];
    A(i,2) = d[i]-a[i];
  }
  tensor inv_A;
  bool is_singular = ! invert_3x3 (A, inv_A);
  check_macro(!is_singular, "inv_piola: singular transformation in tetrahedron");
  point hat_x = inv_A*ax;
  return hat_x; 
}
}// namespace rheolef
#endif // _RHEO_PIOLA_ALGO_H
librheolef-dev 5.93-2 / usr / include / rheolef / piola_algo_v1.h
