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>
25 template <
unsigned int dim >
27 template <
unsigned int dim,
class Field >
35 template <
class Setter>
38 setter.setLocalKeys(localKey_);
44 return localKey_[ i ];
49 return localKey_.size();
53 std::vector< LocalKey > localKey_;
56 template <
unsigned int dim >
61 typedef unsigned int Key;
64 template <
unsigned int dim >
66 public TopologyFactory< RaviartThomasCoefficientsFactoryTraits<dim> >
69 template <
class Topology>
73 if (! supports<Topology>(key) )
75 typename InterpolationFactory::Object *interpol
76 = InterpolationFactory::template create<Topology>(key);
78 InterpolationFactory::release(interpol);
81 template<
class Topology >
100 template <
unsigned int dim,
class Field>
108 typedef FieldVector<Field,dimension>
Normal;
115 TestBasisFactory::release(testBasis_);
116 for (
unsigned int i=0;i<faceStructure_.size();++i)
117 TestFaceBasisFactory::release(faceStructure_[i].basis_);
139 return faceStructure_[f].basis_;
143 return *(faceStructure_[f].normal_);
146 template <
class Topology>
150 topologyId_ = Topology::id;
152 testBasis_ = TestBasisFactory::template create<Topology>(order-1);
157 faceStructure_.reserve( faceSize_ );
158 ForLoop< Creator<Topology>::template GetCodim,0,size-1>::apply(order, faceStructure_ );
159 assert( faceStructure_.size() == faceSize_ );
167 : basis_(tfb), normal_(&n)
171 const Dune::FieldVector<Field,dimension> *normal_;
173 template <
class Topology >
176 template <
int face >
179 typedef typename GenericGeometry::SubTopology<Topology,1,face>::type
FaceTopology;
181 std::vector<FaceStructure> &faceStructure )
183 faceStructure.push_back(
185 TestFaceBasisFactory::template create<FaceTopology>(order),
186 GenericGeometry::ReferenceElement<Topology,Field>::integrationOuterNormal(face)
192 std::vector<FaceStructure> faceStructure_;
194 unsigned int topologyId_, order_, faceSize_;
199 template<
unsigned int dimension,
class F>
215 template<
class Function,
class Fy >
216 void interpolate (
const Function &
function, std::vector< Fy > &coefficients )
const
218 coefficients.resize(
size());
222 template<
class Basis,
class Matrix >
225 matrix.resize(
size(), basis.size() );
238 template <
class Topology>
243 builder_.template build<Topology>(order_);
246 for (
unsigned int f=0;f<builder_.
faceSize();++f )
254 unsigned int row = 0;
255 for (
unsigned int f=0;f<builder_.
faceSize();++f)
262 for (
unsigned int i=0;i<builder_.
testBasis()->
size()*dimension;++i,++row)
264 assert( row ==
size() );
268 template<
class Func,
class Container,
bool type >
271 const Dune::GeometryType geoType( builder_.
topologyId(), dimension );
273 std::vector< Field > testBasisVal;
275 for (
unsigned int i=0;i<
size();++i)
276 for (
unsigned int j=0;j<func.size();++j)
279 unsigned int row = 0;
282 typedef Dune::GenericGeometry::GaussQuadratureProvider< dimension-1,
Field >
283 FaceQuadratureProvider;
285 typedef Dune::GenericReferenceElements< Field, dimension > RefElements;
286 typedef Dune::GenericReferenceElement< Field, dimension > RefElement;
287 typedef typename RefElement::template Codim< 1 >::Mapping Mapping;
289 const RefElement &refElement = RefElements::general( geoType );
291 for (
unsigned int f=0;f<builder_.
faceSize();++f)
297 const Mapping &mapping = refElement.template mapping< 1 >( f );
298 const Dune::GeometryType subGeoType( mapping.type().id(), dimension-1 );
299 const typename FaceQuadratureProvider::Object *faceQuad = FaceQuadratureProvider::create( subGeoType, 2*order_+2 );
301 const unsigned int quadratureSize = faceQuad->size();
302 for(
unsigned int qi = 0; qi < quadratureSize; ++qi )
304 builder_.
testFaceBasis(f)->template evaluate<0>(faceQuad->position(qi),testBasisVal);
305 fillBnd( row, testBasisVal,
306 func.evaluate( mapping.global( faceQuad->position(qi) ) ),
307 builder_.
normal(f), faceQuad->weight(qi),
311 FaceQuadratureProvider::release( faceQuad );
320 typedef Dune::GenericGeometry::GaussQuadratureProvider< dimension, Field > QuadratureProvider;
321 const typename QuadratureProvider::Object *elemQuad = QuadratureProvider::create( geoType, 2*order_+1 );
323 const unsigned int quadratureSize = elemQuad->size();
324 for(
unsigned int qi = 0; qi < quadratureSize; ++qi )
326 builder_.
testBasis()->template evaluate<0>(elemQuad->position(qi),testBasisVal);
327 fillInterior( row, testBasisVal,
328 func.evaluate(elemQuad->position(qi)),
329 elemQuad->weight(qi),
333 QuadratureProvider::release( elemQuad );
342 template <
class MVal,
class RTVal,
class Matrix>
343 void fillBnd (
unsigned int startRow,
346 const FieldVector<Field,dimension> &normal,
348 Matrix &matrix)
const
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)
354 Field cFactor = (*rtiter)*normal;
355 typename MVal::const_iterator miter = mVal.begin();
356 for (
unsigned int row = startRow;
357 row!=endRow; ++miter, ++row )
359 matrix.add(row,col, (weight*cFactor)*(*miter) );
361 assert( miter == mVal.end() );
364 template <
class MVal,
class RTVal,
class Matrix>
365 void fillInterior (
unsigned int startRow,
369 Matrix &matrix)
const
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)
375 typename MVal::const_iterator miter = mVal.begin();
376 for (
unsigned int row = startRow;
377 row!=endRow; ++miter,row+=dimension )
379 for (
unsigned int i=0;i<dimension;++i)
381 matrix.add(row+i,col, (weight*(*miter))*(*rtiter)[i] );
384 assert( miter == mVal.end() );
393 template <
unsigned int dim,
class F >
396 static const unsigned int dimension = dim;
401 template <
unsigned int dim,
class Field >
403 public TopologyFactory< RaviartThomasL2InterpolationFactoryTraits<dim,Field> >
409 template <
class Topology>
412 if ( !supports<Topology>(key) )
415 interpol->template build<Topology>(key);
418 template<
class Topology >
425 #endif // DUNE_RAVIARTTHOMASINTERPOLATION_HH