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
fe.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: 2005-10-06
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 __FiniteElement_H
30 #define __FiniteElement_H 1
31 
32 #include <feel/feelpoly/policy.hpp>
34 
35 namespace Feel
36 {
37 template<typename Poly, template<uint16_type> class PolySetType > class PolynomialSet;
38 namespace detail
39 {
40 template<uint16_type Dim,
41  uint16_type Order,
42  uint16_type RealDim,
43  template<uint16_type> class PolySetType,
44  typename T,
45  uint16_type TheTAG,
46  template<uint16_type,uint16_type,uint16_type> class Convex>
47 class OrthonormalPolynomialSet;
48 }
57 template<typename P,
58  template<class Pr, template<class,uint16_type,class> class Pt> class PDual,
59  template<class,uint16_type,class> class Pts>
61  public mpl::if_<mpl::bool_<P::is_scalar>,
62  mpl::identity<PolynomialSet<P, Scalar> >,
63  mpl::identity<PolynomialSet<P, Vectorial> > >::type::type
64 {
65  typedef typename mpl::if_<mpl::bool_<P::is_scalar>,
66  mpl::identity<PolynomialSet<P, Scalar> >,
67  mpl::identity<PolynomialSet<P, Vectorial> > >::type::type super;
68 
69 public:
70 
74 
76 
77  typedef typename P::value_type value_type;
78 
79  typedef P primal_space_type;
80 
84  typedef typename primal_space_type::polyset_type polyset_type;
85 
86  static const bool is_modal = false;
87 
88  typedef PDual<P, Pts> dual_space_type;
89 
90  typedef typename super::matrix_type matrix_type;
91  typedef typename super::points_type points_type;
92  typedef typename super::self_type polynomialset_type;
93 
94  typedef typename super::polynomial_type polynomial_type;
95  typedef typename super::polynomial_view_type polynomial_view_type;
96 
98  static const uint16_type nLocalDof = dual_space_type::nLocalDof;
100  static const uint16_type nDofPerVertex = dual_space_type::nDofPerVertex;
102  static const uint16_type nDofPerEdge = dual_space_type::nDofPerEdge;
104  static const uint16_type nDofPerFace = dual_space_type::nDofPerFace;
106  static const uint16_type nDofPerVolume = dual_space_type::nDofPerVolume;
107 
108  static const uint16_type nDof = nLocalDof;
109  static const uint16_type nNodes = nDof;
110  static const uint16_type nDofGrad = super::nDim*nDof;
111  static const uint16_type nDofHess = super::nDim*super::nDim*nDof;
112  static const fem::transformation_type trans = ( fem::transformation_type )mpl::if_<mpl::and_<mpl::bool_<P::convex_type::is_simplex>,mpl::equal_to<mpl::int_<super::nOrder>, mpl::int_<1> > >,
113  mpl::int_<fem::LINEAR>,
114  mpl::int_<fem::NONLINEAR> >::type::value;
116 
120 
121  FiniteElement( dual_space_type const& pdual )
122  :
123  super( pdual.primalSpace() ),
124  M_dual( pdual ),
125  M_primal( M_dual.primalSpace() )
126  {
127  DVLOG(2) << "============================================================\n";
128  DVLOG(2) << "New FE \n";
129  ublas::matrix<value_type> A( M_dual( M_primal ) );
130  //std::cout << "[FiniteElement] A = " << A << "\n";
131 
132  ublas::matrix<value_type> D = ublas::identity_matrix<value_type>( A.size1(), A.size2() );
134  ublas::matrix<value_type> C = lu.solve( D );
135  //std::cout << "[FiniteElement] D = " << D << "\n";
136  //std::cout << "[FiniteElement] C = " << C << "\n";
137  DVLOG(2) << "is singular : " << lu.isNonsingular() << "\n"
138  << "det(A) = " << lu.det() << "\n";
139 #if 0
140 
141  if ( !lu.isNonsingular() )
142  {
143  std::cout << "A=" << A << "\n"
144  << "D=" << D << "\n"
145  << "C=" << C << "\n";
146  }
147 
148 #endif
149  FEELPP_ASSERT( lu.isNonsingular() )( A )( D )( C ).error( "vandermonde matrix is singular" );
150 
151  this->setCoefficient( ublas::trans( C ) );
152 
153  //M_pset = polynomialset_type( M_primal, C );
154 
155  //std::cout << "coeff = " << M_pset.coeff() << "\n";
156  //std::cout << "d_x = " << M_pset.derivate(0).coeff() << "\n";
157  //std::cout << "d_x = " << M_pset.derivate(0).coeff() << "\n";
158  //std::cout << "d_x = " << M_pset.derivate(0).coeff() << "\n";
159  DVLOG(2) << "============================================================\n";
160  }
161  FiniteElement( FiniteElement const & fe )
162  :
163  super( fe ),
164  M_dual( fe.M_dual ),
165  M_primal( fe.M_primal )
166  {}
167 
168  ~FiniteElement()
169  {}
170 
172 
176 
177  self_type& operator=( self_type const& fe )
178  {
179  if ( this != &fe )
180  {
181  super::operator=( fe );
182  M_primal = fe.M_primal;
183  M_dual = fe.M_dual;
184  }
185 
186  return *this;
187  }
188 
189  template<typename AE>
190  value_type operator()( uint16_type i, ublas::vector_expression<AE> const& pt ) const
191  {
192  return this->evaluate( i, pt );
193  }
194 
195  template<typename AE>
196  value_type operator()( ublas::vector_expression<AE> const& pt ) const
197  {
198  matrix_type m( pt().size(), 1 );
199  ublas::column( m, 0 ) = pt;
200  ublas::vector<value_type> r( ublas::column( this->evaluate( m ), 0 ) );
201  return ublas::inner_prod( r, ublas::scalar_vector<value_type>( r.size(), 1.0 ) );
202  }
203 
204 
205  matrix_type operator()( points_type const& pts ) const
206  {
207  return this->evaluate( pts );
208  }
209 
211 
215 
216 
220  void domainShape() const {}
221 
225  uint16_type nbPoints() const
226  {
227  return points().size2();
228  }
229 
233  // Que devient cette fonction ??
234  // polynomialset_type const& functionShape() const { return M_pset; }
235 
239  primal_space_type const& primal() const
240  {
241  return M_primal;
242  }
243 
247  dual_space_type const& dual() const
248  {
249  return M_dual;
250  }
251 
255  points_type const& points() const
256  {
257  return M_dual.points();
258  }
259 
269  points_type const& points( uint16_type f ) const
270  {
271  return M_dual.points( f );
272  }
273 
277  virtual std::string familyName() const = 0;
278 
280 
281 
282 private:
283 
284  dual_space_type M_dual;
285  primal_space_type const& M_primal;
286 
287 };
288 
289 template<typename P,
290  template<class Pr, template<class,uint16_type,class> class Pt> class PDual,
291  template<class,uint16_type,class> class Pts>
292 const uint16_type FiniteElement<P,PDual,Pts>::nLocalDof;
293 
294 template<typename P,
295  template<class Pr, template<class,uint16_type,class> class Pt> class PDual,
296  template<class,uint16_type,class> class Pts>
297 const uint16_type FiniteElement<P,PDual,Pts>::nDofPerVertex;
298 
299 template<typename P,
300  template<class Pr, template<class,uint16_type,class> class Pt> class PDual,
301  template<class,uint16_type,class> class Pts>
302 const uint16_type FiniteElement<P,PDual,Pts>::nDofPerEdge;
303 
304 template<typename P,
305  template<class Pr, template<class,uint16_type,class> class Pt> class PDual,
306  template<class,uint16_type,class> class Pts>
307 const uint16_type FiniteElement<P,PDual,Pts>::nDofPerFace;
308 
309 template<typename P,
310  template<class Pr, template<class,uint16_type,class> class Pt> class PDual,
311  template<class,uint16_type,class> class Pts>
312 const uint16_type FiniteElement<P,PDual,Pts>::nDofPerVolume;
313 }
314 #endif /* __FiniteElement_H */

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