initialize.hh 1.35 KB
Newer Older
1 2
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
3 4 5 6
#ifndef __DUNE_GRID_HOWTO_INITIALIZE_HH__
#define __DUNE_GRID_HOWTO_INITIALIZE_HH__

#include <dune/common/fvector.hh>
Peter Bastian's avatar
Peter Bastian committed
7

8 9 10 11 12
//! initialize the vector of unknowns with initial value
template<class G, class M, class V>
void initialize (const G& grid, const M& mapper, V& c)
{
  // first we extract the dimensions of the grid
13
  //const int dim = G::dimension;
14 15 16 17 18
  const int dimworld = G::dimensionworld;

  // type used for coordinates in the grid
  typedef typename G::ctype ct;

19 20 21
  // type of grid view on leaf part
  typedef typename G::LeafGridView GridView;

22
  // leaf iterator type
23 24
  typedef typename GridView::template Codim<0>::Iterator LeafIterator;

25 26 27
  // geometry type
  typedef typename LeafIterator::Entity::Geometry Geometry;

28
  // get grid view on leaf part
29
  GridView gridView = grid.leafGridView();
30 31

  // iterate through leaf grid an evaluate c0 at cell center
32 33
  LeafIterator endit = gridView.template end<0>();
  for (LeafIterator it = gridView.template begin<0>(); it!=endit; ++it)
34
  {
35
    // get geometry
36
    const Geometry geo = it->geometry();
37 38

    // get global coordinate of cell center
39
    Dune::FieldVector<ct,dimworld> global = geo.center();
40

Peter Bastian's avatar
Peter Bastian committed
41
    // initialize cell concentration
42
    c[mapper.index(*it)] = c0(global);
43 44
  }
}
45 46

#endif // __DUNE_GRID_HOWTO_INITIALIZE_HH__