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
boundadapted.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): Gilles Steiner <gilles.steiner@epfl.ch>
6  Date: 2005-12-13
7 
8  Copyright (C) 2005,2006 EPFL
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 __Boundadapted_H
30 #define __Boundadapted_H 1
31 
32 
33 #include <boost/lambda/if.hpp>
34 #include <boost/numeric/ublas/banded.hpp>
36 #include <feel/feelalg/glas.hpp>
37 
40 #include <feel/feelpoly/basis.hpp>
41 
42 namespace Feel
43 {
44 template<class Convex,uint16_type O,typename T> class PointSetWarpBlend;
45 template<uint16_type Dim,uint16_type RealDim,uint16_type Degree,typename NormalizationPolicy,typename T,template<class> class StoragePolicy> class Dubiner;
46 
47 template<uint16_type Dim,
48  uint16_type Degree,
49  typename T,
50  template<class> class StoragePolicy>
52 
53 template<uint16_type Dim,
54  uint16_type Degree,
55  typename T = double,
56  template<class> class StoragePolicy = StorageUBlas>
57 struct BoundaryAdaptedTraits
58 {
59  static const uint16_type nDim = Dim;
60  static const uint16_type nOrder = Degree;
61  static const uint16_type nConvexOrderDiff = nDim+nOrder+1;
62  static const bool is_normalized = false;
63 
64 
65  typedef T value_type;
66 
67  template<uint16_type order, typename V = value_type>
68  struct Convex
69  {
70  typedef Simplex<nDim, order, nDim> type;
71  typedef Reference<Simplex<nDim, order, nDim>, nDim, order, nDim, V> reference_type;
72  };
73 
74  template<typename NewT>
75  struct ChangeValueType
76  {
78  typedef BoundaryAdaptedTraits<Dim, Degree, NewT, StoragePolicy> traits_type;
79 
80  };
81 
82  template<uint16_type NewOrder>
83  struct ChangeOrder
84  {
85  typedef BoundaryAdapted<Dim, NewOrder, T, StoragePolicy> type;
86  typedef BoundaryAdaptedTraits<Dim, NewOrder, T, StoragePolicy> traits_type;
87  };
88 
89  /*
90  * Geometry where the polynomials are defined and constructed
91  */
92 
93  typedef typename Convex<nOrder>::type convex_type;
94  typedef typename Convex<nOrder>::reference_type reference_convex_type;
95  typedef typename Convex<nConvexOrderDiff>::type diff_convex_type;
96  typedef typename Convex<nConvexOrderDiff>::reference_type diff_reference_convex_type;
97  typedef typename mpl::if_<mpl::equal_to<mpl::int_<nDim>, mpl::int_<2> >,
98  mpl::identity<PointSetWarpBlend<diff_convex_type,nConvexOrderDiff,value_type> >,
99  mpl::identity<PointSetEquiSpaced<diff_convex_type, nConvexOrderDiff,value_type> > >::type::type diff_pointset_type;
100 
101  static const uint16_type numVertices = reference_convex_type::numVertices;
102  static const uint16_type numFaces = reference_convex_type::numFaces;
103 
104  /*
105  * storage policy
106  */
107 
108  typedef StoragePolicy<value_type> storage_policy;
109  typedef typename storage_policy::vector_type vector_type;
110  typedef typename storage_policy::matrix_type matrix_type;
111  typedef typename storage_policy::vector_matrix_type vector_matrix_type;
112  typedef typename storage_policy::vector_vector_matrix_type vector_vector_matrix_type;
113  typedef typename storage_policy::matrix_node_type matrix_node_type;
114  typedef typename storage_policy::points_type points_type;
115  typedef typename storage_policy::node_type node_type;
116 };
117 
118 template<int D, int O>
119 struct BoundaryAdaptedTag
120 {
121  static const int Dim = D;
122  static const int Order = O;
123 };
124 
143 template<uint16_type Dim,
144  uint16_type Degree,
145  typename T = double,
146  template<class> class StoragePolicy = StorageUBlas>
147 class BoundaryAdapted : Basis<BoundaryAdaptedTag<Dim, Degree>, T>
148 {
149  typedef Basis<BoundaryAdaptedTag<Dim, Degree>, T > super;
150 
151 public:
152  typedef BoundaryAdaptedTraits<Dim, Degree, T, StoragePolicy> traits_type;
153 
154  static const uint16_type nDim = traits_type::nDim;
155  static const uint16_type nOrder = traits_type::nOrder;
156  static const uint16_type nConvexOrderDiff = traits_type::nConvexOrderDiff;
157  static const bool is_normalized = traits_type::is_normalized;
158  static const uint16_type numVertices = traits_type::numVertices;
159  static const uint16_type numFaces = traits_type::numFaces;
160 
164 
165  typedef BoundaryAdapted<Dim, Degree, T, StoragePolicy> self_type;
166 
167  /*
168  * for now can be used only as a basis but we might be interested
169  * to have then expressed in other basis like in the moment basis
170  */
171  typedef self_type basis_type;
172 
173  /*
174  * numerical type
175  */
176  typedef T value_type;
177  typedef int16_type int_type;
178  typedef uint16_type uint_type;
179 
180  /*
181  * Geometry where the polynomials are defined and constructed
182  */
183  typedef typename traits_type::convex_type convex_type;
184  typedef typename traits_type::reference_convex_type reference_convex_type;
185  typedef typename traits_type::diff_pointset_type diff_pointset_type;
186 
187  /*
188  * storage policy
189  */
190  typedef typename traits_type::storage_policy storage_policy;
191  typedef typename traits_type::vector_type vector_type;
192  typedef typename traits_type::matrix_type matrix_type;
193  typedef typename traits_type::vector_matrix_type vector_matrix_type;
194  typedef typename traits_type::vector_vector_matrix_type vector_vector_matrix_type;
195  typedef typename traits_type::matrix_node_type matrix_node_type;
196  typedef typename traits_type::points_type points_type;
197  typedef typename traits_type::node_type node_type;
198 
199  typedef Principal<Degree, T,StoragePolicy> principal_type;
201 
205 
206  BoundaryAdapted()
207  :
208  super( *this ),
209  M_refconvex(),
210  M_pts( nDim, numVertices ),
211  M_pts_face( reference_convex_type::numVertices )
212  {
213  PointSetEquiSpaced<convex_type,nOrder,value_type> pts;
214 
215  // only the points associated with the vertices
216  M_pts = pts.pointsByEntity( 0 );
217  DVLOG(2) << "[boundaryadapted] pts= " << M_pts << "\n";
218 
219  // get the points for each faces
220  for ( uint16_type e = M_refconvex.entityRange( nDim-1 ).begin();
221  e < M_refconvex.entityRange( nDim-1 ).end();
222  ++e )
223  {
224  M_pts_face[e] = pts.pointsBySubEntity( nDim-1, e, 1 );
225  DVLOG(2) << "[boundaryadapted] face " << e << " pts " << M_pts_face[e] << "\n";
226  }
227  }
228 
229  BoundaryAdapted( BoundaryAdapted const & d )
230  :
231  super( *this ),
232  M_refconvex( d.M_refconvex ),
233  M_pts( d.M_pts ),
234  M_pts_face( d.M_pts_face )
235  {}
236 
237  ~BoundaryAdapted()
238  {}
239 
241 
245  points_type const& points() const
246  {
247  return M_pts;
248  }
249 
253  points_type const& points( int f ) const
254  {
255  return M_pts_face[f];
256  }
257 
261 
262  self_type const& operator=( self_type const& d )
263  {
264  if ( this != &d )
265  {
266  M_refconvex = d.M_refconvex;
267  M_pts = d.M_pts;
268  }
269 
270  return *this;
271  }
272 
273  matrix_type operator()( node_type const& pt ) const
274  {
275  points_type pts( pt.size(), 1 );
276  ublas::column( pts, 0 ) = pt;
277  return evaluate( pts );
278  }
279 
280  matrix_type operator()( points_type const& pts ) const
281  {
282  return evaluate( pts );
283  }
284 
286 
290 
295  uint_type degree() const
296  {
297  return nOrder;
298  }
299 
303  self_type const& basis() const
304  {
305  return *this;
306  }
307 
311  std::string familyName() const
312  {
313  return "dubiner.boundary.adapted";
314  }
315 
317 
321 
322 
324 
328 
329 
330  matrix_type coeff() const
331  {
332  return ublas::identity_matrix<value_type>( reference_convex_type::polyDims( nOrder ), M_pts.size2() );
333  }
334 
335 
341  static matrix_type evaluate( points_type const& __pts )
342  {
343  return evaluate( __pts, mpl::int_<nDim>() );
344  }
345 
346  template<typename AE>
347  static vector_matrix_type derivate( ublas::matrix_expression<AE> const& __pts )
348  {
349  return derivate( __pts, mpl::int_<nDim>() );
350  }
351 
352 
354 
355 private:
356 
367  matrix_type
368  static evaluate( points_type const& __pts, mpl::int_<1> )
369  {
370  matrix_type E = M_pfunc.evaluate_1( ublas::row( __pts,0 ) );
371  matrix_type D;
372  D.resize( E.size1(), E.size2() );
373 
374  ublas::row( D, 0 ) = ublas::row( E, 0 );
375  ublas::row( D, 1 ) = ublas::row( E, E.size2()-1 );
376 
377  for ( unsigned int i=1; i < E.size2()-1 ; ++i )
378  {
379  ublas::row( D, i+1 ) = ublas::row( E, i );
380  }
381 
382  return D;
383  }
384 
388  template<typename AE>
389  vector_matrix_type
390  static derivate( ublas::matrix_expression<AE> const& __pts, mpl::int_<1> )
391  {
392  FEELPP_ASSERT( __pts().size1() == 1 )( __pts().size1() )( __pts().size2() ).error( "invalid points" );
393 
394  vector_matrix_type E( 1 );
395  E[0].resize( nOrder+1, __pts().size2() );
396  E[0] = M_pfunc.derivate_1( ublas::row( __pts(),0 ) );
397 
398  vector_matrix_type D( 1 );
399  D[0].resize( nOrder+1, __pts().size2() );
400 
401  ublas::row( D[0], 0 ) = ublas::row( E[0], 0 );
402  ublas::row( D[0], 1 ) = ublas::row( E[0], E[0].size2()-1 );
403 
404  for ( unsigned int i=1; i < E[0].size2()-1 ; ++i )
405  {
406  ublas::row( D[0], i+1 ) = ublas::row( E[0], i );
407  }
408 
409  return D;
410  }
411 
416  static matrix_type evaluate( points_type const& __pts, mpl::int_<2> );
417 
422  template<typename AE>
423  static vector_matrix_type derivate( ublas::matrix_expression<AE> const& __pts, mpl::int_<2> );
424 
429  static matrix_type evaluate( points_type const& __pts, mpl::int_<3> );
430 
435  template<typename AE>
436  static vector_matrix_type derivate( ublas::matrix_expression<AE> const& __pts, mpl::int_<3> );
437 
438 private:
439  reference_convex_type M_refconvex;
440  points_type M_pts;
441  std::vector<points_type> M_pts_face;
442  static principal_type M_pfunc;
443 
444 };
445 template<uint16_type Dim,
446  uint16_type Degree,
447  typename T,
448  template<class> class StoragePolicy>
449 typename BoundaryAdapted<Dim, Degree, T, StoragePolicy>::principal_type
450 BoundaryAdapted<Dim, Degree, T, StoragePolicy>::M_pfunc;
451 
452 template<uint16_type Dim,
453  uint16_type Degree,
454  typename T,
455  template<class> class StoragePolicy>
456 typename BoundaryAdapted<Dim, Degree, T, StoragePolicy>::matrix_type
457 BoundaryAdapted<Dim, Degree, T, StoragePolicy>::evaluate( points_type const& __pts, mpl::int_<2> )
458 {
459  matrix_type res( convex_type::polyDims( nOrder ), __pts.size2() );
460 
461  details::etas<TRIANGLE, value_type> etas( __pts );
462  vector_type eta1s = ublas::row( etas(), 0 );
463  vector_type eta2s = ublas::row( etas(), 1 );
464 
465  matrix_type psi_1( M_pfunc.evaluate_1( eta1s ) );
466  vector_matrix_type psi_2( M_pfunc.evaluate_2( eta2s ) );
467 
468  /* 3 Vertex */
469 
470  ublas::row( res,0 ) = ublas::element_prod( ublas::row( psi_1,0 ) , ublas::row( psi_2[0],0 ) );
471  ublas::row( res,1 ) = ublas::element_prod( ublas::row( psi_1,nOrder ) , ublas::row( psi_2[nOrder],0 ) );
472 
473  vector_type ones( ublas::scalar_vector<value_type>( __pts.size2(), value_type( 1.0 ) ) );
474 
475  ublas::row( res,2 ) = ( ones + eta2s ) / 2.0;
476 
477  /* Global index */
478  uint_type G_i = 2;
479 
480  /* Edges */
482 
483  /* BD = \Gamma_0 : (nOrder-1)*/
484  for ( uint_type q = 1; q < nOrder; ++q )
485  {
486  ++G_i;
487  ublas::row( res,G_i ) = ublas::element_prod( ublas::row( psi_1,nOrder ) , ublas::row( psi_2[nOrder],q ) );
488  }
489 
490  /* AC = \Gamma_1 : (nOrder-1)*/
491  for ( uint_type q = 1; q < nOrder; ++q )
492  {
493  ++G_i;
494  ublas::row( res,G_i ) = ublas::element_prod( ublas::row( psi_1,0 ) , ublas::row( psi_2[0],q ) );
495 
496  if ( q%2==0 )
497  ublas::row( res,G_i )*= value_type( -1.0 );
498  }
499 
500  /* AB = \Gamma_2 : (nOrder-1) */
501  for ( uint_type p = 1; p < nOrder; ++p )
502  {
503  ++G_i;
504  ublas::row( res,G_i ) = ublas::element_prod( ublas::row( psi_1,p ) , ublas::row( psi_2[p],0 ) );
505  }
506 
508 
509 
510  /* Interior : (nOrder-1)(nOrder-2)/2 */
512  for ( uint_type p = 1; p < nOrder; ++p )
513  for ( uint_type q = 1; q < nOrder-p; ++q )
514  {
515  ++G_i;
516  ublas::row( res,G_i ) = ublas::element_prod( ublas::row( psi_1,p ) , ublas::row( psi_2[p],q ) );
517  }
518 
520 
521  return res;
522 }
523 
524 
525 template<uint16_type Dim,
526  uint16_type Degree,
527  typename T,
528  template<class> class StoragePolicy>
529 template<typename AE>
530 typename BoundaryAdapted<Dim, Degree, T, StoragePolicy>::vector_matrix_type
531 BoundaryAdapted<Dim, Degree, T, StoragePolicy>::derivate( ublas::matrix_expression<AE> const& __pts, mpl::int_<2> )
532 {
533  vector_matrix_type res( 2 );
534  res[0].resize( convex_type::polyDims( nOrder ), __pts().size2() );
535  res[1].resize( convex_type::polyDims( nOrder ), __pts().size2() );
536 
537  details::etas<TRIANGLE, value_type> etas( __pts );
538  vector_type eta1s = ublas::row( etas(), 0 );
539  vector_type eta2s = ublas::row( etas(), 1 );
540 
541  matrix_type psi_1( M_pfunc.evaluate_1( eta1s ) );
542  vector_matrix_type psi_2( M_pfunc.evaluate_2( eta2s ) );
543 
544  matrix_type dpsi_1( M_pfunc.derivate_1( eta1s ) );
545  vector_matrix_type dpsi_2( M_pfunc.derivate_2( eta2s ) );
546 
547  vector_type ones( ublas::scalar_vector<value_type>( __pts().size2(), value_type( 1.0 ) ) );
548  vector_type d1( value_type( 2.0 )*ublas::element_div( ones,ones-eta2s ) ); // 2 / (1 - xi_2)
549  vector_type d2( ublas::element_div( ones+eta1s ,ones-eta2s ) ); // 2(1+xi_1) / (1 - xi_2)^2
550 
551  /* 3 Vertex */
552 
553  ublas::row( res[0],0 ) = ublas::element_prod( ublas::row( dpsi_1,0 ) , ublas::row( psi_2[0],0 ) );
554  ublas::row( res[1],0 ) = ublas::element_prod( ublas::row( res[0],0 ) , d2 );
555 
556  ublas::row( res[0],0 ) = ublas::element_prod( ublas::row( res[0],0 ) , d1 );
557  ublas::row( res[1],0 ) += ublas::element_prod( ublas::row( psi_1,0 ) , ublas::row( dpsi_2[0],0 ) );
558 
559  ublas::row( res[0],1 ) = ublas::element_prod( ublas::row( dpsi_1,nOrder ) , ublas::row( psi_2[nOrder],0 ) );
560  ublas::row( res[1],1 ) = ublas::element_prod( ublas::row( res[0],1 ) , d2 );
561 
562  ublas::row( res[0],1 ) = ublas::element_prod( ublas::row( res[0],1 ) , d1 );
563  ublas::row( res[1],1 ) += ublas::element_prod( ublas::row( psi_1,nOrder ) , ublas::row( dpsi_2[nOrder],0 ) );
564 
565  ublas::row( res[0],2 ) = ublas::scalar_vector<value_type>( __pts().size2(),value_type( 0.0 ) );
566  ublas::row( res[1],2 ) = value_type( 0.5 ) * ones;
567 
568  /* Global index */
569  uint_type G_i = 2;
570 
571  /* Edges */
573 
574  /* BD : (nOrder-1)*/
575  for ( uint_type q = 1; q < nOrder; ++q )
576  {
577  ++G_i;
578  ublas::row( res[0],G_i ) = ublas::element_prod( ublas::row( dpsi_1,nOrder ) , ublas::row( psi_2[nOrder],q ) );
579  ublas::row( res[1],G_i ) = ublas::element_prod( ublas::row( res[0],G_i ) , d2 );
580 
581  ublas::row( res[0],G_i ) = ublas::element_prod( ublas::row( res[0],G_i ) , d1 );
582  ublas::row( res[1],G_i ) += ublas::element_prod( ublas::row( psi_1,nOrder ) , ublas::row( dpsi_2[nOrder],q ) );
583  }
584 
585  /* AC : (nOrder-1)*/
586  for ( uint_type q = 1; q < nOrder; ++q )
587  {
588  ++G_i;
589  ublas::row( res[0],G_i ) = ublas::element_prod( ublas::row( dpsi_1,0 ) , ublas::row( psi_2[0],q ) );
590  ublas::row( res[1],G_i ) = ublas::element_prod( ublas::row( res[0],G_i ) , d2 );
591 
592  ublas::row( res[0],G_i ) = ublas::element_prod( ublas::row( res[0],G_i ) , d1 );
593  ublas::row( res[1],G_i ) += ublas::element_prod( ublas::row( psi_1,0 ) , ublas::row( dpsi_2[0],q ) );
594 
595  if ( q%2==0 )
596  {
597  ublas::row( res[0],G_i )*= value_type( -1.0 );
598  ublas::row( res[1],G_i )*= value_type( -1.0 );
599  }
600  }
601 
602  /* AB : (nOrder-1) */
603  for ( uint_type p = 1; p < nOrder; ++p )
604  {
605  ++G_i;
606  ublas::row( res[0],G_i ) = ublas::element_prod( ublas::row( dpsi_1,p ) , ublas::row( psi_2[p],0 ) );
607  ublas::row( res[1],G_i ) = ublas::element_prod( ublas::row( res[0],G_i ) , d2 );
608 
609  ublas::row( res[0],G_i ) = ublas::element_prod( ublas::row( res[0],G_i ) , d1 );
610  ublas::row( res[1],G_i ) += ublas::element_prod( ublas::row( psi_1,p ) , ublas::row( dpsi_2[p],0 ) );
611  }
612 
613  /* Interior : (nOrder-1)(nOrder-2)/2 */
615  for ( uint_type p = 1; p < nOrder; ++p )
616  for ( uint_type q = 1; q < nOrder-p; ++q )
617  {
618  ++G_i;
619  ublas::row( res[0],G_i ) = ublas::element_prod( ublas::row( dpsi_1,p ) , ublas::row( psi_2[p],q ) );
620  ublas::row( res[1],G_i ) = ublas::element_prod( ublas::row( res[0],G_i ) , d2 );
621 
622  ublas::row( res[0],G_i ) = ublas::element_prod( ublas::row( res[0],G_i ) , d1 );
623  ublas::row( res[1],G_i ) += ublas::element_prod( ublas::row( psi_1,p ) , ublas::row( dpsi_2[p],q ) );
624  }
625 
627 
628 
629  return res;
630 }
631 
632 template<uint16_type Dim,
633  uint16_type Degree,
634  typename T,
635  template<class> class StoragePolicy>
636 typename BoundaryAdapted<Dim, Degree, T, StoragePolicy>::matrix_type
637 BoundaryAdapted<Dim, Degree, T, StoragePolicy>::evaluate( points_type const& __pts, mpl::int_<3> )
638 {
639  matrix_type res( convex_type::polyDims( nOrder ), __pts.size2() );
640 
641  FEELPP_ASSERT( __pts.size1() == 3 )( __pts.size1() ).error( "invalid space dimension" );
642 
643  details::etas<TETRAHEDRON, value_type> etas( __pts );
644  vector_type eta1s = ublas::row( etas(), 0 );
645  vector_type eta2s = ublas::row( etas(), 1 );
646  vector_type eta3s = ublas::row( etas(), 2 );
647 
648  matrix_type psi_1( M_pfunc.evaluate_1( eta1s ) );
649  vector_matrix_type psi_2( M_pfunc.evaluate_2( eta2s ) );
650  vector_vector_matrix_type psi_3( M_pfunc.evaluate_3( eta3s ) );
651 
652  ublas::scalar_vector<value_type> ones( eta3s.size(), value_type( 1.0 ) );
653 
654  /* 4 Vertex */
655 
656  ublas::row( res,0 ) =
657  ublas::element_prod(
658  ublas::element_prod( ublas::row( psi_1,0 ) , ublas::row( psi_2[0],0 ) ),
659  ublas::row( psi_3[0][0],0 )
660  );
661  ublas::row( res,1 ) =
662  ublas::element_prod(
663  ublas::element_prod( ublas::row( psi_1,nOrder ) , ublas::row( psi_2[nOrder],0 ) ),
664  ublas::row( psi_3[nOrder][0],0 )
665  );
666 
667  ublas::row( res,2 ) =
668  ublas::element_prod(
669  ublas::row( psi_2[0],nOrder ), ublas::row( psi_3[0][nOrder],0 )
670  );
671 
672  ublas::row( res,3 ) = ( ones + eta3s ) / 2.0;
673 
674 
675  /* Global index */
676  uint_type G_i = 3;
677 
678  /* Edges */
680  /* AB : (nOrder-1) */
681  for ( uint_type q = 1; q < nOrder; ++q )
682  {
683  ++G_i;
684  ublas::row( res,G_i ) =
685  ublas::element_prod(
686  ublas::element_prod( ublas::row( psi_1,nOrder ) , ublas::row( psi_2[nOrder],q ) ),
687  ublas::row( psi_3[nOrder][q],0 )
688  );
689  }
690 
691  /* AC : (nOrder-1)*/
692  for ( uint_type q = 1; q < nOrder; ++q )
693  {
694  ++G_i;
695  ublas::row( res,G_i ) =
696  ublas::element_prod(
697  ublas::element_prod( ublas::row( psi_1,0 ) , ublas::row( psi_2[0],q ) ),
698  ublas::row( psi_3[0][q],0 )
699  );
700  /* Ensure Consistent anticlockwise orientation of the
701  edges */
702 
703  if ( q%2==0 )
704  {
705  ublas::row( res,G_i )*= value_type( -1.0 );
706  }
707  }
708 
709  /* BC : (nOrder-1)*/
710  for ( uint_type p = 1; p < nOrder; ++p )
711  {
712  ++G_i;
713  ublas::row( res,G_i ) =
714  ublas::element_prod(
715  ublas::element_prod( ublas::row( psi_1,p ) , ublas::row( psi_2[p],0 ) ),
716  ublas::row( psi_3[p][0],0 )
717  );
718  }
719 
720  /* AD : (nOrder-1)*/
721  for ( uint_type r = 1; r < nOrder; ++r )
722  {
723  ++G_i;
724  ublas::row( res,G_i ) =
725  ublas::element_prod(
726  ublas::element_prod( ublas::row( psi_1,0 ) , ublas::row( psi_2[0],0 ) ),
727  ublas::row( psi_3[0][0],r )
728  );
729  }
730 
731  /* BD : (nOrder-1)*/
732  for ( uint_type r = 1; r < nOrder; ++r )
733  {
734  ++G_i;
735  ublas::row( res,G_i ) =
736  ublas::element_prod(
737  ublas::element_prod( ublas::row( psi_1,nOrder ) , ublas::row( psi_2[nOrder],0 ) ),
738  ublas::row( psi_3[nOrder][0],r )
739  );
740  }
741 
742  /* CD : (nOrder-1)*/
743  for ( uint_type r = 1; r < nOrder; ++r )
744  {
745  ++G_i;
746  ublas::row( res,G_i ) =
747  ublas::element_prod( ublas::row( psi_2[0],nOrder ),ublas::row( psi_3[0][nOrder],r ) );
748  }
749 
751 
752 
753  /* Faces */
755 
756  /* BCE : (N-1)(N-2)/2 */
757  for ( uint_type q = 1; q < nOrder; ++q )
758  for ( uint_type r = 1; r < nOrder-q; ++r )
759  {
760  ++G_i;
761  ublas::row( res,G_i ) =
762  ublas::element_prod(
763  ublas::element_prod( ublas::row( psi_1,nOrder ) , ublas::row( psi_2[nOrder],q ) ),
764  ublas::row( psi_3[nOrder][q],r )
765  );
766  }
767 
768  /* ACE : (N-1)(N-2)/2 */
769  for ( uint_type q = 1; q < nOrder; ++q )
770  for ( uint_type r = 1; r < nOrder-q; ++r )
771  {
772  ++G_i;
773  ublas::row( res,G_i ) =
774  ublas::element_prod(
775  ublas::element_prod( ublas::row( psi_1,0 ) , ublas::row( psi_2[0],q ) ),
776  ublas::row( psi_3[0][q],r )
777  );
778  }
779 
780  /* ABE : (N-1)(N-2)/2 */
781  for ( uint_type p = 1; p < nOrder; ++p )
782  for ( uint_type r = 1; r < nOrder-p; ++r )
783  {
784  ++G_i;
785  ublas::row( res,G_i ) =
786  ublas::element_prod(
787  ublas::element_prod( ublas::row( psi_1,p ) , ublas::row( psi_2[p],0 ) ),
788  ublas::row( psi_3[p][0],r )
789  );
790  }
791 
792 
793  /* ABC : (N-1)(N-2)/2 */
794  for ( uint_type p = 1; p < nOrder; ++p )
795  for ( uint_type q = 1; q < nOrder-p; ++q )
796  {
797  ++G_i;
798  ublas::row( res,G_i ) =
799  ublas::element_prod(
800  ublas::element_prod( ublas::row( psi_1,p ) , ublas::row( psi_2[p],q ) ),
801  ublas::row( psi_3[p][q],0 )
802  );
803  }
804 
806 
807 
808  /* Interior : (N-1)(N-2)(N-3)/6 */
809 
810  for ( uint_type p = 1; p < nOrder-2; ++p )
811  for ( uint_type q = 1; q < nOrder-1-p; ++q )
812  for ( uint_type r = 1; r < nOrder-p-q; ++r )
813  {
814  ++G_i;
815  ublas::row( res,G_i ) =
816  ublas::element_prod(
817  ublas::element_prod( ublas::row( psi_1,p ) , ublas::row( psi_2[p],q ) ),
818  ublas::row( psi_3[p][q],r )
819  );
820  }
821 
822  return res;
823 }
824 
825 template<uint16_type Dim,
826  uint16_type Degree,
827  typename T,
828  template<class> class StoragePolicy>
829 template<typename AE>
830 typename BoundaryAdapted<Dim, Degree, T, StoragePolicy>::vector_matrix_type
831 BoundaryAdapted<Dim, Degree, T, StoragePolicy>::derivate( ublas::matrix_expression<AE> const& __pts, mpl::int_<3> )
832 {
833  vector_matrix_type res( 3 );
834  res[0].resize( convex_type::polyDims( nOrder ), __pts().size2() );
835  res[1].resize( convex_type::polyDims( nOrder ), __pts().size2() );
836  res[2].resize( convex_type::polyDims( nOrder ), __pts().size2() );
837 
838  FEELPP_ASSERT( __pts().size1() == 3 )( __pts().size1() ).error( "invalid space dimension" );
839 
840  details::etas<TETRAHEDRON, value_type> etas( __pts );
841  vector_type eta1s = ublas::row( etas(), 0 );
842  vector_type eta2s = ublas::row( etas(), 1 );
843  vector_type eta3s = ublas::row( etas(), 2 );
844 
845  matrix_type psi_1( M_pfunc.evaluate_1( eta1s ) );
846  matrix_type dpsi_1( M_pfunc.derivate_1( eta1s ) );
847 
848  vector_matrix_type psi_2( M_pfunc.evaluate_2( eta2s ) );
849  vector_matrix_type dpsi_2( M_pfunc.derivate_2( eta2s ) );
850 
851  vector_vector_matrix_type psi_3( M_pfunc.evaluate_3( eta3s ) );
852  vector_vector_matrix_type dpsi_3( M_pfunc.derivate_3( eta3s ) );
853 
854  ublas::scalar_vector<value_type> ones( __pts().size2(), value_type( 1.0 ) );
855 
863  matrix_type Jac( 9, __pts().size2() );
864 
867  ublas::row( Jac,0 ) =
868  ublas::element_div( ( -2.0 )*ones, ( ublas::row( __pts(), 1 ) + ublas::row( __pts(), 2 ) ) );
869 
870  ublas::row( Jac,1 ) =
871  ublas::element_div( 2.0*( ones + ublas::row( __pts(), 0 ) ) , ( ublas::row( __pts(), 1 ) + ublas::row( __pts(), 2 ) ) );
872 
873  ublas::row( Jac,1 ) =
874  ublas::element_div( ublas::row( Jac,1 ), ublas::row( __pts(), 1 )+ublas::row( __pts(), 2 ) );
875 
876  ublas::row( Jac,2 ) =
877  ublas::row( Jac,1 );
878 
881  ublas::row( Jac,3 ) =
882  ublas::scalar_vector<value_type>( __pts().size2(), value_type( 0.0 ) );
883 
884  ublas::row( Jac,4 ) =
885  ublas::element_div( 2.0*ones, ( ones - ublas::row( __pts(), 2 ) ) );
886 
887  ublas::row( Jac,5 ) =
888  ublas::element_div( 2.0*( ones + ublas::row( __pts(), 1 ) ) , ( ones - ublas::row( __pts(), 2 ) ) );
889 
890  ublas::row( Jac,5 ) =
891  ublas::element_div( ublas::row( Jac,5 ), ones - ublas::row( __pts(), 2 ) );
892 
895  ublas::row( Jac,6 ) =
896  ublas::row( Jac,3 );
897 
898  ublas::row( Jac,7 ) =
899  ublas::row( Jac,3 );
900 
901  ublas::row( Jac,8 ) = ones;
906  matrix_type d1 ( res[1].size1(), res[1].size2() );
907  matrix_type d2 ( res[1].size1(), res[1].size2() );
908  matrix_type d3 ( res[1].size1(), res[1].size2() );
909 
910  /* 4 Vertex */
911 
912  ublas::row( d1,0 ) =
913  ublas::element_prod(
914  ublas::element_prod( ublas::row( dpsi_1,0 ) , ublas::row( psi_2[0],0 ) ),
915  ublas::row( psi_3[0][0],0 )
916  );
917  ublas::row( d2,0 ) =
918  ublas::element_prod(
919  ublas::element_prod( ublas::row( psi_1,0 ) , ublas::row( dpsi_2[0],0 ) ),
920  ublas::row( psi_3[0][0],0 )
921  );
922  ublas::row( d3,0 ) =
923  ublas::element_prod(
924  ublas::element_prod( ublas::row( psi_1,0 ) , ublas::row( psi_2[0],0 ) ),
925  ublas::row( dpsi_3[0][0],0 )
926  );
927 
928  ublas::row( d1,1 ) =
929  ublas::element_prod(
930  ublas::element_prod( ublas::row( dpsi_1,nOrder ) , ublas::row( psi_2[nOrder],0 ) ),
931  ublas::row( psi_3[nOrder][0],0 )
932  );
933  ublas::row( d2,1 ) =
934  ublas::element_prod(
935  ublas::element_prod( ublas::row( psi_1,nOrder ) , ublas::row( dpsi_2[nOrder],0 ) ),
936  ublas::row( psi_3[nOrder][0],0 )
937  );
938  ublas::row( d3,1 ) =
939  ublas::element_prod(
940  ublas::element_prod( ublas::row( psi_1,nOrder ) , ublas::row( psi_2[nOrder],0 ) ),
941  ublas::row( dpsi_3[nOrder][0],0 )
942  );
943 
944  ublas::row( d1,2 ) = ublas::scalar_vector<value_type>( __pts().size2(), value_type( 0.0 ) );
945 
946  ublas::row( d2,2 ) =
947  ublas::element_prod(
948  ublas::row( dpsi_2[0],nOrder ), ublas::row( psi_3[0][nOrder],0 )
949  );
950 
951  ublas::row( d3,2 ) =
952  ublas::element_prod(
953  ublas::row( psi_2[0],nOrder ), ublas::row( dpsi_3[0][nOrder],0 )
954  );
955 
956  ublas::row( d1,3 ) =
957  ublas::scalar_vector<value_type>( __pts().size2(), value_type( 0.0 ) );
958 
959  ublas::row( d2,3 ) =
960  ublas::scalar_vector<value_type>( __pts().size2(), value_type( 0.0 ) );
961 
962  ublas::row( d3,3 ) =
963  ublas::scalar_vector<value_type>( __pts().size2(), value_type( 0.5 ) );
964 
965  /* Global index */
966  uint_type G_i = 3;
967 
968  /* Edges */
970  /* AB : (nOrder-1) */
971  for ( uint_type q = 1; q < nOrder; ++q )
972  {
973  ++G_i;
974  ublas::row( d1,G_i ) =
975  ublas::element_prod(
976  ublas::element_prod( ublas::row( dpsi_1,nOrder ) , ublas::row( psi_2[nOrder],q ) ),
977  ublas::row( psi_3[nOrder][q],0 )
978  );
979  ublas::row( d2,G_i ) =
980  ublas::element_prod(
981  ublas::element_prod( ublas::row( psi_1,nOrder ) , ublas::row( dpsi_2[nOrder],q ) ),
982  ublas::row( psi_3[nOrder][q],0 )
983  );
984  ublas::row( d3,G_i ) =
985  ublas::element_prod(
986  ublas::element_prod( ublas::row( psi_1,nOrder ) , ublas::row( psi_2[nOrder],q ) ),
987  ublas::row( dpsi_3[nOrder][q],0 )
988  );
989  }
990 
991  /* AC : (nOrder-1)*/
992  for ( uint_type q = 1; q < nOrder; ++q )
993  {
994  ++G_i;
995  ublas::row( d1,G_i ) =
996  ublas::element_prod(
997  ublas::element_prod( ublas::row( dpsi_1,0 ) , ublas::row( psi_2[0],q ) ),
998  ublas::row( psi_3[0][q],0 )
999  );
1000  ublas::row( d2,G_i ) =
1001  ublas::element_prod(
1002  ublas::element_prod( ublas::row( psi_1,0 ) , ublas::row( dpsi_2[0],q ) ),
1003  ublas::row( psi_3[0][q],0 )
1004  );
1005  ublas::row( d3,G_i ) =
1006  ublas::element_prod(
1007  ublas::element_prod( ublas::row( psi_1,0 ) , ublas::row( psi_2[0],q ) ),
1008  ublas::row( dpsi_3[0][q],0 )
1009  );
1010 
1011  /* Ensure Consistent anticlockwise orientation of the
1012  edges */
1013  if ( q%2==0 )
1014  {
1015  ublas::row( d1,G_i )*= value_type( -1.0 );
1016  ublas::row( d2,G_i )*= value_type( -1.0 );
1017  ublas::row( d3,G_i )*= value_type( -1.0 );
1018  }
1019  }
1020 
1021  /* BC : (nOrder-1)*/
1022  for ( uint_type p = 1; p < nOrder; ++p )
1023  {
1024  ++G_i;
1025  ublas::row( d1,G_i ) =
1026  ublas::element_prod(
1027  ublas::element_prod( ublas::row( dpsi_1,p ) , ublas::row( psi_2[p],0 ) ),
1028  ublas::row( psi_3[p][0],0 )
1029  );
1030  ublas::row( d2,G_i ) =
1031  ublas::element_prod(
1032  ublas::element_prod( ublas::row( psi_1,p ) , ublas::row( dpsi_2[p],0 ) ),
1033  ublas::row( psi_3[p][0],0 )
1034  );
1035  ublas::row( d3,G_i ) =
1036  ublas::element_prod(
1037  ublas::element_prod( ublas::row( psi_1,p ) , ublas::row( psi_2[p],0 ) ),
1038  ublas::row( dpsi_3[p][0],0 )
1039  );
1040  }
1041 
1042  /* AD : (nOrder-1)*/
1043  for ( uint_type r = 1; r < nOrder; ++r )
1044  {
1045  ++G_i;
1046  ublas::row( d1,G_i ) =
1047  ublas::element_prod(
1048  ublas::element_prod( ublas::row( dpsi_1,0 ) , ublas::row( psi_2[0],0 ) ),
1049  ublas::row( psi_3[0][0],r )
1050  );
1051  ublas::row( d2,G_i ) =
1052  ublas::element_prod(
1053  ublas::element_prod( ublas::row( psi_1,0 ) , ublas::row( dpsi_2[0],0 ) ),
1054  ublas::row( psi_3[0][0],r )
1055  );
1056  ublas::row( d3,G_i ) =
1057  ublas::element_prod(
1058  ublas::element_prod( ublas::row( psi_1,0 ) , ublas::row( psi_2[0],0 ) ),
1059  ublas::row( dpsi_3[0][0],r )
1060  );
1061  }
1062 
1063  /* BD : (nOrder-1)*/
1064  for ( uint_type r = 1; r < nOrder; ++r )
1065  {
1066  ++G_i;
1067  ublas::row( d1,G_i ) =
1068  ublas::element_prod(
1069  ublas::element_prod( ublas::row( dpsi_1,nOrder ) , ublas::row( psi_2[nOrder],0 ) ),
1070  ublas::row( psi_3[nOrder][0],r )
1071  );
1072  ublas::row( d2,G_i ) =
1073  ublas::element_prod(
1074  ublas::element_prod( ublas::row( psi_1,nOrder ) , ublas::row( dpsi_2[nOrder],0 ) ),
1075  ublas::row( psi_3[nOrder][0],r )
1076  );
1077  ublas::row( d3,G_i ) =
1078  ublas::element_prod(
1079  ublas::element_prod( ublas::row( psi_1,nOrder ) , ublas::row( psi_2[nOrder],0 ) ),
1080  ublas::row( dpsi_3[nOrder][0],r )
1081  );
1082  }
1083 
1084  /* CD : (nOrder-1)*/
1085  for ( uint_type r = 1; r < nOrder; ++r )
1086  {
1087  ++G_i;
1088  ublas::row( d1,G_i ) = ublas::scalar_vector<value_type>( __pts().size2(), value_type( 0.0 ) );
1089  ublas::row( d2,G_i ) =
1090  ublas::element_prod( ublas::row( dpsi_2[0],nOrder ),ublas::row( psi_3[0][nOrder],r ) );
1091  ublas::row( d3,G_i ) =
1092  ublas::element_prod( ublas::row( psi_2[0],nOrder ),ublas::row( dpsi_3[0][nOrder],r ) );
1093  }
1094 
1096 
1097 
1098  /* Faces */
1100 
1101  /* BCE : (N-1)(N-2)/2 */
1102  for ( uint_type q = 1; q < nOrder; ++q )
1103  for ( uint_type r = 1; r < nOrder-q; ++r )
1104  {
1105  ++G_i;
1106  ublas::row( d1,G_i ) =
1107  ublas::element_prod(
1108  ublas::element_prod( ublas::row( dpsi_1,nOrder ) , ublas::row( psi_2[nOrder],q ) ),
1109  ublas::row( psi_3[nOrder][q],r )
1110  );
1111  ublas::row( d2,G_i ) =
1112  ublas::element_prod(
1113  ublas::element_prod( ublas::row( psi_1,nOrder ) , ublas::row( dpsi_2[nOrder],q ) ),
1114  ublas::row( psi_3[nOrder][q],r )
1115  );
1116  ublas::row( d3,G_i ) =
1117  ublas::element_prod(
1118  ublas::element_prod( ublas::row( psi_1,nOrder ) , ublas::row( psi_2[nOrder],q ) ),
1119  ublas::row( dpsi_3[nOrder][q],r )
1120  );
1121  }
1122 
1123  /* ACE : (N-1)(N-2)/2 */
1124  for ( uint_type q = 1; q < nOrder; ++q )
1125  for ( uint_type r = 1; r < nOrder-q; ++r )
1126  {
1127  ++G_i;
1128  ublas::row( d1,G_i ) =
1129  ublas::element_prod(
1130  ublas::element_prod( ublas::row( dpsi_1,0 ) , ublas::row( psi_2[0],q ) ),
1131  ublas::row( psi_3[0][q],r )
1132  );
1133  ublas::row( d2,G_i ) =
1134  ublas::element_prod(
1135  ublas::element_prod( ublas::row( psi_1,0 ) , ublas::row( dpsi_2[0],q ) ),
1136  ublas::row( psi_3[0][q],r )
1137  );
1138  ublas::row( d3,G_i ) =
1139  ublas::element_prod(
1140  ublas::element_prod( ublas::row( psi_1,0 ) , ublas::row( psi_2[0],q ) ),
1141  ublas::row( dpsi_3[0][q],r )
1142  );
1143  }
1144 
1145  /* ABE : (N-1)(N-2)/2 */
1146  for ( uint_type p = 1; p < nOrder; ++p )
1147  for ( uint_type r = 1; r < nOrder-p; ++r )
1148  {
1149  ++G_i;
1150  ublas::row( d1,G_i ) =
1151  ublas::element_prod(
1152  ublas::element_prod( ublas::row( dpsi_1,p ) , ublas::row( psi_2[p],0 ) ),
1153  ublas::row( psi_3[p][0],r )
1154  );
1155  ublas::row( d2,G_i ) =
1156  ublas::element_prod(
1157  ublas::element_prod( ublas::row( psi_1,p ) , ublas::row( dpsi_2[p],0 ) ),
1158  ublas::row( psi_3[p][0],r )
1159  );
1160  ublas::row( d3,G_i ) =
1161  ublas::element_prod(
1162  ublas::element_prod( ublas::row( psi_1,p ) , ublas::row( psi_2[p],0 ) ),
1163  ublas::row( dpsi_3[p][0],r )
1164  );
1165  }
1166 
1167 
1168  /* ABC : (N-1)(N-2)/2 */
1169  for ( uint_type p = 1; p < nOrder; ++p )
1170  for ( uint_type q = 1; q < nOrder-p; ++q )
1171  {
1172  ++G_i;
1173  ublas::row( d1,G_i ) =
1174  ublas::element_prod(
1175  ublas::element_prod( ublas::row( dpsi_1,p ) , ublas::row( psi_2[p],q ) ),
1176  ublas::row( psi_3[p][q],0 )
1177  );
1178  ublas::row( d2,G_i ) =
1179  ublas::element_prod(
1180  ublas::element_prod( ublas::row( psi_1,p ) , ublas::row( dpsi_2[p],q ) ),
1181  ublas::row( psi_3[p][q],0 )
1182  );
1183  ublas::row( d3,G_i ) =
1184  ublas::element_prod(
1185  ublas::element_prod( ublas::row( psi_1,p ) , ublas::row( psi_2[p],q ) ),
1186  ublas::row( dpsi_3[p][q],0 )
1187  );
1188  }
1189 
1191 
1192 
1193  /* Interior : (N-1)(N-2)(N-3)/6 */
1194 
1195  for ( uint_type p = 1; p < nOrder-2; ++p )
1196  for ( uint_type q = 1; q < nOrder-1-p; ++q )
1197  for ( uint_type r = 1; r < nOrder-p-q; ++r )
1198  {
1199  ++G_i;
1200  ublas::row( d1,G_i ) =
1201  ublas::element_prod( ublas::element_prod( ublas::row( dpsi_1,p ) , ublas::row( psi_2[p],q ) ), ublas::row( psi_3[p][q],r ) );
1202  ublas::row( d2,G_i ) =
1203  ublas::element_prod( ublas::element_prod( ublas::row( psi_1,p ) , ublas::row( dpsi_2[p],q ) ), ublas::row( psi_3[p][q],r ) );
1204  ublas::row( d3,G_i ) =
1205  ublas::element_prod( ublas::element_prod( ublas::row( psi_1,p ) , ublas::row( psi_2[p],q ) ), ublas::row( dpsi_3[p][q],r ) );
1206  }
1207 
1208 
1211 #if defined (FEELPP_HAS_QD_REAL)
1212 
1213  /* Warning : if qd are enable we can not use
1214  * the ublas::diagonal_matrix type
1215  */
1216 
1217  matrix_type j0( __pts().size2(), __pts().size2() );
1218  matrix_type j1( __pts().size2(), __pts().size2() );
1219  matrix_type j2( __pts().size2(), __pts().size2() );
1220  matrix_type j4( __pts().size2(), __pts().size2() );
1221  matrix_type j5( __pts().size2(), __pts().size2() );
1222 
1223  for ( uint_type i=0; i < __pts().size2(); ++i )
1224  for ( uint_type j=i; j < __pts().size2(); ++j )
1225  {
1226  if ( j==i )
1227  {
1228  j0( i,i ) = ublas::row( Jac,0 )( i );
1229  j1( i,i ) = ublas::row( Jac,1 )( i );
1230  j2( i,i ) = ublas::row( Jac,2 )( i );
1231  j4( i,i ) = ublas::row( Jac,4 )( i );
1232  j5( i,i ) = ublas::row( Jac,5 )( i );
1233  }
1234 
1235  else
1236  {
1237  j0( i,j ) = value_type( 0.0 );
1238  j1( i,j ) = value_type( 0.0 );
1239  j2( i,j ) = value_type( 0.0 );
1240  j4( i,j ) = value_type( 0.0 );
1241  j5( i,j ) = value_type( 0.0 );
1242  j0( j,i ) = value_type( 0.0 );
1243  j1( j,i ) = value_type( 0.0 );
1244  j2( j,i ) = value_type( 0.0 );
1245  j4( j,i ) = value_type( 0.0 );
1246  j5( j,i ) = value_type( 0.0 );
1247  }
1248  }
1249 
1250 #else
1251  ublas::diagonal_matrix<value_type> j0( __pts().size2() );
1252  ublas::diagonal_matrix<value_type> j1( __pts().size2() );
1253  ublas::diagonal_matrix<value_type> j2( __pts().size2() );
1254  ublas::diagonal_matrix<value_type> j4( __pts().size2() );
1255  ublas::diagonal_matrix<value_type> j5( __pts().size2() );
1256 
1257  for ( uint_type i=0; i < __pts().size2(); ++i )
1258  {
1259  j0( i,i ) = ublas::row( Jac,0 )( i );
1260  j1( i,i ) = ublas::row( Jac,1 )( i );
1261  j2( i,i ) = ublas::row( Jac,2 )( i );
1262  j4( i,i ) = ublas::row( Jac,4 )( i );
1263  j5( i,i ) = ublas::row( Jac,5 )( i );
1264  }
1265 
1266 #endif
1267 
1268  res[0] = ublas::prod( d1, j0 );
1269  res[1] = ublas::prod( d1, j1 ) + ublas::prod( d2, j4 );
1270  res[2] = ublas::prod( d1, j2 ) + ublas::prod( d2, j5 ) + d3;
1271 
1272  return res;
1273 }
1274 
1275 
1276 }
1277 #endif /* __BoundAdapted_H */

Generated on Sun Oct 20 2013 08:24:54 for Feel++ by doxygen 1.8.4