1 #ifndef DUNE_LOBATTOBASIS_HH
2 #define DUNE_LOBATTOBASIS_HH
6 #include <dune/common/forloop.hh>
8 #include <dune/geometry/topologyfactory.hh>
9 #include <dune/geometry/quadraturerules/lobattoquadrature.hh>
10 #include <dune/geometry/referenceelements.hh>
11 #include <dune/geometry/genericgeometry/referenceelements.hh>
21 template<
class Field >
33 typedef Field MPField;
35 GenericGeometry::LobattoPoints<MPField> lobatto(order+1);
37 for (
unsigned int i=1;i<
order;++i) {
57 template <
class Field,
class Topology>
58 struct LobattoInnerPoints;
60 template <
class Field>
61 struct LobattoInnerPoints<Field,GenericGeometry::Point>
63 static const unsigned int dimension = 0;
64 static unsigned int size(
const unsigned int order)
68 template <
unsigned int dim>
69 static unsigned int setupSimplex(
70 const unsigned int iup,
71 const unsigned int dimStart,
72 unsigned int startPos,
73 const std::vector<Field> &points1D,
76 const unsigned int order = points1D.size()+1;
77 unsigned int i = order+1-iup;
78 if (i-2>=points1D.size())
82 for (
unsigned int d=0;d<dimStart;++d )
84 points[startPos].
point_[d] = -points1D[i-2];
89 template <
unsigned int dim>
90 static void setup(
const std::vector<Field> &points1D,
96 template <
class Field,
class Base>
97 struct LobattoInnerPoints<Field,GenericGeometry::Pyramid<Base> >
100 static const unsigned int dimension = Base::dimension+1;
101 static unsigned int size(
const unsigned int order)
103 if (order<=dimension)
return 0;
105 for (
unsigned int o=0;o<order-dimension;++o)
106 size += LobattoBase::size(o+dimension);
109 template <
unsigned int dim>
110 static unsigned int setupSimplex(
111 const unsigned int iup,
112 const unsigned int dimStart,
113 unsigned int startPos,
114 const std::vector<Field> &points1D,
117 const unsigned int order = points1D.size()+1;
118 unsigned int endPos = startPos;
119 for (
unsigned int i=2;i<=order-iup;++i)
121 endPos = LobattoBase::template setupSimplex<dim>(iup+i-1,dimStart,startPos,points1D,points);
122 for (
unsigned int j=startPos;j<endPos;++j)
124 for (
unsigned int d=0;d<dimStart;++d)
126 if ( d==dimStart-dimension )
127 points[j].
point_[d] += dimStart*points1D[i-2];
129 points[j].
point_[d] -= points1D[i-2];
136 template <
unsigned int dim>
137 static void setup(
const std::vector<Field> &points1D,
140 const unsigned int order = points1D.size()+1;
141 unsigned int startPos=0,endPos=0;
142 for (
unsigned int i=2;i<=order;++i)
144 endPos = LobattoBase::template setupSimplex<dim>(i-1,dimension,startPos,points1D,points);
146 for (
unsigned int j=startPos;j<endPos;++j)
148 for (
unsigned int d=0;d<dimension;++d)
151 points[j].
point_[d] += 1.+int(dimension)*points1D[i-2];
153 points[j].
point_[d] += 1.-points1D[i-2];
154 points[j].
point_[d] /= (dimension+1);
161 template <
class Field,
class Base>
162 struct LobattoInnerPoints<Field,GenericGeometry::Prism<Base> >
165 static const unsigned int dimension = Base::dimension+1;
166 static unsigned int size(
const unsigned int order)
168 return LobattoBase::size(order)*(order-1);
170 template <
unsigned int dim>
171 static unsigned int setupSimplex(
172 const unsigned int iup,
173 const unsigned int dimStart,
174 unsigned int startPos,
175 const std::vector<Field> &points1D,
178 const unsigned int order = points1D.size()+1;
179 unsigned int endPos = startPos;
180 for (
unsigned int i=2;i<=order-iup;++i)
182 endPos = LobattoBase::template setupSimplex<dim>(iup+i-1,dimStart,startPos,points1D,points);
183 for (
unsigned int j=startPos;j<endPos;++j)
185 for (
unsigned int d=0;d<dimStart;++d)
187 if ( d==dimStart-dimension )
188 points[j].
point_[d] += dimStart*points1D[i-2];
190 points[j].
point_[d] -= points1D[i-2];
197 template <
unsigned int dim>
198 static void setup(
const std::vector<Field> &points1D,
201 const unsigned int order = points1D.size()+1;
202 assert(dim>=dimension);
203 assert(points1D.size()==order-1);
204 LobattoBase::template setup<dim>(points1D,points);
205 const unsigned int baseSize = LobattoBase::size(order);
206 for (
unsigned int q=0;q<points1D.size();q++)
208 for (
unsigned int i=0;i<baseSize;++i)
210 const unsigned int pos = q*baseSize+i;
211 for (
unsigned int d=0;d<dimension-1;++d)
212 points[pos].point_[d] = points[i].point_[d];
213 points[pos].
point_[dimension-1]=points1D[q];
219 template<
class F,
unsigned int dim >
232 template<
class Topology >
237 ForLoop<Setup<Topology>::template InitCodim,0,
dimension>::
241 template<
class Topology >
253 template <
class Topology>
261 const std::vector<Field> &points1D,
262 std::vector<Point> &points)
270 typedef typename GenericGeometry::SubTopology<Topology,codim,i>::type
SubTopology;
272 const std::vector<Field> &points1D,
273 std::vector<Point> &points)
279 template <
class SubTopology>
285 typedef GenericReferenceElement< Field, dimension >
RefElement;
286 typedef typename RefElement::template Codim< codimension >::Mapping
Mapping;
289 template <
unsigned int subEntity>
291 const std::vector<Field> &points1D,
292 std::vector<Point> &points)
294 unsigned int oldSize = points.size();
298 points.resize(oldSize+size);
303 const GeometryType geoType( Topology::id,
dimension );
304 const RefElement &refElement = RefElements::general( geoType );
305 const Mapping &mapping = refElement.template mapping< codimension >( subEntity );
308 for (
unsigned int nr = 0; nr<
size; ++nr, ++p)
310 p->
point_ = mapping.global( subPoints[nr].point_ );
311 p->localKey_ =
LocalKey( subEntity, codimension, nr );
313 bool test = GenericGeometry::ReferenceElement<Topology,Field>::checkInside(p->point_);
315 std::cout <<
"not inside" << std::endl;
324 #endif // DUNE_LOBATTOBASIS_HH