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
crouzeixraviart.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-12
7 
8  Copyright (C) 2006 EPFL
9  Copyright (C) 2011 Universite 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 __CrouzeixRaviart_H
31 #define __CrouzeixRaviart_H 1
32 
33 #include <boost/ptr_container/ptr_vector.hpp>
34 
35 #include <boost/numeric/ublas/vector.hpp>
36 #include <boost/numeric/ublas/io.hpp>
37 #include <boost/numeric/ublas/matrix.hpp>
38 #include <boost/numeric/ublas/matrix_proxy.hpp>
39 #include <boost/numeric/ublas/lu.hpp>
40 
41 #include <feel/feelcore/feel.hpp>
42 #include <feel/feelcore/traits.hpp>
43 #include <feel/feelalg/lu.hpp>
44 
47 
48 #include <feel/feelpoly/fe.hpp>
52 #include <feel/feelpoly/moment.hpp>
53 
54 namespace Feel
55 {
56 namespace fem
57 {
58 namespace detail
59 {
60 template<uint16_type N,
61  template<uint16_type Dim> class PolySetType = Scalar,
62  typename T = double>
63 class RannacherTurekPolynomialSet
64  :
65 public MomentPolynomialSet<N, 2, N, PolySetType, T, Hypercube>
66 {
67  typedef MomentPolynomialSet<N, 2, N, PolySetType, T, Hypercube> super;
68 
69 public:
70 
71  typedef RannacherTurekPolynomialSet<N, PolySetType, T> self_type;
72 
73  typedef typename super::value_type value_type;
74  typedef typename super::convex_type convex_type;
75  typedef typename super::matrix_type matrix_type;
76  typedef typename super::points_type points_type;
77 
78  static const uint16_type nDim = super::nDim;
79  static const uint16_type nOrder = super::nOrder;
80  static const uint16_type nComponents = super::nComponents;
81  static const bool is_product = true;
82 
83  RannacherTurekPolynomialSet()
84  :
85  super()
86  {
87  Moment<2,2,Hypercube<2> > m;
88 
89  for ( int c = 0; c < nComponents; ++c )
90  {
91  // 1
92  this->insert( m.template pick<PolySetType>( 0, c ).toSet( true ), c==0 );
93  // x
94  this->insert( m.template pick<PolySetType>( 1, c ).toSet( true ) );
95  // y
96  this->insert( m.template pick<PolySetType>( 3, c ).toSet( true ) );
97  // x^2 - y^2
98  Polynomial<Moment<2,2,Hypercube<2> >, PolySetType> p( m );
99  p = m.template pick<PolySetType>( 2, c );
100  p -= m.template pick<PolySetType>( 6, c );
101  //std::cout << "x^2-y^2:" << p.coeff() << "\n";
102  this->insert( p.toSet( true ) );
103  }
104 
105  }
106 
107 
108 };
109 
110 template<typename Basis, template<class, uint16_type, class> class PointSetType>
111 class CrouzeixRaviartDual
112  :
113 public DualBasis<Basis>
114 {
115  typedef DualBasis<Basis> super;
116 public:
117 
118  static const uint16_type nDim = super::nDim;
119  static const uint16_type nOrder= super::nOrder;
120 
121  typedef typename super::primal_space_type primal_space_type;
122  typedef typename primal_space_type::value_type value_type;
123  typedef typename primal_space_type::points_type points_type;
124  typedef typename primal_space_type::matrix_type matrix_type;
125  typedef typename primal_space_type::template convex<2>::type convex_type;
126  typedef Reference<convex_type, nDim, 2, nDim, value_type> reference_convex_type;
127  typedef typename reference_convex_type::node_type node_type;
128 
129  // point set type associated with the functionals
130  typedef PointSet<convex_type, value_type> pointset_type;
131  typedef PointSetType<convex_type, 2, value_type> equispaced_pointset_type;
132 
133  static const uint16_type nVertices = reference_convex_type::numVertices;
134  static const uint16_type nFaces = reference_convex_type::numFaces;
135  static const uint16_type nGeometricFaces = reference_convex_type::numFaces;
136  static const uint16_type nEdges = reference_convex_type::numEdges;
137  static const uint16_type nNormals = reference_convex_type::numNormals;
138 
139  static const uint16_type nbPtsPerVertex = 0;
140  static const uint16_type nbPtsPerEdge = mpl::if_<mpl::equal_to<mpl::int_<nDim>,mpl::int_<2> >,
141  mpl::int_<reference_convex_type::nbPtsPerEdge>,
142  mpl::int_<0> >::type::value;
143  static const uint16_type nbPtsPerFace = mpl::if_<mpl::equal_to<mpl::int_<nDim>,mpl::int_<3> >,
144  mpl::int_<reference_convex_type::nbPtsPerFace>,
145  mpl::int_<0> >::type::value;
146  static const uint16_type nbPtsPerVolume = 0;
147  static const uint16_type numPoints = ( reference_convex_type::numGeometricFaces*nbPtsPerFace+
148  reference_convex_type::numEdges*nbPtsPerEdge );
149 
151  static const uint16_type nDofPerVertex = nbPtsPerVertex;
152 
154  static const uint16_type nDofPerEdge = nbPtsPerEdge;
155 
157  static const uint16_type nDofPerFace = nbPtsPerFace;
158 
160  static const uint16_type nDofPerVolume = nbPtsPerVolume;
161 
163  static const uint16_type nLocalDof = numPoints;
164 
165  static const uint16_type nFacesInConvex = mpl::if_< mpl::equal_to<mpl::int_<nDim>, mpl::int_<1> >,
166  mpl::int_<nVertices>,
167  typename mpl::if_<mpl::equal_to<mpl::int_<nDim>, mpl::int_<2> >,
168  mpl::int_<nEdges>,
169  mpl::int_<nFaces> >::type >::type::value;
170 
171 
172  CrouzeixRaviartDual( primal_space_type const& primal )
173  :
174  super( primal ),
175  M_convex_ref(),
176  M_eid( M_convex_ref.topologicalDimension()+1 ),
177  M_pts( nDim, numPoints ),
178  M_points_face( nFacesInConvex ),
179  M_fset( primal )
180  {
181 #if 1
182  std::cout << "Lagrange finite element: \n";
183  std::cout << " o- dim = " << nDim << "\n";
184  std::cout << " o- order = " << nOrder << "\n";
185  std::cout << " o- numPoints = " << numPoints << "\n";
186  std::cout << " o- nbPtsPerVertex = " << ( int )nbPtsPerVertex << "\n";
187  std::cout << " o- nbPtsPerEdge = " << ( int )nbPtsPerEdge << "\n";
188  std::cout << " o- nbPtsPerFace = " << ( int )nbPtsPerFace << "\n";
189  std::cout << " o- nbPtsPerVolume = " << ( int )nbPtsPerVolume << "\n";
190 #endif
191  equispaced_pointset_type epts;
192  // in d-dimension, consider only the d-1 entity mid-points
193  int d = M_convex_ref.topologicalDimension()-1;
194  int p = 0;
195 
196  // loop on each entity forming the convex of topological
197  // dimension d
198  for ( int e = M_convex_ref.entityRange( d ).begin();
199  e < M_convex_ref.entityRange( d ).end();
200  ++e )
201  {
202  M_points_face[e] = epts.pointsBySubEntity( nDim-1, e, 0 );
203  ublas::subrange( M_pts, 0, nDim, p, p+M_points_face[e].size2() ) = M_points_face[e];
204  p+=M_points_face[e].size2();
205  }
206 
207  //std::cout << "[CrouzeixRaviartDual] points= " << M_pts << "\n";
208  setFset( primal, M_pts, mpl::bool_<primal_space_type::is_scalar>() );
209 
210 
211  }
212 
213  points_type const& points() const
214  {
215  return M_pts;
216  }
217  points_type const& points( uint16_type f ) const
218  {
219  return M_points_face[f];
220  }
221 
222 
223  matrix_type operator()( primal_space_type const& pset ) const
224  {
225  return M_fset( pset );
226  }
227 private:
228 
229  void setFset( primal_space_type const& primal, points_type const& __pts, mpl::bool_<true> )
230  {
231  M_fset.setFunctionalSet( functional::PointsEvaluation<primal_space_type>( primal,
232  __pts ) );
233  }
234 
235  void setFset( primal_space_type const& primal, points_type const& __pts, mpl::bool_<false> )
236  {
237  M_fset.setFunctionalSet( functional::ComponentsPointsEvaluation<primal_space_type>( primal,
238  __pts ) );
239  }
240 
241 
242 private:
243  reference_convex_type M_convex_ref;
244  std::vector<std::vector<uint16_type> > M_eid;
245  points_type M_pts;
246  std::vector<points_type> M_points_face;
247  FunctionalSet<primal_space_type> M_fset;
248 
249 
250 };
251 }// detail
252 
253 
260 template<uint16_type N,
261  uint16_type RealDim,
262  template<uint16_type Dim> class PolySetType,
263  typename T = double,
264  template<uint16_type, uint16_type, uint16_type> class Convex = Simplex,
265  uint16_type TheTAG=0 >
267  :
268 public FiniteElement<typename mpl::if_<mpl::bool_<Convex<N,1,N>::is_simplex>,
269  mpl::identity<Feel::detail::OrthonormalPolynomialSet<N, 1, RealDim, PolySetType, T, TheTAG, Convex> >,
270  mpl::identity<fem::detail::RannacherTurekPolynomialSet<N, PolySetType, T> > >::type::type,
271  detail::CrouzeixRaviartDual,
272  PointSetEquiSpaced >
273 {
275  mpl::identity<Feel::detail::OrthonormalPolynomialSet<N, 1, RealDim, PolySetType, T, TheTAG, Convex> >,
276  mpl::identity<fem::detail::RannacherTurekPolynomialSet<N, PolySetType, T> > >::type::type,
277  detail::CrouzeixRaviartDual,
278  PointSetEquiSpaced > super;
279 public:
280 
281  BOOST_STATIC_ASSERT( N > 1 );
282 
286 
287  static const uint16_type nDim = N;
288  static const uint16_type nOrder = super::nOrder;
289  static const bool isTransformationEquivalent = true;
290  static const bool isContinuous = true;
291 
292  typedef typename super::value_type value_type;
293  typedef typename super::primal_space_type primal_space_type;
294  typedef typename super::dual_space_type dual_space_type;
295  typedef Continuous continuity_type;
296  static const uint16_type TAG = TheTAG;
297 
302  static const bool is_vectorial = polyset_type::is_vectorial;
303  static const bool is_scalar = polyset_type::is_scalar;
304  static const uint16_type nComponents = polyset_type::nComponents;
305  static const bool is_product = true;
306 
308 
309  typedef typename dual_space_type::convex_type convex_type;
310  typedef typename dual_space_type::pointset_type pointset_type;
311  typedef typename dual_space_type::reference_convex_type reference_convex_type;
312  typedef typename reference_convex_type::node_type node_type;
313  typedef typename reference_convex_type::points_type points_type;
314 
315 
316  static const uint16_type nbPtsPerVertex = 0;
317  static const uint16_type nbPtsPerEdge = mpl::if_<mpl::equal_to<mpl::int_<nDim>,mpl::int_<2> >,
318  mpl::int_<reference_convex_type::nbPtsPerEdge>,
319  mpl::int_<0> >::type::value;
320  static const uint16_type nbPtsPerFace = mpl::if_<mpl::equal_to<mpl::int_<nDim>,mpl::int_<3> >,
321  mpl::int_<reference_convex_type::nbPtsPerFace>,
322  mpl::int_<0> >::type::value;
323  static const uint16_type nbPtsPerVolume = 0;
324  static const uint16_type numPoints = ( reference_convex_type::numGeometricFaces*nbPtsPerFace+
325  reference_convex_type::numEdges*nbPtsPerEdge );
327 
331 
333  :
334  super( dual_space_type( primal_space_type() ) ),
335  M_refconvex()
336  {}
337  CrouzeixRaviart( CrouzeixRaviart const & cr )
338  :
339  super( cr ),
340  M_refconvex()
341  {}
342  ~CrouzeixRaviart()
343  {}
344 
346 
350 
351 
353 
357 
361  reference_convex_type const& referenceConvex() const
362  {
363  return M_refconvex;
364  }
365 
367 
371 
372 
374 
378 
382  std::string familyName() const
383  {
384  return "CrouzeixRaviart";
385  }
386 
387  template<typename ExprType>
388  static auto
389  isomorphism( ExprType& expr ) -> decltype( expr )
390  {
391  return expr;
392  //return expr;
393  }
395 
396 
397 
398 protected:
399  reference_convex_type M_refconvex;
400 private:
401 
402 };
403 template<uint16_type N,
404  uint16_type RealDim,
405  template<uint16_type Dim> class PolySetType,
406  typename T,
407  template<uint16_type, uint16_type, uint16_type> class Convex,
408  uint16_type TheTAG >
409 const uint16_type CrouzeixRaviart<N,RealDim,PolySetType,T,Convex,TheTAG>::nDim;
410 template<uint16_type N,
411  uint16_type RealDim,
412  template<uint16_type Dim> class PolySetType,
413  typename T,
414  template<uint16_type, uint16_type, uint16_type> class Convex,
415  uint16_type TheTAG >
416 const uint16_type CrouzeixRaviart<N,RealDim,PolySetType,T,Convex,TheTAG>::nOrder;
417 
418 } // fem
419 
420 template<uint16_type Order,
421  template<uint16_type Dim> class PolySetType = Scalar,
422  template<class, uint16_type, class> class Pts = PointSetEquiSpaced,
423  uint16_type TheTAG=0 >
424 class CrouzeixRaviart
425 {
426 public:
427  template<uint16_type N,
428  uint16_type RealDim,
429  typename T = double,
430  typename Convex = Simplex<N> >
431  struct apply
432  {
433  typedef typename mpl::if_<mpl::bool_<Convex::is_simplex>,
434  mpl::identity<fem::CrouzeixRaviart<N,RealDim,PolySetType,T,Simplex,TheTAG> >,
435  mpl::identity<fem::CrouzeixRaviart<N,RealDim,PolySetType,T,Hypercube,TheTAG> > >::type::type result_type;
436  typedef result_type type;
437  };
438 
439  template<uint16_type TheNewTAG>
440  struct ChangeTag
441  {
442  typedef CrouzeixRaviart<Order,PolySetType,Pts,TheNewTAG> type;
443  };
444 
445  typedef CrouzeixRaviart<Order,Scalar,Pts> component_basis_type;
446 
447  static const uint16_type nOrder = Order;
448  static const uint16_type TAG = TheTAG;
449 
450 };
451 
452 template<uint16_type Order,
453  template<uint16_type Dim> class PolySetType,
454  template<class, uint16_type, class> class Pts,
455  uint16_type TheTAG >
456 const uint16_type CrouzeixRaviart<Order,PolySetType,Pts,TheTAG>::nOrder;
457 
458 
459 } // Feel
460 #endif /* __CrouzeixRaviart_H */

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