1 #ifndef DUNE_RAVIARTTHOMASINTERPOLATION_HH
2 #define DUNE_RAVIARTTHOMASINTERPOLATION_HH
7 #include <dune/common/exceptions.hh>
8 #include <dune/common/forloop.hh>
10 #include <dune/geometry/topologyfactory.hh>
11 #include <dune/geometry/referenceelements.hh>
12 #include <dune/geometry/genericgeometry/referenceelements.hh>
13 #include <dune/geometry/quadraturerules/gaussquadrature.hh>
26 template <
unsigned int dim >
27 struct RaviartThomasCoefficientsFactory;
28 template <
unsigned int dim,
class Field >
29 struct RaviartThomasL2InterpolationFactory;
36 template <
class Setter>
39 setter.setLocalKeys(localKey_);
45 return localKey_[ i ];
50 return localKey_.size();
54 std::vector< LocalKey > localKey_;
57 template <
unsigned int dim >
62 typedef unsigned int Key;
65 template <
unsigned int dim >
70 template <
class Topology>
74 if (! supports<Topology>(key) )
76 typename InterpolationFactory::Object *interpol
77 = InterpolationFactory::template create<Topology>(key);
79 InterpolationFactory::release(interpol);
82 template<
class Topology >
101 template <
unsigned int dim,
class Field>
109 typedef FieldVector<Field,dimension>
Normal;
116 TestBasisFactory::release(testBasis_);
117 for (
unsigned int i=0;i<faceStructure_.size();++i)
118 TestFaceBasisFactory::release(faceStructure_[i].basis_);
140 return faceStructure_[f].basis_;
144 return *(faceStructure_[f].normal_);
147 template <
class Topology>
151 topologyId_ = Topology::id;
153 testBasis_ = TestBasisFactory::template create<Topology>(order-1);
158 faceStructure_.reserve( faceSize_ );
159 ForLoop< Creator<Topology>::template GetCodim,0,size-1>::apply(order, faceStructure_ );
160 assert( faceStructure_.size() == faceSize_ );
168 : basis_(tfb), normal_(&n)
172 const Dune::FieldVector<Field,dimension> *normal_;
174 template <
class Topology >
177 template <
int face >
180 typedef typename GenericGeometry::SubTopology<Topology,1,face>::type
FaceTopology;
182 std::vector<FaceStructure> &faceStructure )
184 faceStructure.push_back(
186 TestFaceBasisFactory::template create<FaceTopology>(order),
187 GenericGeometry::ReferenceElement<Topology,Field>::integrationOuterNormal(face)
193 std::vector<FaceStructure> faceStructure_;
195 unsigned int topologyId_, order_, faceSize_;
200 template<
unsigned int dimension,
class F>
216 template<
class Function,
class Fy >
217 void interpolate (
const Function &
function, std::vector< Fy > &coefficients )
const
219 coefficients.resize(
size());
220 typename Base::template Helper<Function,std::vector<Fy>,
true> func(
function,coefficients );
223 template<
class Basis,
class Matrix >
226 matrix.resize(
size(), basis.size() );
239 template <
class Topology>
244 builder_.template build<Topology>(order_);
247 for (
unsigned int f=0;f<builder_.
faceSize();++f )
255 unsigned int row = 0;
256 for (
unsigned int f=0;f<builder_.
faceSize();++f)
263 for (
unsigned int i=0;i<builder_.
testBasis()->
size()*dimension;++i,++row)
265 assert( row ==
size() );
269 template<
class Func,
class Container,
bool type >
270 void interpolate (
typename Base::template Helper<Func,Container,type> &func )
const
272 const Dune::GeometryType geoType( builder_.
topologyId(), dimension );
274 std::vector< Field > testBasisVal;
276 for (
unsigned int i=0;i<
size();++i)
277 for (
unsigned int j=0;j<func.size();++j)
280 unsigned int row = 0;
283 typedef Dune::GenericGeometry::GaussQuadratureProvider< dimension-1,
Field >
284 FaceQuadratureProvider;
286 typedef Dune::GenericReferenceElements< Field, dimension > RefElements;
287 typedef Dune::GenericReferenceElement< Field, dimension > RefElement;
288 typedef typename RefElement::template Codim< 1 >::Mapping Mapping;
290 const RefElement &refElement = RefElements::general( geoType );
292 for (
unsigned int f=0;f<builder_.
faceSize();++f)
298 const Mapping &mapping = refElement.template mapping< 1 >( f );
299 const Dune::GeometryType subGeoType( mapping.type().id(), dimension-1 );
300 const typename FaceQuadratureProvider::Object *faceQuad = FaceQuadratureProvider::create( subGeoType, 2*order_+2 );
302 const unsigned int quadratureSize = faceQuad->size();
303 for(
unsigned int qi = 0; qi < quadratureSize; ++qi )
305 builder_.
testFaceBasis(f)->template evaluate<0>(faceQuad->position(qi),testBasisVal);
306 fillBnd( row, testBasisVal,
307 func.evaluate( mapping.global( faceQuad->position(qi) ) ),
308 builder_.
normal(f), faceQuad->weight(qi),
312 FaceQuadratureProvider::release( faceQuad );
321 typedef Dune::GenericGeometry::GaussQuadratureProvider< dimension, Field > QuadratureProvider;
322 const typename QuadratureProvider::Object *elemQuad = QuadratureProvider::create( geoType, 2*order_+1 );
324 const unsigned int quadratureSize = elemQuad->size();
325 for(
unsigned int qi = 0; qi < quadratureSize; ++qi )
327 builder_.
testBasis()->template evaluate<0>(elemQuad->position(qi),testBasisVal);
328 fillInterior( row, testBasisVal,
329 func.evaluate(elemQuad->position(qi)),
330 elemQuad->weight(qi),
334 QuadratureProvider::release( elemQuad );
343 template <
class MVal,
class RTVal,
class Matrix>
344 void fillBnd (
unsigned int startRow,
347 const FieldVector<Field,dimension> &normal,
349 Matrix &matrix)
const
351 const unsigned int endRow = startRow+mVal.size();
352 typename RTVal::const_iterator rtiter = rtVal.begin();
353 for (
unsigned int col = 0; col < rtVal.size() ; ++rtiter,++col)
355 Field cFactor = (*rtiter)*normal;
356 typename MVal::const_iterator miter = mVal.begin();
357 for (
unsigned int row = startRow;
358 row!=endRow; ++miter, ++row )
360 matrix.add(row,col, (weight*cFactor)*(*miter) );
362 assert( miter == mVal.end() );
365 template <
class MVal,
class RTVal,
class Matrix>
366 void fillInterior (
unsigned int startRow,
370 Matrix &matrix)
const
372 const unsigned int endRow = startRow+mVal.size()*dimension;
373 typename RTVal::const_iterator rtiter = rtVal.begin();
374 for (
unsigned int col = 0; col < rtVal.size() ; ++rtiter,++col)
376 typename MVal::const_iterator miter = mVal.begin();
377 for (
unsigned int row = startRow;
378 row!=endRow; ++miter,row+=dimension )
380 for (
unsigned int i=0;i<dimension;++i)
382 matrix.add(row+i,col, (weight*(*miter))*(*rtiter)[i] );
385 assert( miter == mVal.end() );
394 template <
unsigned int dim,
class F >
397 static const unsigned int dimension = dim;
402 template <
unsigned int dim,
class Field >
404 public TopologyFactory< RaviartThomasL2InterpolationFactoryTraits<dim,Field> >
410 template <
class Topology>
413 if ( !supports<Topology>(key) )
416 interpol->template build<Topology>(key);
419 template<
class Topology >
426 #endif // DUNE_RAVIARTTHOMASINTERPOLATION_HH