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
raviartthomas.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: 2006-01-14
7 
8  Copyright (C) 2006 EPFL
9  Copyright (C) 2008, 2009 Université de Grenoble 1 (Joseph Fourier)
10 
11  This library is free software; you can redistribute it and/or
12  modify it under the terms of the GNU Lesser General Public
13  License as published by the Free Software Foundation; either
14  version 3.0 of the License, or (at your option) any later version.
15 
16  This library is distributed in the hope that it will be useful,
17  but WITHOUT ANY WARRANTY; without even the implied warranty of
18  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19  Lesser General Public License for more details.
20 
21  You should have received a copy of the GNU Lesser General Public
22  License along with this library; if not, write to the Free Software
23  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
24 */
30 #ifndef __RaviartThomas_H
31 #define __RaviartThomas_H 1
32 
33 #include <boost/ptr_container/ptr_vector.hpp>
34 #include <boost/assign/std/vector.hpp> // for 'operator+=()'
35 #include <boost/numeric/ublas/vector.hpp>
36 #include <boost/numeric/ublas/io.hpp>
37 #include <boost/numeric/ublas/matrix.hpp>
38 #include <boost/numeric/ublas/matrix_proxy.hpp>
39 #include <boost/numeric/ublas/lu.hpp>
40 
41 #include <boost/assign/list_of.hpp>
42 #include <boost/assign/std/vector.hpp>
43 
44 #include <feel/feelcore/feel.hpp>
45 #include <feel/feelcore/traits.hpp>
46 #include <feel/feelalg/lu.hpp>
47 
50 
51 
58 #include <feel/feelpoly/quadpoint.hpp>
59 #include <feel/feelpoly/fe.hpp>
60 
61 #include <feel/feelvf/vf.hpp>
62 
63 namespace Feel
64 {
65 namespace detail
66 {
67 template<typename P>
68 struct times_x
69 {
70  typedef typename P::value_type value_type;
71  typedef typename P::points_type points_type;
72  times_x ( P const& p, int c )
73  :
74  M_p ( p ),
75  M_c( c )
76  {
77  //std::cout << "component : " << c << std::endl;
78  }
79  typename ublas::vector<value_type> operator() ( points_type const& __pts ) const
80  {
81 #if 0
82  std::cout << "times_x(pts) : " << __pts << std::endl;
83  std::cout << "times_x(pts) : " << M_p.evaluate( __pts ) << std::endl;
84  std::cout << "times_x(coeff) : " << M_p.coefficients() << std::endl;
85 #endif
86  // __pts[c] * p( __pts )
87  return ublas::element_prod( ublas::row( __pts, M_c ),
88  ublas::row( M_p.evaluate( __pts ), 0 ) );
89  }
90  //P const& M_p;
91  P M_p;
92  int M_c;
93 };
94 
95 template< class T >
96 struct extract_all_poly_indices
97 {
98  T start;
99  extract_all_poly_indices( T comp, T N )
100  :
101  start( comp*N )
102  { }
103 
104  T operator()()
105  {
106  return start++;
107  }
108 };
109 }// detail
110 
111 template<uint16_type N,
112  uint16_type O,
113  typename T = double,
114  template<uint16_type, uint16_type, uint16_type> class Convex = Simplex>
115 class RaviartThomasPolynomialSet
116  :
117  public Feel::detail::OrthonormalPolynomialSet<N, N, O+1, Vectorial, T, Convex>
118 {
119  typedef Feel::detail::OrthonormalPolynomialSet<N, N, O+1, Vectorial, T, Convex> super;
120 
121 public:
122  static const uint16_type Om1 = (O==0)?0:O-1;
123  typedef Feel::detail::OrthonormalPolynomialSet<N, N, O, Vectorial, T, Convex> Pk_v_type;
124  typedef Feel::detail::OrthonormalPolynomialSet<N, N, O+1, Vectorial, T, Convex> Pkp1_v_type;
125  typedef Feel::detail::OrthonormalPolynomialSet<N, N, Om1, Vectorial, T, Convex> Pkm1_v_type;
126  typedef Feel::detail::OrthonormalPolynomialSet<N, N, O, Scalar, T, Convex> Pk_s_type;
127  typedef Feel::detail::OrthonormalPolynomialSet<N, N, O+1, Scalar, T, Convex> Pkp1_s_type;
128 
129  typedef PolynomialSet<typename super::basis_type,Vectorial> vectorial_polynomialset_type;
130  typedef typename vectorial_polynomialset_type::polynomial_type vectorial_polynomial_type;
131  typedef PolynomialSet<typename super::basis_type,Scalar> scalar_polynomialset_type;
132  typedef typename scalar_polynomialset_type::polynomial_type scalar_polynomial_type;
133 
134  typedef RaviartThomasPolynomialSet<N, O, T> self_type;
135 
136  typedef typename super::value_type value_type;
137  typedef typename super::convex_type convex_type;
138  typedef typename super::matrix_type matrix_type;
139  typedef typename super::points_type points_type;
140 
141  static const uint16_type nDim = super::nDim;
142  static const uint16_type nOrder = super::nOrder;
143  static const uint16_type nComponents = super::nComponents;
144  static const bool is_product = false;
145  RaviartThomasPolynomialSet()
146  :
147  super()
148  {
149  uint16_type dim_Pkp1 = convex_type::polyDims( nOrder );
150  uint16_type dim_Pk = convex_type::polyDims( nOrder-1 );
151  uint16_type dim_Pkm1 = ( nOrder==1 )?0:convex_type::polyDims( nOrder-2 );
152 #if 0
153  std::cout << "[RTPset] dim_Pkp1 = " << dim_Pkp1 << "\n";
154  std::cout << "[RTPset] dim_Pk = " << dim_Pk << "\n";
155  std::cout << "[RTPset] dim_Pkm1 = " << dim_Pkm1 << "\n";
156 #endif
157  // (P_k)^d
158  Pkp1_v_type Pkp1_v;
159  vectorial_polynomialset_type Pk_v( Pkp1_v.polynomialsUpToDimension( dim_Pk ) );
160 #if 0
161  std::cout << "[RTPset] Pk_v =" << Pk_v.coeff() << "\n";
162 #endif
163  // P_k
164  Pkp1_s_type Pkp1;
165  scalar_polynomialset_type Pk ( Pkp1.polynomialsUpToDimension( dim_Pk ) );
166 #if 0
167  std::cout << "[RTPset] Pk =" << Pk.coeff() << "\n";
168  std::cout << "[RTPset] Pk(0) =" << Pk.polynomial( 0 ).coefficients() << "\n";
169 #endif
170 
171  // x P_k \ P_{k-1}
172  IMGeneral<convex_type::nDim, 2*nOrder,value_type> im;
173  //std::cout << "[RTPset] im.points() = " << im.points() << std::endl;
174  ublas::matrix<value_type> xPkc( nComponents*( dim_Pk-dim_Pkm1 ),Pk.coeff().size2() );
175 
176  //std::cout << "[RTPset] before xPkc = " << xPkc << "\n";
177  for ( int l = dim_Pkm1, i = 0; l < dim_Pk; ++l, ++i )
178  {
179  for ( int j = 0; j < convex_type::nDim; ++j )
180  {
181  Feel::detail::times_x<scalar_polynomial_type> xp( Pk.polynomial( l ), j );
182  ublas::row( xPkc,i*nComponents+j )=
183  ublas::row( Feel::project( Pkp1,
184  xp,
185  im ).coeff(), 0 );
186  }
187  }
188 
189 
190  //std::cout << "[RTPset] after xPkc = " << xPkc << "\n";
191  vectorial_polynomialset_type xPk( typename super::basis_type(), xPkc, true );
192  //std::cout << "[RTPset] here 1\n";
193  // (P_k)^d + x P_k
194  //std::cout << "[RTPset] RT Poly coeff = " << unite( Pk_v, xPk ).coeff() << "\n";
195  this->setCoefficient( unite( Pk_v, xPk ).coeff(), true );
196  //std::cout << "[RTPset] here 2\n";
197  }
198 
199 
200 };
201 
202 namespace fem
203 {
204 
205 namespace detail
206 {
207 
208 
209 
210 template<typename Basis,
211  template<class, uint16_type, class> class PointSetType>
212 class RaviartThomasDual
213  :
214 public DualBasis<Basis>
215 {
216  typedef DualBasis<Basis> super;
217 public:
218 
219  static const uint16_type nDim = super::nDim;
220  static const uint16_type nOrder= super::nOrder;
221 
222  typedef typename super::primal_space_type primal_space_type;
223  typedef typename primal_space_type::value_type value_type;
224  typedef typename primal_space_type::points_type points_type;
225  typedef typename primal_space_type::matrix_type matrix_type;
226  typedef typename primal_space_type::template convex<nDim+nOrder-1>::type convex_type;
227  typedef Reference<convex_type, nDim, nDim+nOrder-1, nDim, value_type> reference_convex_type;
228  typedef typename reference_convex_type::node_type node_type;
229 
230  typedef typename primal_space_type::Pkp1_v_type Pkp1_v_type;
231  typedef typename primal_space_type::vectorial_polynomialset_type vectorial_polynomialset_type;
232 
233  // point set type associated with the functionals
234  typedef PointSetType<convex_type, nOrder, value_type> pointset_type;
235 
236  static const uint16_type nbPtsPerVertex = 0;
237  static const uint16_type nbPtsPerEdge = mpl::if_<mpl::equal_to<mpl::int_<nDim>,mpl::int_<2> >,mpl::int_<reference_convex_type::nbPtsPerEdge>,mpl::int_<0> >::type::value;
238  static const uint16_type nbPtsPerFace =mpl::if_<mpl::equal_to<mpl::int_<nDim>,mpl::int_<3> >,mpl::int_<reference_convex_type::nbPtsPerFace>,mpl::int_<0> >::type::value;
239  static const uint16_type nbPtsPerVolume = 0;
240  static const uint16_type numPoints = ( reference_convex_type::numGeometricFaces*nbPtsPerFace+reference_convex_type::numEdges*nbPtsPerEdge );
241 
243  static const uint16_type nDofPerVertex = 0;
244 
246  static const uint16_type nDofPerEdge = nbPtsPerEdge;
247 
249  static const uint16_type nDofPerFace = nbPtsPerFace;
250 
252  static const uint16_type nDofPerVolume = 0;
253 
255  static const uint16_type nLocalDof = numPoints;
256 
257  RaviartThomasDual( primal_space_type const& primal )
258  :
259  super( primal ),
260  M_convex_ref(),
261  M_eid( M_convex_ref.topologicalDimension()+1 ),
262  M_pts( nDim, numPoints ),
263  M_pts_per_face( convex_type::numTopologicalFaces ),
264  M_fset( primal )
265  {
266 #if 1
267  std::cout << "Raviart-Thomas finite element(dual): \n";
268  std::cout << " o- dim = " << nDim << "\n";
269  std::cout << " o- order = " << nOrder << "\n";
270  std::cout << " o- numPoints = " << numPoints << "\n";
271  std::cout << " o- nbPtsPerVertex = " << ( int )nbPtsPerVertex << "\n";
272  std::cout << " o- nbPtsPerEdge = " << ( int )nbPtsPerEdge << "\n";
273  std::cout << " o- nbPtsPerFace = " << ( int )nbPtsPerFace << "\n";
274  std::cout << " o- nbPtsPerVolume = " << ( int )nbPtsPerVolume << "\n";
275  std::cout << " o- nLocalDof = " << nLocalDof << "\n";
276 #endif
277 
278  // loop on each entity forming the convex of topological
279  // dimension nDim-1 ( the faces)
280  for ( int p = 0, e = M_convex_ref.entityRange( nDim-1 ).begin();
281  e < M_convex_ref.entityRange( nDim-1 ).end();
282  ++e )
283  {
284  points_type Gt ( M_convex_ref.makePoints( nDim-1, e ) );
285  M_pts_per_face[e] = Gt ;
286 
287  if ( Gt.size2() )
288  {
289  //VLOG(1) << "Gt = " << Gt << "\n";
290  //VLOG(1) << "p = " << p << "\n";
291  ublas::subrange( M_pts, 0, nDim, p, p+Gt.size2() ) = Gt;
292  //for ( size_type j = 0; j < Gt.size2(); ++j )
293  //M_eid[d].push_back( p+j );
294  p+=Gt.size2();
295  }
296  }
297 
298  //std::cout << "[RT Dual] done 1\n";
299  // compute \f$ \ell_e( U ) = (U * n[e]) (edge_pts(e)) \f$
300  typedef Functional<primal_space_type> functional_type;
301  std::vector<functional_type> fset;
302 
303  // jacobian of the transformation from reference face to the face in the
304  // reference element
305  std::vector<double> j;
306  {
307  // bring 'operator+=()' into scope
308  using namespace boost::assign;
309 
310  if ( nDim == 2 )
311  j += 2.8284271247461903,2.0,2.0;
312 
313  if ( nDim == 3 )
314  j+= 3.464101615137754, 2, 2, 2;
315  }
316 
317  //for( int k = 0; k < nDim; ++k )
318  {
319  // loopover the each edge entities and add the correponding functionals
320  for ( int e = M_convex_ref.entityRange( nDim-1 ).begin();
321  e < M_convex_ref.entityRange( nDim-1 ).end();
322  ++e )
323  {
324  typedef Feel::functional::DirectionalComponentPointsEvaluation<primal_space_type> dcpe_type;
325  node_type dir = M_convex_ref.normal( e )*j[e];
326  //dcpe_type __dcpe( primal, 1, dir, pts_per_face[e] );
327  dcpe_type __dcpe( primal, dir, M_pts_per_face[e] );
328  std::copy( __dcpe.begin(), __dcpe.end(), std::back_inserter( fset ) );
329  }
330  }
331 
332  //std::cout << "[RT Dual] done 2" << std::endl;
333  if ( nOrder-1 > 0 )
334  {
335  // we need more equations : add interior moment
336  // indeed the space is orthogonal to Pk-1
337  uint16_type dim_Pkp1 = convex_type::polyDims( nOrder );
338  uint16_type dim_Pk = convex_type::polyDims( nOrder-1 );
339  uint16_type dim_Pm1 = convex_type::polyDims( nOrder-2 );
340 
341  Pkp1_v_type Pkp1;
342 
343  vectorial_polynomialset_type Pkm1 ( Pkp1.polynomialsUpToDimension( dim_Pm1 ) );
344 
345  //std::cout << "Pkm1 = " << Pkm1.coeff() << "\n";
346  //std::cout << "Primal = " << primal.coeff() << "\n";
347  for ( int i = 0; i < Pkm1.polynomialDimension(); ++i )
348  {
349  typedef functional::IntegralMoment<primal_space_type, vectorial_polynomialset_type> fim_type;
350  //typedef functional::IntegralMoment<Pkp1_v_type, vectorial_polynomialset_type> fim_type;
351  //std::cout << "P(" << i << ")=" << Pkm1.polynomial( i ).coeff() << "\n";
352  fset.push_back( fim_type( primal, Pkm1.polynomial( i ) ) );
353  }
354  }
355 
356  //std::cout << "[RT Dual] done 3, n fset = " << fset.size() << std::endl;
357  M_fset.setFunctionalSet( fset );
358  // std::cout << "[RT DUAL matrix] mat = " << M_fset.rep() << "\n";
359  //std::cout << "[RT Dual] done 4\n";
360 
361  }
362 
366  //pointset_type const& pointSet() const { return M_pset; }
367 
368  points_type const& points() const
369  {
370  return M_pts;
371  }
372 
373 
374  matrix_type operator()( primal_space_type const& pset ) const
375  {
376  //std::cout << "RT matrix = " << M_fset( pset ) << std::endl;
377  return M_fset( pset );
378  }
379 
380  points_type const& points( uint16_type f ) const
381  {
382  return M_pts_per_face[f];
383  }
384  ublas::matrix_column<points_type const> point( uint16_type f, uint32_type __i ) const
385  {
386  return ublas::column( M_pts_per_face[f], __i );
387  }
388  ublas::matrix_column<points_type> point( uint16_type f, uint32_type __i )
389  {
390  return ublas::column( M_pts_per_face[f], __i );
391  }
392 
393 private:
397  void setPoints( uint16_type f, points_type const& n )
398  {
399  M_pts_per_face[f].resize( n.size1(), n.size2(), false );
400  M_pts_per_face[f] = n;
401  }
402 
403 private:
404  reference_convex_type M_convex_ref;
405  std::vector<std::vector<uint16_type> > M_eid;
406  points_type M_pts;
407  std::vector<points_type> M_pts_per_face;
408  FunctionalSet<primal_space_type> M_fset;
409 
410 
411 };
412 }// detail
413 
414 
423 template<uint16_type N,
424  uint16_type O,
425  typename T = double,
426  template<uint16_type, uint16_type, uint16_type> class Convex = Simplex,
427  uint16_type TheTAG=0 >
429  :
430 public FiniteElement<RaviartThomasPolynomialSet<N, O, T, Convex>,
431  fem::detail::RaviartThomasDual,
432  PointSetEquiSpaced >,
433 public boost::enable_shared_from_this<RaviartThomas<N,O,T,Convex> >
434 {
436  fem::detail::RaviartThomasDual,
437  PointSetEquiSpaced > super;
438 public:
439 
440  BOOST_STATIC_ASSERT( N > 1 );
441 
445 
446  static const uint16_type nDim = N;
447  //static const bool isTransformationEquivalent = false;
448  static const bool isTransformationEquivalent = true;
449  static const bool isContinuous = true;
450  typedef Continuous continuity_type;
451  static const uint16_type TAG = TheTAG;
452 
453  //static const polynomial_transformation_type transformation = POLYNOMIAL_CONTEXT_NEEDS_1ST_PIOLA_TRANSFORMATION;
454  typedef typename super::value_type value_type;
455  typedef typename super::primal_space_type primal_space_type;
456  typedef typename super::dual_space_type dual_space_type;
457 
462  static const bool is_vectorial = polyset_type::is_vectorial;
463  static const bool is_scalar = polyset_type::is_scalar;
464  static const uint16_type nComponents = polyset_type::nComponents;
465  static const bool is_product = false;
466 
467 
468  typedef typename dual_space_type::convex_type convex_type;
469  typedef typename dual_space_type::pointset_type pointset_type;
470  typedef typename dual_space_type::reference_convex_type reference_convex_type;
471  typedef typename reference_convex_type::node_type node_type;
472  typedef typename reference_convex_type::points_type points_type;
473 
474 
475  static const uint16_type nOrder = dual_space_type::nOrder;
476  static const uint16_type nbPtsPerVertex = 0;
477  static const uint16_type nbPtsPerEdge = mpl::if_<mpl::equal_to<mpl::int_<nDim>,mpl::int_<2> >,
478  mpl::int_<reference_convex_type::nbPtsPerEdge>,
479  mpl::int_<0> >::type::value;
480  static const uint16_type nbPtsPerFace = mpl::if_<mpl::equal_to<mpl::int_<nDim>,mpl::int_<3> >,
481  mpl::int_<reference_convex_type::nbPtsPerFace>,
482  mpl::int_<0> >::type::value;
483  static const uint16_type nbPtsPerVolume = 0;
484  static const uint16_type numPoints = ( reference_convex_type::numGeometricFaces*nbPtsPerFace+
485  reference_convex_type::numEdges*nbPtsPerEdge );
486 
487  static const uint16_type nLocalDof = dual_space_type::nLocalDof;
489 
493 
494  RaviartThomas()
495  :
496  super( dual_space_type( primal_space_type() ) ),
497  M_refconvex()
498  {
499 #if 0
500  std::cout << "[RT] nPtsPerEdge = " << nbPtsPerEdge << "\n";
501  std::cout << "[RT] nPtsPerFace = " << nbPtsPerFace << "\n";
502  std::cout << "[RT] numPoints = " << numPoints << "\n";
503 
504  std::cout << "[RT] nDof = " << super::nDof << "\n";
505 
506  std::cout << "[RT] coeff : " << this->coeff() << "\n";
507  std::cout << "[RT] pts : " << this->points() << "\n";
508  std::cout << "[RT] eval at pts : " << this->evaluate( this->points() ) << "\n";
509  std::cout << "[RT] is_product : " << is_product << "\n";
510 #endif
511  }
512  RaviartThomas( RaviartThomas const & cr )
513  :
514  super( cr ),
515  M_refconvex()
516  {}
517  ~RaviartThomas()
518  {}
519 
521 
525 
526 
528 
532 
536  reference_convex_type const& referenceConvex() const
537  {
538  return M_refconvex;
539  }
540 
541  std::string familyName() const
542  {
543  return "raviartthomas";
544  }
545 
547 
551 
552 
554 
558 
559  template<typename ExprType>
560  static auto
561  isomorphism( ExprType expr ) -> decltype( Feel::vf::detJ()*( trans( Feel::vf::JinvT() )*expr )*Feel::vf::Nref() )
562  //isomorphism( ExprType& expr ) -> decltype( expr )
563  {
564  using namespace Feel::vf;
565  return detJ()*( trans( JinvT() )*expr )*Nref();
566  //return expr;
567  }
568 #if 0
569 
573  template<typename ExprType, typename ContextType>
574  std::vector<value_type>
575  interpolate( boost::shared_ptr<ContextType>& ctx, ExprType & expr )
576  {
577  using namespace Feel::vf;
578  typedef boost::shared_ptr<ContextType> gmc_ptrtype;
579  typedef fusion::map<fusion::pair<vf::detail::gmc<0>, gmc_ptrtype> > map_gmc_type;
580 
581  std::vector<value_type> v( nLocalDof );
582 
583  // First deal with the face dof
584  for ( int face = 0; face < numTopologicalFaces; ++face )
585  {
586  // update the geomap at dof on face
587  ctx->update( _face=face, _element=ctx->id() );
588 
589  map_gmc_type mapgmc( fusion::make_pair<vf::detail::gmc<0> >( ctx ) );
590  expr.update( mapgmc, face );
591 
592  for ( int q = 0; q < nDofPerFace; ++q )
593  {
594  int ldof = nDofPerFace*face+i;
595  v[ldof] = expr.evalq( 0,0,i );
596  }
597  }
598 
599  // evaluate expr \cdot n on each face
600 
601  // evaluate moments of the expression
602  }
603 #endif
604 
605 #if 0
606  template<typename GMContext, typename PC, typename Phi, typename GPhi, typename HPhi >
607  static void transform( boost::shared_ptr<GMContext> gmc, boost::shared_ptr<PC> const& pc,
608  Phi& phi_t,
609  GPhi& g_phi_t, const bool do_gradient,
610  HPhi& h_phi_t, const bool do_hessian
611 
612  )
613  {
614  transform ( *gmc, *pc, phi_t, g_phi_t, do_gradient, h_phi_t, do_hessian );
615  }
616  template<typename GMContext, typename PC, typename Phi, typename GPhi, typename HPhi >
617  static void transform( GMContext const& gmc,
618  PC const& pc,
619  Phi& phi_t,
620  GPhi& g_phi_t, const bool do_gradient,
621  HPhi& h_phi_t, const bool do_hessian
622 
623  )
624  {
625  //phi_t = phi; return ;
626  typename GMContext::gm_type::matrix_type const B = gmc.B( 0 );
627  typename GMContext::gm_type::matrix_type const K = gmc.K( 0 );
628  typename GMContext::gm_type::matrix_type JB( K/gmc.J( 0 ) );
629 #if 0
630  std::cout << "K= " << gmc.K( 0 ) << "\n";
631  std::cout << "B= " << B << "\n";
632  std::cout << "J= " << gmc.J( 0 ) << "\n";
633  std::cout << "JB= " << JB << "\n";
634 #endif
635  std::fill( phi_t.data(), phi_t.data()+phi_t.num_elements(), value_type( 0 ) );
636 
637  if ( do_gradient )
638  {
639  //std::cout << "compute gradient\n";
640  std::fill( g_phi_t.data(), g_phi_t.data()+g_phi_t.num_elements(), value_type( 0 ) );
641  }
642 
643  if ( do_hessian )
644  std::fill( h_phi_t.data(), h_phi_t.data()+h_phi_t.num_elements(), value_type( 0 ) );
645 
646  const uint16_type Q = gmc.nPoints();//M_grad.size2();
647 
648  // transform
649  for ( uint16_type i = 0; i < nLocalDof; ++i )
650  {
651  for ( uint16_type l = 0; l < nDim; ++l )
652  {
653 
654  for ( uint16_type p = 0; p < nDim; ++p )
655  {
656  for ( uint16_type q = 0; q < Q; ++q )
657  {
658  // \warning : here the transformation depends on the
659  // numbering of the degrees of freedom of the finite
660  // element
661  //phi_t[i][l][0][q] = pc.phi(i,l,0,q);
662  phi_t[i][l][0][q] += JB( l, p ) * pc.phi( i,p,0,q );
663  //phi_t[i][l][0][q] = gmc.J( 0 ) * B( p, l ) * pc.phi(i,p,0,q);
664  //std::cout << "pc[" << i << "][" << l << "][" << q << "]=" << pc.phi(i,l,0,q) << "\n";
665  //std::cout << "phi_t[" << i << "][" << l << "][" << q << "]=" << phi_t[i][l][0][q] << "\n";
666  }
667  }
668  }
669 
670  //if ( do_gradient )
671  {
672 
673  for ( uint16_type p = 0; p < nDim; ++p )
674  {
675  for ( uint16_type r = 0; r < nDim; ++r )
676  {
677  for ( uint16_type s = 0; s < Q; ++s )
678  {
679  g_phi_t[i][p][r][s] = 0;
680 
681  for ( uint16_type q = 0; q < nDim; ++q )
682  {
683  g_phi_t[i][p][r][s] += JB( p, q ) * pc.grad( i,q,r,s );
684  //g_phi_t[i][p][r][s] = pc.grad(i,p,r,s);
685 
686  //std::cout << "J G[" << i << "][" << q << "][" << r << "][" << s << "=" << JB( p, q ) * pc.grad(i,q,r,s) << "\n";
687  }
688  }
689 
690  //std::cout << "g_phi_t[" << i << "][" << p << "][" << r << "][" << 0 << "=" << g_phi_t[i][p][r][0] << "\n";
691 
692  }
693  }
694  }
695 
696  if ( do_hessian )
697  {
698  }
699  }
700  }
701 #endif
702 
704 
705 
706 
707 protected:
708  reference_convex_type M_refconvex;
709 private:
710 
711 };
712 template<uint16_type N,
713  uint16_type O,
714  typename T,
715  template<uint16_type, uint16_type, uint16_type> class Convex,
716  uint16_type TheTAG>
717 const uint16_type RaviartThomas<N,O,T,Convex,TheTAG>::nDim;
718 template<uint16_type N,
719  uint16_type O,
720  typename T,
721  template<uint16_type, uint16_type, uint16_type> class Convex,
722  uint16_type TheTAG >
723 const uint16_type RaviartThomas<N,O,T,Convex,TheTAG>::nOrder;
724 template<uint16_type N,
725  uint16_type O,
726  typename T,
727  template<uint16_type, uint16_type, uint16_type> class Convex,
728  uint16_type TheTAG>
729 const uint16_type RaviartThomas<N,O,T,Convex,TheTAG>::nLocalDof;
730 
731 } // fem
732 template<uint16_type Order,
733  uint16_type TheTAG=0 >
734 class RaviartThomas
735 {
736 public:
737  template<uint16_type N,
738  uint16_type R = N,
739  typename T = double,
740  typename Convex = Simplex<N> >
741  struct apply
742  {
743  typedef typename mpl::if_<mpl::bool_<Convex::is_simplex>,
744  mpl::identity<fem::RaviartThomas<N,Order,T,Simplex,TheTAG> >,
745  mpl::identity<fem::RaviartThomas<N,Order,T,Hypercube,TheTAG> > >::type::type result_type;
746  typedef result_type type;
747  };
748 
749  template<uint16_type TheNewTAG>
750  struct ChangeTag
751  {
752  typedef RaviartThomas<Order,TheNewTAG> type;
753  };
754 
755  typedef Lagrange<Order,Scalar> component_basis_type;
756 
757  static const uint16_type TAG = TheTAG;
758 
759 };
760 
761 } // Feel
762 #endif /* __RaviartThomas_H */

Generated on Sun Oct 20 2013 08:25:04 for Feel++ by doxygen 1.8.4