evolve.hh 4.51 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_EVOLVE_HH__
#define __DUNE_GRID_HOWTO_EVOLVE_HH__

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

8 9 10 11 12 13 14 15 16
template<class G, class M, class V>
void evolve (const G& grid, const M& mapper, V& c, double t, double& dt)
{
  // first we extract the dimensions of the grid
  const int dimworld = G::dimensionworld;

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

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

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

23 24 25
  // leaf entity geometry
  typedef typename LeafIterator::Entity::Geometry LeafGeometry;

26
  // intersection iterator type
27
  typedef typename GridView::IntersectionIterator IntersectionIterator;
28

29 30 31
  // intersection geometry
  typedef typename IntersectionIterator::Intersection::Geometry IntersectionGeometry;

32 33
  // entity type
  typedef typename G::template Codim<0>::Entity Entity;
34

35
  // get grid view on leaf part
36
  GridView gridView = grid.leafGridView();
37

38
  // allocate a temporary vector for the update
39
  V update(c.size());                                  /*@\label{evh:update}@*/
Christian Engwer's avatar
Christian Engwer committed
40
  for (typename V::size_type i=0; i<c.size(); i++) update[i] = 0;
41 42 43 44 45

  // initialize dt very large
  dt = 1E100;

  // compute update vector and optimum dt in one grid traversal
46 47
  LeafIterator endit = gridView.template end<0>();     /*@\label{evh:loop0}@*/
  for (LeafIterator it = gridView.template begin<0>(); it!=endit; ++it)
48
  {
49
    // cell geometry
50
    const LeafGeometry geo = it->geometry();
51 52 53


    // cell volume, assume linear map here
54
    double volume = geo.volume();
55 56

    // cell index
57
    int indexi = mapper.index(*it);
58 59 60 61 62

    // variable to compute sum of positive factors
    double sumfactor = 0.0;

    // run through all intersections with neighbors and boundary
63 64
    IntersectionIterator isend = gridView.iend(*it);       /*@\label{evh:flux0}@*/
    for (IntersectionIterator is = gridView.ibegin(*it); is!=isend; ++is)
65 66
    {
      // get geometry type of face
67
      const IntersectionGeometry igeo = is->geometry();
68

69
      // get normal vector scaled with volume
70
      Dune::FieldVector<ct,dimworld> integrationOuterNormal
71
        = is->centerUnitOuterNormal();
72
      integrationOuterNormal *= igeo.volume();
73 74

      // center of face in global coordinates
75
      Dune::FieldVector<ct,dimworld> faceglobal = igeo.center();
76 77

      // evaluate velocity at face center
78
      Dune::FieldVector<double,dimworld> velocity = u(faceglobal,t);
79 80 81 82 83 84 85 86

      // compute factor occuring in flux formula
      double factor = velocity*integrationOuterNormal/volume;

      // for time step calculation
      if (factor>=0) sumfactor += factor;

      // handle interior face
87
      if (is->neighbor())             // "correct" version /*@\label{evh:neighbor}@*/
88
      {
Peter Bastian's avatar
Peter Bastian committed
89
        // access neighbor
90 91
        Entity outside = is->outside();
        int indexj = mapper.index(outside);
Peter Bastian's avatar
Peter Bastian committed
92

93
        // compute flux from one side only
94
        if (indexi<indexj)
Peter Bastian's avatar
Peter Bastian committed
95
        {
96
          // compute factor in neighbor
97
          const LeafGeometry nbgeo = outside.geometry();
98
          double nbvolume = nbgeo.volume();
99 100 101
          double nbfactor = velocity*integrationOuterNormal/nbvolume;

          if (factor<0)                         // inflow
Peter Bastian's avatar
Peter Bastian committed
102
          {
103 104 105 106 107 108 109
            update[indexi] -= c[indexj]*factor;
            update[indexj] += c[indexj]*nbfactor;
          }
          else                         // outflow
          {
            update[indexi] -= c[indexi]*factor;
            update[indexj] += c[indexi]*nbfactor;
Peter Bastian's avatar
Peter Bastian committed
110 111 112
          }
        }
      }
113 114

      // handle boundary face
115 116
      if (is->boundary())                               /*@\label{evh:bndry}@*/
      {
117 118 119 120
        if (factor<0)                 // inflow, apply boundary condition
          update[indexi] -= b(faceglobal,t)*factor;
        else                 // outflow
          update[indexi] -= c[indexi]*factor;
121
      }
122
    }             // end all intersections             /*@\label{evh:flux1}@*/
123 124

    // compute dt restriction
125
    dt = std::min(dt,1.0/sumfactor);                   /*@\label{evh:dt}@*/
126

127
  }       // end grid traversal                        /*@\label{evh:loop1}@*/
128 129

  // scale dt with safety factor
130
  dt *= 0.99;                                          /*@\label{evh:.99}@*/
131 132 133

  // update the concentration vector
  for (unsigned int i=0; i<c.size(); ++i)
134
    c[i] += dt*update[i];                              /*@\label{evh:updc}@*/
135 136 137

  return;
}
138 139

#endif //__DUNE_GRID_HOWTO_EVOLVE_HH__