33 #include <feel/feelcore/feel.hpp>
44 namespace ublas = boost::numeric::ublas;
50 template<u
int16_type N>
53 static const uint16_type jacobi_degree = N;
54 static const uint16_type integration_degree = 2*jacobi_degree-1;
60 template<u
int16_type N>
63 static const uint16_type integration_degree = N;
64 static const uint16_type jacobi_degree = ( integration_degree+1 )/2+1;
87 template<u
int16_type, u
int16_type, u
int16_type>
class Entity =
Simplex,
88 template<
class Convex, u
int16_type O,
typename T2>
class QPS =
Gauss,
92 public QPS<Entity<Dim,1,Dim>, DegreePolicy<Order>::integration_degree, T>
94 typedef QPS<Entity<Dim,1,Dim>, DegreePolicy<Order>::integration_degree, T> super;
96 static const bool is_exact =
false;
97 typedef DegreePolicy<Order> degree_policy_type;
98 static const uint16_type nDim = Dim;
99 static const uint16_type nNodes = degree_policy_type::jacobi_degree;
100 static const uint16_type nOrder = degree_policy_type::integration_degree;
103 typedef Entity<Dim,nOrder,Dim> convex_type;
104 typedef typename super::value_type value_type;
105 typedef typename super::node_type node_type;
106 typedef typename super::weights_type weights_type;
107 typedef typename super::nodes_type nodes_type;
109 typedef boost::tuple<nodes_type, weights_type> quadrature_data_type;
110 typedef typename super::face_quadrature_type face_quadrature_type;
121 quadrature_data_type data()
const
123 return boost::make_tuple( this->
points(), this->weights() );
131 template<u
int16_type, u
int16_type, u
int16_type>
class Entity =
Simplex>
134 public mpl::if_<mpl::and_<mpl::or_<mpl::and_<mpl::less_equal<mpl::int_<IMORDER>,mpl::int_<20> >,
135 mpl::equal_to<mpl::int_<DIM>,mpl::int_<2> > >,
136 mpl::and_<mpl::less_equal<mpl::int_<IMORDER>,mpl::int_<20> >,
137 mpl::equal_to<mpl::int_<DIM>,mpl::int_<3> > > >,
138 mpl::bool_<Entity<DIM,1,DIM>::is_simplex> >,
139 mpl::identity<IMSimplex<DIM, IMORDER, T> >,
140 mpl::identity<IMGeneral<DIM, IMORDER, T, Entity> > >::type::type
144 template<u
int16_type, u
int16_type, u
int16_type>
class Entity1>
147 typedef typename mpl::if_<mpl::and_<mpl::or_<mpl::and_<mpl::less_equal<mpl::int_<IMORDER>,mpl::int_<20> >,
148 mpl::equal_to<mpl::int_<DIM1>,mpl::int_<2> > >,
149 mpl::and_<mpl::less_equal<mpl::int_<IMORDER>,mpl::int_<20> >,
150 mpl::equal_to<mpl::int_<DIM1>,mpl::int_<3> > > >,
151 mpl::bool_<Entity1<DIM1,1,DIM1>::is_simplex> >,
152 mpl::identity<IMSimplex<DIM1, IMORDER, T1> >,
153 mpl::identity<IMGeneral<DIM1, IMORDER, T1, Entity1> > >::type::type type;
157 template<
int IMORDER>
162 template<u
int16_type, u
int16_type, u
int16_type>
class Entity>
165 typedef typename mpl::if_<mpl::and_<mpl::or_<mpl::and_<mpl::less_equal<mpl::int_<IMORDER>,mpl::int_<20> >,
166 mpl::equal_to<mpl::int_<DIM>,mpl::int_<2> > >,
167 mpl::and_<mpl::less_equal<mpl::int_<IMORDER>,mpl::int_<20> >,
168 mpl::equal_to<mpl::int_<DIM>,mpl::int_<3> > > >,
169 mpl::bool_<Entity<DIM,1,DIM>::is_simplex> >,
170 mpl::identity<IMSimplex<DIM, IMORDER, T> >,
171 mpl::identity<IMGeneral<DIM, IMORDER, T, Entity> > >::type::type type;
174 template<
typename ContextType>
177 static const int DIM = ContextType::PDim;
178 typedef typename ContextType::value_type T;
179 typedef typename mpl::if_<mpl::and_<mpl::or_<mpl::and_<mpl::less_equal<mpl::int_<IMORDER>,mpl::int_<20> >,
180 mpl::equal_to<mpl::int_<DIM>,mpl::int_<2> > >,
181 mpl::and_<mpl::less_equal<mpl::int_<IMORDER>,mpl::int_<20> >,
182 mpl::equal_to<mpl::int_<DIM>,mpl::int_<3> > > >,
183 mpl::bool_<ContextType::element_type::is_simplex> >,
184 mpl::identity<IMSimplex<DIM, IMORDER, T> >,
185 typename mpl::if_<mpl::bool_<ContextType::element_type::is_simplex>,
186 mpl::identity<IMGeneral<DIM, IMORDER, T, Simplex> >,
187 mpl::identity<IMGeneral<DIM, IMORDER, T, Hypercube> > >::type>::type::type type;
196 template<u
int16_type,u
int16_type,u
int16_type>
class Entity,
197 template<
class Convex, u
int16_type O,
typename T2>
class QPS,
198 template<u
int16_type N>
class DegreePolicy>
200 operator<<( std::ostream& os,
201 Feel::IM<Dim, Order, T, Entity,QPS, DegreePolicy>
const& qps )
203 os <<
"quadrature point set:\n"
204 <<
"number of points: " << qps.nPoints() <<
"\n"
205 <<
"points : " << qps.points() <<
"\n"
206 <<
"weights : " << qps.weights() <<
"\n";
208 for (
size_type f = 0; f < qps.nFaces(); ++f )
210 os <<
" o Face " << f <<
"\n";
211 os <<
" number of points: " << qps.nPointsOnFace( f ) <<
"\n"
212 <<
" points : " << qps.fpoints( f ) <<
"\n"
213 <<
" weights : " << qps.weights( f ) <<
"\n";
219 template<
int Dim,
int IMORDER,
typename T>
struct ImBestSimplex
221 public mpl::if_<mpl::less_equal<mpl::int_<IMORDER>,mpl::int_<6> >,
222 mpl::identity<IMSimplex<Dim, IMORDER, T> >,
223 mpl::identity<IM<Dim, IMORDER, T, Simplex> > >::type::type