31 #define __lagrange_H 1
33 #include <boost/ptr_container/ptr_vector.hpp>
35 #include <feel/feelcore/feel.hpp>
36 #include <feel/feelcore/traits.hpp>
64 template<
typename Basis,
template<
class, u
int16_type,
class>
class PointSetType>
67 public DualBasis<Basis>
69 typedef DualBasis<Basis> super;
72 static const uint16_type nDim = super::nDim;
73 static const uint16_type nOrder= super::nOrder;
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;
84 typedef PointSetType<convex_type, nOrder, value_type> pointset_type;
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;
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;
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;
111 static const uint16_type nDofPerVertex = nbPtsPerVertex;
114 static const uint16_type nDofPerEdge = nbPtsPerEdge;
117 static const uint16_type nDofPerFace = nbPtsPerFace;
120 static const uint16_type nDofPerVolume = nbPtsPerVolume;
123 static const uint16_type nLocalDof = numPoints;
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> >,
129 mpl::int_<nFaces> >::type >::type::value;
131 LagrangeDual( primal_space_type
const& primal )
135 M_eid( M_convex_ref.topologicalDimension()+1 ),
136 M_pts( nDim, numPoints ),
137 M_points_face( nFacesInConvex ),
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";
151 M_pts = pts.points();
155 for ( uint16_type e = M_convex_ref.entityRange( nDim-1 ).begin();
156 e < M_convex_ref.entityRange( nDim-1 ).end();
159 M_points_face[e] = pts.pointsBySubEntity( nDim-1, e, 1 );
160 DVLOG(2) <<
"face " << e <<
" pts " << M_points_face[e] <<
"\n";
164 setFset( primal, M_pts, mpl::bool_<primal_space_type::is_scalar>() );
167 LagrangeDual( primal_space_type
const& primal, pointset_type
const& pts )
171 M_eid( M_convex_ref.topologicalDimension()+1 ),
173 M_points_face( nFacesInConvex ),
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";
187 for ( uint16_type e = M_convex_ref.entityRange( nDim-1 ).begin();
188 e < M_convex_ref.entityRange( nDim-1 ).end();
191 M_points_face[e] = pts.pointsBySubEntity( nDim-1, e, 1 );
192 DVLOG(2) <<
"face " << e <<
" pts " << M_points_face[e] <<
"\n";
196 setFset( primal, M_pts, mpl::bool_<primal_space_type::is_scalar>() );
203 points_type
const&
points()
const
208 points_type
const&
points( uint16_type f )
const
210 return M_points_face[f];
212 ublas::matrix_column<points_type const> point( uint16_type f, uint32_type __i )
const
214 return ublas::column( M_points_face[f], __i );
216 ublas::matrix_column<points_type> point( uint16_type f, uint32_type __i )
218 return ublas::column( M_points_face[f], __i );
221 matrix_type operator()( primal_space_type
const& pset )
const
223 return M_fset( pset );
227 void setFset( primal_space_type
const& primal, points_type
const& __pts, mpl::bool_<true> )
229 M_fset.setFunctionalSet( functional::PointsEvaluation<primal_space_type>( primal,
233 void setFset( primal_space_type
const& primal, points_type
const& __pts, mpl::bool_<false> )
235 M_fset.setFunctionalSet( functional::ComponentsPointsEvaluation<primal_space_type>( primal,
242 void setPoints( uint16_type f, points_type
const& n )
244 M_points_face[f].resize( n.size1(), n.size2(), false );
245 M_points_face[f] = n;
249 reference_convex_type M_convex_ref;
250 std::vector<std::vector<uint16_type> > M_eid;
252 std::vector<points_type> M_points_face;
253 FunctionalSet<primal_space_type> M_fset;
275 template<uint16_type N,
278 template<u
int16_type Dim>
class PolySetType,
279 typename ContinuityType = Continuous,
281 template<u
int16_type, u
int16_type, u
int16_type>
class Convex = Simplex,
282 template<
class, u
int16_type,
class>
class Pts = PointSetEquiSpaced,
283 uint16_type TheTAG = 0 >
286 public FiniteElement<detail::OrthonormalPolynomialSet<N, O, RealDim, PolySetType, T, TheTAG, Convex>, details::LagrangeDual, Pts >
291 BOOST_STATIC_ASSERT( ( boost::is_same<PolySetType<N>,
Scalar<N> >::value ||
293 boost::is_same<PolySetType<N>,
Tensor2<N> >::value ) );
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;
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;
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;
326 typedef boost::shared_ptr<face_basis_type> face_basis_ptrtype;
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;
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;
343 typedef Lagrange<N-1, RealDim, O, PolySetType, continuity_type, T,
Convex, Pts, TheTAG> type;
354 super( dual_space_type( primal_space_type() ) ),
364 virtual ~Lagrange() {}
408 template<
typename ExprType>
410 isomorphism( ExprType& expr ) -> decltype( expr )
419 reference_convex_type M_refconvex;
420 face_basis_ptrtype M_bdylag;
422 template<uint16_type N,
425 template<u
int16_type Dim>
class PolySetType,
426 typename ContinuityType,
428 template<u
int16_type, u
int16_type, u
int16_type>
class Convex,
429 template<
class, u
int16_type,
class>
class Pts,
431 const uint16_type Lagrange<N,RealDim,O,PolySetType,ContinuityType,T,Convex,Pts,TheTAG>::nDim;
433 template<uint16_type N,
436 template<u
int16_type Dim>
class PolySetType,
437 typename ContinuityType,
439 template<u
int16_type, u
int16_type, u
int16_type>
class Convex,
440 template<
class, u
int16_type,
class>
class Pts,
442 const uint16_type Lagrange<N,RealDim,O,PolySetType,ContinuityType,T,Convex,Pts,TheTAG>::nOrder;
445 template<uint16_type Order,
446 template<u
int16_type Dim>
class PolySetType = Scalar,
447 typename ContinuityType = Continuous,
448 template<
class, u
int16_type,
class>
class Pts = PointSetEquiSpaced,
449 uint16_type TheTAG=0 >
453 template<uint16_type N,
456 typename Convex = Simplex<N> >
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;
465 template<u
int16_type TheNewTAG>
468 typedef Lagrange<Order,PolySetType,ContinuityType,Pts,TheNewTAG> type;
471 typedef Lagrange<Order,Scalar,ContinuityType,Pts,TheTAG> component_basis_type;
473 static const uint16_type nOrder = Order;
474 static const uint16_type TAG = TheTAG;
477 template<uint16_type Order,
478 template<u
int16_type Dim>
class PolySetType,
479 typename ContinuityType,
480 template<
class, u
int16_type,
class>
class Pts,
482 const uint16_type Lagrange<Order,PolySetType,ContinuityType,Pts,TheTAG>::nOrder;
super::polyset_type polyset_type
Definition: lagrange.hpp:313
std::string familyName() const
Definition: lagrange.hpp:390
Definition: feelpoly/policy.hpp:98
boost::tuple< mpl::size_t< MESH_POINTS >, typename MeshTraits< MeshType >::point_const_iterator, typename MeshTraits< MeshType >::point_const_iterator > points(MeshType const &mesh)
Definition: filters.hpp:1296
Lagrange polynomial set.
Definition: lagrange.hpp:284
Convex base class.
Definition: convex.hpp:43
reference_convex_type const & referenceConvex() const
Definition: lagrange.hpp:382
Definition: feelpoly/policy.hpp:59
primal_space_type::polyset_type polyset_type
Definition: fe.hpp:84
Finite element following Ciarlet framework.
Definition: fe.hpp:60
Definition: feelpoly/policy.hpp:242