elementdata.hh 2.72 KB
Newer Older
Peter Bastian's avatar
Peter Bastian committed
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
#ifndef __DUNE_GRID_HOWTO_ELEMENT_DATA_HH
#define __DUNE_GRID_HOWTO_ELEMENT_DATA_HH

6 7
#include <dune/grid/common/mcmgmapper.hh>
#include <dune/grid/io/file/vtk/vtkwriter.hh>
Peter Bastian's avatar
Peter Bastian committed
8
#if HAVE_GRAPE
9
#include <dune/grid/io/visual/grapedatadisplay.hh>
Peter Bastian's avatar
Peter Bastian committed
10 11 12 13 14 15 16
#endif

// demonstrate attaching data to elements
template<class G, class F>
void elementdata (const G& grid, const F& f)
{
  // the usual stuff
17
  //const int dim = G::dimension;
Peter Bastian's avatar
Peter Bastian committed
18 19
  const int dimworld = G::dimensionworld;
  typedef typename G::ctype ct;
20 21
  typedef typename G::LeafGridView GridView;
  typedef typename GridView::template Codim<0>::Iterator ElementLeafIterator;
22
  typedef typename ElementLeafIterator::Entity::Geometry LeafGeometry;
23 24

  // get grid view on leaf part
25
  GridView gridView = grid.leafGridView();
Peter Bastian's avatar
Peter Bastian committed
26 27

  // make a mapper for codim 0 entities in the leaf grid
28 29
  Dune::LeafMultipleCodimMultipleGeomTypeMapper<G>
  mapper(grid, Dune::mcmgElementLayout());             /*@\label{edh:mapper}@*/
Peter Bastian's avatar
Peter Bastian committed
30 31

  // allocate a vector for the data
32
  std::vector<double> c(mapper.size());                /*@\label{edh:c}@*/
Peter Bastian's avatar
Peter Bastian committed
33

34
  // iterate through all entities of codim 0 at the leaves
35 36
  for (ElementLeafIterator it = gridView.template begin<0>(); /*@\label{edh:loop0}@*/
       it!=gridView.template end<0>(); ++it)
Peter Bastian's avatar
Peter Bastian committed
37
  {
38 39
    // cell geometry
    const LeafGeometry geo = it->geometry();
Peter Bastian's avatar
Peter Bastian committed
40 41

    // get global coordinate of cell center
42
    Dune::FieldVector<ct,dimworld> global = geo.center();
Peter Bastian's avatar
Peter Bastian committed
43 44

    // evaluate functor and store value
45
    c[mapper.index(*it)] = f(global);                  /*@\label{edh:feval}@*/
46
  }                                                    /*@\label{edh:loop1}@*/
Peter Bastian's avatar
Peter Bastian committed
47 48

  // generate a VTK file
49
  // Dune::LeafP0Function<G,double> cc(grid,c);
50
  Dune::VTKWriter<typename G::LeafGridView> vtkwriter(gridView); /*@\label{edh:vtk0}@*/
51
  vtkwriter.addCellData(c,"data");
52
  vtkwriter.write( "elementdata", Dune::VTK::appendedraw ); /*@\label{edh:vtk1}@*/
Peter Bastian's avatar
Peter Bastian committed
53 54

  // online visualization with Grape
55
#if HAVE_GRAPE                                         /*@\label{edh:grape0}@*/
56 57 58 59 60 61 62 63
  {
    const int polynomialOrder = 0; // we piecewise constant data
    const int dimRange = 1; // we have scalar data here
    // create instance of data display
    Dune::GrapeDataDisplay<G> grape(grid);
    // display data
    grape.displayVector("concentration", // name of data that appears in grape
                        c,  // data vector
64
                        gridView.indexSet(), // used index set
65 66 67
                        polynomialOrder, // polynomial order of data
                        dimRange); // dimRange of data
  }
68
#endif                                                 /*@\label{edh:grape1}@*/
Peter Bastian's avatar
Peter Bastian committed
69
}
70 71

#endif //__DUNE_GRID_HOWTO_ELEMENT_DATA_HH