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>
20 template<
class Field >
32 typedef Field MPField;
34 GenericGeometry::LobattoPoints<MPField> lobatto(order+1);
36 for (
unsigned int i=1;i<
order;++i) {
56 template <
class Field,
class Topology>
59 template <
class Field>
62 static const unsigned int dimension = 0;
63 static unsigned int size(
const unsigned int order)
67 template <
unsigned int dim>
68 static unsigned int setupSimplex(
69 const unsigned int iup,
70 const unsigned int dimStart,
71 unsigned int startPos,
72 const std::vector<Field> &points1D,
75 const unsigned int order = points1D.size()+1;
76 unsigned int i = order+1-iup;
77 if (i-2>=points1D.size())
81 for (
unsigned int d=0;d<dimStart;++d )
83 points[startPos].
point_[d] = -points1D[i-2];
88 template <
unsigned int dim>
89 static void setup(
const std::vector<Field> &points1D,
95 template <
class Field,
class Base>
99 static const unsigned int dimension = Base::dimension+1;
100 static unsigned int size(
const unsigned int order)
102 if (order<=dimension)
return 0;
104 for (
unsigned int o=0;o<order-dimension;++o)
105 size += LobattoBase::size(o+dimension);
108 template <
unsigned int dim>
109 static unsigned int setupSimplex(
110 const unsigned int iup,
111 const unsigned int dimStart,
112 unsigned int startPos,
113 const std::vector<Field> &points1D,
116 const unsigned int order = points1D.size()+1;
117 unsigned int endPos = startPos;
118 for (
unsigned int i=2;i<=order-iup;++i)
120 endPos = LobattoBase::template setupSimplex<dim>(iup+i-1,dimStart,startPos,points1D,points);
121 for (
unsigned int j=startPos;j<endPos;++j)
123 for (
unsigned int d=0;d<dimStart;++d)
125 if ( d==dimStart-dimension )
126 points[j].
point_[d] += dimStart*points1D[i-2];
128 points[j].
point_[d] -= points1D[i-2];
135 template <
unsigned int dim>
136 static void setup(
const std::vector<Field> &points1D,
139 const unsigned int order = points1D.size()+1;
140 unsigned int startPos=0,endPos=0;
141 for (
unsigned int i=2;i<=order;++i)
143 endPos = LobattoBase::template setupSimplex<dim>(i-1,dimension,startPos,points1D,points);
145 for (
unsigned int j=startPos;j<endPos;++j)
147 for (
unsigned int d=0;d<dimension;++d)
150 points[j].
point_[d] += 1.+int(dimension)*points1D[i-2];
152 points[j].
point_[d] += 1.-points1D[i-2];
153 points[j].
point_[d] /= (dimension+1);
160 template <
class Field,
class Base>
164 static const unsigned int dimension = Base::dimension+1;
165 static unsigned int size(
const unsigned int order)
167 return LobattoBase::size(order)*(order-1);
169 template <
unsigned int dim>
170 static unsigned int setupSimplex(
171 const unsigned int iup,
172 const unsigned int dimStart,
173 unsigned int startPos,
174 const std::vector<Field> &points1D,
177 const unsigned int order = points1D.size()+1;
178 unsigned int endPos = startPos;
179 for (
unsigned int i=2;i<=order-iup;++i)
181 endPos = LobattoBase::template setupSimplex<dim>(iup+i-1,dimStart,startPos,points1D,points);
182 for (
unsigned int j=startPos;j<endPos;++j)
184 for (
unsigned int d=0;d<dimStart;++d)
186 if ( d==dimStart-dimension )
187 points[j].
point_[d] += dimStart*points1D[i-2];
189 points[j].
point_[d] -= points1D[i-2];
196 template <
unsigned int dim>
197 static void setup(
const std::vector<Field> &points1D,
200 const unsigned int order = points1D.size()+1;
201 assert(dim>=dimension);
202 assert(points1D.size()==order-1);
203 LobattoBase::template setup<dim>(points1D,points);
204 const unsigned int baseSize = LobattoBase::size(order);
205 for (
unsigned int q=0;q<points1D.size();q++)
207 for (
unsigned int i=0;i<baseSize;++i)
209 const unsigned int pos = q*baseSize+i;
210 for (
unsigned int d=0;d<dimension-1;++d)
211 points[pos].point_[d] = points[i].point_[d];
212 points[pos].
point_[dimension-1]=points1D[q];
218 template<
class F,
unsigned int dim >
231 template<
class Topology >
236 ForLoop<Setup<Topology>::template InitCodim,0,
dimension>::
240 template<
class Topology >
252 template <
class Topology>
260 const std::vector<Field> &points1D,
261 std::vector<Point> &points)
264 ForLoop<InitSub,0,size-1>::apply(order,points1D,points);
269 typedef typename GenericGeometry::SubTopology<Topology,codim,i>::type
SubTopology;
271 const std::vector<Field> &points1D,
272 std::vector<Point> &points)
278 template <
class SubTopology>
284 typedef GenericReferenceElement< Field, dimension >
RefElement;
285 typedef typename RefElement::template Codim< codimension >::Mapping
Mapping;
288 template <
unsigned int subEntity>
290 const std::vector<Field> &points1D,
291 std::vector<Point> &points)
293 unsigned int oldSize = points.size();
294 unsigned int size = InnerPoints::size(order);
297 points.resize(oldSize+size);
302 const GeometryType geoType( Topology::id,
dimension );
303 const RefElement &refElement = RefElements::general( geoType );
304 const Mapping &mapping = refElement.template mapping< codimension >( subEntity );
307 for (
unsigned int nr = 0; nr<
size; ++nr, ++p)
309 p->
point_ = mapping.global( subPoints[nr].point_ );
310 p->localKey_ =
LocalKey( subEntity, codimension, nr );
312 bool test = GenericGeometry::ReferenceElement<Topology,Field>::checkInside(p->point_);
314 std::cout <<
"not inside" << std::endl;
323 #endif // DUNE_LOBATTOBASIS_HH