30 #ifndef __CrouzeixRaviart_H
31 #define __CrouzeixRaviart_H 1
33 #include <boost/ptr_container/ptr_vector.hpp>
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>
41 #include <feel/feelcore/feel.hpp>
42 #include <feel/feelcore/traits.hpp>
60 template<uint16_type N,
61 template<u
int16_type Dim>
class PolySetType = Scalar,
63 class RannacherTurekPolynomialSet
65 public MomentPolynomialSet<N, 2, N, PolySetType, T, Hypercube>
67 typedef MomentPolynomialSet<N, 2, N, PolySetType, T, Hypercube> super;
71 typedef RannacherTurekPolynomialSet<N, PolySetType, T> self_type;
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;
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;
83 RannacherTurekPolynomialSet()
87 Moment<2,2,Hypercube<2> > m;
89 for (
int c = 0; c < nComponents; ++c )
92 this->insert( m.template pick<PolySetType>( 0, c ).toSet(
true ), c==0 );
94 this->insert( m.template pick<PolySetType>( 1, c ).toSet(
true ) );
96 this->insert( m.template pick<PolySetType>( 3, c ).toSet(
true ) );
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 );
102 this->insert( p.toSet(
true ) );
110 template<
typename Basis,
template<
class, u
int16_type,
class>
class PointSetType>
111 class CrouzeixRaviartDual
113 public DualBasis<Basis>
115 typedef DualBasis<Basis> super;
118 static const uint16_type nDim = super::nDim;
119 static const uint16_type nOrder= super::nOrder;
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;
130 typedef PointSet<convex_type, value_type> pointset_type;
131 typedef PointSetType<convex_type, 2, value_type> equispaced_pointset_type;
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;
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 );
151 static const uint16_type nDofPerVertex = nbPtsPerVertex;
154 static const uint16_type nDofPerEdge = nbPtsPerEdge;
157 static const uint16_type nDofPerFace = nbPtsPerFace;
160 static const uint16_type nDofPerVolume = nbPtsPerVolume;
163 static const uint16_type nLocalDof = numPoints;
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> >,
169 mpl::int_<nFaces> >::type >::type::value;
172 CrouzeixRaviartDual( primal_space_type
const& primal )
176 M_eid( M_convex_ref.topologicalDimension()+1 ),
177 M_pts( nDim, numPoints ),
178 M_points_face( nFacesInConvex ),
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";
191 equispaced_pointset_type epts;
193 int d = M_convex_ref.topologicalDimension()-1;
198 for (
int e = M_convex_ref.entityRange( d ).begin();
199 e < M_convex_ref.entityRange( d ).end();
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();
208 setFset( primal, M_pts, mpl::bool_<primal_space_type::is_scalar>() );
213 points_type
const&
points()
const
217 points_type
const&
points( uint16_type f )
const
219 return M_points_face[f];
223 matrix_type operator()( primal_space_type
const& pset )
const
225 return M_fset( pset );
229 void setFset( primal_space_type
const& primal, points_type
const& __pts, mpl::bool_<true> )
231 M_fset.setFunctionalSet( functional::PointsEvaluation<primal_space_type>( primal,
235 void setFset( primal_space_type
const& primal, points_type
const& __pts, mpl::bool_<false> )
237 M_fset.setFunctionalSet( functional::ComponentsPointsEvaluation<primal_space_type>( primal,
243 reference_convex_type M_convex_ref;
244 std::vector<std::vector<uint16_type> > M_eid;
246 std::vector<points_type> M_points_face;
247 FunctionalSet<primal_space_type> M_fset;
260 template<uint16_type N,
262 template<u
int16_type Dim>
class PolySetType,
264 template<u
int16_type, u
int16_type, u
int16_type>
class Convex = Simplex,
265 uint16_type TheTAG=0 >
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,
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;
281 BOOST_STATIC_ASSERT( N > 1 );
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;
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;
296 static const uint16_type TAG = TheTAG;
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;
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;
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 );
334 super( dual_space_type( primal_space_type() ) ),
384 return "CrouzeixRaviart";
387 template<
typename ExprType>
389 isomorphism( ExprType& expr ) -> decltype( expr )
399 reference_convex_type M_refconvex;
403 template<uint16_type N,
405 template<u
int16_type Dim>
class PolySetType,
407 template<u
int16_type, u
int16_type, u
int16_type>
class Convex,
409 const uint16_type CrouzeixRaviart<N,RealDim,PolySetType,T,Convex,TheTAG>::nDim;
410 template<uint16_type N,
412 template<u
int16_type Dim>
class PolySetType,
414 template<u
int16_type, u
int16_type, u
int16_type>
class Convex,
416 const uint16_type CrouzeixRaviart<N,RealDim,PolySetType,T,Convex,TheTAG>::nOrder;
420 template<uint16_type Order,
421 template<u
int16_type Dim>
class PolySetType = Scalar,
422 template<
class, u
int16_type,
class>
class Pts = PointSetEquiSpaced,
423 uint16_type TheTAG=0 >
424 class CrouzeixRaviart
427 template<uint16_type N,
430 typename Convex = Simplex<N> >
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;
439 template<u
int16_type TheNewTAG>
442 typedef CrouzeixRaviart<Order,PolySetType,Pts,TheNewTAG> type;
445 typedef CrouzeixRaviart<Order,Scalar,Pts> component_basis_type;
447 static const uint16_type nOrder = Order;
448 static const uint16_type TAG = TheTAG;
452 template<uint16_type Order,
453 template<u
int16_type Dim>
class PolySetType,
454 template<
class, u
int16_type,
class>
class Pts,
456 const uint16_type CrouzeixRaviart<Order,PolySetType,Pts,TheTAG>::nOrder;