29 #ifndef __FiniteElement_H
30 #define __FiniteElement_H 1
32 #include <feel/feelpoly/policy.hpp>
37 template<
typename Poly,
template<u
int16_type>
class PolySetType >
class PolynomialSet;
40 template<uint16_type Dim,
43 template<u
int16_type>
class PolySetType,
46 template<u
int16_type,u
int16_type,u
int16_type>
class Convex>
47 class OrthonormalPolynomialSet;
58 template<
class Pr,
template<
class,u
int16_type,
class>
class Pt>
class PDual,
59 template<class,uint16_type,class> class Pts>
61 public mpl::if_<mpl::bool_<P::is_scalar>,
62 mpl::identity<PolynomialSet<P, Scalar> >,
63 mpl::identity<PolynomialSet<P, Vectorial> > >::type::type
65 typedef typename mpl::if_<mpl::bool_<P::is_scalar>,
66 mpl::identity<PolynomialSet<P, Scalar> >,
67 mpl::identity<PolynomialSet<P, Vectorial> > >::type::type super;
77 typedef typename P::value_type value_type;
79 typedef P primal_space_type;
86 static const bool is_modal =
false;
88 typedef PDual<P, Pts> dual_space_type;
90 typedef typename super::matrix_type matrix_type;
91 typedef typename super::points_type points_type;
92 typedef typename super::self_type polynomialset_type;
94 typedef typename super::polynomial_type polynomial_type;
98 static const uint16_type nLocalDof = dual_space_type::nLocalDof;
100 static const uint16_type nDofPerVertex = dual_space_type::nDofPerVertex;
102 static const uint16_type nDofPerEdge = dual_space_type::nDofPerEdge;
104 static const uint16_type nDofPerFace = dual_space_type::nDofPerFace;
106 static const uint16_type nDofPerVolume = dual_space_type::nDofPerVolume;
108 static const uint16_type nDof = nLocalDof;
109 static const uint16_type nNodes = nDof;
110 static const uint16_type nDofGrad = super::nDim*nDof;
111 static const uint16_type nDofHess = super::nDim*super::nDim*nDof;
112 static const fem::transformation_type trans = ( fem::transformation_type )mpl::if_<mpl::and_<mpl::bool_<P::convex_type::is_simplex>,mpl::equal_to<mpl::int_<super::nOrder>, mpl::int_<1> > >,
113 mpl::int_<fem::LINEAR>,
114 mpl::int_<fem::NONLINEAR> >::type::value;
123 super( pdual.primalSpace() ),
125 M_primal( M_dual.primalSpace() )
127 DVLOG(2) <<
"============================================================\n";
128 DVLOG(2) <<
"New FE \n";
129 ublas::matrix<value_type> A( M_dual( M_primal ) );
132 ublas::matrix<value_type> D = ublas::identity_matrix<value_type>( A.size1(), A.size2() );
134 ublas::matrix<value_type> C = lu.solve( D );
137 DVLOG(2) <<
"is singular : " << lu.isNonsingular() <<
"\n"
138 <<
"det(A) = " << lu.det() <<
"\n";
141 if ( !lu.isNonsingular() )
143 std::cout <<
"A=" << A <<
"\n"
145 <<
"C=" << C <<
"\n";
149 FEELPP_ASSERT( lu.isNonsingular() )( A )( D )( C ).error(
"vandermonde matrix is singular" );
151 this->setCoefficient( ublas::trans( C ) );
159 DVLOG(2) <<
"============================================================\n";
165 M_primal( fe.M_primal )
177 self_type&
operator=( self_type
const& fe )
182 M_primal = fe.M_primal;
189 template<
typename AE>
190 value_type operator()( uint16_type i, ublas::vector_expression<AE>
const& pt )
const
192 return this->evaluate( i, pt );
195 template<
typename AE>
196 value_type operator()( ublas::vector_expression<AE>
const& pt )
const
198 matrix_type m( pt().size(), 1 );
199 ublas::column( m, 0 ) = pt;
200 ublas::vector<value_type> r( ublas::column( this->evaluate( m ), 0 ) );
201 return ublas::inner_prod( r, ublas::scalar_vector<value_type>( r.size(), 1.0 ) );
205 matrix_type operator()( points_type
const& pts )
const
207 return this->evaluate( pts );
247 dual_space_type
const&
dual()
const
257 return M_dual.points();
269 points_type
const&
points( uint16_type f )
const
271 return M_dual.points( f );
277 virtual std::string familyName()
const = 0;
284 dual_space_type M_dual;
285 primal_space_type
const& M_primal;
290 template<
class Pr,
template<
class,u
int16_type,
class>
class Pt>
class PDual,
291 template<class,uint16_type,class> class Pts>
292 const uint16_type FiniteElement<P,PDual,Pts>::nLocalDof;
295 template<class Pr, template<class,uint16_type,class> class Pt> class PDual,
296 template<class,uint16_type,class> class Pts>
297 const uint16_type FiniteElement<P,PDual,Pts>::nDofPerVertex;
300 template<class Pr, template<class,uint16_type,class> class Pt> class PDual,
301 template<class,uint16_type,class> class Pts>
302 const uint16_type FiniteElement<P,PDual,Pts>::nDofPerEdge;
305 template<class Pr, template<class,uint16_type,class> class Pt> class PDual,
306 template<class,uint16_type,class> class Pts>
307 const uint16_type FiniteElement<P,PDual,Pts>::nDofPerFace;
310 template<class Pr, template<class,uint16_type,class> class Pt> class PDual,
311 template<class,uint16_type,class> class Pts>
312 const uint16_type FiniteElement<P,PDual,Pts>::nDofPerVolume;