librheolef-dev 6.7-6 / usr / include / rheolef / geo_element_curved_ball.h

#ifndef _RHEOLEF_GEO_ELEMENT_CURVED_BALL_H
#define _RHEOLEF_GEO_ELEMENT_CURVED_BALL_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
///
/// =========================================================================
//
// transfinite piola transformaton from
// * a triangle (a,b,c) or a quadrangle(a,b,c,d)
//   to a curved triangle or quadrangle
//   where edge (a,b) is an arc of a circle.
// * a tetra (a,b,c,d) or an hexa(a,b,c,d,e,f,g,h)
//   to a curved tetra or hexa
//   where face (a,b,c), or (a,b,c,d) for hexa,
//   is a part of sphere boundary
//
// Remark: can be generalized from a circle/ball
//   to any curved boundary, by sending a function
//   that project or replace internal nodes onto the boundary.
//   See SherKar-2005, page 224.
//
// author: Pierre.Saramito@imag.fr
//
// date: 22 november 2011
// 
// TODO:
//  * tetra, hexa
//  * ellipse(r1,r2,r3,c)
//  * parametrized by a function that project on the bdry ?
//
#include "point.h"

namespace rheolef {

// ====================================================================================
// quadrangle (a,b,c,d) with curved edge (a,b)
// ====================================================================================
template<class T>
class curved_ball_q {
public:
// allocator:
  curved_ball_q (const point_basic<T>& a0, const point_basic<T>& b0, const point_basic<T>& c0, const point_basic<T>& d0,
	const point_basic<T>& center0 = point_basic<T>(0,0), const T& radius0 = 1)
   : a(a0), b(b0), c(c0), d(d0), center(center0), radius(radius0) {}
// accessor:
  point_basic<T> operator() (const point_basic<T>& hat_x) const {
     // transform large square [-1:1]^2 into small one [0:1]^2
     T x = 0.5*(hat_x[0]+1);
     T y = 0.5*(hat_x[1]+1);
     return f_small_carre (x, y);
  }
// internals:
protected:
  point_basic<T> f_ab (const T& s) const {
	point_basic<T> x = (1-s)*a + s*b;
	return center + radius*(x-center)/norm(x-center);
  }
  point_basic<T> f_bc (const T& s) const { return (1-s)*b + s*c; }
  point_basic<T> f_dc (const T& s) const { return (1-s)*d + s*c; }
  point_basic<T> f_ad (const T& s) const { return (1-s)*a + s*d; }
  /*
    Transformation for a square: from eqn (17) in:
    @article{GorHal-1972,
     author = {W. J. Gordon and C. A. Hall},
     title = {Transfinite element methods:
            blending-function interpolation over
            arbitrary curved element domains},
     journal = {Numer. Math.},
     volume = 21,
     pages = {109--129},
     year = 1973,
     url = {http://www.springerlink.com/content/p5517206551228j6/fulltext.pdf},
     keywords = {method!fem!curved}
    }
  */
  point_basic<T> f_small_carre (const T& x, const T& y) const {
    // small carre: 0 <= x,y <= 1
    return (1-x)*f_ad(y) + x*f_bc(y)
         + (1-y)*f_ab(x) + y*f_dc(x)
           - (1-x)*(1-y)*a
           -    x *(1-y)*b
           -    x *   y *c
           - (1-x)*   y *d
  	;
  }
protected:
// data:
  point_basic<T> a, b, c, d, center;
  T radius;
};
#ifdef TO_CLEAN
// ====================================================================================
// triangle (a,b,c) with curved edge (a,b)
// ====================================================================================
template<class T>
class curved_ball_t_old : public curved_ball_q {
public:
// allocator:
  curved_ball_t_old (const point_basic<T>& a0, const point_basic<T>& b0, const point_basic<T>& c0,
	const point_basic<T>& center0 = point_basic<T>(0,0), const T& radius0 = 1)
   : curved_ball_q (a0, b0, c0, c0, center0, radius0) {}
// accessor:
  point_basic<T> operator() (const point_basic<T>& hat_x) const {
    /* Transform for a triangle, from eqn (27) in:
	@article{DeySheFla-1997,
	 title = {Geometry representation issues associated with p-version finite element computations},
	 author = {S. Dey and M. S. Shephard and J. E. Flaherty},
	 journal = {Comput. Meth. Appl. Mech. Engrg.},
	 volume = 150, 
	 year = 1997, 
	 pages = {39--55}
	}
      Remark : SheKar-1999, p. 162-163 : cite GorHal-1973 for a quadrangle
      but says that the triangle-to-square transform + the quadrangle formulae
      does not work here (its true, I have checked) because of a singular jacobian.
      Triangle requires thus a specific formulae as it is done here.
      Also: similar ref: eqn (3.17), page 169 in SolSegDol-2004 is incorrect.
    */
    T la = 1 - hat_x[0] - hat_x[1]; // barycentric coords
    T lb = hat_x[0];
    T lc = hat_x[1];
    return 0.5*( la/(1-lb)*f_ab(lb) + lb/(1-la)*f_ab(1-la)
               + lb/(1-lc)*f_bc(lc) + lc/(1-lb)*f_bc(1-lb) 
	       + lc/(1-la)*f_ca(la) + la/(1-lc)*f_ca(1-lc)
	       - la*a - lb*b - lc*c);
  }
// internals:
protected:
  point_basic<T> f_ca (const T& s) const { return f_ad(1-s); }
};
#endif // TO_CLEAN
// ====================================================================================
// triangle (a,b,c) with i-th curved edge
// ====================================================================================
template<class T>
class curved_ball_t {
public:
// allocator:
  curved_ball_t (const point_basic<T>& a0, const point_basic<T>& b0, const point_basic<T>& c0,
	size_t loc_curved_iedg, const point_basic<T>& center0 = point_basic<T>(0,0), const T& radius0 = 1)
   : node(), center(center0), radius(radius0), is_bdry_edg()
  {
    node[0] = a0;
    node[1] = b0;
    node[2] = c0;
    is_bdry_edg.fill (false);
    is_bdry_edg [loc_curved_iedg] = true;
  }
// accessor:
  point_basic<T> operator() (const point_basic<T>& hat_x) const {
    /* Transform for a triangle, from eqn (27) in:
	@article{DeySheFla-1997,
	 title = {Geometry representation issues associated with p-version finite element computations},
	 author = {S. Dey and M. S. Shephard and J. E. Flaherty},
	 journal = {Comput. Meth. Appl. Mech. Engrg.},
	 volume = 150, 
	 year = 1997, 
	 pages = {39--55}
	}
      Remark : SheKar-1999, p. 162-163 : cite GorHal-1973 for a quadrangle
      but says that the triangle-to-square transform + the quadrangle formulae
      does not work here (its true, I have checked) because of a singular jacobian.
      Triangle requires thus a specific formulae as it is done here.
      Also: similar ref: eqn (3.17), page 169 in SolSegDol-2004 is incorrect.
    */
    T la = 1 - hat_x[0] - hat_x[1]; // barycentric coords
    T lb = hat_x[0];
    T lc = hat_x[1];
    return 0.5*( la/(1-lb)*edge(0,lb) + lb/(1-la)*edge(0,1-la)
               + lb/(1-lc)*edge(1,lc) + lc/(1-lb)*edge(1,1-lb) 
	       + lc/(1-la)*edge(2,la) + la/(1-lc)*edge(2,1-lc)
	       - la*node[0] - lb*node[1] - lc*node[2]);
  }
protected:
// internals:
  point_basic<T> project_on_boundary (const point_basic<T>& x) const {
    return center + radius*(x-center)/norm(x-center);
  }
  point_basic<T> edge (size_t loc_iedg, const T& hat_x) const {
    const point_basic<T>& a = node [loc_iedg];
    const point_basic<T>& b = node [(loc_iedg+1)%3];
    point_basic<T> x = (1-hat_x)*a + hat_x*b;
    if (is_bdry_edg [loc_iedg]) {
      return project_on_boundary (x);
    } else {
      return x;
    }
  }
// data:
  std::array<point_basic<T>,3> node;
  point_basic<T> center;
  T radius;
  std::array<bool,3> is_bdry_edg;
};
#ifdef TO_CLEAN
// ====================================================================================
// tetrahedron (a,b,c,d) with curved face (a,b,c)
// ====================================================================================
template<class T>
class curved_ball_T_old {
public:
// allocator:
  curved_ball_T_old (const point_basic<T>& a0, const point_basic<T>& b0, const point_basic<T>& c0, const point_basic<T>& d0,
	const point_basic<T>& center0 = point_basic<T>(0,0), const T& radius0 = 1)
   : a(a0), b(b0), c(c0), d(d0), center(center0), radius(radius0) {}
// accessor:
  point_basic<T> operator() (const point_basic<T>& hat_x) const {
    T la = 1 - hat_x[0] - hat_x[1] - hat_x[2]; // barycentric coords
    T lb = hat_x[0];
    T lc = hat_x[1];
    T ld = hat_x[2];
 
    curved_ball_t_old G_abd (a, b, d, center, radius);
    curved_ball_t_old G_acd (a, c, d, center, radius);
    curved_ball_t_old G_bcd (b, c, d, center, radius);

    point_basic<T> s_abc = f_abc((la*a + lb*b + lc*c)/(la+lb+lc));
    point_basic<T> s_abd = G_abd(point_basic<T>(lb, ld)/(la+lb+ld));
    point_basic<T> s_acd = G_acd(point_basic<T>(lc, ld)/(la+lc+ld));
    point_basic<T> s_bcd = G_bcd(point_basic<T>(lc, ld)/(lb+lc+ld));
    point_basic<T> e_ab  = f_ab ((la*a + lb*b)/(la+lb));
    point_basic<T> e_ac  = f_ac ((la*a + lc*c)/(la+lc));
    point_basic<T> e_ad  = f_ad ((la*a + ld*d)/(la+ld));
    point_basic<T> e_bc  = f_bc ((lb*b + lc*c)/(lb+lc));
    point_basic<T> e_cd  = f_cd ((lc*c + ld*d)/(lc+ld));
    point_basic<T> e_bd  = f_bd ((lb*b + ld*d)/(lb+ld));
    return (1-ld)*s_abc
         + (1-lc)*s_abd
         + (1-lb)*s_acd
         + (1-la)*s_bcd
         - (1-lc-ld)*e_ab
         - (1-ld-lb)*e_ac
         - (1-lb-lc)*e_ad
         - (1-la-ld)*e_bc
         - (1-la-lb)*e_cd
         - (1-la-lc)*e_bd
         + la*a + lb*b + lc*c + ld*d;
  }
// internals:
protected:
    /* Transform for a tetra, from eqn (30) in:
	@article{DeySheFla-1997,
	 title = {Geometry representation issues associated with p-version finite element computations},
	 author = {S. Dey and M. S. Shephard and J. E. Flaherty},
	 journal = {Comput. Meth. Appl. Mech. Engrg.},
	 volume = 150, 
	 pages = {39--55},
	 year = 1997
	}
    */
  point_basic<T> project_on_boundary (const point_basic<T>& x) const {
    return center + radius*(x-center)/norm(x-center);
  }
  point_basic<T> f_abc (const point_basic<T>& x) const { return project_on_boundary (x); }
  point_basic<T> f_abd (const point_basic<T>& x) const { return x; }
  point_basic<T> f_acd (const point_basic<T>& x) const { return x; }
  point_basic<T> f_bcd (const point_basic<T>& x) const { return x; }
  point_basic<T> f_ab  (const point_basic<T>& x) const { return project_on_boundary (x); }
  point_basic<T> f_ac  (const point_basic<T>& x) const { return project_on_boundary (x); }
  point_basic<T> f_bc  (const point_basic<T>& x) const { return project_on_boundary (x); }
  point_basic<T> f_ad  (const point_basic<T>& x) const { return x; }
  point_basic<T> f_bd  (const point_basic<T>& x) const { return x; }
  point_basic<T> f_cd  (const point_basic<T>& x) const { return x; }
// data:
  point_basic<T> a, b, c, d, center;
  T radius;
};
#endif // TO_CLEAN
// ====================================================================================
// tetrahedron (a,b,c,d) with i-th curved face
// ====================================================================================
template<class T>
class curved_ball_T {
public:
// allocator:
  curved_ball_T (const point_basic<T>& a0, const point_basic<T>& b0, const point_basic<T>& c0, const point_basic<T>& d0,
	const point_basic<T>& center0 = point_basic<T>(0,0), const T& radius0 = 1)
   : node(), center(center0), radius(radius0), is_bdry_edg(), is_bdry_fac(), is_curved_fac()
  {
    node[0] = a0;
    node[1] = b0;
    node[2] = c0;
    node[3] = d0;
    is_bdry_fac.fill(false);
    is_bdry_edg.fill(false);
    is_curved_fac.fill(false);
  }
// modifiers:
  void set_boundary_face (size_t loc_ifac_curved)
  {
    is_bdry_fac[loc_ifac_curved] = true;
    // then, all faces are either on boundary or have a curved edge:
    for (size_t loc_ifac = 0; loc_ifac < 4; loc_ifac++) {
      if (!is_bdry_fac [loc_ifac]) {
        is_curved_fac[loc_ifac] = true;
      }
    }
    // the three edges of the bdry face are also on the boundary
    for (size_t loc_ifac_jedg = 0; loc_ifac_jedg < 3; loc_ifac_jedg++) {
      size_t loc_jedg = reference_element_T::face2edge (loc_ifac_curved, loc_ifac_jedg);
      is_bdry_edg [loc_jedg] = true;
    }
  }
  void set_boundary_edge (size_t loc_iedg_curved)
  {
    is_bdry_edg [loc_iedg_curved] = true;
    // also, the two faces that contains this edge are then curved (if not yet on the boundary)
    for (size_t loc_ifac = 0; loc_ifac < 4; loc_ifac++) {
      for (size_t loc_ifac_jedg = 0; loc_ifac_jedg < 3; loc_ifac_jedg++) {
        size_t loc_jedg = reference_element_T::face2edge (loc_ifac, loc_ifac_jedg);
        if (loc_jedg != loc_iedg_curved) continue;
	// this face contains the boundary edge:
        if (!is_bdry_fac [loc_ifac]) {
          is_curved_fac[loc_ifac] = true;
          break;
        }
      }
    }
  }
// accessor:
  point_basic<T> operator() (const point_basic<T>& hat_x) const {
    /* Transform for a tetra, from eqn (30) in:
	@article{DeySheFla-1997,
	 title = {Geometry representation issues associated with p-version finite element computations},
	 author = {S. Dey and M. S. Shephard and J. E. Flaherty},
	 journal = {Comput. Meth. Appl. Mech. Engrg.},
	 volume = 150, 
	 pages = {39--55},
	 year = 1997
	}
    */
    // 1) compute the baycentric coordinates: lambda_i(xj) = delta_ij
    std::array<T,4> l;
    l[0] = 1 - hat_x[0] - hat_x[1] - hat_x[2];
    l[1] = hat_x[0];
    l[2] = hat_x[1];
    l[3] = hat_x[2];

    // 2) compute the contribution of all faces
    std::array<size_t,3> S;       // table of local vertices indexes
    std::array<bool,4>   is_on_S; // when vertex is on S
    point_basic<T> x_all_faces (0,0,0);
    for (size_t loc_ifac = 0; loc_ifac < 4; loc_ifac++) {
      is_on_S.fill (false);
      for (size_t loc_ifac_jv = 0; loc_ifac_jv < 3; loc_ifac_jv++) {
        S[loc_ifac_jv] = reference_element_T::subgeo_local_node (1, 2, loc_ifac, loc_ifac_jv);
        is_on_S [S[loc_ifac_jv]] = true;
      }
      // jv_3 is the only index in 0..3 that is not on the face ifac:
      size_t jv_3 = size_t(-1);
      for (size_t loc_jv = 0; loc_jv < 4; loc_jv++) {
        if (!is_on_S [loc_jv]) jv_3 = loc_jv;
      }
      T sum = l[S[0]] + l[S[1]] + l[S[2]];
      point_basic<T> hat_xi (l[S[1]]/sum, l[S[2]]/sum);
      x_all_faces += (1 - l[jv_3])*x_face (loc_ifac, hat_xi);
    }
    // 2) compute the contribution of all edges
    point_basic<T> x_all_edges (0,0,0);
    for (size_t loc_iedg = 0; loc_iedg < 6; loc_iedg++) {
      // edge = [a,b] and c & d are the two others verices of the tetra
      size_t ia = reference_element_T::subgeo_local_node (1, 1, loc_iedg, 0);
      size_t ib = reference_element_T::subgeo_local_node (1, 1, loc_iedg, 1);
      size_t ic = size_t(-1);
      size_t id = size_t(-1);
      for (size_t loc_iv = 0; loc_iv < 4; loc_iv++) {
        if (loc_iv == ia || loc_iv == ib) continue;
        if (ic == size_t(-1)) { ic = loc_iv; continue; }
        id = loc_iv;
      }
      assert_macro (ic != size_t(-1) && id != size_t(-1), "unexpected");
      T hat_xi = l[ib]/(l[ia] + l[ib]);
      x_all_edges += (1 - l[ic] - l[id])*x_edge (loc_iedg, hat_xi);
    }
    // 3) compute the contribution of all vertices
    point_basic<T> x_all_vertices (0,0,0);
    for (size_t loc_iv = 0; loc_iv < 4; loc_iv++) {
      x_all_vertices += l[loc_iv]*node[loc_iv];
    }
    return x_all_faces - x_all_edges + x_all_vertices;
  }
// internals:
protected:
  point_basic<T> x_face (size_t loc_ifac, const point_basic<T>& hat_x) const {
    size_t ia = reference_element_T::subgeo_local_node (1, 2, loc_ifac, 0);
    size_t ib = reference_element_T::subgeo_local_node (1, 2, loc_ifac, 1);
    size_t ic = reference_element_T::subgeo_local_node (1, 2, loc_ifac, 2);
    const point_basic<T>& a = node[ia];
    const point_basic<T>& b = node[ib];
    const point_basic<T>& c = node[ic];
    if (is_curved_fac [loc_ifac]) {
      // search the edge who lives on the boundary
      size_t loc_ifac_jedg_bdry = size_t(-1);
      for (size_t loc_ifac_jedg = 0; loc_ifac_jedg < 3; loc_ifac_jedg++) {
        size_t loc_jedg = reference_element_T::face2edge (loc_ifac, loc_ifac_jedg);
        if (is_bdry_edg [loc_jedg]) {
          loc_ifac_jedg_bdry = loc_ifac_jedg;
          break;
        }
      }
      assert_macro (loc_ifac_jedg_bdry != size_t(-1), "unexpected");
      curved_ball_t<T> F (a, b, c, loc_ifac_jedg_bdry, center, radius);
      return F (hat_x);
    } else {
      point_basic<T> x = (1 - hat_x[0] - hat_x[1])*a + hat_x[0]*b + hat_x[1]*c;
      if (is_bdry_fac [loc_ifac]) {
        return project_on_boundary (x);
      } else { // strait face
        return x;
      }
    }
  }
  point_basic<T> x_edge (size_t loc_iedg, const T& hat_x) const {
    size_t ia = reference_element_T::subgeo_local_node (1, 1, loc_iedg, 0);
    size_t ib = reference_element_T::subgeo_local_node (1, 1, loc_iedg, 1);
    const point_basic<T>& a = node[ia];
    const point_basic<T>& b = node[ib];
    point_basic<T> x = (1 - hat_x)*a + hat_x*b;
    if (is_bdry_edg [loc_iedg]) {
      return project_on_boundary (x);
    } else {
      return x;
    }
  }
  point_basic<T> project_on_boundary (const point_basic<T>& x) const {
	return center + radius*(x-center)/norm(x-center);
  }
  point_basic<T> f_abc (const point_basic<T>& x) const { return project_on_boundary (x); }
  point_basic<T> f_abd (const point_basic<T>& x) const { return x; }
  point_basic<T> f_acd (const point_basic<T>& x) const { return x; }
  point_basic<T> f_bcd (const point_basic<T>& x) const { return x; }
  point_basic<T> f_ab  (const point_basic<T>& x) const { return project_on_boundary (x); }
  point_basic<T> f_ac  (const point_basic<T>& x) const { return project_on_boundary (x); }
  point_basic<T> f_bc  (const point_basic<T>& x) const { return project_on_boundary (x); }
  point_basic<T> f_ad  (const point_basic<T>& x) const { return x; }
  point_basic<T> f_bd  (const point_basic<T>& x) const { return x; }
  point_basic<T> f_cd  (const point_basic<T>& x) const { return x; }
// data:
  std::array<point_basic<T>,4> node;
  point_basic<T> center;
  T radius;
  std::array<bool,6> is_bdry_edg;
  std::array<bool,4> is_bdry_fac;
  std::array<bool,4> is_curved_fac;
};
// ====================================================================================
// TODO: hexahedron with a curved face
// ====================================================================================
template<class T>
class curved_ball_H {
public:
// allocator:
  curved_ball_H (const point_basic<T>& a0, const point_basic<T>& b0, const point_basic<T>& c0, const point_basic<T>& d0,
                 const point_basic<T>& e0, const point_basic<T>& f0, const point_basic<T>& g0, const point_basic<T>& h0,
	const point_basic<T>& center0 = point_basic<T>(0,0), const T& radius0 = 1)
   : a(a0), b(b0), c(c0), d(d0), e(e0), f(f0), g(g0), h(h0), center(center0), radius(radius0) {}
// accessor:
  point_basic<T> operator() (const point_basic<T>& z) const {
     error_macro ("curved H: not yet");
     return point_basic<T>();
  }
// internals:
protected:
    /* Transform for an hexa, from annex A in:
	@article{KirSza-1997, 
	 author = {G. Kir\'alyfalvi and B. A. Szab\'o},
	 title = {Quasi-regional mapping for the $p$-version
		of the finite element method},
	 journal = {Finite Elements in Analysis and Design},
	 volume = 27, 
	 pages = {85--97},
	 year = 1997, 
	 keywords = {method!fem!curved!blending}
	}
    */
// data:
  point_basic<T> a, b, c, d, e, f, g, h, center;
  T radius;
};

} // namespace rheolef
#endif // _RHEOLEF_GEO_ELEMENT_CURVED_BALL_H