dune-localfunctions 3.0-git
raviartthomassimplexinterpolation.hh
Go to the documentation of this file.
1// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2// vi: set et ts=4 sw=2 sts=2:
3#ifndef DUNE_RAVIARTTHOMASINTERPOLATION_HH
4#define DUNE_RAVIARTTHOMASINTERPOLATION_HH
5
6#include <fstream>
7#include <utility>
8
9#include <dune/common/exceptions.hh>
10#include <dune/common/forloop.hh>
11
12#include <dune/geometry/topologyfactory.hh>
13#include <dune/geometry/referenceelements.hh>
14#include <dune/geometry/quadraturerules.hh>
15
20
21namespace Dune
22{
23
24 // LocalCoefficientsContainer
25 // -------------------
26 template < unsigned int dim >
27 struct RaviartThomasCoefficientsFactory;
28 template < unsigned int dim, class Field >
29 struct RaviartThomasL2InterpolationFactory;
30
32 {
34
35 public:
36 template <class Setter>
37 LocalCoefficientsContainer ( const Setter &setter )
38 {
39 setter.setLocalKeys(localKey_);
40 }
41
42 const LocalKey &localKey ( const unsigned int i ) const
43 {
44 assert( i < size() );
45 return localKey_[ i ];
46 }
47
48 unsigned int size () const
49 {
50 return localKey_.size();
51 }
52
53 private:
54 std::vector< LocalKey > localKey_;
55 };
56
57 template < unsigned int dim >
59 {
60 static const unsigned int dimension = dim;
62 typedef unsigned int Key;
64 };
65 template < unsigned int dim >
67 public TopologyFactory< RaviartThomasCoefficientsFactoryTraits<dim> >
68 {
70 template <class Topology>
71 static typename Traits::Object *createObject( const typename Traits::Key &key )
72 {
73 typedef RaviartThomasL2InterpolationFactory<dim,double> InterpolationFactory;
74 if (! supports<Topology>(key) )
75 return 0;
76 typename InterpolationFactory::Object *interpol
77 = InterpolationFactory::template create<Topology>(key);
78 typename Traits::Object *localKeys = new typename Traits::Object(*interpol);
79 InterpolationFactory::release(interpol);
80 return localKeys;
81 }
82 template< class Topology >
83 static bool supports ( const typename Traits::Key &key )
84 {
85 return GenericGeometry::IsSimplex<Topology>::value;
86 }
87 };
88
89 // LocalInterpolation
90 // -------------------
91
92 // -----------------------------------------
93 // RTL2InterpolationBuilder
94 // -----------------------------------------
95 // L2 Interpolation requires:
96 // - for element
97 // - test basis
98 // - for each face (dynamic)
99 // - test basis
100 // - normal
101 template <unsigned int dim, class Field>
103 {
104 static const unsigned int dimension = dim;
109 typedef FieldVector<Field,dimension> Normal;
110
113
115 {
116 TestBasisFactory::release(testBasis_);
117 for (unsigned int i=0; i<faceStructure_.size(); ++i)
118 TestFaceBasisFactory::release(faceStructure_[i].basis_);
119 }
120
121 unsigned int topologyId() const
122 {
123 return topologyId_;
124 }
125 unsigned int order() const
126 {
127 return order_;
128 }
129 unsigned int faceSize() const
130 {
131 return faceSize_;
132 }
134 {
135 return testBasis_;
136 }
137 TestFaceBasis *testFaceBasis( unsigned int f ) const
138 {
139 assert( f < faceSize() );
140 return faceStructure_[f].basis_;
141 }
142 const Normal &normal( unsigned int f ) const
143 {
144 return *(faceStructure_[f].normal_);
145 }
146
147 template <class Topology>
148 void build(unsigned int order)
149 {
150 order_ = order;
151 topologyId_ = Topology::id;
152 if (order>0)
153 testBasis_ = TestBasisFactory::template create<Topology>(order-1);
154 else
155 testBasis_ = 0;
156 const unsigned int size = GenericGeometry::Size<Topology,1>::value;
157 faceSize_ = size;
158 faceStructure_.reserve( faceSize_ );
159 ForLoop< Creator<Topology>::template GetCodim,0,size-1>::apply(order, faceStructure_ );
160 assert( faceStructure_.size() == faceSize_ );
161 }
162
163 private:
164 struct FaceStructure
165 {
166 FaceStructure( TestFaceBasis *tfb,
167 const Normal &n )
168 : basis_(tfb), normal_(&n)
169 {}
170 TestFaceBasis *basis_;
171 const Dune::FieldVector<Field,dimension> *normal_;
172 };
173 template < class Topology >
174 struct Creator
175 {
176 template < int face >
177 struct GetCodim
178 {
179 typedef typename GenericGeometry::SubTopology<Topology,1,face>::type FaceTopology;
180 static void apply( const unsigned int order,
181 std::vector<FaceStructure> &faceStructure )
182 {
183 static const int dimTopo = Topology::dimension;
184 faceStructure.push_back(
185 FaceStructure(
186 TestFaceBasisFactory::template create<FaceTopology>(order),
187 ReferenceElements<Field,dimTopo>::general(GeometryType( Topology() )).integrationOuterNormal(face)
188 ) );
189 }
190 };
191 };
192
193 std::vector<FaceStructure> faceStructure_;
194 TestBasis *testBasis_;
195 unsigned int topologyId_, order_, faceSize_;
196 };
197
198 // A L2 based interpolation for Raviart Thomas
199 // --------------------------------------------------
200 template< unsigned int dimension, class F>
202 : public InterpolationHelper<F,dimension>
203 {
206
207 public:
208 typedef F Field;
211 : order_(0),
212 size_(0)
213 {}
214
215 template< class Function, class Fy >
216 void interpolate ( const Function &function, std::vector< Fy > &coefficients ) const
217 {
218 coefficients.resize(size());
219 typename Base::template Helper<Function,std::vector<Fy>,true> func( function,coefficients );
220 interpolate(func);
221 }
222 template< class Basis, class Matrix >
223 void interpolate ( const Basis &basis, Matrix &matrix ) const
224 {
225 matrix.resize( size(), basis.size() );
226 typename Base::template Helper<Basis,Matrix,false> func( basis,matrix );
227 interpolate(func);
228 }
229
230 unsigned int order() const
231 {
232 return order_;
233 }
234 unsigned int size() const
235 {
236 return size_;
237 }
238 template <class Topology>
239 void build( unsigned int order )
240 {
241 size_ = 0;
242 order_ = order;
243 builder_.template build<Topology>(order_);
244 if (builder_.testBasis())
245 size_ += dimension*builder_.testBasis()->size();
246 for ( unsigned int f=0; f<builder_.faceSize(); ++f )
247 if (builder_.testFaceBasis(f))
248 size_ += builder_.testFaceBasis(f)->size();
249 }
250
251 void setLocalKeys(std::vector< LocalKey > &keys) const
252 {
253 keys.resize(size());
254 unsigned int row = 0;
255 for (unsigned int f=0; f<builder_.faceSize(); ++f)
256 {
257 if (builder_.faceSize())
258 for (unsigned int i=0; i<builder_.testFaceBasis(f)->size(); ++i,++row)
259 keys[row] = LocalKey(f,1,i);
260 }
261 if (builder_.testBasis())
262 for (unsigned int i=0; i<builder_.testBasis()->size()*dimension; ++i,++row)
263 keys[row] = LocalKey(0,0,i);
264 assert( row == size() );
265 }
266
267 protected:
268 template< class Func, class Container, bool type >
269 void interpolate ( typename Base::template Helper<Func,Container,type> &func ) const
270 {
271 const Dune::GeometryType geoType( builder_.topologyId(), dimension );
272
273 std::vector< Field > testBasisVal;
274
275 for (unsigned int i=0; i<size(); ++i)
276 for (unsigned int j=0; j<func.size(); ++j)
277 func.set(i,j,0);
278
279 unsigned int row = 0;
280
281 // boundary dofs:
282 typedef Dune::QuadratureRule<Field, dimension-1> FaceQuadrature;
283 typedef Dune::QuadratureRules<Field, dimension-1> FaceQuadratureRules;
284
285 typedef Dune::ReferenceElements< Field, dimension > RefElements;
286 typedef Dune::ReferenceElement< Field, dimension > RefElement;
287 typedef typename RefElement::template Codim< 1 >::Geometry Geometry;
288
289 const RefElement &refElement = RefElements::general( geoType );
290
291 for (unsigned int f=0; f<builder_.faceSize(); ++f)
292 {
293 if (!builder_.testFaceBasis(f))
294 continue;
295 testBasisVal.resize(builder_.testFaceBasis(f)->size());
296
297 const Geometry &geometry = refElement.template geometry< 1 >( f );
298 const Dune::GeometryType subGeoType( geometry.type().id(), dimension-1 );
299 const FaceQuadrature &faceQuad = FaceQuadratureRules::rule( subGeoType, 2*order_+2 );
300
301 const unsigned int quadratureSize = faceQuad.size();
302 for( unsigned int qi = 0; qi < quadratureSize; ++qi )
303 {
304 if (dimension>1)
305 builder_.testFaceBasis(f)->template evaluate<0>(faceQuad[qi].position(),testBasisVal);
306 else
307 testBasisVal[0] = 1.;
308 fillBnd( row, testBasisVal,
309 func.evaluate( geometry.global( faceQuad[qi].position() ) ),
310 builder_.normal(f), faceQuad[qi].weight(),
311 func);
312 }
313
314 row += builder_.testFaceBasis(f)->size();
315 }
316 // element dofs
317 if (builder_.testBasis())
318 {
319 testBasisVal.resize(builder_.testBasis()->size());
320
321 typedef Dune::QuadratureRule<Field, dimension> Quadrature;
322 typedef Dune::QuadratureRules<Field, dimension> QuadratureRules;
323 const Quadrature &elemQuad = QuadratureRules::rule( geoType, 2*order_+1 );
324
325 const unsigned int quadratureSize = elemQuad.size();
326 for( unsigned int qi = 0; qi < quadratureSize; ++qi )
327 {
328 builder_.testBasis()->template evaluate<0>(elemQuad[qi].position(),testBasisVal);
329 fillInterior( row, testBasisVal,
330 func.evaluate(elemQuad[qi].position()),
331 elemQuad[qi].weight(),
332 func );
333 }
334
335 row += builder_.testBasis()->size()*dimension;
336 }
337 assert(row==size());
338 }
339
340 private:
342 template <class MVal, class RTVal,class Matrix>
343 void fillBnd (unsigned int startRow,
344 const MVal &mVal,
345 const RTVal &rtVal,
346 const FieldVector<Field,dimension> &normal,
347 const Field &weight,
348 Matrix &matrix) const
349 {
350 const unsigned int endRow = startRow+mVal.size();
351 typename RTVal::const_iterator rtiter = rtVal.begin();
352 for ( unsigned int col = 0; col < rtVal.size() ; ++rtiter,++col)
353 {
354 Field cFactor = (*rtiter)*normal;
355 typename MVal::const_iterator miter = mVal.begin();
356 for (unsigned int row = startRow;
357 row!=endRow; ++miter, ++row )
358 {
359 matrix.add(row,col, (weight*cFactor)*(*miter) );
360 }
361 assert( miter == mVal.end() );
362 }
363 }
364 template <class MVal, class RTVal,class Matrix>
365 void fillInterior (unsigned int startRow,
366 const MVal &mVal,
367 const RTVal &rtVal,
368 Field weight,
369 Matrix &matrix) const
370 {
371 const unsigned int endRow = startRow+mVal.size()*dimension;
372 typename RTVal::const_iterator rtiter = rtVal.begin();
373 for ( unsigned int col = 0; col < rtVal.size() ; ++rtiter,++col)
374 {
375 typename MVal::const_iterator miter = mVal.begin();
376 for (unsigned int row = startRow;
377 row!=endRow; ++miter,row+=dimension )
378 {
379 for (unsigned int i=0; i<dimension; ++i)
380 {
381 matrix.add(row+i,col, (weight*(*miter))*(*rtiter)[i] );
382 }
383 }
384 assert( miter == mVal.end() );
385 }
386 }
387
388 Builder builder_;
389 unsigned int order_;
390 unsigned int size_;
391 };
392
393 template < unsigned int dim, class F >
401 template < unsigned int dim, class Field >
403 public TopologyFactory< RaviartThomasL2InterpolationFactoryTraits<dim,Field> >
404 {
407 typedef typename Traits::Object Object;
408 typedef typename std::remove_const<Object>::type NonConstObject;
409 template <class Topology>
410 static typename Traits::Object *createObject( const typename Traits::Key &key )
411 {
412 if ( !supports<Topology>(key) )
413 return 0;
414 NonConstObject *interpol = new NonConstObject();
415 interpol->template build<Topology>(key);
416 return interpol;
417 }
418 template< class Topology >
419 static bool supports ( const typename Traits::Key &key )
420 {
421 return GenericGeometry::IsSimplex<Topology>::value;
422 }
423 };
424}
425#endif // DUNE_RAVIARTTHOMASINTERPOLATION_HH
Definition brezzidouglasmarini1cube2d.hh:14
Describe position of one degree of freedom.
Definition localkey.hh:21
Definition orthonormalbasis.hh:38
Definition raviartthomassimplexinterpolation.hh:68
static bool supports(const typename Traits::Key &key)
Definition raviartthomassimplexinterpolation.hh:83
RaviartThomasCoefficientsFactoryTraits< dim > Traits
Definition raviartthomassimplexinterpolation.hh:69
static Traits::Object * createObject(const typename Traits::Key &key)
Definition raviartthomassimplexinterpolation.hh:71
Definition raviartthomassimplexinterpolation.hh:404
std::remove_const< Object >::type NonConstObject
Definition raviartthomassimplexinterpolation.hh:408
static Traits::Object * createObject(const typename Traits::Key &key)
Definition raviartthomassimplexinterpolation.hh:410
RaviartThomasL2InterpolationFactoryTraits< dim, Field > Traits
Definition raviartthomassimplexinterpolation.hh:405
RTL2InterpolationBuilder< dim, Field > Builder
Definition raviartthomassimplexinterpolation.hh:406
static bool supports(const typename Traits::Key &key)
Definition raviartthomassimplexinterpolation.hh:419
Traits::Object Object
Definition raviartthomassimplexinterpolation.hh:407
Definition raviartthomassimplexinterpolation.hh:32
unsigned int size() const
Definition raviartthomassimplexinterpolation.hh:48
LocalCoefficientsContainer(const Setter &setter)
Definition raviartthomassimplexinterpolation.hh:37
const LocalKey & localKey(const unsigned int i) const
Definition raviartthomassimplexinterpolation.hh:42
Definition raviartthomassimplexinterpolation.hh:59
RaviartThomasCoefficientsFactory< dim > Factory
Definition raviartthomassimplexinterpolation.hh:63
const LocalCoefficientsContainer Object
Definition raviartthomassimplexinterpolation.hh:61
static const unsigned int dimension
Definition raviartthomassimplexinterpolation.hh:60
unsigned int Key
Definition raviartthomassimplexinterpolation.hh:62
Definition raviartthomassimplexinterpolation.hh:103
TestBasis * testBasis() const
Definition raviartthomassimplexinterpolation.hh:133
RTL2InterpolationBuilder()
Definition raviartthomassimplexinterpolation.hh:111
TestBasisFactory::Object TestBasis
Definition raviartthomassimplexinterpolation.hh:106
TestFaceBasisFactory::Object TestFaceBasis
Definition raviartthomassimplexinterpolation.hh:108
unsigned int faceSize() const
Definition raviartthomassimplexinterpolation.hh:129
void build(unsigned int order)
Definition raviartthomassimplexinterpolation.hh:148
OrthonormalBasisFactory< dimension, Field > TestBasisFactory
Definition raviartthomassimplexinterpolation.hh:105
TestFaceBasis * testFaceBasis(unsigned int f) const
Definition raviartthomassimplexinterpolation.hh:137
const Normal & normal(unsigned int f) const
Definition raviartthomassimplexinterpolation.hh:142
static const unsigned int dimension
Definition raviartthomassimplexinterpolation.hh:104
unsigned int order() const
Definition raviartthomassimplexinterpolation.hh:125
OrthonormalBasisFactory< dimension-1, Field > TestFaceBasisFactory
Definition raviartthomassimplexinterpolation.hh:107
unsigned int topologyId() const
Definition raviartthomassimplexinterpolation.hh:121
~RTL2InterpolationBuilder()
Definition raviartthomassimplexinterpolation.hh:114
FieldVector< Field, dimension > Normal
Definition raviartthomassimplexinterpolation.hh:109
Definition raviartthomassimplexinterpolation.hh:178
static void apply(const unsigned int order, std::vector< FaceStructure > &faceStructure)
Definition raviartthomassimplexinterpolation.hh:180
GenericGeometry::SubTopology< Topology, 1, face >::type FaceTopology
Definition raviartthomassimplexinterpolation.hh:179
Definition raviartthomassimplexinterpolation.hh:203
RaviartThomasL2Interpolation()
Definition raviartthomassimplexinterpolation.hh:210
void interpolate(typename Base::template Helper< Func, Container, type > &func) const
Definition raviartthomassimplexinterpolation.hh:269
unsigned int order() const
Definition raviartthomassimplexinterpolation.hh:230
void interpolate(const Function &function, std::vector< Fy > &coefficients) const
Definition raviartthomassimplexinterpolation.hh:216
RTL2InterpolationBuilder< dimension, Field > Builder
Definition raviartthomassimplexinterpolation.hh:209
F Field
Definition raviartthomassimplexinterpolation.hh:208
unsigned int size() const
Definition raviartthomassimplexinterpolation.hh:234
void interpolate(const Basis &basis, Matrix &matrix) const
Definition raviartthomassimplexinterpolation.hh:223
void build(unsigned int order)
Definition raviartthomassimplexinterpolation.hh:239
void setLocalKeys(std::vector< LocalKey > &keys) const
Definition raviartthomassimplexinterpolation.hh:251
Definition raviartthomassimplexinterpolation.hh:395
static const unsigned int dimension
Definition raviartthomassimplexinterpolation.hh:396
unsigned int Key
Definition raviartthomassimplexinterpolation.hh:397
const RaviartThomasL2Interpolation< dim, F > Object
Definition raviartthomassimplexinterpolation.hh:398
RaviartThomasL2InterpolationFactory< dim, F > Factory
Definition raviartthomassimplexinterpolation.hh:399
Definition interpolationhelper.hh:18
Definition interpolationhelper.hh:20
Definition polynomialbasis.hh:63
unsigned int size() const
Definition polynomialbasis.hh:108