30 #ifndef __RaviartThomas_H
31 #define __RaviartThomas_H 1
33 #include <boost/ptr_container/ptr_vector.hpp>
34 #include <boost/assign/std/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 <boost/assign/list_of.hpp>
42 #include <boost/assign/std/vector.hpp>
44 #include <feel/feelcore/feel.hpp>
45 #include <feel/feelcore/traits.hpp>
58 #include <feel/feelpoly/quadpoint.hpp>
70 typedef typename P::value_type value_type;
71 typedef typename P::points_type points_type;
72 times_x ( P
const& p,
int c )
79 typename ublas::vector<value_type> operator() ( points_type
const& __pts )
const
82 std::cout <<
"times_x(pts) : " << __pts << std::endl;
83 std::cout <<
"times_x(pts) : " << M_p.evaluate( __pts ) << std::endl;
84 std::cout <<
"times_x(coeff) : " << M_p.coefficients() << std::endl;
87 return ublas::element_prod( ublas::row( __pts, M_c ),
88 ublas::row( M_p.evaluate( __pts ), 0 ) );
96 struct extract_all_poly_indices
99 extract_all_poly_indices( T comp, T N )
111 template<uint16_type N,
114 template<u
int16_type, u
int16_type, u
int16_type>
class Convex = Simplex>
115 class RaviartThomasPolynomialSet
117 public Feel::detail::OrthonormalPolynomialSet<N, N, O+1, Vectorial, T, Convex>
119 typedef Feel::detail::OrthonormalPolynomialSet<N, N, O+1, Vectorial, T, Convex> super;
122 static const uint16_type Om1 = (O==0)?0:O-1;
123 typedef Feel::detail::OrthonormalPolynomialSet<N, N, O, Vectorial, T, Convex> Pk_v_type;
124 typedef Feel::detail::OrthonormalPolynomialSet<N, N, O+1, Vectorial, T, Convex> Pkp1_v_type;
125 typedef Feel::detail::OrthonormalPolynomialSet<N, N, Om1, Vectorial, T, Convex> Pkm1_v_type;
126 typedef Feel::detail::OrthonormalPolynomialSet<N, N, O, Scalar, T, Convex> Pk_s_type;
127 typedef Feel::detail::OrthonormalPolynomialSet<N, N, O+1, Scalar, T, Convex> Pkp1_s_type;
129 typedef PolynomialSet<typename super::basis_type,Vectorial> vectorial_polynomialset_type;
130 typedef typename vectorial_polynomialset_type::polynomial_type vectorial_polynomial_type;
131 typedef PolynomialSet<typename super::basis_type,Scalar> scalar_polynomialset_type;
132 typedef typename scalar_polynomialset_type::polynomial_type scalar_polynomial_type;
134 typedef RaviartThomasPolynomialSet<N, O, T> self_type;
136 typedef typename super::value_type value_type;
137 typedef typename super::convex_type convex_type;
138 typedef typename super::matrix_type matrix_type;
139 typedef typename super::points_type points_type;
141 static const uint16_type nDim = super::nDim;
142 static const uint16_type nOrder = super::nOrder;
143 static const uint16_type nComponents = super::nComponents;
144 static const bool is_product =
false;
145 RaviartThomasPolynomialSet()
149 uint16_type dim_Pkp1 = convex_type::polyDims( nOrder );
150 uint16_type dim_Pk = convex_type::polyDims( nOrder-1 );
151 uint16_type dim_Pkm1 = ( nOrder==1 )?0:convex_type::polyDims( nOrder-2 );
153 std::cout <<
"[RTPset] dim_Pkp1 = " << dim_Pkp1 <<
"\n";
154 std::cout <<
"[RTPset] dim_Pk = " << dim_Pk <<
"\n";
155 std::cout <<
"[RTPset] dim_Pkm1 = " << dim_Pkm1 <<
"\n";
159 vectorial_polynomialset_type Pk_v( Pkp1_v.polynomialsUpToDimension( dim_Pk ) );
161 std::cout <<
"[RTPset] Pk_v =" << Pk_v.coeff() <<
"\n";
165 scalar_polynomialset_type Pk ( Pkp1.polynomialsUpToDimension( dim_Pk ) );
167 std::cout <<
"[RTPset] Pk =" << Pk.coeff() <<
"\n";
168 std::cout <<
"[RTPset] Pk(0) =" << Pk.polynomial( 0 ).coefficients() <<
"\n";
172 IMGeneral<convex_type::nDim, 2*nOrder,value_type> im;
174 ublas::matrix<value_type> xPkc( nComponents*( dim_Pk-dim_Pkm1 ),Pk.coeff().size2() );
177 for (
int l = dim_Pkm1, i = 0; l < dim_Pk; ++l, ++i )
179 for (
int j = 0; j < convex_type::nDim; ++j )
181 Feel::detail::times_x<scalar_polynomial_type> xp( Pk.polynomial( l ), j );
182 ublas::row( xPkc,i*nComponents+j )=
191 vectorial_polynomialset_type xPk(
typename super::basis_type(), xPkc,
true );
195 this->setCoefficient(
unite( Pk_v, xPk ).coeff(),
true );
210 template<
typename Basis,
211 template<
class, u
int16_type,
class>
class PointSetType>
212 class RaviartThomasDual
214 public DualBasis<Basis>
216 typedef DualBasis<Basis> super;
219 static const uint16_type nDim = super::nDim;
220 static const uint16_type nOrder= super::nOrder;
222 typedef typename super::primal_space_type primal_space_type;
223 typedef typename primal_space_type::value_type value_type;
224 typedef typename primal_space_type::points_type points_type;
225 typedef typename primal_space_type::matrix_type matrix_type;
226 typedef typename primal_space_type::template convex<nDim+nOrder-1>::type convex_type;
227 typedef Reference<convex_type, nDim, nDim+nOrder-1, nDim, value_type> reference_convex_type;
228 typedef typename reference_convex_type::node_type node_type;
230 typedef typename primal_space_type::Pkp1_v_type Pkp1_v_type;
231 typedef typename primal_space_type::vectorial_polynomialset_type vectorial_polynomialset_type;
234 typedef PointSetType<convex_type, nOrder, value_type> pointset_type;
236 static const uint16_type nbPtsPerVertex = 0;
237 static const uint16_type nbPtsPerEdge = mpl::if_<mpl::equal_to<mpl::int_<nDim>,mpl::int_<2> >,mpl::int_<reference_convex_type::nbPtsPerEdge>,mpl::int_<0> >::type::value;
238 static const uint16_type nbPtsPerFace =mpl::if_<mpl::equal_to<mpl::int_<nDim>,mpl::int_<3> >,mpl::int_<reference_convex_type::nbPtsPerFace>,mpl::int_<0> >::type::value;
239 static const uint16_type nbPtsPerVolume = 0;
240 static const uint16_type numPoints = ( reference_convex_type::numGeometricFaces*nbPtsPerFace+reference_convex_type::numEdges*nbPtsPerEdge );
243 static const uint16_type nDofPerVertex = 0;
246 static const uint16_type nDofPerEdge = nbPtsPerEdge;
249 static const uint16_type nDofPerFace = nbPtsPerFace;
252 static const uint16_type nDofPerVolume = 0;
255 static const uint16_type nLocalDof = numPoints;
257 RaviartThomasDual( primal_space_type
const& primal )
261 M_eid( M_convex_ref.topologicalDimension()+1 ),
262 M_pts( nDim, numPoints ),
263 M_pts_per_face( convex_type::numTopologicalFaces ),
267 std::cout <<
"Raviart-Thomas finite element(dual): \n";
268 std::cout <<
" o- dim = " << nDim <<
"\n";
269 std::cout <<
" o- order = " << nOrder <<
"\n";
270 std::cout <<
" o- numPoints = " << numPoints <<
"\n";
271 std::cout <<
" o- nbPtsPerVertex = " << ( int )nbPtsPerVertex <<
"\n";
272 std::cout <<
" o- nbPtsPerEdge = " << ( int )nbPtsPerEdge <<
"\n";
273 std::cout <<
" o- nbPtsPerFace = " << ( int )nbPtsPerFace <<
"\n";
274 std::cout <<
" o- nbPtsPerVolume = " << ( int )nbPtsPerVolume <<
"\n";
275 std::cout <<
" o- nLocalDof = " << nLocalDof <<
"\n";
280 for (
int p = 0, e = M_convex_ref.entityRange( nDim-1 ).begin();
281 e < M_convex_ref.entityRange( nDim-1 ).end();
284 points_type Gt ( M_convex_ref.makePoints( nDim-1, e ) );
285 M_pts_per_face[e] = Gt ;
291 ublas::subrange( M_pts, 0, nDim, p, p+Gt.size2() ) = Gt;
300 typedef Functional<primal_space_type> functional_type;
301 std::vector<functional_type> fset;
305 std::vector<double> j;
308 using namespace boost::assign;
311 j += 2.8284271247461903,2.0,2.0;
314 j+= 3.464101615137754, 2, 2, 2;
320 for (
int e = M_convex_ref.entityRange( nDim-1 ).begin();
321 e < M_convex_ref.entityRange( nDim-1 ).end();
325 node_type dir = M_convex_ref.normal( e )*j[e];
327 dcpe_type __dcpe( primal, dir, M_pts_per_face[e] );
328 std::copy( __dcpe.begin(), __dcpe.end(), std::back_inserter( fset ) );
337 uint16_type dim_Pkp1 = convex_type::polyDims( nOrder );
338 uint16_type dim_Pk = convex_type::polyDims( nOrder-1 );
339 uint16_type dim_Pm1 = convex_type::polyDims( nOrder-2 );
343 vectorial_polynomialset_type Pkm1 ( Pkp1.polynomialsUpToDimension( dim_Pm1 ) );
347 for (
int i = 0; i < Pkm1.polynomialDimension(); ++i )
349 typedef functional::IntegralMoment<primal_space_type, vectorial_polynomialset_type> fim_type;
352 fset.push_back( fim_type( primal, Pkm1.polynomial( i ) ) );
357 M_fset.setFunctionalSet( fset );
368 points_type
const&
points()
const
374 matrix_type operator()( primal_space_type
const& pset )
const
377 return M_fset( pset );
380 points_type
const&
points( uint16_type f )
const
382 return M_pts_per_face[f];
384 ublas::matrix_column<points_type const> point( uint16_type f, uint32_type __i )
const
386 return ublas::column( M_pts_per_face[f], __i );
388 ublas::matrix_column<points_type> point( uint16_type f, uint32_type __i )
390 return ublas::column( M_pts_per_face[f], __i );
397 void setPoints( uint16_type f, points_type
const& n )
399 M_pts_per_face[f].resize( n.size1(), n.size2(), false );
400 M_pts_per_face[f] = n;
404 reference_convex_type M_convex_ref;
405 std::vector<std::vector<uint16_type> > M_eid;
407 std::vector<points_type> M_pts_per_face;
408 FunctionalSet<primal_space_type> M_fset;
423 template<uint16_type N,
426 template<u
int16_type, u
int16_type, u
int16_type>
class Convex = Simplex,
427 uint16_type TheTAG=0 >
430 public FiniteElement<RaviartThomasPolynomialSet<N, O, T, Convex>,
431 fem::detail::RaviartThomasDual,
432 PointSetEquiSpaced >,
433 public boost::enable_shared_from_this<RaviartThomas<N,O,T,Convex> >
436 fem::detail::RaviartThomasDual,
437 PointSetEquiSpaced >
super;
440 BOOST_STATIC_ASSERT( N > 1 );
446 static const uint16_type nDim = N;
448 static const bool isTransformationEquivalent =
true;
449 static const bool isContinuous =
true;
451 static const uint16_type TAG = TheTAG;
454 typedef typename super::value_type value_type;
455 typedef typename super::primal_space_type primal_space_type;
456 typedef typename super::dual_space_type dual_space_type;
462 static const bool is_vectorial = polyset_type::is_vectorial;
463 static const bool is_scalar = polyset_type::is_scalar;
464 static const uint16_type nComponents = polyset_type::nComponents;
465 static const bool is_product =
false;
468 typedef typename dual_space_type::convex_type convex_type;
469 typedef typename dual_space_type::pointset_type pointset_type;
470 typedef typename dual_space_type::reference_convex_type reference_convex_type;
471 typedef typename reference_convex_type::node_type node_type;
472 typedef typename reference_convex_type::points_type points_type;
475 static const uint16_type nOrder = dual_space_type::nOrder;
476 static const uint16_type nbPtsPerVertex = 0;
477 static const uint16_type nbPtsPerEdge = mpl::if_<mpl::equal_to<mpl::int_<nDim>,mpl::int_<2> >,
478 mpl::int_<reference_convex_type::nbPtsPerEdge>,
479 mpl::int_<0> >::type::value;
480 static const uint16_type nbPtsPerFace = mpl::if_<mpl::equal_to<mpl::int_<nDim>,mpl::int_<3> >,
481 mpl::int_<reference_convex_type::nbPtsPerFace>,
482 mpl::int_<0> >::type::value;
483 static const uint16_type nbPtsPerVolume = 0;
484 static const uint16_type numPoints = ( reference_convex_type::numGeometricFaces*nbPtsPerFace+
485 reference_convex_type::numEdges*nbPtsPerEdge );
487 static const uint16_type nLocalDof = dual_space_type::nLocalDof;
496 super( dual_space_type( primal_space_type() ) ),
500 std::cout <<
"[RT] nPtsPerEdge = " << nbPtsPerEdge <<
"\n";
501 std::cout <<
"[RT] nPtsPerFace = " << nbPtsPerFace <<
"\n";
502 std::cout <<
"[RT] numPoints = " << numPoints <<
"\n";
504 std::cout <<
"[RT] nDof = " << super::nDof <<
"\n";
506 std::cout <<
"[RT] coeff : " << this->coeff() <<
"\n";
507 std::cout <<
"[RT] pts : " << this->
points() <<
"\n";
508 std::cout <<
"[RT] eval at pts : " << this->evaluate( this->
points() ) <<
"\n";
509 std::cout <<
"[RT] is_product : " << is_product <<
"\n";
543 return "raviartthomas";
559 template<
typename ExprType>
561 isomorphism( ExprType expr ) -> decltype( Feel::vf::detJ()*( trans( Feel::vf::JinvT() )*expr )*Feel::vf::Nref() )
564 using namespace Feel::vf;
565 return detJ()*( trans( JinvT() )*expr )*Nref();
573 template<
typename ExprType,
typename ContextType>
574 std::vector<value_type>
575 interpolate( boost::shared_ptr<ContextType>& ctx, ExprType & expr )
577 using namespace Feel::vf;
578 typedef boost::shared_ptr<ContextType> gmc_ptrtype;
579 typedef fusion::map<fusion::pair<vf::detail::gmc<0>, gmc_ptrtype> > map_gmc_type;
581 std::vector<value_type> v( nLocalDof );
584 for (
int face = 0; face < numTopologicalFaces; ++face )
587 ctx->update( _face=face, _element=ctx->id() );
589 map_gmc_type mapgmc( fusion::make_pair<vf::detail::gmc<0> >( ctx ) );
590 expr.update( mapgmc, face );
592 for (
int q = 0; q < nDofPerFace; ++q )
594 int ldof = nDofPerFace*face+i;
595 v[ldof] = expr.evalq( 0,0,i );
606 template<
typename GMContext,
typename PC,
typename Phi,
typename GPhi,
typename HPhi >
607 static void transform( boost::shared_ptr<GMContext> gmc, boost::shared_ptr<PC>
const& pc,
609 GPhi& g_phi_t,
const bool do_gradient,
610 HPhi& h_phi_t,
const bool do_hessian
614 transform ( *gmc, *pc, phi_t, g_phi_t, do_gradient, h_phi_t, do_hessian );
616 template<
typename GMContext,
typename PC,
typename Phi,
typename GPhi,
typename HPhi >
617 static void transform( GMContext
const& gmc,
620 GPhi& g_phi_t,
const bool do_gradient,
621 HPhi& h_phi_t,
const bool do_hessian
626 typename GMContext::gm_type::matrix_type
const B = gmc.B( 0 );
627 typename GMContext::gm_type::matrix_type
const K = gmc.K( 0 );
628 typename GMContext::gm_type::matrix_type JB( K/gmc.J( 0 ) );
630 std::cout <<
"K= " << gmc.K( 0 ) <<
"\n";
631 std::cout <<
"B= " << B <<
"\n";
632 std::cout <<
"J= " << gmc.J( 0 ) <<
"\n";
633 std::cout <<
"JB= " << JB <<
"\n";
635 std::fill( phi_t.data(), phi_t.data()+phi_t.num_elements(), value_type( 0 ) );
640 std::fill( g_phi_t.data(), g_phi_t.data()+g_phi_t.num_elements(), value_type( 0 ) );
644 std::fill( h_phi_t.data(), h_phi_t.data()+h_phi_t.num_elements(), value_type( 0 ) );
646 const uint16_type Q = gmc.nPoints();
649 for ( uint16_type i = 0; i < nLocalDof; ++i )
651 for ( uint16_type l = 0; l < nDim; ++l )
654 for ( uint16_type p = 0; p < nDim; ++p )
656 for ( uint16_type q = 0; q < Q; ++q )
662 phi_t[i][l][0][q] += JB( l, p ) * pc.phi( i,p,0,q );
673 for ( uint16_type p = 0; p < nDim; ++p )
675 for ( uint16_type r = 0; r < nDim; ++r )
677 for ( uint16_type s = 0; s < Q; ++s )
679 g_phi_t[i][p][r][s] = 0;
681 for ( uint16_type q = 0; q < nDim; ++q )
683 g_phi_t[i][p][r][s] += JB( p, q ) * pc.grad( i,q,r,s );
708 reference_convex_type M_refconvex;
712 template<uint16_type N,
715 template<u
int16_type, u
int16_type, u
int16_type>
class Convex,
717 const uint16_type RaviartThomas<N,O,T,Convex,TheTAG>::nDim;
718 template<uint16_type N,
721 template<u
int16_type, u
int16_type, u
int16_type>
class Convex,
723 const uint16_type RaviartThomas<N,O,T,Convex,TheTAG>::nOrder;
724 template<uint16_type N,
727 template<u
int16_type, u
int16_type, u
int16_type>
class Convex,
729 const uint16_type RaviartThomas<N,O,T,Convex,TheTAG>::nLocalDof;
732 template<uint16_type Order,
733 uint16_type TheTAG=0 >
737 template<uint16_type N,
740 typename Convex = Simplex<N> >
743 typedef typename mpl::if_<mpl::bool_<Convex::is_simplex>,
744 mpl::identity<fem::RaviartThomas<N,Order,T,Simplex,TheTAG> >,
745 mpl::identity<fem::RaviartThomas<N,Order,T,Hypercube,TheTAG> > >::type::type result_type;
746 typedef result_type type;
749 template<u
int16_type TheNewTAG>
752 typedef RaviartThomas<Order,TheNewTAG> type;
755 typedef Lagrange<Order,Scalar> component_basis_type;
757 static const uint16_type TAG = TheTAG;
functional associate with directional component point evaluation
Definition: functionals.hpp:405
reference_convex_type const & referenceConvex() const
Definition: raviartthomas.hpp:536
Polynomial< Pset, Scalar > project(Pset const &pset, Func const &f, IM const &im)
Definition: operations.hpp:436
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
describe continuous functions
Definition: continuity.hpp:51
PolynomialSet< P, Type > unite(Poly1< P, Type > const &pset1, Poly2< P, Type > const &pset2)
union of two polynomial sets P1 and P2 :
Definition: operations.hpp:461
super::polyset_type polyset_type
Definition: raviartthomas.hpp:461
primal_space_type::polyset_type polyset_type
Definition: fe.hpp:84
Finite element following Ciarlet framework.
Definition: fe.hpp:60
RaviartThomas Finite Element.
Definition: raviartthomas.hpp:428
void interpolate(boost::shared_ptr< SpaceType > const &space, FunctionType const &f, typename SpaceType::element_type &interp, int same_mesh=INTERPOLATE_DIFFERENT_MESH)
Definition: interpolate.hpp:51
std::string familyName() const
Definition: raviartthomas.hpp:541