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
functionals2.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-22
7 
8  Copyright (C) 2006 EPFL
9  Copyright (C) 2009 Université de Grenoble 1 (Joseph Fourier)
10 
11  This library is free software; you can redistribute it and/or
12  modify it under the terms of the GNU Lesser General Public
13  License as published by the Free Software Foundation; either
14  version 3.0 of the License, or (at your option) any later version.
15 
16  This library is distributed in the hope that it will be useful,
17  but WITHOUT ANY WARRANTY; without even the implied warranty of
18  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19  Lesser General Public License for more details.
20 
21  You should have received a copy of the GNU Lesser General Public
22  License along with this library; if not, write to the Free Software
23  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
24 */
30 #ifndef __FEELPP_FUNCTIONALS2_HPP
31 #define __FEELPP_FUNCTIONALS2_HPP 1
32 
34 
35 #include <feel/feelpoly/quadpoint.hpp>
36 #include <feel/feelpoly/im.hpp>
37 
38 namespace Feel
39 {
40 
41 namespace functional
42 {
53 template<typename Space, typename Poly = Space>
55  :
56 public Functional<Space>
57 {
58  typedef Functional<Space> super;
59 public:
60 
62  typedef typename super::space_type space_type;
63  typedef typename Poly::polynomial_type polynomial_type;
64  typedef typename space_type::value_type value_type;
65  typedef typename node<value_type>::type node_type;
66 
67  //BOOST_STATIC_ASSERT( ( boost::is_same<typename space_type::basis_type, typename Poly::basis_type>::value ) );
68 
70  :
71  super()
72  {}
73 
74  IntegralMoment( space_type const& p, polynomial_type const& q )
75  :
76  super( p )
77  {
78 #if 0
79  std::cout << "p.coeff = " << p.coeff() << "\n";
80  std::cout << "p.basis.coeff = " << p.basis().coeff() << "\n";
81  std::cout << "q.coeff = " << q.coeff() << "\n";
82 #endif
83  //ublas::matrix<value_type> m ( ublas::prod( q.coeff(), ublas::trans( p.basis().coeff() ) ) );
84  typename space_type::Pkp1_v_type l;
85  ublas::matrix<value_type> m ( ublas::prod( q.coeff(), ublas::trans( l.coeff() ) ) );
86 
87  //ublas::matrix<value_type> m ( ublas::prod( q.coeff(), p.coeff() ) );
88  //std::cout << "[IntegralMoment] m = " << m << "\n";
89  this->setCoefficient( m );
90  }
91 
92  IntegralMoment( IntegralMoment const& im )
93  :
94  super( im )
95  {}
96 
98  {
99  if ( this != &im )
100  {
101  super::operator=( im );
102  }
103 
104  return *this;
105  }
106 
107 };
108 
109 namespace detail
110 {
111 template<typename P1, typename P2>
112 struct prod
113 {
114  typedef typename P1::value_type value_type;
115  prod ( P1 const& p1, P2 const& p2 )
116  :
117  M_p1 ( p1 ),
118  M_p2 ( p2 )
119  {}
120 #if 0
121  value_type operator() ( typename node<value_type>::type const& n )
122  {
123  return M_p1.evaluate( n )( 0,0 ) * M_p2.evaluate( n )( 0,0 );
124  }
125 #endif
126  typename node<value_type>::type operator() ( typename node<value_type>::type const& n )
127  {
128  return ublas::column( ublas::element_prod( M_p1.evaluate( n ), M_p2.evaluate( n ) ), 0 );
129  }
130  P1 M_p1;
131  P2 M_p2;
132 };
133 } // detail
141 template<typename Space>
143  :
144 public Functional<Space>
145 {
146  typedef Functional<Space> super;
147 public:
148 
150  typedef typename super::space_type space_type;
151  typedef typename super::polynomial_type polynomial_type;
152  typedef typename space_type::value_type value_type;
153  typedef typename node<value_type>::type node_type;
154 
155  const static uint16_type nComponents = space_type::nComponents;
156 
164  IntegralMomentOnFace( space_type const& p, uint16_type k, IntegrationFaceEnum face = ALL_FACES )
165  :
166  super( p ),
167  M_k ( k ),
168  M_q()
169  {
170  ublas::matrix<value_type> __rep( nComponents, p.polynomialDimensionPerComponent() );
171  typedef detail::prod<typename space_type::polynomial_type, polynomial_type> prod_fun;
172 
173  for ( uint16_type i = 0; i < p.polynomialDimensionPerComponent(); ++i )
174  {
175  //if ( p.is_scalar )
176  //__rep( 0, i ) = M_q.integrate( face, prod_fun( p.polynomial(i), p.polynomial( k ) ) );
177  //else
178  typedef typename node<value_type>::type node_type;
179  ublas::column( __rep, i ) = M_q.integrate( face, prod_fun( p.polynomial( i ), p.polynomial( k ) ) );
180  }
181 
182  this->setCoefficient( __rep );
183  }
184 
193  IntegralMomentOnFace( space_type const& p, uint16_type k, uint16_type c, IntegrationFaceEnum face = ALL_FACES )
194  :
195  super( p ),
196  M_k ( k ),
197  M_q()
198  {
199  ublas::matrix<value_type> __rep( ublas::zero_matrix<value_type>( nComponents, p.polynomialDimensionPerComponent() ) );
200  typedef detail::prod<typename space_type::polynomial_type::component_type,
201  typename polynomial_type::component_type> prod_fun;
202 
203  int nc = p.polynomialDimensionPerComponent()*c;
204  int ind_p2 = nc + k;
205 
206  for ( uint16_type i = 0; i < p.polynomialDimensionPerComponent(); ++i )
207  {
208  int ind_p1 = nc + i;
209  typedef typename node<value_type>::type node_type;
210  __rep( c, i ) = M_q.integrate( face, prod_fun( p.polynomial( ind_p1 )[c],
211  p.polynomial( ind_p2 )[c] ) )( 0 );
212  }
213 
214  this->setCoefficient( __rep );
215  }
216 
217 private:
218 
219  // disabled
221 
222 private:
223 
224  // polynomial degree integrate against
225  uint16_type M_k;
226 
227  // quadrature rule on the element and faces of the element
228  IM<Space::nDim,2*Space::nOrder+1, value_type> M_q;
229 };
230 
231 
239 template<typename Space,typename BasisType>
241  :
242 public std::vector<Functional<Space> >
243 {
244  typedef std::vector<Functional<Space> > super;
245 public:
246 
249  typedef BasisType basis_type;
250  typedef Space space_type;
251  typedef typename space_type::reference_convex_type reference_convex_type;
252  typedef typename super::polynomial_type polynomial_type;
253  typedef typename space_type::value_type value_type;
254  typedef typename node<value_type>::type node_type;
255 
256  const static uint16_type nComponents = space_type::nComponents;
257 
265  IntegralMomentsOnFace( space_type const& p,
266  basis_type const& l,
267  IntegrationFaceEnum face = ALL_FACES )
268  :
269  super()
270  {
271  reference_convex_type ref_convex;
272  typedef typename reference_convex_type::topological_face_type element_type;
273  element_type ref_convex_face = ref_convex.topologicalFace( face );
274 
276  typedef typename gm_type::face_gm_type::precompute_type face_pc_type;
277  typedef typename gm_type::face_gm_type::precompute_ptrtype face_pc_ptrtype;
278  gm_type __gm;
279  IM<reference_convex_type::nDim-1,2*space_type::nOrder-1> __qr_face;
280  face_pc_ptrtype __geopc( new face_pc_type( __gm->boundaryMap(),__qr_face.points() ) );
281 
282  DVLOG(2) << "[nc] ref_convex_face " << face << "=" << ref_convex_face.points() << "\n";
283 
284 
285  typename gm_type::template Context<vm::POINT,element_type> __c( __gm->boundaryMap(),
286  ref_convex_face,
287  __geopc );
288 
289  __c.update( ref_convex_face, __geopc );
290  DVLOG(2) << "[nc] ref_convex_face " << face << " xref" << __c.xRefs() << "\n";
291  DVLOG(2) << "[nc] ref_convex_face " << face << " xreal" << __c.xReal() << "\n";
292 
293  for ( uint16_type k = 0; k < l.polynomialDimensionPerComponent(); ++k )
294  {
295  ublas::matrix<value_type> __rep( nComponents, p.polynomialDimensionPerComponent() );
296 
297  for ( uint16_type i = 0; i < p.polynomialDimensionPerComponent(); ++i )
298  {
299  /*
300  typedef typename node<value_type>::type node_type;
301  double __res = 0;
302  double __len = 0;
303  for ( uint16_type __ip = 0; __ip < __qr_face.nPoints();++__ip )
304  {
305  __res += ( __qr_face.weight( __ip )*
306  p.polynomial(i).evaluate( __c.xReal(__ip) )*
307  l.polynomial(k).evaluate( __c.xRef( __ip ) ));
308 
309  __len += __qr_face.weight( __ip );
310  }
311  DVLOG(2) << "[nc] length = " << __len << "\n";
312  DVLOG(2) << "[nc] res = " << __res << "\n";
313  ublas::column( __rep, i ) = __res/__len;
314  }
315  this->push_back( functional_type( p, __rep ) );
316  */
317  }
318 
319  }
320  }
321 
322 private:
323 
324  // disabled
326 
327 private:
328 
329 
330 };
331 
338 template<typename Space>
340  :
341 public Functional<Space>
342 {
343  typedef Functional<Space> super;
344 public:
345 
347  typedef typename super::space_type space_type;
348  typedef typename space_type::value_type value_type;
349  typedef typename node<value_type>::type node_type;
350 
352  :
353  super()
354  {}
355  template<typename P>
356  IntegralMomentOfDerivative( space_type const& b, uint16_type i )
357  :
358  super( b, b.d( i ) )
359  {}
360 };
361 
368 template<typename Space, typename Poly = Space>
370  :
371 public Functional<Space>
372 {
373  typedef Functional<Space> super;
374 public:
375 
376  BOOST_STATIC_ASSERT( Space::is_vectorial );
377 
379  typedef typename Poly::polynomial_type polynomial_type;
380  typedef typename super::space_type space_type;
381  typedef typename space_type::value_type value_type;
382  typedef typename node<value_type>::type node_type;
383 
385  :
386  super()
387  {}
388 
389  IntegralMomentOfDivergence( space_type const& p, polynomial_type const& q )
390  :
391  super( p )
392  {
393  ublas::matrix<value_type> __rep( space_type::nComponents, p.polynomialDimensionPerComponent() );
394 
395  //std::cout << "q = " << q.coeff() << "\n";
396  for ( int i = 0; i < space_type::nComponents; ++i )
397  {
398 
399  std::cout << "p.d("<< i << ") = " << p.basis().d( i ) << "\n"
400  << "prod = " << ublas::prod( q.coeff(), p.basis().d( i ) ) << "\n";
401 
402  ublas::row( __rep,i ) = ublas::row( ublas::prod( q.coeff(), p.basis().d( i ) ), 0 ) ;
403  }
404 
405  //std::cout << "[IntegralMomentOfDivergence] rep = " << __rep << "\n";
406  this->setCoefficient( __rep );
407  }
408 };
409 
410 } // functional
411 
412 } // Feel
413 
414 #endif // __FEELPP_FUNCTIONALS2_HPP

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