32 #include <boost/ptr_container/ptr_vector.hpp>
34 #include <feel/feelcore/feel.hpp>
35 #include <feel/feelcore/traits.hpp>
57 template<
typename Basis,
template<
class, u
int16_type,
class>
class PointSetType>
60 public DualBasis<Basis>
62 typedef DualBasis<Basis> super;
65 static const uint16_type nDim = super::nDim;
66 static const uint16_type nOrder= super::nOrder;
68 typedef typename super::primal_space_type primal_space_type;
69 typedef typename primal_space_type::value_type value_type;
70 typedef typename primal_space_type::points_type points_type;
71 typedef typename primal_space_type::matrix_type matrix_type;
72 typedef typename primal_space_type::convex_type convex_type;
73 typedef typename primal_space_type::reference_convex_type reference_convex_type;
74 typedef typename reference_convex_type::node_type node_type;
77 typedef PointSetType<convex_type, 3, value_type> pointset_type;
85 typedef mpl::vector_c<uint16_type, 0, 1, 3, ( 4 ) + ( 4 ) - 2> edges_t;
86 typedef mpl::vector_c<uint16_type, 0, 0, 1, 4> geo_faces_t;
87 typedef mpl::vector_c<uint16_type, 0, 2, 3, 4> faces_t;
88 typedef mpl::vector_c<uint16_type, 0, 2, 3, 4> normals_t;
91 static const uint16_type nVertices = reference_convex_type::numVertices;
92 static const uint16_type nFaces = reference_convex_type::numFaces;
93 static const uint16_type nGeometricFaces = reference_convex_type::numFaces;
94 static const uint16_type nEdges = reference_convex_type::numEdges;
95 static const uint16_type nNormals = reference_convex_type::numNormals;
97 static const uint16_type nbPtsPerVertex = ( nDim+1 );
98 static const uint16_type nbPtsPerEdge = 0;
99 static const uint16_type nbPtsPerFace = ( nDim==2 )?1:0;
100 static const uint16_type nbPtsPerVolume = 0;
101 static const uint16_type numPoints = nVertices*nbPtsPerVertex+nGeometricFaces*nbPtsPerFace;
105 static const uint16_type nDofPerVertex = nbPtsPerVertex;
108 static const uint16_type nDofPerEdge = 0;
111 static const uint16_type nDofPerFace = nbPtsPerFace;
114 static const uint16_type nDofPerVolume = 0;
117 static const uint16_type nLocalDof = nDofPerVertex*nVertices + nDofPerFace*nGeometricFaces;
119 static const uint16_type nFacesInConvex = mpl::if_< mpl::equal_to<mpl::int_<nDim>, mpl::int_<1> >,
120 mpl::int_<nVertices>,
121 typename mpl::if_<mpl::equal_to<mpl::int_<nDim>, mpl::int_<2> >,
123 mpl::int_<nFaces> >::type >::type::value;
125 HermiteDual( primal_space_type
const& primal )
129 M_eid( M_convex_ref.topologicalDimension()+1 ),
130 M_pts( nDim, nLocalDof ),
131 M_points_face( nFacesInConvex ),
134 DVLOG(2) <<
"Hermite finite element: \n";
135 DVLOG(2) <<
" o- dim = " << nDim <<
"\n";
136 DVLOG(2) <<
" o- order = " << nOrder <<
"\n";
137 DVLOG(2) <<
" o- nVertices = " << nVertices <<
"\n";
138 DVLOG(2) <<
" o- nGeometricFaces= " << nGeometricFaces <<
"\n";
139 DVLOG(2) <<
" o- numPoints = " << numPoints <<
"\n";
140 DVLOG(2) <<
" o- nbPtsPerVertex = " << nbPtsPerVertex <<
"\n";
141 DVLOG(2) <<
" o- nbPtsPerEdge = " << nbPtsPerEdge <<
"\n";
142 DVLOG(2) <<
" o- nbPtsPerFace = " << nbPtsPerFace <<
"\n";
143 DVLOG(2) <<
" o- nbPtsPerVolume = " << nbPtsPerVolume <<
"\n";
145 DVLOG(2) <<
" o- nbDofPerVertex = " << nDofPerVertex <<
"\n";
146 DVLOG(2) <<
" o- nbDofPerEdge = " << nDofPerEdge <<
"\n";
147 DVLOG(2) <<
" o- nbDofPerFace = " << nDofPerFace <<
"\n";
148 DVLOG(2) <<
" o- nbDofPerVolume = " << nDofPerVolume <<
"\n";
149 DVLOG(2) <<
" o- nLocalDof = " << nLocalDof <<
"\n";
157 for ( uint16_type e = M_convex_ref.entityRange( 0 ).begin();
158 e < M_convex_ref.entityRange( 0 ).end();
161 M_points_face[e] = pts.pointsBySubEntity( 0, e, 1 );
162 points_type _pts = pts.pointsBySubEntity( 0, e, 1 );
164 for (
int j = 0; j < _pts.size2(); ++j, ++i )
166 ublas::column( M_pts, i ) = ublas::column( _pts, j );
167 DVLOG(2) <<
"pts " << i <<
" = " << ublas::column( M_pts, i ) <<
"\n";
171 for (
int j = 0; j < nDim; ++j )
173 ublas::subrange( M_pts, 0, nDim, i, i+nDim+1 ) = ublas::subrange( M_pts, 0, nDim, 0, nDim+1 );
179 points_type _pts = pts.pointsBySubEntity( 2, 0, 0 );
180 ublas::column( M_pts, i ) = ublas::column( _pts, 0 );
181 DVLOG(2) <<
"pts " << i <<
" = " << ublas::column( M_pts, i ) <<
"\n";
185 setFset( primal, M_pts, mpl::bool_<primal_space_type::is_scalar>() );
191 points_type
const&
points()
const
196 points_type
const&
points( uint16_type f )
const
198 return M_points_face[f];
200 ublas::matrix_column<points_type const> point( uint16_type f, uint32_type __i )
const
202 return ublas::column( M_points_face[f], __i );
204 ublas::matrix_column<points_type> point( uint16_type f, uint32_type __i )
206 return ublas::column( M_points_face[f], __i );
209 matrix_type operator()( primal_space_type
const& pset )
const
211 matrix_type m = M_fset( pset );
217 void setFset( primal_space_type
const& primal, points_type
const& __pts, mpl::bool_<true> )
219 int nfs = 1+nDim+( ( nDim==2 )?1:0 );
220 std::vector<std::vector<Functional<primal_space_type> > > pd( nfs );
222 ublas::range( 0,nDim ),
223 ublas::range( 0, nDim+1 ) );
224 pd[0] = functional::PointsEvaluation<primal_space_type>( primal, pts1 );
227 for (
int d = 0; d < nDim; ++d )
229 pd[d+1] = functional::PointsDerivative<primal_space_type>( primal, d, pts1 );
236 ublas::range( 0,nDim ),
237 ublas::range( ( nDim+1 )*( nDim+1 ), numPoints ) );
239 pd[3] = functional::PointsEvaluation<primal_space_type>( primal, pts3 );
243 std::vector<Functional<primal_space_type> > fs( nLocalDof );
244 typename std::vector<Functional<primal_space_type> >::iterator it = fs.begin();
246 for (
int i = 0; i < nfs; ++i )
249 it = std::copy( pd[i].begin(), pd[i].end(), it );
252 M_fset.setFunctionalSet( fs );
255 void setFset( primal_space_type
const& primal, points_type
const& __pts, mpl::bool_<false> )
264 void setPoints( uint16_type f, points_type
const& n )
266 M_points_face[f].resize( n.size1(), n.size2(), false );
267 M_points_face[f] = n;
271 reference_convex_type M_convex_ref;
272 std::vector<std::vector<uint16_type> > M_eid;
274 std::vector<points_type> M_points_face;
275 FunctionalSet<primal_space_type> M_fset;
297 template<uint16_type N,
299 template<u
int16_type Dim>
class PolySetType,
301 template<u
int16_type, u
int16_type, u
int16_type>
class Convex = Simplex,
302 template<
class, u
int16_type,
class>
class Pts = PointSetEquiSpaced >
305 public FiniteElement<detail::OrthonormalPolynomialSet<N, N, O, PolySetType, T, Convex>, details::HermiteDual, Pts >
310 BOOST_STATIC_ASSERT( ( boost::is_same<PolySetType<N>,
Scalar<N> >::value ||
312 boost::is_same<PolySetType<N>,
Tensor2<N> >::value ) );
318 static const uint16_type nDim = N;
319 static const uint16_type nOrder = O;
320 static const bool isTransformationEquivalent =
false;
321 static const bool isContinuous =
true;
323 typedef typename super::value_type value_type;
324 typedef typename super::primal_space_type primal_space_type;
325 typedef typename super::dual_space_type dual_space_type;
331 static const bool is_tensor2 = polyset_type::is_tensor2;
332 static const bool is_vectorial = polyset_type::is_vectorial;
333 static const bool is_scalar = polyset_type::is_scalar;
334 static const uint16_type nComponents = polyset_type::nComponents;
339 typedef typename mpl::if_<mpl::equal_to<mpl::int_<nDim>, mpl::int_<1> >,
340 mpl::identity<boost::none_t>,
343 typedef boost::shared_ptr<face_basis_type> face_basis_ptrtype;
345 typedef typename dual_space_type::convex_type convex_type;
346 typedef typename dual_space_type::pointset_type pointset_type;
347 typedef typename dual_space_type::reference_convex_type reference_convex_type;
348 typedef typename reference_convex_type::node_type node_type;
349 typedef typename reference_convex_type::points_type points_type;
351 static const uint16_type numPoints = reference_convex_type::numPoints;
352 static const uint16_type nbPtsPerVertex = reference_convex_type::nbPtsPerVertex;
353 static const uint16_type nbPtsPerEdge = reference_convex_type::nbPtsPerEdge;
354 static const uint16_type nbPtsPerFace = reference_convex_type::nbPtsPerFace;
355 static const uint16_type nbPtsPerVolume = reference_convex_type::nbPtsPerVolume;
366 super( dual_space_type( primal_space_type() ) ),
418 template<
typename GMContext,
typename PC,
typename Phi,
typename GPhi,
typename HPhi >
419 static void transform( boost::shared_ptr<GMContext> gmc, boost::shared_ptr<PC>
const& pc,
421 GPhi& g_phi_t,
const bool do_gradient,
422 HPhi& h_phi_t,
const bool do_hessian
426 transform ( *gmc, *pc, phi_t, g_phi_t, do_gradient, h_phi_t, do_hessian );
428 template<
typename GMContext,
typename PC,
typename Phi,
typename GPhi,
typename HPhi >
429 static void transform( GMContext
const& gmc,
432 GPhi& g_phi_t,
const bool do_gradient,
433 HPhi& h_phi_t,
const bool do_hessian
438 typename GMContext::gm_type::matrix_type
const& B = gmc.B( 0 );
439 std::fill( phi_t.data(), phi_t.data()+phi_t.num_elements(), value_type( 0 ) );
442 std::fill( g_phi_t.data(), g_phi_t.data()+g_phi_t.num_elements(), value_type( 0 ) );
445 std::fill( h_phi_t.data(), h_phi_t.data()+g_phi_t.num_elements(), value_type( 0 ) );
447 const uint16_type Q = gmc.nPoints();
450 for ( uint16_type i = numPoints; i < nLocalDof; i+=nDim )
452 for ( uint16_type l = 0; l < nDim; ++l )
454 for ( uint16_type p = 0; p < nDim; ++p )
456 for ( uint16_type q = 0; q < Q; ++q )
461 phi_t[i+l][0][0][q] += B( l, p ) * pc.phi( i+p,0,0,q );
465 for ( uint16_type j = 0; j < nDim; ++j )
467 g_phi_t[i+l][0][p][q] += B( l, j ) * pc.grad( i+j,0,p,q );
483 reference_convex_type M_refconvex;
484 face_basis_ptrtype M_bdylag;
486 template<uint16_type N,
488 template<u
int16_type Dim>
class PolySetType,
490 template<u
int16_type, u
int16_type, u
int16_type>
class Convex,
491 template<
class, u
int16_type,
class>
class Pts >
492 const uint16_type Hermite<N,O,PolySetType,T,Convex,Pts>::nDim;
494 template<uint16_type N,
496 template<u
int16_type Dim>
class PolySetType,
498 template<u
int16_type, u
int16_type, u
int16_type>
class Convex,
499 template<
class, u
int16_type,
class>
class Pts >
500 const uint16_type Hermite<N,O,PolySetType,T,Convex,Pts>::nOrder;
502 template<uint16_type N,
504 template<u
int16_type Dim>
class PolySetType,
506 template<u
int16_type, u
int16_type, u
int16_type>
class Convex,
507 template<
class, u
int16_type,
class>
class Pts >
508 const uint16_type Hermite<N,O,PolySetType,T,Convex,Pts>::numPoints;
512 template<uint16_type Order,
513 template<u
int16_type Dim>
class PolySetType = Scalar,
514 template<
class, u
int16_type,
class>
class Pts = PointSetEquiSpaced>
518 template<uint16_type N,
520 typename Convex = Simplex<N> >
523 typedef typename mpl::if_<mpl::bool_<Convex::is_simplex>,
524 mpl::identity<fem::Hermite<N,Order,PolySetType,T,Simplex, Pts> >,
525 mpl::identity<fem::Hermite<N,Order,PolySetType,T,Hypercube, Pts> > >::type::type result_type;
526 typedef result_type type;
529 typedef Hermite<Order,Scalar,Pts> component_basis_type;
531 static const uint16_type nOrder = Order;