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
lagrange.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-08-18
7 
8  Copyright (C) 2005,2006 EPFL
9  Copyright (C) 2008-2012 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 __lagrange_H
31 #define __lagrange_H 1
32 
33 #include <boost/ptr_container/ptr_vector.hpp>
34 
35 #include <feel/feelcore/feel.hpp>
36 #include <feel/feelcore/traits.hpp>
37 #include <feel/feelalg/lu.hpp>
38 
42 
43 
48 #include <feel/feelpoly/fe.hpp>
49 #include <feel/feelpoly/mortar.hpp>
50 
51 
52 
53 
54 
55 namespace Feel
56 {
57 
58 namespace fem
59 {
60 
62 namespace details
63 {
64 template<typename Basis, template<class, uint16_type, class> class PointSetType>
65 class LagrangeDual
66  :
67 public DualBasis<Basis>
68 {
69  typedef DualBasis<Basis> super;
70 public:
71 
72  static const uint16_type nDim = super::nDim;
73  static const uint16_type nOrder= super::nOrder;
74 
75  typedef typename super::primal_space_type primal_space_type;
76  typedef typename primal_space_type::value_type value_type;
77  typedef typename primal_space_type::points_type points_type;
78  typedef typename primal_space_type::matrix_type matrix_type;
79  typedef typename primal_space_type::convex_type convex_type;
80  typedef typename primal_space_type::reference_convex_type reference_convex_type;
81  typedef typename reference_convex_type::node_type node_type;
82 
83  // point set type associated with the functionals
84  typedef PointSetType<convex_type, nOrder, value_type> pointset_type;
85 
86  static const uint16_type numPoints = reference_convex_type::numPoints;
87  static const uint16_type nbPtsPerVertex = reference_convex_type::nbPtsPerVertex;
88  static const uint16_type nbPtsPerEdge = reference_convex_type::nbPtsPerEdge;
89  static const uint16_type nbPtsPerFace = reference_convex_type::nbPtsPerFace;
90  static const uint16_type nbPtsPerVolume = reference_convex_type::nbPtsPerVolume;
91 
92 #if 0
93 
97  typedef mpl::vector_c<uint16_type, 0, 1, 3, ( 4 ) + ( 4 ) - 2> edges_t;
98  typedef mpl::vector_c<uint16_type, 0, 0, 1, 4> geo_faces_t;
99  typedef mpl::vector_c<uint16_type, 0, 2, 3, 4> faces_t;
100  typedef mpl::vector_c<uint16_type, 0, 2, 3, 4> normals_t;
101 #endif
102 
103  static const uint16_type nVertices = reference_convex_type::numVertices;
104  static const uint16_type nFaces = reference_convex_type::numFaces;
105  static const uint16_type nGeometricFaces = reference_convex_type::numFaces;
106  static const uint16_type nEdges = reference_convex_type::numEdges;
107  static const uint16_type nNormals = reference_convex_type::numNormals;
108 
109 
111  static const uint16_type nDofPerVertex = nbPtsPerVertex;
112 
114  static const uint16_type nDofPerEdge = nbPtsPerEdge;
115 
117  static const uint16_type nDofPerFace = nbPtsPerFace;
118 
120  static const uint16_type nDofPerVolume = nbPtsPerVolume;
121 
123  static const uint16_type nLocalDof = numPoints;
124 
125  static const uint16_type nFacesInConvex = mpl::if_< mpl::equal_to<mpl::int_<nDim>, mpl::int_<1> >,
126  mpl::int_<nVertices>,
127  typename mpl::if_<mpl::equal_to<mpl::int_<nDim>, mpl::int_<2> >,
128  mpl::int_<nEdges>,
129  mpl::int_<nFaces> >::type >::type::value;
130 
131  LagrangeDual( primal_space_type const& primal )
132  :
133  super( primal ),
134  M_convex_ref(),
135  M_eid( M_convex_ref.topologicalDimension()+1 ),
136  M_pts( nDim, numPoints ),
137  M_points_face( nFacesInConvex ),
138  M_fset( primal )
139  {
140  DVLOG(2) << "Lagrange finite element: \n";
141  DVLOG(2) << " o- dim = " << nDim << "\n";
142  DVLOG(2) << " o- order = " << nOrder << "\n";
143  DVLOG(2) << " o- numPoints = " << numPoints << "\n";
144  DVLOG(2) << " o- nbPtsPerVertex = " << nbPtsPerVertex << "\n";
145  DVLOG(2) << " o- nbPtsPerEdge = " << nbPtsPerEdge << "\n";
146  DVLOG(2) << " o- nbPtsPerFace = " << nbPtsPerFace << "\n";
147  DVLOG(2) << " o- nbPtsPerVolume = " << nbPtsPerVolume << "\n";
148 
149  pointset_type pts;
150 
151  M_pts = pts.points();
152 
153  if ( nOrder > 0 )
154  {
155  for ( uint16_type e = M_convex_ref.entityRange( nDim-1 ).begin();
156  e < M_convex_ref.entityRange( nDim-1 ).end();
157  ++e )
158  {
159  M_points_face[e] = pts.pointsBySubEntity( nDim-1, e, 1 );
160  DVLOG(2) << "face " << e << " pts " << M_points_face[e] << "\n";
161  }
162  }
163 
164  setFset( primal, M_pts, mpl::bool_<primal_space_type::is_scalar>() );
165  }
166 
167  LagrangeDual( primal_space_type const& primal, pointset_type const& pts )
168  :
169  super( primal ),
170  M_convex_ref(),
171  M_eid( M_convex_ref.topologicalDimension()+1 ),
172  M_pts( pts.points() ),
173  M_points_face( nFacesInConvex ),
174  M_fset( primal )
175  {
176  DVLOG(2) << "Lagrange finite element: \n";
177  DVLOG(2) << " o- dim = " << nDim << "\n";
178  DVLOG(2) << " o- order = " << nOrder << "\n";
179  DVLOG(2) << " o- numPoints = " << numPoints << "\n";
180  DVLOG(2) << " o- nbPtsPerVertex = " << nbPtsPerVertex << "\n";
181  DVLOG(2) << " o- nbPtsPerEdge = " << nbPtsPerEdge << "\n";
182  DVLOG(2) << " o- nbPtsPerFace = " << nbPtsPerFace << "\n";
183  DVLOG(2) << " o- nbPtsPerVolume = " << nbPtsPerVolume << "\n";
184 
185  if ( nOrder > 0 )
186  {
187  for ( uint16_type e = M_convex_ref.entityRange( nDim-1 ).begin();
188  e < M_convex_ref.entityRange( nDim-1 ).end();
189  ++e )
190  {
191  M_points_face[e] = pts.pointsBySubEntity( nDim-1, e, 1 );
192  DVLOG(2) << "face " << e << " pts " << M_points_face[e] << "\n";
193  }
194  }
195 
196  setFset( primal, M_pts, mpl::bool_<primal_space_type::is_scalar>() );
197  }
198 
199  ~LagrangeDual()
200  {
201 
202  }
203  points_type const& points() const
204  {
205  return M_pts;
206  }
207 
208  points_type const& points( uint16_type f ) const
209  {
210  return M_points_face[f];
211  }
212  ublas::matrix_column<points_type const> point( uint16_type f, uint32_type __i ) const
213  {
214  return ublas::column( M_points_face[f], __i );
215  }
216  ublas::matrix_column<points_type> point( uint16_type f, uint32_type __i )
217  {
218  return ublas::column( M_points_face[f], __i );
219  }
220 
221  matrix_type operator()( primal_space_type const& pset ) const
222  {
223  return M_fset( pset );
224  }
225 private:
226 
227  void setFset( primal_space_type const& primal, points_type const& __pts, mpl::bool_<true> )
228  {
229  M_fset.setFunctionalSet( functional::PointsEvaluation<primal_space_type>( primal,
230  __pts ) );
231  }
232 
233  void setFset( primal_space_type const& primal, points_type const& __pts, mpl::bool_<false> )
234  {
235  M_fset.setFunctionalSet( functional::ComponentsPointsEvaluation<primal_space_type>( primal,
236  __pts ) );
237  }
238 
242  void setPoints( uint16_type f, points_type const& n )
243  {
244  M_points_face[f].resize( n.size1(), n.size2(), false );
245  M_points_face[f] = n;
246  }
247 
248 private:
249  reference_convex_type M_convex_ref;
250  std::vector<std::vector<uint16_type> > M_eid;
251  points_type M_pts;
252  std::vector<points_type> M_points_face;
253  FunctionalSet<primal_space_type> M_fset;
254 
255 
256 };
257 }// details
259 
275 template<uint16_type N,
276  uint16_type RealDim,
277  uint16_type O,
278  template<uint16_type Dim> class PolySetType,
279  typename ContinuityType = Continuous,
280  typename T = double,
281  template<uint16_type, uint16_type, uint16_type> class Convex = Simplex,
282  template<class, uint16_type, class> class Pts = PointSetEquiSpaced,
283  uint16_type TheTAG = 0 >
284 class Lagrange
285  :
286  public FiniteElement<detail::OrthonormalPolynomialSet<N, O, RealDim, PolySetType, T, TheTAG, Convex>, details::LagrangeDual, Pts >
287 {
289 public:
290 
291  BOOST_STATIC_ASSERT( ( boost::is_same<PolySetType<N>, Scalar<N> >::value ||
292  boost::is_same<PolySetType<N>, Vectorial<N> >::value ||
293  boost::is_same<PolySetType<N>, Tensor2<N> >::value ) );
294 
298 
299  static const uint16_type nDim = N;
300  static const uint16_type nRealDim = RealDim;
301  static const uint16_type nOrder = O;
302  static const bool isTransformationEquivalent = true;
303  static const bool isContinuous = ContinuityType::is_continuous;
304  typedef typename super::value_type value_type;
305  typedef typename super::primal_space_type primal_space_type;
306  typedef typename super::dual_space_type dual_space_type;
307  typedef ContinuityType continuity_type;
308  static const uint16_type TAG = TheTAG;
309 
314  static const bool is_tensor2 = polyset_type::is_tensor2;
315  static const bool is_vectorial = polyset_type::is_vectorial;
316  static const bool is_scalar = polyset_type::is_scalar;
317  static const uint16_type nComponents = polyset_type::nComponents;
318  static const bool is_product = true;
319 
321 
322  typedef typename mpl::if_<mpl::equal_to<mpl::int_<nDim>, mpl::int_<1> >,
323  mpl::identity<boost::none_t>,
324  mpl::identity< Lagrange<N-1, RealDim, O, Scalar, continuity_type, T, Convex, Pts, TheTAG> > >::type::type face_basis_type;
325 
326  typedef boost::shared_ptr<face_basis_type> face_basis_ptrtype;
327 
328  typedef typename dual_space_type::convex_type convex_type;
329  typedef typename dual_space_type::pointset_type pointset_type;
330  typedef typename dual_space_type::reference_convex_type reference_convex_type;
331  typedef typename reference_convex_type::node_type node_type;
332  typedef typename reference_convex_type::points_type points_type;
333 
334  static const uint16_type numPoints = reference_convex_type::numPoints;
335  static const uint16_type nbPtsPerVertex = reference_convex_type::nbPtsPerVertex;
336  static const uint16_type nbPtsPerEdge = reference_convex_type::nbPtsPerEdge;
337  static const uint16_type nbPtsPerFace = reference_convex_type::nbPtsPerFace;
338  static const uint16_type nbPtsPerVolume = reference_convex_type::nbPtsPerVolume;
339 
340  template<int subN>
341  struct SubSpace
342  {
343  typedef Lagrange<N-1, RealDim, O, PolySetType, continuity_type, T, Convex, Pts, TheTAG> type;
344  };
345 
347 
351 
352  Lagrange()
353  :
354  super( dual_space_type( primal_space_type() ) ),
355  M_refconvex()
356  // M_bdylag( new face_basis_type )
357  {
358 
359 
360 
361  // std::cout << "[LagrangeDual] points= " << M_pts << "\n";
362  }
363 
364  virtual ~Lagrange() {}
365 
367 
371 
372 
374 
378 
382  reference_convex_type const& referenceConvex() const
383  {
384  return M_refconvex;
385  }
386 
390  std::string familyName() const
391  {
392  return "lagrange";
393  }
394 
396 
400 
401 
403 
407 
408  template<typename ExprType>
409  static auto
410  isomorphism( ExprType& expr ) -> decltype( expr )
411  {
412  return expr;
413  //return expr;
414  }
416 
417 private:
418 
419  reference_convex_type M_refconvex;
420  face_basis_ptrtype M_bdylag;
421 };
422 template<uint16_type N,
423  uint16_type RealDim,
424  uint16_type O,
425  template<uint16_type Dim> class PolySetType,
426  typename ContinuityType,
427  typename T,
428  template<uint16_type, uint16_type, uint16_type> class Convex,
429  template<class, uint16_type, class> class Pts,
430  uint16_type TheTAG >
431 const uint16_type Lagrange<N,RealDim,O,PolySetType,ContinuityType,T,Convex,Pts,TheTAG>::nDim;
432 
433 template<uint16_type N,
434  uint16_type RealDim,
435  uint16_type O,
436  template<uint16_type Dim> class PolySetType,
437  typename ContinuityType,
438  typename T,
439  template<uint16_type, uint16_type, uint16_type> class Convex,
440  template<class, uint16_type, class> class Pts,
441  uint16_type TheTAG >
442 const uint16_type Lagrange<N,RealDim,O,PolySetType,ContinuityType,T,Convex,Pts,TheTAG>::nOrder;
443 
444 } // namespace fem
445 template<uint16_type Order,
446  template<uint16_type Dim> class PolySetType = Scalar,
447  typename ContinuityType = Continuous,
448  template<class, uint16_type, class> class Pts = PointSetEquiSpaced,
449  uint16_type TheTAG=0 >
450 class Lagrange
451 {
452 public:
453  template<uint16_type N,
454  uint16_type RealDim,
455  typename T = double,
456  typename Convex = Simplex<N> >
457  struct apply
458  {
459  typedef typename mpl::if_<mpl::bool_<Convex::is_simplex>,
460  mpl::identity<fem::Lagrange<N,RealDim,Order,PolySetType,ContinuityType,T,Simplex, Pts, TheTAG > >,
461  mpl::identity<fem::Lagrange<N,RealDim,Order,PolySetType,ContinuityType,T,Hypercube, Pts, TheTAG> > >::type::type result_type;
462  typedef result_type type;
463  };
464 
465  template<uint16_type TheNewTAG>
466  struct ChangeTag
467  {
468  typedef Lagrange<Order,PolySetType,ContinuityType,Pts,TheNewTAG> type;
469  };
470 
471  typedef Lagrange<Order,Scalar,ContinuityType,Pts,TheTAG> component_basis_type;
472 
473  static const uint16_type nOrder = Order;
474  static const uint16_type TAG = TheTAG;
475 
476 };
477 template<uint16_type Order,
478  template<uint16_type Dim> class PolySetType,
479  typename ContinuityType,
480  template<class, uint16_type, class> class Pts,
481  uint16_type TheTAG>
482 const uint16_type Lagrange<Order,PolySetType,ContinuityType,Pts,TheTAG>::nOrder;
483 
484 } // namespace Feel
485 #endif /* __lagrange_H */

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