Logo  0.95.0-final
Finite Element Embedded Library and Language in C++
Feel++ Feel++ on Github Feel++ community
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
hermite.hpp
Go to the documentation of this file.
1 /* -*- mode: c++; coding: utf-8; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; show-trailing-whitespace: t -*- vim:fenc=utf-8:ft=tcl:et:sw=4:ts=4:sts=4
2 
3  This file is part of the Feel library
4 
5  Author(s): Christophe Prud'homme <christophe.prudhomme@feelpp.org>
6  Date: 2010-03-04
7 
8  Copyright (C) 2010 Université Joseph Fourier (Grenoble I)
9 
10  This library is free software; you can redistribute it and/or
11  modify it under the terms of the GNU Lesser General Public
12  License as published by the Free Software Foundation; either
13  version 2.1 of the License, or (at your option) any later version.
14 
15  This library is distributed in the hope that it will be useful,
16  but WITHOUT ANY WARRANTY; without even the implied warranty of
17  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18  Lesser General Public License for more details.
19 
20  You should have received a copy of the GNU Lesser General Public
21  License along with this library; if not, write to the Free Software
22  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
23 */
29 #ifndef __hermite_H
30 #define __hermite_H 1
31 
32 #include <boost/ptr_container/ptr_vector.hpp>
33 
34 #include <feel/feelcore/feel.hpp>
35 #include <feel/feelcore/traits.hpp>
36 #include <feel/feelalg/lu.hpp>
37 
41 
42 
47 #include <feel/feelpoly/fe.hpp>
48 namespace Feel
49 {
50 
51 namespace fem
52 {
53 
55 namespace details
56 {
57 template<typename Basis, template<class, uint16_type, class> class PointSetType>
58 class HermiteDual
59  :
60 public DualBasis<Basis>
61 {
62  typedef DualBasis<Basis> super;
63 public:
64 
65  static const uint16_type nDim = super::nDim;
66  static const uint16_type nOrder= super::nOrder;
67 
68  typedef typename super::primal_space_type primal_space_type;
69  typedef typename primal_space_type::value_type value_type;
70  typedef typename primal_space_type::points_type points_type;
71  typedef typename primal_space_type::matrix_type matrix_type;
72  typedef typename primal_space_type::convex_type convex_type;
73  typedef typename primal_space_type::reference_convex_type reference_convex_type;
74  typedef typename reference_convex_type::node_type node_type;
75 
76  // point set type associated with the functionals
77  typedef PointSetType<convex_type, 3, value_type> pointset_type;
78 
79 
80 #if 0
81 
85  typedef mpl::vector_c<uint16_type, 0, 1, 3, ( 4 ) + ( 4 ) - 2> edges_t;
86  typedef mpl::vector_c<uint16_type, 0, 0, 1, 4> geo_faces_t;
87  typedef mpl::vector_c<uint16_type, 0, 2, 3, 4> faces_t;
88  typedef mpl::vector_c<uint16_type, 0, 2, 3, 4> normals_t;
89 #endif
90 
91  static const uint16_type nVertices = reference_convex_type::numVertices;
92  static const uint16_type nFaces = reference_convex_type::numFaces;
93  static const uint16_type nGeometricFaces = reference_convex_type::numFaces;
94  static const uint16_type nEdges = reference_convex_type::numEdges;
95  static const uint16_type nNormals = reference_convex_type::numNormals;
96 
97  static const uint16_type nbPtsPerVertex = ( nDim+1 ); //reference_convex_type::nbPtsPerVertex;
98  static const uint16_type nbPtsPerEdge = 0;//reference_convex_type::nbPtsPerEdge;
99  static const uint16_type nbPtsPerFace = ( nDim==2 )?1:0; //reference_convex_type::nbPtsPerFace;
100  static const uint16_type nbPtsPerVolume = 0;//reference_convex_type::nbPtsPerVolume;
101  static const uint16_type numPoints = nVertices*nbPtsPerVertex+nGeometricFaces*nbPtsPerFace;
102 
103 
105  static const uint16_type nDofPerVertex = nbPtsPerVertex;
106 
108  static const uint16_type nDofPerEdge = 0;//(nDim+1)*nbPtsPerEdge;
109 
111  static const uint16_type nDofPerFace = nbPtsPerFace;
112 
114  static const uint16_type nDofPerVolume = 0;//(nDim+1)*nbPtsPerVolume;
115 
117  static const uint16_type nLocalDof = nDofPerVertex*nVertices + nDofPerFace*nGeometricFaces;
118 
119  static const uint16_type nFacesInConvex = mpl::if_< mpl::equal_to<mpl::int_<nDim>, mpl::int_<1> >,
120  mpl::int_<nVertices>,
121  typename mpl::if_<mpl::equal_to<mpl::int_<nDim>, mpl::int_<2> >,
122  mpl::int_<nEdges>,
123  mpl::int_<nFaces> >::type >::type::value;
124 
125  HermiteDual( primal_space_type const& primal )
126  :
127  super( primal ),
128  M_convex_ref(),
129  M_eid( M_convex_ref.topologicalDimension()+1 ),
130  M_pts( nDim, nLocalDof ),
131  M_points_face( nFacesInConvex ),
132  M_fset( primal )
133  {
134  DVLOG(2) << "Hermite finite element: \n";
135  DVLOG(2) << " o- dim = " << nDim << "\n";
136  DVLOG(2) << " o- order = " << nOrder << "\n";
137  DVLOG(2) << " o- nVertices = " << nVertices << "\n";
138  DVLOG(2) << " o- nGeometricFaces= " << nGeometricFaces << "\n";
139  DVLOG(2) << " o- numPoints = " << numPoints << "\n";
140  DVLOG(2) << " o- nbPtsPerVertex = " << nbPtsPerVertex << "\n";
141  DVLOG(2) << " o- nbPtsPerEdge = " << nbPtsPerEdge << "\n";
142  DVLOG(2) << " o- nbPtsPerFace = " << nbPtsPerFace << "\n";
143  DVLOG(2) << " o- nbPtsPerVolume = " << nbPtsPerVolume << "\n";
144 
145  DVLOG(2) << " o- nbDofPerVertex = " << nDofPerVertex << "\n";
146  DVLOG(2) << " o- nbDofPerEdge = " << nDofPerEdge << "\n";
147  DVLOG(2) << " o- nbDofPerFace = " << nDofPerFace << "\n";
148  DVLOG(2) << " o- nbDofPerVolume = " << nDofPerVolume << "\n";
149  DVLOG(2) << " o- nLocalDof = " << nLocalDof << "\n";
150 
151  pointset_type pts;
152 
153  //M_pts = pts.points();
154 
155  int i = 0;
156 
157  for ( uint16_type e = M_convex_ref.entityRange( 0 ).begin();
158  e < M_convex_ref.entityRange( 0 ).end();
159  ++e )
160  {
161  M_points_face[e] = pts.pointsBySubEntity( 0, e, 1 );
162  points_type _pts = pts.pointsBySubEntity( 0, e, 1 );
163 
164  for ( int j = 0; j < _pts.size2(); ++j, ++i )
165  {
166  ublas::column( M_pts, i ) = ublas::column( _pts, j );
167  DVLOG(2) << "pts " << i << " = " << ublas::column( M_pts, i ) << "\n";
168  }
169  }
170 
171  for ( int j = 0; j < nDim; ++j )
172  {
173  ublas::subrange( M_pts, 0, nDim, i, i+nDim+1 ) = ublas::subrange( M_pts, 0, nDim, 0, nDim+1 );
174  i+=nDim+1;
175  }
176 
177  if ( nDim == 2 )
178  {
179  points_type _pts = pts.pointsBySubEntity( 2, 0, 0 );
180  ublas::column( M_pts, i ) = ublas::column( _pts, 0 );
181  DVLOG(2) << "pts " << i << " = " << ublas::column( M_pts, i ) << "\n";
182  }
183 
184  //std::cout << "pts = " << M_pts << "\n";
185  setFset( primal, M_pts, mpl::bool_<primal_space_type::is_scalar>() );
186  }
187  ~HermiteDual()
188  {
189 
190  }
191  points_type const& points() const
192  {
193  return M_pts;
194  }
195 
196  points_type const& points( uint16_type f ) const
197  {
198  return M_points_face[f];
199  }
200  ublas::matrix_column<points_type const> point( uint16_type f, uint32_type __i ) const
201  {
202  return ublas::column( M_points_face[f], __i );
203  }
204  ublas::matrix_column<points_type> point( uint16_type f, uint32_type __i )
205  {
206  return ublas::column( M_points_face[f], __i );
207  }
208 
209  matrix_type operator()( primal_space_type const& pset ) const
210  {
211  matrix_type m = M_fset( pset );
212  //std::cout << "m=" << m << "\n";
213  return m;
214  }
215 private:
216 
217  void setFset( primal_space_type const& primal, points_type const& __pts, mpl::bool_<true> )
218  {
219  int nfs = 1+nDim+( ( nDim==2 )?1:0 );
220  std::vector<std::vector<Functional<primal_space_type> > > pd( nfs );
221  points_type pts1 = ublas::project( __pts,
222  ublas::range( 0,nDim ),
223  ublas::range( 0, nDim+1 ) );
224  pd[0] = functional::PointsEvaluation<primal_space_type>( primal, pts1 );
225 
226  //std::cout << "pd[" << 0 << "].size()=" << pd[0].size() << "\n";
227  for ( int d = 0; d < nDim; ++d )
228  {
229  pd[d+1] = functional::PointsDerivative<primal_space_type>( primal, d, pts1 );
230  //std::cout << "pd[" << d+1 << "].size()=" << pd[d+1].size() << "\n";
231  }
232 
233  if ( nDim == 2 )
234  {
235  points_type pts3 = ublas::project( __pts,
236  ublas::range( 0,nDim ),
237  ublas::range( ( nDim+1 )*( nDim+1 ), numPoints ) );
238  //std::cout << "pts3 = " << pts3 << "\n";
239  pd[3] = functional::PointsEvaluation<primal_space_type>( primal, pts3 );
240  //std::cout << "pd[2]=" << pd[3][0].coeff() << "\n";
241  }
242 
243  std::vector<Functional<primal_space_type> > fs( nLocalDof );
244  typename std::vector<Functional<primal_space_type> >::iterator it = fs.begin();
245 
246  for ( int i = 0; i < nfs; ++i )
247  {
248  //std::cout << "pd[" << i << "].size()=" << pd[i].size() << "\n";
249  it = std::copy( pd[i].begin(), pd[i].end(), it );
250  }
251 
252  M_fset.setFunctionalSet( fs );
253  }
254 
255  void setFset( primal_space_type const& primal, points_type const& __pts, mpl::bool_<false> )
256  {
257  //M_fset.setFunctionalSet( functional::ComponentsPointsEvaluation<primal_space_type>( primal,
258  //__pts ) );
259  }
260 
264  void setPoints( uint16_type f, points_type const& n )
265  {
266  M_points_face[f].resize( n.size1(), n.size2(), false );
267  M_points_face[f] = n;
268  }
269 
270 private:
271  reference_convex_type M_convex_ref;
272  std::vector<std::vector<uint16_type> > M_eid;
273  points_type M_pts;
274  std::vector<points_type> M_points_face;
275  FunctionalSet<primal_space_type> M_fset;
276 
277 
278 };
279 }// details
281 
297 template<uint16_type N,
298  uint16_type O,
299  template<uint16_type Dim> class PolySetType,
300  typename T = double,
301  template<uint16_type, uint16_type, uint16_type> class Convex = Simplex,
302  template<class, uint16_type, class> class Pts = PointSetEquiSpaced >
303 class Hermite
304  :
305 public FiniteElement<detail::OrthonormalPolynomialSet<N, N, O, PolySetType, T, Convex>, details::HermiteDual, Pts >
306 {
308 public:
309 
310  BOOST_STATIC_ASSERT( ( boost::is_same<PolySetType<N>, Scalar<N> >::value ||
311  boost::is_same<PolySetType<N>, Vectorial<N> >::value ||
312  boost::is_same<PolySetType<N>, Tensor2<N> >::value ) );
313 
317 
318  static const uint16_type nDim = N;
319  static const uint16_type nOrder = O;
320  static const bool isTransformationEquivalent = false;
321  static const bool isContinuous = true;
322 
323  typedef typename super::value_type value_type;
324  typedef typename super::primal_space_type primal_space_type;
325  typedef typename super::dual_space_type dual_space_type;
326  typedef Continuous continuity_type;
331  static const bool is_tensor2 = polyset_type::is_tensor2;
332  static const bool is_vectorial = polyset_type::is_vectorial;
333  static const bool is_scalar = polyset_type::is_scalar;
334  static const uint16_type nComponents = polyset_type::nComponents;
335 
336 
338 
339  typedef typename mpl::if_<mpl::equal_to<mpl::int_<nDim>, mpl::int_<1> >,
340  mpl::identity<boost::none_t>,
341  mpl::identity< Hermite<N-1, O, Scalar, T, Convex, Pts> > >::type::type face_basis_type;
342 
343  typedef boost::shared_ptr<face_basis_type> face_basis_ptrtype;
344 
345  typedef typename dual_space_type::convex_type convex_type;
346  typedef typename dual_space_type::pointset_type pointset_type;
347  typedef typename dual_space_type::reference_convex_type reference_convex_type;
348  typedef typename reference_convex_type::node_type node_type;
349  typedef typename reference_convex_type::points_type points_type;
350 
351  static const uint16_type numPoints = reference_convex_type::numPoints;
352  static const uint16_type nbPtsPerVertex = reference_convex_type::nbPtsPerVertex;
353  static const uint16_type nbPtsPerEdge = reference_convex_type::nbPtsPerEdge;
354  static const uint16_type nbPtsPerFace = reference_convex_type::nbPtsPerFace;
355  static const uint16_type nbPtsPerVolume = reference_convex_type::nbPtsPerVolume;
356 
357  static const uint16_type nLocalDof = super::nLocalDof;
359 
363 
364  Hermite()
365  :
366  super( dual_space_type( primal_space_type() ) ),
367  M_refconvex()
368  // M_bdylag( new face_basis_type )
369  {
370 
371 
372 
373  }
374 
375  virtual ~Hermite() {}
376 
378 
382 
383 
385 
389 
393  reference_convex_type const& referenceConvex() const
394  {
395  return M_refconvex;
396  }
397 
401  std::string familyName() const
402  {
403  return "hermite";
404  }
405 
407 
411 
412 
414 
418  template<typename GMContext, typename PC, typename Phi, typename GPhi, typename HPhi >
419  static void transform( boost::shared_ptr<GMContext> gmc, boost::shared_ptr<PC> const& pc,
420  Phi& phi_t,
421  GPhi& g_phi_t, const bool do_gradient,
422  HPhi& h_phi_t, const bool do_hessian
423 
424  )
425  {
426  transform ( *gmc, *pc, phi_t, g_phi_t, do_gradient, h_phi_t, do_hessian );
427  }
428  template<typename GMContext, typename PC, typename Phi, typename GPhi, typename HPhi >
429  static void transform( GMContext const& gmc,
430  PC const& pc,
431  Phi& phi_t,
432  GPhi& g_phi_t, const bool do_gradient,
433  HPhi& h_phi_t, const bool do_hessian
434 
435  )
436  {
437  //phi_t = phi; return ;
438  typename GMContext::gm_type::matrix_type const& B = gmc.B( 0 );
439  std::fill( phi_t.data(), phi_t.data()+phi_t.num_elements(), value_type( 0 ) );
440 
441  if ( do_gradient )
442  std::fill( g_phi_t.data(), g_phi_t.data()+g_phi_t.num_elements(), value_type( 0 ) );
443 
444  if ( do_hessian )
445  std::fill( h_phi_t.data(), h_phi_t.data()+g_phi_t.num_elements(), value_type( 0 ) );
446 
447  const uint16_type Q = gmc.nPoints();//M_grad.size2();
448 
449  // transform
450  for ( uint16_type i = numPoints; i < nLocalDof; i+=nDim )
451  {
452  for ( uint16_type l = 0; l < nDim; ++l )
453  {
454  for ( uint16_type p = 0; p < nDim; ++p )
455  {
456  for ( uint16_type q = 0; q < Q; ++q )
457  {
458  // \warning : here the transformation depends on the
459  // numbering of the degrees of freedom of the finite
460  // element
461  phi_t[i+l][0][0][q] += B( l, p ) * pc.phi( i+p,0,0,q );
462 
463  if ( do_gradient )
464  {
465  for ( uint16_type j = 0; j < nDim; ++j )
466  {
467  g_phi_t[i+l][0][p][q] += B( l, j ) * pc.grad( i+j,0,p,q );
468  }
469  }
470 
471  if ( do_hessian )
472  {
473  }
474  }
475  }
476  }
477  }
478  }
480 
481 private:
482 
483  reference_convex_type M_refconvex;
484  face_basis_ptrtype M_bdylag;
485 };
486 template<uint16_type N,
487  uint16_type O,
488  template<uint16_type Dim> class PolySetType,
489  typename T,
490  template<uint16_type, uint16_type, uint16_type> class Convex,
491  template<class, uint16_type, class> class Pts >
492 const uint16_type Hermite<N,O,PolySetType,T,Convex,Pts>::nDim;
493 
494 template<uint16_type N,
495  uint16_type O,
496  template<uint16_type Dim> class PolySetType,
497  typename T,
498  template<uint16_type, uint16_type, uint16_type> class Convex,
499  template<class, uint16_type, class> class Pts >
500 const uint16_type Hermite<N,O,PolySetType,T,Convex,Pts>::nOrder;
501 
502 template<uint16_type N,
503  uint16_type O,
504  template<uint16_type Dim> class PolySetType,
505  typename T,
506  template<uint16_type, uint16_type, uint16_type> class Convex,
507  template<class, uint16_type, class> class Pts >
508 const uint16_type Hermite<N,O,PolySetType,T,Convex,Pts>::numPoints;
509 
510 } // namespace fem
511 
512 template<uint16_type Order,
513  template<uint16_type Dim> class PolySetType = Scalar,
514  template<class, uint16_type, class> class Pts = PointSetEquiSpaced>
515 class Hermite
516 {
517 public:
518  template<uint16_type N,
519  typename T = double,
520  typename Convex = Simplex<N> >
521  struct apply
522  {
523  typedef typename mpl::if_<mpl::bool_<Convex::is_simplex>,
524  mpl::identity<fem::Hermite<N,Order,PolySetType,T,Simplex, Pts> >,
525  mpl::identity<fem::Hermite<N,Order,PolySetType,T,Hypercube, Pts> > >::type::type result_type;
526  typedef result_type type;
527  };
528 
529  typedef Hermite<Order,Scalar,Pts> component_basis_type;
530 
531  static const uint16_type nOrder = Order;
532 
533 };
534 
535 } // namespace Feel
536 #endif /* __hermite_H */

Generated on Fri Oct 25 2013 14:24:15 for Feel++ by doxygen 1.8.4