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
tensorisedboundadapted.hpp
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 __TensorisedBoundadapted_H
30 #define __TensorisedBoundadapted_H 1
31 
32 
33 #include <boost/lambda/if.hpp>
35 #include <feel/feelalg/glas.hpp>
36 
38 
39 namespace Feel
40 {
41 
60 template<uint16_type Dim,
61  uint16_type Degree,
62  typename T = double,
63  template<class> class StoragePolicy = StorageUBlas>
65 {
66 public:
67 
68  static const uint16_type nDim = Dim;
69  static const uint16_type nOrder = Degree;
70  static const uint16_type nConvexOrder = nDim+nOrder+1;
71 
75 
77 
78  /*
79  * for now can be used only as a basis but we might be interested
80  * to have then expressed in other basis like in the moment basis
81  */
82  typedef self_type basis_type;
83 
84  /*
85  * numerical type
86  */
87  typedef T value_type;
88  typedef int16_type int_type;
89  typedef uint16_type uint_type;
90 
91  /*
92  * Geometry where the polynomials are defined and constructed
93  */
94  typedef Hypercube<nDim, nConvexOrder, nDim> convex_type;
96  typedef typename reference_convex_type::points_type points_type;
97 
98  static const uint16_type numVertices = reference_convex_type::numVertices;
99  static const uint16_type numFaces = reference_convex_type::numFaces;
100 
101  template<uint16_type order>
102  struct Convex
103  {
104  typedef Hypercube<nDim, order, nDim> type;
105  };
106 
107  /*
108  * storage policy
109  */
110  typedef StoragePolicy<value_type> storage_policy;
111  typedef typename storage_policy::vector_type vector_type;
112  typedef typename storage_policy::matrix_type matrix_type;
113  typedef typename storage_policy::vector_matrix_type vector_matrix_type;
114  typedef typename storage_policy::vector_vector_matrix_type vector_vector_matrix_type;
115  typedef typename storage_policy::matrix_node_type matrix_node_type;
116  typedef typename storage_policy::node_type node_type;
117 
119 
123 
125  :
126  M_refconvex(),
127  M_pts( nDim, numVertices ),
128  M_pts_face( numVertices )
129  {
130  PointSetEquiSpaced<convex_type, nOrder, value_type> pts;
131 
132  // only the points associated with the vertices
133  M_pts = pts.pointsByEntity( 0 );
134 
135  // get the points for each faces
136  for ( uint16_type e = M_refconvex.entityRange( nDim-1 ).begin();
137  e < M_refconvex.entityRange( nDim-1 ).end();
138  ++e )
139  {
140  M_pts_face[e] = pts.pointsBySubEntity( nDim-1, e, 1 );
141  }
142 
143  matrix_type A( ublas::trans( evaluate( pts.points() ) ) );
144  matrix_type D = ublas::identity_matrix<value_type>( A.size1(), A.size2() );
145  LU<matrix_type> lu( A );
146  matrix_type C = lu.solve( D );
147  vector_matrix_type d ( derivate( pts.points() ) );
148  M_D.resize( d.size() );
149 
150  for ( size_type i = 0; i < d.size(); ++i )
151  {
152  M_D[i] = ublas::prod( d[i], C );
153  glas::clean( M_D[i] );
154  }
155  }
156 
158  :
159  M_refconvex( d.M_refconvex ),
160  M_pts( d.M_pts ),
161  M_pts_face( d.M_pts_face )
162  {
163  }
164 
166  {
167  }
168 
170 
175  points_type points()
176  {
177  return M_pts;
178  }
179 
183  points_type const& points( int f ) const
184  {
185  return M_pts_face[f];
186  }
187 
191 
192  self_type const& operator=( self_type const& d )
193  {
194  if ( this != &d )
195  {
196  M_pts = d.M_pts;
197  M_D = d.M_D;
198  }
199 
200  return *this;
201  }
202 
203  matrix_type operator()( node_type const& pt ) const
204  {
205  points_type pts( pt.size(), 1 );
206  ublas::column( pts, 0 ) = pt;
207  return evaluate( pts );
208  }
209 
210  matrix_type operator()( points_type const& pts ) const
211  {
212  return evaluate( pts );
213  }
214 
216 
220 
225  uint_type degree() const
226  {
227  return nOrder;
228  }
229 
233  self_type const& basis() const
234  {
235  return *this;
236  }
237 
241  std::string familyName() const
242  {
243  return "tensorizedboundaryadapted";
244  }
245 
246 
248 
252 
253 
255 
259 
260 
261  matrix_type coeff() const
262  {
263  return ublas::identity_matrix<value_type>( reference_convex_type::polyDims( nOrder ), M_pts.size2() );
264  }
265 
266 
272  matrix_type evaluate( points_type const& __pts ) const
273  {
274  return evaluate( __pts, mpl::int_<nDim>() );
275  }
276 
277  template<typename AE>
278  vector_matrix_type derivate( ublas::matrix_expression<AE> const& __pts ) const
279  {
280  return derivate( __pts, mpl::int_<nDim>() );
281  }
282 
289  matrix_type const& d( uint16_type i )
290  {
291  return M_D[i];
292  }
293 
300  matrix_type const& derivate( uint16_type i )
301  {
302  return M_D[i];
303  }
304 
306 
307 private:
308 
319  matrix_type
320  evaluate( points_type const& __pts, mpl::int_<1> ) const
321  {
322  return M_pfunc.evaluate_1( ublas::row( __pts,0 ) );
323  }
324 
328  template<typename AE>
329  vector_matrix_type
330  derivate( ublas::matrix_expression<AE> const& __pts, mpl::int_<1> ) const
331  {
332  FEELPP_ASSERT( __pts().size1() == 1 )( __pts().size1() )( __pts().size2() ).error( "invalid points" );
333 
334  vector_matrix_type D( 1 );
335  D[0].resize( nOrder+1, __pts().size2() );
336  D[0] = M_pfunc.derivate_1( ublas::row( __pts(),0 ) );
337 
338  return D;
339  }
340 
345  matrix_type evaluate( points_type const& __pts, mpl::int_<2> ) const;
346 
351  template<typename AE>
352  vector_matrix_type derivate( ublas::matrix_expression<AE> const& __pts, mpl::int_<2> ) const;
353 
358  matrix_type evaluate( points_type const& __pts, mpl::int_<3> ) const;
359 
364  template<typename AE>
365  vector_matrix_type derivate( ublas::matrix_expression<AE> const& __pts, mpl::int_<3> ) const;
366 
367 private:
368  reference_convex_type M_refconvex;
369  points_type M_pts;
370  std::vector<points_type> M_pts_face;
371 
372  Principal<Degree, T, StoragePolicy> M_pfunc;
373 
374  vector_matrix_type M_D;
375 };
376 
377 template<uint16_type Dim,
378  uint16_type Degree,
379  typename T,
380  template<class> class StoragePolicy>
381 typename TensorisedBoundaryAdapted<Dim, Degree, T, StoragePolicy>::matrix_type
382 TensorisedBoundaryAdapted<Dim, Degree, T, StoragePolicy>::evaluate( points_type const& __pts, mpl::int_<2> ) const
383 {
384  matrix_type res( convex_type::polyDims( nOrder ), __pts.size2() );
385 
386  vector_type eta1s = ublas::row( __pts, 0 );
387  vector_type eta2s = ublas::row( __pts, 1 );
388 
389  matrix_type psi_1( M_pfunc.evaluate_1( eta1s ) );
390  matrix_type psi_2( M_pfunc.evaluate_1( eta2s ) );
391 
392  /*
393 
394  C ------- D y
395  | | ^
396  | | |
397  | | |
398  | | |
399  | | |-----> x
400  A ------- B
401 
402  */
403 
404  /* 4 Vertex */
405 
406  /* A */
407  ublas::row( res,0 ) = ublas::element_prod( ublas::row( psi_1,0 ) , ublas::row( psi_2,0 ) );
408  /* B */
409  ublas::row( res,1 ) = ublas::element_prod( ublas::row( psi_1,nOrder ) , ublas::row( psi_2,0 ) );
410  /* C */
411  ublas::row( res,2 ) = ublas::element_prod( ublas::row( psi_1,0 ) , ublas::row( psi_2,nOrder ) );
412  /* D */
413  ublas::row( res,3 ) = ublas::element_prod( ublas::row( psi_1,nOrder ) , ublas::row( psi_2,nOrder ) );
414 
415  /* Global index */
416  uint_type G_i = 3;
417 
418  /* Edges */
420  /* AB = (nOrder-1) */
421  for ( uint_type p = 1; p < nOrder; ++p )
422  {
423  ++G_i;
424  ublas::row( res,G_i ) = ublas::element_prod( ublas::row( psi_1,p ) , ublas::row( psi_2,0 ) );
425  }
426 
427  /* AC = (nOrder-1)*/
428  for ( uint_type q = 1; q < nOrder; ++q )
429  {
430  ++G_i;
431  ublas::row( res,G_i ) = ublas::element_prod( ublas::row( psi_1,0 ) , ublas::row( psi_2,q ) );
432 
433  if ( q%2==0 )
434  ublas::row( res,G_i )*= value_type( -1.0 );
435  }
436 
437  /* CD = (nOrder-1) */
438  for ( uint_type p = 1; p < nOrder; ++p )
439  {
440  ++G_i;
441  ublas::row( res,G_i ) = ublas::element_prod( ublas::row( psi_1,p ) , ublas::row( psi_2,nOrder ) );
442 
443  if ( p%2==0 )
444  ublas::row( res,G_i )*= value_type( -1.0 );
445  }
446 
447  /* BD = (nOrder-1)*/
448  for ( uint_type q = 1; q < nOrder; ++q )
449  {
450  ++G_i;
451  ublas::row( res,G_i ) = ublas::element_prod( ublas::row( psi_1,nOrder ) , ublas::row( psi_2,q ) );
452  }
453 
455 
456  /* Interior : (N-1)*(N-1) */
458  for ( uint_type p = 1; p < nOrder; ++p )
459  for ( uint_type q = 1; q < nOrder; ++q )
460  {
461  ++G_i;
462  ublas::row( res,G_i ) = ublas::element_prod( ublas::row( psi_1,p ) , ublas::row( psi_2,q ) );
463  }
464 
466 
467  return res;
468 }
469 
470 
471 template<uint16_type Dim,
472  uint16_type Degree,
473  typename T,
474  template<class> class StoragePolicy>
475 template<typename AE>
476 typename TensorisedBoundaryAdapted<Dim, Degree, T, StoragePolicy>::vector_matrix_type
477 TensorisedBoundaryAdapted<Dim, Degree, T, StoragePolicy>::derivate( ublas::matrix_expression<AE> const& __pts, mpl::int_<2> ) const
478 {
479  vector_matrix_type res( 2 );
480  res[0].resize( convex_type::polyDims( nOrder ), __pts().size2() );
481  res[1].resize( convex_type::polyDims( nOrder ), __pts().size2() );
482 
483  vector_type eta1s = ublas::row( __pts(), 0 );
484  vector_type eta2s = ublas::row( __pts(), 1 );
485 
486  matrix_type psi_1( M_pfunc.evaluate_1( eta1s ) );
487  matrix_type psi_2( M_pfunc.evaluate_1( eta2s ) );
488 
489  matrix_type dpsi_1( M_pfunc.derivate_1( eta1s ) );
490  matrix_type dpsi_2( M_pfunc.derivate_1( eta2s ) );
491 
492  /* 4 Vertex */
493 
494  /* A */
495  ublas::row( res[0],0 ) = ublas::element_prod( ublas::row( dpsi_1,0 ) , ublas::row( psi_2,0 ) );
496  ublas::row( res[1],0 ) = ublas::element_prod( ublas::row( psi_1,0 ) , ublas::row( dpsi_2,0 ) );
497  /* B */
498  ublas::row( res[0],1 ) = ublas::element_prod( ublas::row( dpsi_1,nOrder ) , ublas::row( psi_2,0 ) );
499  ublas::row( res[1],1 ) = ublas::element_prod( ublas::row( psi_1,nOrder ) , ublas::row( dpsi_2,0 ) );
500  /* C */
501  ublas::row( res[0],2 ) = ublas::element_prod( ublas::row( dpsi_1,0 ) , ublas::row( psi_2,nOrder ) );
502  ublas::row( res[1],2 ) = ublas::element_prod( ublas::row( psi_1,0 ) , ublas::row( dpsi_2,nOrder ) );
503  /* D */
504  ublas::row( res[0],3 ) = ublas::element_prod( ublas::row( dpsi_1,nOrder ) , ublas::row( psi_2,nOrder ) );
505  ublas::row( res[1],3 ) = ublas::element_prod( ublas::row( psi_1,nOrder ) , ublas::row( dpsi_2,nOrder ) );
506 
507  /* Global index */
508  uint_type G_i = 3;
509 
510  /* Edges */
512  /* AB = (nOrder-1) */
513  for ( uint_type p = 1; p < nOrder; ++p )
514  {
515  ++G_i;
516  ublas::row( res[0],G_i ) = ublas::element_prod( ublas::row( dpsi_1,p ) , ublas::row( psi_2,0 ) );
517  ublas::row( res[1],G_i ) = ublas::element_prod( ublas::row( psi_1,p ) , ublas::row( dpsi_2,0 ) );
518  }
519 
520  /* AC = (nOrder-1)*/
521  for ( uint_type q = 1; q < nOrder; ++q )
522  {
523  ++G_i;
524  ublas::row( res[0],G_i ) = ublas::element_prod( ublas::row( dpsi_1,0 ) , ublas::row( psi_2,q ) );
525  ublas::row( res[1],G_i ) = ublas::element_prod( ublas::row( psi_1,0 ) , ublas::row( dpsi_2,q ) );
526  }
527 
528  /* CD = (nOrder-1) */
529  for ( uint_type p = 1; p < nOrder; ++p )
530  {
531  ++G_i;
532  ublas::row( res[0],G_i ) = ublas::element_prod( ublas::row( dpsi_1,p ) , ublas::row( psi_2,nOrder ) );
533  ublas::row( res[1],G_i ) = ublas::element_prod( ublas::row( psi_1,p ) , ublas::row( dpsi_2,nOrder ) );
534  }
535 
536  /* BD = (nOrder-1)*/
537  for ( uint_type q = 1; q < nOrder; ++q )
538  {
539  ++G_i;
540  ublas::row( res[0],G_i ) = ublas::element_prod( ublas::row( dpsi_1,nOrder ) , ublas::row( psi_2,q ) );
541  ublas::row( res[1],G_i ) = ublas::element_prod( ublas::row( psi_1,nOrder ) , ublas::row( dpsi_2,q ) );
542  }
543 
545 
546  /* Interior : (N-1)*(N-1) */
548  for ( uint_type p = 1; p < nOrder; ++p )
549  for ( uint_type q = 1; q < nOrder; ++q )
550  {
551  ++G_i;
552  ublas::row( res[0],G_i ) = ublas::element_prod( ublas::row( dpsi_1,p ) , ublas::row( psi_2,q ) );
553  ublas::row( res[1],G_i ) = ublas::element_prod( ublas::row( psi_1,p ) , ublas::row( dpsi_2,q ) );
554  }
555 
557 
558  return res;
559 }
560 
561 template<uint16_type Dim,
562  uint16_type Degree,
563  typename T,
564  template<class> class StoragePolicy>
565 typename TensorisedBoundaryAdapted<Dim, Degree, T, StoragePolicy>::matrix_type
566 TensorisedBoundaryAdapted<Dim, Degree, T, StoragePolicy>::evaluate( points_type const& __pts, mpl::int_<3> ) const
567 {
568  matrix_type res( convex_type::polyDims( nOrder ), __pts.size2() );
569 
570  ublas::vector<value_type> eta1s = ublas::row( __pts, 0 );
571  ublas::vector<value_type> eta2s = ublas::row( __pts, 1 );
572  ublas::vector<value_type> eta3s = ublas::row( __pts, 2 );
573 
574  matrix_type psi_1( M_pfunc.evaluate_1( eta1s ) );
575  matrix_type psi_2( M_pfunc.evaluate_1( eta2s ) );
576  matrix_type psi_3( M_pfunc.evaluate_1( eta3s ) );
577 
578  /*
579 
580  H---------G
581  /| /|
582  / | / |
583  / | / |
584  E---------F | z y
585  | | | | ^ ^
586  | D-----|---C | /
587  | / | / | /
588  | / | / |/
589  |/ |/ |-----> x
590  A---------B
591 
592  */
593 
594 
595  /* 8 Vertex */
596 
597  /* A */
598  ublas::row( res,0 ) =
599  ublas::element_prod( ublas::element_prod( ublas::row( psi_1,0 ), ublas::row( psi_2,0 ) ), ublas::row( psi_3,0 ) );
600  /* B */
601  ublas::row( res,1 ) =
602  ublas::element_prod( ublas::element_prod( ublas::row( psi_1,nOrder ), ublas::row( psi_2,0 ) ), ublas::row( psi_3,0 ) );
603  /* C */
604  ublas::row( res,2 ) =
605  ublas::element_prod( ublas::element_prod( ublas::row( psi_1,nOrder ), ublas::row( psi_2,nOrder ) ), ublas::row( psi_3,0 ) );
606  /* D */
607  ublas::row( res,3 ) =
608  ublas::element_prod( ublas::element_prod( ublas::row( psi_1,0 ), ublas::row( psi_2,nOrder ) ), ublas::row( psi_3,0 ) );
609  /* E */
610  ublas::row( res,4 ) =
611  ublas::element_prod( ublas::element_prod( ublas::row( psi_1,0 ), ublas::row( psi_2,0 ) ), ublas::row( psi_3,nOrder ) );
612  /* F */
613  ublas::row( res,5 ) =
614  ublas::element_prod( ublas::element_prod( ublas::row( psi_1,nOrder ), ublas::row( psi_2,0 ) ), ublas::row( psi_3,nOrder ) );
615  /* G */
616  ublas::row( res,6 ) =
617  ublas::element_prod( ublas::element_prod( ublas::row( psi_1,nOrder ), ublas::row( psi_2,nOrder ) ), ublas::row( psi_3,nOrder ) );
618  /* H */
619  ublas::row( res,7 ) =
620  ublas::element_prod( ublas::element_prod( ublas::row( psi_1,0 ), ublas::row( psi_2,nOrder ) ), ublas::row( psi_3,nOrder ) );
621 
622 
623  /* Global index */
624  uint_type G_i = 7;
625 
626  /* 12 Edges */
628  /* AB : (nOrder-1) */
629 
630  for ( uint_type p = 1; p < nOrder; ++p )
631  {
632  ++G_i;
633  ublas::row( res,G_i ) =
634  ublas::element_prod( ublas::element_prod( ublas::row( psi_1,p ), ublas::row( psi_2,0 ) ), ublas::row( psi_3,0 ) );
635  }
636 
637  /* BC : (nOrder-1) */
638 
639  for ( uint_type p = 1; p < nOrder; ++p )
640  {
641  ++G_i;
642  ublas::row( res,G_i ) =
643  ublas::element_prod( ublas::element_prod( ublas::row( psi_1,nOrder ), ublas::row( psi_2,p ) ), ublas::row( psi_3,0 ) );
644  }
645 
646  /* DC : (nOrder-1) */
647 
648  for ( uint_type p = 1; p < nOrder; ++p )
649  {
650  ++G_i;
651  ublas::row( res,G_i ) =
652  ublas::element_prod( ublas::element_prod( ublas::row( psi_1,p ), ublas::row( psi_2,nOrder ) ), ublas::row( psi_3,0 ) );
653  }
654 
655  /* AD : (nOrder-1) */
656 
657  for ( uint_type p = 1; p < nOrder; ++p )
658  {
659  ++G_i;
660  ublas::row( res,G_i ) =
661  ublas::element_prod( ublas::element_prod( ublas::row( psi_1,0 ), ublas::row( psi_2,p ) ), ublas::row( psi_3,0 ) );
662  }
663 
664  /* AE : (nOrder-1) */
665 
666  for ( uint_type p = 1; p < nOrder; ++p )
667  {
668  ++G_i;
669  ublas::row( res,G_i ) =
670  ublas::element_prod( ublas::element_prod( ublas::row( psi_1,0 ), ublas::row( psi_2,0 ) ), ublas::row( psi_3,p ) );
671  }
672 
673  /* BF : (nOrder-1) */
674 
675  for ( uint_type p = 1; p < nOrder; ++p )
676  {
677  ++G_i;
678  ublas::row( res,G_i ) =
679  ublas::element_prod( ublas::element_prod( ublas::row( psi_1,nOrder ), ublas::row( psi_2,0 ) ), ublas::row( psi_3,p ) );
680  }
681 
682  /* CG : (nOrder-1) */
683 
684  for ( uint_type p = 1; p < nOrder; ++p )
685  {
686  ++G_i;
687  ublas::row( res,G_i ) =
688  ublas::element_prod( ublas::element_prod( ublas::row( psi_1,nOrder ), ublas::row( psi_2,nOrder ) ), ublas::row( psi_3,p ) );
689  }
690 
691  /* DH : (nOrder-1) */
692 
693  for ( uint_type p = 1; p < nOrder; ++p )
694  {
695  ++G_i;
696  ublas::row( res,G_i ) =
697  ublas::element_prod( ublas::element_prod( ublas::row( psi_1,0 ), ublas::row( psi_2,nOrder ) ), ublas::row( psi_3,p ) );
698  }
699 
700  /* EF : (nOrder-1) */
701 
702  for ( uint_type p = 1; p < nOrder; ++p )
703  {
704  ++G_i;
705  ublas::row( res,G_i ) =
706  ublas::element_prod( ublas::element_prod( ublas::row( psi_1,p ), ublas::row( psi_2,0 ) ), ublas::row( psi_3,nOrder ) );
707  }
708 
709  /* FG : (nOrder-1) */
710 
711  for ( uint_type p = 1; p < nOrder; ++p )
712  {
713  ++G_i;
714  ublas::row( res,G_i ) =
715  ublas::element_prod( ublas::element_prod( ublas::row( psi_1,nOrder ), ublas::row( psi_2,p ) ), ublas::row( psi_3,nOrder ) );
716  }
717 
718  /* HG : (nOrder-1) */
719 
720  for ( uint_type p = 1; p < nOrder; ++p )
721  {
722  ++G_i;
723  ublas::row( res,G_i ) =
724  ublas::element_prod( ublas::element_prod( ublas::row( psi_1,p ), ublas::row( psi_2,nOrder ) ), ublas::row( psi_3,nOrder ) );
725  }
726 
727  /* EH : (nOrder-1) */
728 
729  for ( uint_type p = 1; p < nOrder; ++p )
730  {
731  ++G_i;
732  ublas::row( res,G_i ) =
733  ublas::element_prod( ublas::element_prod( ublas::row( psi_1,0 ), ublas::row( psi_2,p ) ), ublas::row( psi_3,nOrder ) );
734  }
735 
737 
738 
739  /* 6 Faces : (N-1)(N-1) modes */
741 
742  /* ABCD */
743 
744  for ( uint_type p = 1; p < nOrder; ++p )
745  for ( uint_type q = 1; q < nOrder; ++q )
746  {
747  ++G_i;
748  ublas::row( res,G_i ) =
749  ublas::element_prod( ublas::element_prod( ublas::row( psi_1,p ) , ublas::row( psi_2,q ) ), ublas::row( psi_3,0 ) );
750  }
751 
752  /* ABFE */
753 
754  for ( uint_type p = 1; p < nOrder; ++p )
755  for ( uint_type q = 1; q < nOrder; ++q )
756  {
757  ++G_i;
758  ublas::row( res,G_i ) =
759  ublas::element_prod( ublas::element_prod( ublas::row( psi_1,p ) , ublas::row( psi_2,0 ) ), ublas::row( psi_3,q ) );
760  }
761 
762  /* BCGF */
763 
764  for ( uint_type p = 1; p < nOrder; ++p )
765  for ( uint_type q = 1; q < nOrder; ++q )
766  {
767  ++G_i;
768  ublas::row( res,G_i ) =
769  ublas::element_prod( ublas::element_prod( ublas::row( psi_1,nOrder ) , ublas::row( psi_2,p ) ), ublas::row( psi_3,q ) );
770  }
771 
772  /* CGHD */
773 
774  for ( uint_type p = 1; p < nOrder; ++p )
775  for ( uint_type q = 1; q < nOrder; ++q )
776  {
777  ++G_i;
778  ublas::row( res,G_i ) =
779  ublas::element_prod( ublas::element_prod( ublas::row( psi_1,p ) , ublas::row( psi_2,nOrder ) ), ublas::row( psi_3,q ) );
780  }
781 
782  /* ADHE */
783 
784  for ( uint_type p = 1; p < nOrder; ++p )
785  for ( uint_type q = 1; q < nOrder; ++q )
786  {
787  ++G_i;
788  ublas::row( res,G_i ) =
789  ublas::element_prod( ublas::element_prod( ublas::row( psi_1,0 ) , ublas::row( psi_2,p ) ), ublas::row( psi_3,q ) );
790  }
791 
792  /* EFGH */
793 
794  for ( uint_type p = 1; p < nOrder; ++p )
795  for ( uint_type q = 1; q < nOrder; ++q )
796  {
797  ++G_i;
798  ublas::row( res,G_i ) =
799  ublas::element_prod( ublas::element_prod( ublas::row( psi_1,p ) , ublas::row( psi_2,q ) ), ublas::row( psi_3,nOrder ) );
800  }
801 
803 
804 
805  /* Interior : (N-1)^3 */
806 
807  for ( uint_type p = 1; p < nOrder; ++p )
808  for ( uint_type q = 1; q < nOrder; ++q )
809  for ( uint_type r = 1; r < nOrder; ++r )
810  {
811  ++G_i;
812  ublas::row( res,G_i ) =
813  ublas::element_prod( ublas::element_prod( ublas::row( psi_1,p ) , ublas::row( psi_2,q ) ), ublas::row( psi_3,r ) );
814  }
815 
816  return res;
817 }
818 
819 template<uint16_type Dim,
820  uint16_type Degree,
821  typename T,
822  template<class> class StoragePolicy>
823 template<typename AE>
824 typename TensorisedBoundaryAdapted<Dim, Degree, T, StoragePolicy>::vector_matrix_type
825 TensorisedBoundaryAdapted<Dim, Degree, T, StoragePolicy>::derivate( ublas::matrix_expression<AE> const& __pts, mpl::int_<3> ) const
826 {
827  vector_matrix_type res( 3 );
828  res[0].resize( convex_type::polyDims( nOrder ), __pts().size2() );
829  res[1].resize( convex_type::polyDims( nOrder ), __pts().size2() );
830  res[2].resize( convex_type::polyDims( nOrder ), __pts().size2() );
831 
832 
833  vector_type eta1s = ublas::row( __pts(), 0 );
834  vector_type eta2s = ublas::row( __pts(), 1 );
835  vector_type eta3s = ublas::row( __pts(), 2 );
836 
837  matrix_type psi_1( M_pfunc.evaluate_1( eta1s ) );
838  matrix_type dpsi_1( M_pfunc.derivate_1( eta1s ) );
839 
840  matrix_type psi_2( M_pfunc.evaluate_1( eta2s ) );
841  matrix_type dpsi_2( M_pfunc.derivate_1( eta2s ) );
842 
843  matrix_type psi_3( M_pfunc.evaluate_1( eta3s ) );
844  matrix_type dpsi_3( M_pfunc.derivate_1( eta3s ) );
845 
846  /* 8 Vertex */
847 
848  ublas::row( res[0],0 ) =
849  ublas::element_prod( ublas::element_prod( ublas::row( dpsi_1,0 ), ublas::row( psi_2,0 ) ), ublas::row( psi_3,0 ) );
850  ublas::row( res[1],0 ) =
851  ublas::element_prod( ublas::element_prod( ublas::row( psi_1,0 ), ublas::row( dpsi_2,0 ) ), ublas::row( psi_3,0 ) );
852  ublas::row( res[2],0 ) =
853  ublas::element_prod( ublas::element_prod( ublas::row( psi_1,0 ), ublas::row( psi_2,0 ) ), ublas::row( dpsi_3,0 ) );
854 
855  ublas::row( res[0],1 ) =
856  ublas::element_prod( ublas::element_prod( ublas::row( dpsi_1,nOrder ), ublas::row( psi_2,0 ) ), ublas::row( psi_3,0 ) );
857  ublas::row( res[1],1 ) =
858  ublas::element_prod( ublas::element_prod( ublas::row( psi_1,nOrder ), ublas::row( dpsi_2,0 ) ), ublas::row( psi_3,0 ) );
859  ublas::row( res[2],1 ) =
860  ublas::element_prod( ublas::element_prod( ublas::row( psi_1,nOrder ), ublas::row( psi_2,0 ) ), ublas::row( dpsi_3,0 ) );
861 
862  ublas::row( res[0],2 ) =
863  ublas::element_prod( ublas::element_prod( ublas::row( dpsi_1,nOrder ), ublas::row( psi_2,nOrder ) ), ublas::row( psi_3,0 ) );
864  ublas::row( res[1],2 ) =
865  ublas::element_prod( ublas::element_prod( ublas::row( psi_1,nOrder ), ublas::row( dpsi_2,nOrder ) ), ublas::row( psi_3,0 ) );
866  ublas::row( res[2],2 ) =
867  ublas::element_prod( ublas::element_prod( ublas::row( psi_1,nOrder ), ublas::row( psi_2,nOrder ) ), ublas::row( dpsi_3,0 ) );
868 
869  ublas::row( res[0],3 ) =
870  ublas::element_prod( ublas::element_prod( ublas::row( dpsi_1,0 ), ublas::row( psi_2,nOrder ) ), ublas::row( psi_3,0 ) );
871  ublas::row( res[1],3 ) =
872  ublas::element_prod( ublas::element_prod( ublas::row( psi_1,0 ), ublas::row( dpsi_2,nOrder ) ), ublas::row( psi_3,0 ) );
873  ublas::row( res[2],3 ) =
874  ublas::element_prod( ublas::element_prod( ublas::row( psi_1,0 ), ublas::row( psi_2,nOrder ) ), ublas::row( dpsi_3,0 ) );
875 
876  ublas::row( res[0],4 ) =
877  ublas::element_prod( ublas::element_prod( ublas::row( dpsi_1,0 ), ublas::row( psi_2,0 ) ), ublas::row( psi_3,nOrder ) );
878  ublas::row( res[1],4 ) =
879  ublas::element_prod( ublas::element_prod( ublas::row( psi_1,0 ), ublas::row( dpsi_2,0 ) ), ublas::row( psi_3,nOrder ) );
880  ublas::row( res[2],4 ) =
881  ublas::element_prod( ublas::element_prod( ublas::row( psi_1,0 ), ublas::row( psi_2,0 ) ), ublas::row( dpsi_3,nOrder ) );
882 
883  ublas::row( res[0],5 ) =
884  ublas::element_prod( ublas::element_prod( ublas::row( dpsi_1,nOrder ), ublas::row( psi_2,0 ) ), ublas::row( psi_3,nOrder ) );
885  ublas::row( res[1],5 ) =
886  ublas::element_prod( ublas::element_prod( ublas::row( psi_1,nOrder ), ublas::row( dpsi_2,0 ) ), ublas::row( psi_3,nOrder ) );
887  ublas::row( res[2],5 ) =
888  ublas::element_prod( ublas::element_prod( ublas::row( psi_1,nOrder ), ublas::row( psi_2,0 ) ), ublas::row( dpsi_3,nOrder ) );
889 
890  ublas::row( res[0],6 ) =
891  ublas::element_prod( ublas::element_prod( ublas::row( dpsi_1,nOrder ), ublas::row( psi_2,nOrder ) ), ublas::row( psi_3,nOrder ) );
892  ublas::row( res[1],6 ) =
893  ublas::element_prod( ublas::element_prod( ublas::row( psi_1,nOrder ), ublas::row( dpsi_2,nOrder ) ), ublas::row( psi_3,nOrder ) );
894  ublas::row( res[2],6 ) =
895  ublas::element_prod( ublas::element_prod( ublas::row( psi_1,nOrder ), ublas::row( psi_2,nOrder ) ), ublas::row( dpsi_3,nOrder ) );
896 
897  ublas::row( res[0],7 ) =
898  ublas::element_prod( ublas::element_prod( ublas::row( dpsi_1,0 ), ublas::row( psi_2,nOrder ) ), ublas::row( psi_3,nOrder ) );
899  ublas::row( res[1],7 ) =
900  ublas::element_prod( ublas::element_prod( ublas::row( psi_1,0 ), ublas::row( dpsi_2,nOrder ) ), ublas::row( psi_3,nOrder ) );
901  ublas::row( res[2],7 ) =
902  ublas::element_prod( ublas::element_prod( ublas::row( psi_1,0 ), ublas::row( psi_2,nOrder ) ), ublas::row( dpsi_3,nOrder ) );
903 
904 
905  /* Global index */
906  uint_type G_i = 7;
907 
908  /* 12 Edges */
910  /* AB : (nOrder-1) */
911 
912  for ( uint_type p = 1; p < nOrder; ++p )
913  {
914  ++G_i;
915  ublas::row( res[0],G_i ) =
916  ublas::element_prod( ublas::element_prod( ublas::row( dpsi_1,p ), ublas::row( psi_2,0 ) ), ublas::row( psi_3,0 ) );
917  ublas::row( res[1],G_i ) =
918  ublas::element_prod( ublas::element_prod( ublas::row( psi_1,p ), ublas::row( dpsi_2,0 ) ), ublas::row( psi_3,0 ) );
919  ublas::row( res[2],G_i ) =
920  ublas::element_prod( ublas::element_prod( ublas::row( psi_1,p ), ublas::row( psi_2,0 ) ), ublas::row( dpsi_3,0 ) );
921  }
922 
923  /* BC : (nOrder-1) */
924 
925  for ( uint_type p = 1; p < nOrder; ++p )
926  {
927  ++G_i;
928  ublas::row( res[0],G_i ) =
929  ublas::element_prod( ublas::element_prod( ublas::row( dpsi_1,nOrder ), ublas::row( psi_2,p ) ), ublas::row( psi_3,0 ) );
930  ublas::row( res[1],G_i ) =
931  ublas::element_prod( ublas::element_prod( ublas::row( psi_1,nOrder ), ublas::row( dpsi_2,p ) ), ublas::row( psi_3,0 ) );
932  ublas::row( res[2],G_i ) =
933  ublas::element_prod( ublas::element_prod( ublas::row( psi_1,nOrder ), ublas::row( psi_2,p ) ), ublas::row( dpsi_3,0 ) );
934  }
935 
936  /* DC : (nOrder-1) */
937 
938  for ( uint_type p = 1; p < nOrder; ++p )
939  {
940  ++G_i;
941  ublas::row( res[0],G_i ) =
942  ublas::element_prod( ublas::element_prod( ublas::row( dpsi_1,p ), ublas::row( psi_2,nOrder ) ), ublas::row( psi_3,0 ) );
943  ublas::row( res[1],G_i ) =
944  ublas::element_prod( ublas::element_prod( ublas::row( psi_1,p ), ublas::row( dpsi_2,nOrder ) ), ublas::row( psi_3,0 ) );
945  ublas::row( res[2],G_i ) =
946  ublas::element_prod( ublas::element_prod( ublas::row( psi_1,p ), ublas::row( psi_2,nOrder ) ), ublas::row( dpsi_3,0 ) );
947  }
948 
949  /* AD : (nOrder-1) */
950 
951  for ( uint_type p = 1; p < nOrder; ++p )
952  {
953  ++G_i;
954  ublas::row( res[0],G_i ) =
955  ublas::element_prod( ublas::element_prod( ublas::row( dpsi_1,0 ), ublas::row( psi_2,p ) ), ublas::row( psi_3,0 ) );
956  ublas::row( res[1],G_i ) =
957  ublas::element_prod( ublas::element_prod( ublas::row( psi_1,0 ), ublas::row( dpsi_2,p ) ), ublas::row( psi_3,0 ) );
958  ublas::row( res[2],G_i ) =
959  ublas::element_prod( ublas::element_prod( ublas::row( psi_1,0 ), ublas::row( psi_2,p ) ), ublas::row( dpsi_3,0 ) );
960  }
961 
962  /* AE : (nOrder-1) */
963 
964  for ( uint_type p = 1; p < nOrder; ++p )
965  {
966  ++G_i;
967  ublas::row( res[0],G_i ) =
968  ublas::element_prod( ublas::element_prod( ublas::row( dpsi_1,0 ), ublas::row( psi_2,0 ) ), ublas::row( psi_3,p ) );
969  ublas::row( res[1],G_i ) =
970  ublas::element_prod( ublas::element_prod( ublas::row( psi_1,0 ), ublas::row( dpsi_2,0 ) ), ublas::row( psi_3,p ) );
971  ublas::row( res[2],G_i ) =
972  ublas::element_prod( ublas::element_prod( ublas::row( psi_1,0 ), ublas::row( psi_2,0 ) ), ublas::row( dpsi_3,p ) );
973  }
974 
975  /* BF : (nOrder-1) */
976 
977  for ( uint_type p = 1; p < nOrder; ++p )
978  {
979  ++G_i;
980  ublas::row( res[0],G_i ) =
981  ublas::element_prod( ublas::element_prod( ublas::row( dpsi_1,nOrder ), ublas::row( psi_2,0 ) ), ublas::row( psi_3,p ) );
982  ublas::row( res[1],G_i ) =
983  ublas::element_prod( ublas::element_prod( ublas::row( psi_1,nOrder ), ublas::row( dpsi_2,0 ) ), ublas::row( psi_3,p ) );
984  ublas::row( res[2],G_i ) =
985  ublas::element_prod( ublas::element_prod( ublas::row( psi_1,nOrder ), ublas::row( psi_2,0 ) ), ublas::row( dpsi_3,p ) );
986  }
987 
988  /* CG : (nOrder-1) */
989 
990  for ( uint_type p = 1; p < nOrder; ++p )
991  {
992  ++G_i;
993  ublas::row( res[0],G_i ) =
994  ublas::element_prod( ublas::element_prod( ublas::row( dpsi_1,nOrder ), ublas::row( psi_2,nOrder ) ), ublas::row( psi_3,p ) );
995  ublas::row( res[1],G_i ) =
996  ublas::element_prod( ublas::element_prod( ublas::row( psi_1,nOrder ), ublas::row( dpsi_2,nOrder ) ), ublas::row( psi_3,p ) );
997  ublas::row( res[2],G_i ) =
998  ublas::element_prod( ublas::element_prod( ublas::row( psi_1,nOrder ), ublas::row( psi_2,nOrder ) ), ublas::row( dpsi_3,p ) );
999  }
1000 
1001  /* DH : (nOrder-1) */
1002 
1003  for ( uint_type p = 1; p < nOrder; ++p )
1004  {
1005  ++G_i;
1006  ublas::row( res[0],G_i ) =
1007  ublas::element_prod( ublas::element_prod( ublas::row( dpsi_1,0 ), ublas::row( psi_2,nOrder ) ), ublas::row( psi_3,p ) );
1008  ublas::row( res[1],G_i ) =
1009  ublas::element_prod( ublas::element_prod( ublas::row( psi_1,0 ), ublas::row( dpsi_2,nOrder ) ), ublas::row( psi_3,p ) );
1010  ublas::row( res[2],G_i ) =
1011  ublas::element_prod( ublas::element_prod( ublas::row( psi_1,0 ), ublas::row( psi_2,nOrder ) ), ublas::row( dpsi_3,p ) );
1012  }
1013 
1014  /* EF : (nOrder-1) */
1015 
1016  for ( uint_type p = 1; p < nOrder; ++p )
1017  {
1018  ++G_i;
1019  ublas::row( res[0],G_i ) =
1020  ublas::element_prod( ublas::element_prod( ublas::row( dpsi_1,p ), ublas::row( psi_2,0 ) ), ublas::row( psi_3,nOrder ) );
1021  ublas::row( res[1],G_i ) =
1022  ublas::element_prod( ublas::element_prod( ublas::row( psi_1,p ), ublas::row( dpsi_2,0 ) ), ublas::row( psi_3,nOrder ) );
1023  ublas::row( res[2],G_i ) =
1024  ublas::element_prod( ublas::element_prod( ublas::row( psi_1,p ), ublas::row( psi_2,0 ) ), ublas::row( dpsi_3,nOrder ) );
1025  }
1026 
1027  /* FG : (nOrder-1) */
1028 
1029  for ( uint_type p = 1; p < nOrder; ++p )
1030  {
1031  ++G_i;
1032  ublas::row( res[0],G_i ) =
1033  ublas::element_prod( ublas::element_prod( ublas::row( dpsi_1,nOrder ), ublas::row( psi_2,p ) ), ublas::row( psi_3,nOrder ) );
1034  ublas::row( res[1],G_i ) =
1035  ublas::element_prod( ublas::element_prod( ublas::row( psi_1,nOrder ), ublas::row( dpsi_2,p ) ), ublas::row( psi_3,nOrder ) );
1036  ublas::row( res[2],G_i ) =
1037  ublas::element_prod( ublas::element_prod( ublas::row( psi_1,nOrder ), ublas::row( psi_2,p ) ), ublas::row( dpsi_3,nOrder ) );
1038  }
1039 
1040  /* HG : (nOrder-1) */
1041 
1042  for ( uint_type p = 1; p < nOrder; ++p )
1043  {
1044  ++G_i;
1045  ublas::row( res[0],G_i ) =
1046  ublas::element_prod( ublas::element_prod( ublas::row( dpsi_1,p ), ublas::row( psi_2,nOrder ) ), ublas::row( psi_3,nOrder ) );
1047  ublas::row( res[1],G_i ) =
1048  ublas::element_prod( ublas::element_prod( ublas::row( psi_1,p ), ublas::row( dpsi_2,nOrder ) ), ublas::row( psi_3,nOrder ) );
1049  ublas::row( res[2],G_i ) =
1050  ublas::element_prod( ublas::element_prod( ublas::row( psi_1,p ), ublas::row( psi_2,nOrder ) ), ublas::row( dpsi_3,nOrder ) );
1051  }
1052 
1053  /* EH : (nOrder-1) */
1054 
1055  for ( uint_type p = 1; p < nOrder; ++p )
1056  {
1057  ++G_i;
1058  ublas::row( res[0],G_i ) =
1059  ublas::element_prod( ublas::element_prod( ublas::row( dpsi_1,0 ), ublas::row( psi_2,p ) ), ublas::row( psi_3,nOrder ) );
1060  ublas::row( res[1],G_i ) =
1061  ublas::element_prod( ublas::element_prod( ublas::row( psi_1,0 ), ublas::row( dpsi_2,p ) ), ublas::row( psi_3,nOrder ) );
1062  ublas::row( res[2],G_i ) =
1063  ublas::element_prod( ublas::element_prod( ublas::row( psi_1,0 ), ublas::row( psi_2,p ) ), ublas::row( dpsi_3,nOrder ) );
1064  }
1065 
1067 
1068 
1069  /* 6 Faces : (N-1)(N-1) modes */
1071 
1072  /* ABCD */
1073 
1074  for ( uint_type p = 1; p < nOrder; ++p )
1075  for ( uint_type q = 1; q < nOrder; ++q )
1076  {
1077  ++G_i;
1078  ublas::row( res[0],G_i ) =
1079  ublas::element_prod( ublas::element_prod( ublas::row( dpsi_1,p ) , ublas::row( psi_2,q ) ), ublas::row( psi_3,0 ) );
1080  ublas::row( res[1],G_i ) =
1081  ublas::element_prod( ublas::element_prod( ublas::row( psi_1,p ) , ublas::row( dpsi_2,q ) ), ublas::row( psi_3,0 ) );
1082  ublas::row( res[2],G_i ) =
1083  ublas::element_prod( ublas::element_prod( ublas::row( psi_1,p ) , ublas::row( psi_2,q ) ), ublas::row( dpsi_3,0 ) );
1084  }
1085 
1086  /* ABFE */
1087 
1088  for ( uint_type p = 1; p < nOrder; ++p )
1089  for ( uint_type q = 1; q < nOrder; ++q )
1090  {
1091  ++G_i;
1092  ublas::row( res[0],G_i ) =
1093  ublas::element_prod( ublas::element_prod( ublas::row( dpsi_1,p ) , ublas::row( psi_2,0 ) ), ublas::row( psi_3,q ) );
1094  ublas::row( res[1],G_i ) =
1095  ublas::element_prod( ublas::element_prod( ublas::row( psi_1,p ) , ublas::row( dpsi_2,0 ) ), ublas::row( psi_3,q ) );
1096  ublas::row( res[2],G_i ) =
1097  ublas::element_prod( ublas::element_prod( ublas::row( psi_1,p ) , ublas::row( psi_2,0 ) ), ublas::row( dpsi_3,q ) );
1098  }
1099 
1100  /* BCGF */
1101 
1102  for ( uint_type p = 1; p < nOrder; ++p )
1103  for ( uint_type q = 1; q < nOrder; ++q )
1104  {
1105  ++G_i;
1106  ublas::row( res[0],G_i ) =
1107  ublas::element_prod( ublas::element_prod( ublas::row( dpsi_1,nOrder ) , ublas::row( psi_2,p ) ), ublas::row( psi_3,q ) );
1108  ublas::row( res[1],G_i ) =
1109  ublas::element_prod( ublas::element_prod( ublas::row( psi_1,nOrder ) , ublas::row( dpsi_2,p ) ), ublas::row( psi_3,q ) );
1110  ublas::row( res[2],G_i ) =
1111  ublas::element_prod( ublas::element_prod( ublas::row( psi_1,nOrder ) , ublas::row( psi_2,p ) ), ublas::row( dpsi_3,q ) );
1112  }
1113 
1114  /* CGHD */
1115 
1116  for ( uint_type p = 1; p < nOrder; ++p )
1117  for ( uint_type q = 1; q < nOrder; ++q )
1118  {
1119  ++G_i;
1120  ublas::row( res[0],G_i ) =
1121  ublas::element_prod( ublas::element_prod( ublas::row( dpsi_1,p ) , ublas::row( psi_2,nOrder ) ), ublas::row( psi_3,q ) );
1122  ublas::row( res[1],G_i ) =
1123  ublas::element_prod( ublas::element_prod( ublas::row( psi_1,p ) , ublas::row( dpsi_2,nOrder ) ), ublas::row( psi_3,q ) );
1124  ublas::row( res[2],G_i ) =
1125  ublas::element_prod( ublas::element_prod( ublas::row( psi_1,p ) , ublas::row( psi_2,nOrder ) ), ublas::row( dpsi_3,q ) );
1126  }
1127 
1128  /* ADHE */
1129 
1130  for ( uint_type p = 1; p < nOrder; ++p )
1131  for ( uint_type q = 1; q < nOrder; ++q )
1132  {
1133  ++G_i;
1134  ublas::row( res[0],G_i ) =
1135  ublas::element_prod( ublas::element_prod( ublas::row( dpsi_1,0 ) , ublas::row( psi_2,p ) ), ublas::row( psi_3,q ) );
1136  ublas::row( res[1],G_i ) =
1137  ublas::element_prod( ublas::element_prod( ublas::row( psi_1,0 ) , ublas::row( dpsi_2,p ) ), ublas::row( psi_3,q ) );
1138  ublas::row( res[2],G_i ) =
1139  ublas::element_prod( ublas::element_prod( ublas::row( psi_1,0 ) , ublas::row( psi_2,p ) ), ublas::row( dpsi_3,q ) );
1140  }
1141 
1142  /* EFGH */
1143 
1144  for ( uint_type p = 1; p < nOrder; ++p )
1145  for ( uint_type q = 1; q < nOrder; ++q )
1146  {
1147  ++G_i;
1148  ublas::row( res[0],G_i ) =
1149  ublas::element_prod( ublas::element_prod( ublas::row( dpsi_1,p ) , ublas::row( psi_2,q ) ), ublas::row( psi_3,nOrder ) );
1150  ublas::row( res[1],G_i ) =
1151  ublas::element_prod( ublas::element_prod( ublas::row( psi_1,p ) , ublas::row( dpsi_2,q ) ), ublas::row( psi_3,nOrder ) );
1152  ublas::row( res[2],G_i ) =
1153  ublas::element_prod( ublas::element_prod( ublas::row( psi_1,p ) , ublas::row( psi_2,q ) ), ublas::row( dpsi_3,nOrder ) );
1154  }
1155 
1157 
1158 
1159  /* Interior : (N-1)^3 */
1160 
1161  for ( uint_type p = 1; p < nOrder; ++p )
1162  for ( uint_type q = 1; q < nOrder; ++q )
1163  for ( uint_type r = 1; r < nOrder; ++r )
1164  {
1165  ++G_i;
1166  ublas::row( res[0],G_i ) =
1167  ublas::element_prod( ublas::element_prod( ublas::row( dpsi_1,p ) , ublas::row( psi_2,q ) ), ublas::row( psi_3,r ) );
1168  ublas::row( res[1],G_i ) =
1169  ublas::element_prod( ublas::element_prod( ublas::row( psi_1,p ) , ublas::row( dpsi_2,q ) ), ublas::row( psi_3,r ) );
1170  ublas::row( res[2],G_i ) =
1171  ublas::element_prod( ublas::element_prod( ublas::row( psi_1,p ) , ublas::row( psi_2,q ) ), ublas::row( dpsi_3,r ) );
1172  }
1173 
1174 
1175  return res;
1176 }
1177 
1178 
1179 
1180 
1181 }
1182 #endif /* __BoundAdapted_H */

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