32 #include <boost/ptr_container/ptr_vector.hpp>
33 #include <boost/assign/std/vector.hpp>
34 #include <boost/numeric/ublas/vector.hpp>
35 #include <boost/numeric/ublas/io.hpp>
36 #include <boost/numeric/ublas/matrix.hpp>
37 #include <boost/numeric/ublas/matrix_proxy.hpp>
38 #include <boost/numeric/ublas/lu.hpp>
40 #include <boost/assign/list_of.hpp>
41 #include <boost/assign/std/vector.hpp>
43 #include <feel/feelcore/feel.hpp>
44 #include <feel/feelcore/traits.hpp>
57 #include <feel/feelpoly/quadpoint.hpp>
69 typedef typename P::value_type value_type;
70 typedef typename P::points_type points_type;
71 times_rotx ( P
const& p,
int c )
78 typename ublas::vector<value_type> operator() ( points_type
const& __pts )
const
81 std::cout <<
"times_x(pts) : " << __pts << std::endl;
82 std::cout <<
"times_x(pts) : " << M_p.evaluate( __pts ) << std::endl;
83 std::cout <<
"times_x(coeff) : " << M_p.coefficients() << std::endl;
88 return ublas::element_prod( ublas::row( __pts, 1 ),
89 ublas::row( M_p.evaluate( __pts ), 0 ) );
92 return ublas::element_prod( -ublas::row( __pts, 0 ),
93 ublas::row( M_p.evaluate( __pts ), 0 ) );
101 struct extract_all_poly_indices
104 extract_all_poly_indices( T comp, T N )
116 template<uint16_type N,
119 template<u
int16_type, u
int16_type, u
int16_type>
class Convex = Simplex>
120 class NedelecPolynomialSet
122 public Feel::detail::OrthonormalPolynomialSet<N, O+1, N, Vectorial, T, Convex>
124 typedef Feel::detail::OrthonormalPolynomialSet<N, O+1, N, Vectorial, T, Convex> super;
127 static const uint16_type Om1 = (O==0)?0:O-1;
128 typedef Feel::detail::OrthonormalPolynomialSet<N, O, N, Vectorial, T, Convex> Pk_v_type;
129 typedef Feel::detail::OrthonormalPolynomialSet<N, O+1, N, Vectorial, T, Convex> Pkp1_v_type;
130 typedef Feel::detail::OrthonormalPolynomialSet<N, Om1, N, Vectorial, T, Convex> Pkm1_v_type;
131 typedef Feel::detail::OrthonormalPolynomialSet<N, O, N, Scalar, T, Convex> Pk_s_type;
132 typedef Feel::detail::OrthonormalPolynomialSet<N, O+1, N, Scalar, T, Convex> Pkp1_s_type;
134 typedef PolynomialSet<typename super::basis_type,Vectorial> vectorial_polynomialset_type;
135 typedef typename vectorial_polynomialset_type::polynomial_type vectorial_polynomial_type;
136 typedef PolynomialSet<typename super::basis_type,Scalar> scalar_polynomialset_type;
137 typedef typename scalar_polynomialset_type::polynomial_type scalar_polynomial_type;
139 typedef NedelecPolynomialSet<N, O, T> self_type;
141 typedef typename super::value_type value_type;
142 typedef typename super::convex_type convex_type;
143 typedef typename super::matrix_type matrix_type;
144 typedef typename super::points_type points_type;
146 static const uint16_type nDim = super::nDim;
147 static const uint16_type nOrder = super::nOrder;
148 static const uint16_type nComponents = super::nComponents;
149 static const bool is_product =
false;
150 NedelecPolynomialSet()
154 std::cout <<
"[NPset] nOrder = " << nOrder <<
"\n";
155 std::cout <<
"[NPset] O = " << O <<
"\n";
156 uint16_type dim_Pkp1 = convex_type::polyDims( nOrder );
157 uint16_type dim_Pk = convex_type::polyDims( nOrder-1 );
158 uint16_type dim_Pkm1 = ( nOrder==1 )?0:convex_type::polyDims( nOrder-2 );
160 std::cout <<
"[NPset] dim_Pkp1 = " << dim_Pkp1 <<
"\n";
161 std::cout <<
"[NPset] dim_Pk = " << dim_Pk <<
"\n";
162 std::cout <<
"[NPset] dim_Pkm1 = " << dim_Pkm1 <<
"\n";
166 vectorial_polynomialset_type Pk_v( Pkp1_v.polynomialsUpToDimension( dim_Pk ) );
168 std::cout <<
"[NPset] Pk_v =" << Pk_v.coeff() <<
"\n";
172 scalar_polynomialset_type Pk ( Pkp1.polynomialsUpToDimension( dim_Pk ) );
174 std::cout <<
"[NPset] Pk =" << Pk.coeff() <<
"\n";
175 std::cout <<
"[NPset] Pk(0) =" << Pk.polynomial( 0 ).coefficients() <<
"\n";
179 IMGeneral<convex_type::nDim, 2*nOrder,value_type> im;
181 ublas::matrix<value_type> xPkc( nComponents*( dim_Pk-dim_Pkm1 ),Pk.coeff().size2() );
184 for (
int l = dim_Pkm1, i = 0; l < dim_Pk; ++l, ++i )
186 for (
int j = 0; j < convex_type::nDim; ++j )
188 Feel::detail::times_rotx<scalar_polynomial_type> xp( Pk.polynomial( l ), j );
189 ublas::row( xPkc,i*nComponents+j )=
198 vectorial_polynomialset_type xPk(
typename super::basis_type(), xPkc,
true );
202 this->setCoefficient(
unite( Pk_v, xPk ).coeff(),
true );
217 template<
typename Basis,
218 template<
class, u
int16_type,
class>
class PointSetType>
221 public DualBasis<Basis>
223 typedef DualBasis<Basis> super;
226 static const uint16_type nDim = super::nDim;
227 static const uint16_type nOrder= super::nOrder;
229 typedef typename super::primal_space_type primal_space_type;
230 typedef typename primal_space_type::value_type value_type;
231 typedef typename primal_space_type::points_type points_type;
232 typedef typename primal_space_type::matrix_type matrix_type;
233 typedef typename primal_space_type::template convex<nDim+nOrder-1>::type convex_type;
234 typedef Reference<convex_type, nDim, nDim+nOrder-1, nDim, value_type> reference_convex_type;
235 typedef typename reference_convex_type::node_type node_type;
237 typedef typename primal_space_type::Pkp1_v_type Pkp1_v_type;
238 typedef typename primal_space_type::vectorial_polynomialset_type vectorial_polynomialset_type;
241 typedef PointSetType<convex_type, nOrder, value_type> pointset_type;
243 static const uint16_type nbPtsPerVertex = 0;
244 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;
245 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;
246 static const uint16_type nbPtsPerVolume = 0;
247 static const uint16_type numPoints = ( reference_convex_type::numGeometricFaces*nbPtsPerFace+reference_convex_type::numEdges*nbPtsPerEdge );
250 static const uint16_type nDofPerVertex = 0;
253 static const uint16_type nDofPerEdge = nbPtsPerEdge;
256 static const uint16_type nDofPerFace = nbPtsPerFace;
259 static const uint16_type nDofPerVolume = 0;
262 static const uint16_type nLocalDof = numPoints;
264 NedelecDual( primal_space_type
const& primal )
268 M_eid( M_convex_ref.topologicalDimension()+1 ),
269 M_pts( nDim, numPoints ),
270 M_pts_per_face( convex_type::numTopologicalFaces ),
274 std::cout <<
"Nedelec finite element(dual): \n";
275 std::cout <<
" o- dim = " << nDim <<
"\n";
276 std::cout <<
" o- order = " << nOrder <<
"\n";
277 std::cout <<
" o- numPoints = " << numPoints <<
"\n";
278 std::cout <<
" o- nbPtsPerVertex = " << ( int )nbPtsPerVertex <<
"\n";
279 std::cout <<
" o- nbPtsPerEdge = " << ( int )nbPtsPerEdge <<
"\n";
280 std::cout <<
" o- nbPtsPerFace = " << ( int )nbPtsPerFace <<
"\n";
281 std::cout <<
" o- nbPtsPerVolume = " << ( int )nbPtsPerVolume <<
"\n";
282 std::cout <<
" o- nLocalDof = " << nLocalDof <<
"\n";
287 for (
int p = 0, e = M_convex_ref.entityRange( nDim-1 ).begin();
288 e < M_convex_ref.entityRange( nDim-1 ).end();
291 points_type Gt ( M_convex_ref.makePoints( nDim-1, e ) );
292 M_pts_per_face[e] = Gt ;
298 ublas::subrange( M_pts, 0, nDim, p, p+Gt.size2() ) = Gt;
307 typedef Functional<primal_space_type> functional_type;
308 std::vector<functional_type> fset;
312 std::vector<double> j;
315 using namespace boost::assign;
318 j += 2.8284271247461903,2.0,2.0;
321 j+= 3.464101615137754, 2, 2, 2;
327 for (
int e = M_convex_ref.entityRange( nDim-1 ).begin();
328 e < M_convex_ref.entityRange( nDim-1 ).end();
331 typedef Feel::functional::DirectionalComponentPointsEvaluation<primal_space_type> dcpe_type;
332 std::cout <<
"tangent " << e <<
":" << M_convex_ref.tangent( e ) <<
"\n";
333 node_type dir= M_convex_ref.tangent( e )*j[e];
337 dcpe_type __dcpe( primal, dir, M_pts_per_face[e] );
338 std::copy( __dcpe.begin(), __dcpe.end(), std::back_inserter( fset ) );
347 uint16_type dim_Pkp1 = convex_type::polyDims( nOrder );
348 uint16_type dim_Pk = convex_type::polyDims( nOrder-1 );
349 uint16_type dim_Pm1 = convex_type::polyDims( nOrder-2 );
353 vectorial_polynomialset_type Pkm1 ( Pkp1.polynomialsUpToDimension( dim_Pm1 ) );
357 for (
int i = 0; i < Pkm1.polynomialDimension(); ++i )
359 typedef functional::IntegralMoment<primal_space_type, vectorial_polynomialset_type> fim_type;
362 fset.push_back( fim_type( primal, Pkm1.polynomial( i ) ) );
367 M_fset.setFunctionalSet( fset );
378 points_type
const&
points()
const
384 matrix_type operator()( primal_space_type
const& pset )
const
387 return M_fset( pset );
390 points_type
const&
points( uint16_type f )
const
392 return M_pts_per_face[f];
394 ublas::matrix_column<points_type const> point( uint16_type f, uint32_type __i )
const
396 return ublas::column( M_pts_per_face[f], __i );
398 ublas::matrix_column<points_type> point( uint16_type f, uint32_type __i )
400 return ublas::column( M_pts_per_face[f], __i );
407 void setPoints( uint16_type f, points_type
const& n )
409 M_pts_per_face[f].resize( n.size1(), n.size2(), false );
410 M_pts_per_face[f] = n;
414 reference_convex_type M_convex_ref;
415 std::vector<std::vector<uint16_type> > M_eid;
417 std::vector<points_type> M_pts_per_face;
418 FunctionalSet<primal_space_type> M_fset;
433 template<uint16_type N,
436 template<u
int16_type, u
int16_type, u
int16_type>
class Convex = Simplex,
437 uint16_type TheTAG = 0 >
441 fem::detail::NedelecDual,
442 PointSetEquiSpaced >,
443 public boost::enable_shared_from_this<Nedelec<N,O,T,Convex> >
446 fem::detail::NedelecDual,
447 PointSetEquiSpaced >
super;
450 BOOST_STATIC_ASSERT( N > 1 );
456 static const uint16_type nDim = N;
458 static const bool isTransformationEquivalent =
true;
459 static const bool isContinuous =
true;
461 static const uint16_type TAG = TheTAG;
464 typedef typename super::value_type value_type;
465 typedef typename super::primal_space_type primal_space_type;
466 typedef typename super::dual_space_type dual_space_type;
472 static const bool is_vectorial = polyset_type::is_vectorial;
473 static const bool is_scalar = polyset_type::is_scalar;
474 static const uint16_type nComponents = polyset_type::nComponents;
475 static const bool is_product =
false;
478 typedef typename dual_space_type::convex_type convex_type;
479 typedef typename dual_space_type::pointset_type pointset_type;
480 typedef typename dual_space_type::reference_convex_type reference_convex_type;
481 typedef typename reference_convex_type::node_type node_type;
482 typedef typename reference_convex_type::points_type points_type;
485 static const uint16_type nOrder = dual_space_type::nOrder;
486 static const uint16_type nbPtsPerVertex = 0;
487 static const uint16_type nbPtsPerEdge = mpl::if_<mpl::equal_to<mpl::int_<nDim>,mpl::int_<2> >,
488 mpl::int_<reference_convex_type::nbPtsPerEdge>,
489 mpl::int_<0> >::type::value;
490 static const uint16_type nbPtsPerFace = mpl::if_<mpl::equal_to<mpl::int_<nDim>,mpl::int_<3> >,
491 mpl::int_<reference_convex_type::nbPtsPerFace>,
492 mpl::int_<0> >::type::value;
493 static const uint16_type nbPtsPerVolume = 0;
494 static const uint16_type numPoints = ( reference_convex_type::numGeometricFaces*nbPtsPerFace+
495 reference_convex_type::numEdges*nbPtsPerEdge );
497 static const uint16_type nLocalDof = dual_space_type::nLocalDof;
506 super( dual_space_type( primal_space_type() ) ),
510 std::cout <<
"[N] nPtsPerEdge = " << nbPtsPerEdge <<
"\n";
511 std::cout <<
"[N] nPtsPerFace = " << nbPtsPerFace <<
"\n";
512 std::cout <<
"[N] numPoints = " << numPoints <<
"\n";
514 std::cout <<
"[N] nDof = " << super::nDof <<
"\n";
516 std::cout <<
"[N] coeff : " << this->coeff() <<
"\n";
517 std::cout <<
"[N] pts : " << this->
points() <<
"\n";
518 std::cout <<
"[N] eval at pts : " << this->evaluate( this->
points() ) <<
"\n";
519 std::cout <<
"[N] is_product : " << is_product <<
"\n";
569 template<
typename ExprType>
571 isomorphism( ExprType expr ) -> decltype( Feel::vf::detJ()*( trans( Feel::vf::JinvT() )*expr )*Feel::vf::Nref() )
574 using namespace Feel::vf;
575 return detJ()*( trans( JinvT() )*expr )*Nref();
583 template<
typename ExprType,
typename ContextType>
584 std::vector<value_type>
585 interpolate( boost::shared_ptr<ContextType>& ctx, ExprType & expr )
587 using namespace Feel::vf;
588 typedef boost::shared_ptr<ContextType> gmc_ptrtype;
589 typedef fusion::map<fusion::pair<Feel::detail::gmc<0>, gmc_ptrtype> > map_gmc_type;
591 std::vector<value_type> v( nLocalDof );
594 for (
int face = 0; face < numTopologicalFaces; ++face )
597 ctx->update( _face=face, _element=ctx->id() );
599 map_gmc_type mapgmc( fusion::make_pair<Feel::detail::gmc<0> >( ctx ) );
600 expr.update( mapgmc, face );
602 for (
int q = 0; q < nDofPerFace; ++q )
604 int ldof = nDofPerFace*face+i;
605 v[ldof] = expr.evalq( 0,0,i );
614 template<
typename GMContext,
typename PC,
typename Phi,
typename GPhi,
typename HPhi >
615 static void transform( boost::shared_ptr<GMContext> gmc, boost::shared_ptr<PC>
const& pc,
617 GPhi& g_phi_t,
const bool do_gradient,
618 HPhi& h_phi_t,
const bool do_hessian
622 transform ( *gmc, *pc, phi_t, g_phi_t, do_gradient, h_phi_t, do_hessian );
624 template<
typename GMContext,
typename PC,
typename Phi,
typename GPhi,
typename HPhi >
625 static void transform( GMContext
const& gmc,
628 GPhi& g_phi_t,
const bool do_gradient,
629 HPhi& h_phi_t,
const bool do_hessian
634 typename GMContext::gm_type::matrix_type
const B = gmc.B( 0 );
635 typename GMContext::gm_type::matrix_type
const K = gmc.K( 0 );
636 typename GMContext::gm_type::matrix_type JB( K/gmc.J( 0 ) );
638 std::cout <<
"K= " << gmc.K( 0 ) <<
"\n";
639 std::cout <<
"B= " << B <<
"\n";
640 std::cout <<
"J= " << gmc.J( 0 ) <<
"\n";
641 std::cout <<
"JB= " << JB <<
"\n";
643 std::fill( phi_t.data(), phi_t.data()+phi_t.num_elements(), value_type( 0 ) );
648 std::fill( g_phi_t.data(), g_phi_t.data()+g_phi_t.num_elements(), value_type( 0 ) );
652 std::fill( h_phi_t.data(), h_phi_t.data()+h_phi_t.num_elements(), value_type( 0 ) );
654 const uint16_type Q = gmc.nPoints();
657 for ( uint16_type i = 0; i < nLocalDof; ++i )
659 for ( uint16_type l = 0; l < nDim; ++l )
662 for ( uint16_type p = 0; p < nDim; ++p )
664 for ( uint16_type q = 0; q < Q; ++q )
670 phi_t[i][l][0][q] += JB( l, p ) * pc.phi( i,p,0,q );
681 for ( uint16_type p = 0; p < nDim; ++p )
683 for ( uint16_type r = 0; r < nDim; ++r )
685 for ( uint16_type s = 0; s < Q; ++s )
687 g_phi_t[i][p][r][s] = 0;
689 for ( uint16_type q = 0; q < nDim; ++q )
691 g_phi_t[i][p][r][s] += JB( p, q ) * pc.grad( i,q,r,s );
716 reference_convex_type M_refconvex;
720 template<uint16_type N,
723 template<u
int16_type, u
int16_type, u
int16_type>
class Convex,
725 const uint16_type Nedelec<N,O,T,Convex,TheTAG>::nDim;
726 template<uint16_type N,
729 template<u
int16_type, u
int16_type, u
int16_type>
class Convex,
731 const uint16_type Nedelec<N,O,T,Convex,TheTAG>::nOrder;
732 template<uint16_type N,
735 template<u
int16_type, u
int16_type, u
int16_type>
class Convex,
737 const uint16_type Nedelec<N,O,T,Convex,TheTAG>::nLocalDof;
740 template<uint16_type Order,
741 uint16_type TheTAG=0>
745 template<uint16_type N,
748 typename Convex = Simplex<N> >
751 typedef typename mpl::if_<mpl::bool_<Convex::is_simplex>,
752 mpl::identity<fem::Nedelec<N,Order,T,Simplex,TheTAG> >,
753 mpl::identity<fem::Nedelec<N,Order,T,Hypercube,TheTAG> > >::type::type result_type;
754 typedef result_type type;
757 template<u
int16_type TheNewTAG>
760 typedef Nedelec<Order,TheNewTAG> type;
763 typedef Lagrange<Order,Scalar> component_basis_type;
765 static const uint16_type TAG = TheTAG;