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