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

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