parevolve.hh 5.86 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:
Christian Engwer's avatar
Christian Engwer committed
3
#include <dune/grid/common/referenceelements.hh>
4

Peter Bastian's avatar
Peter Bastian committed
5 6 7 8
template<class G, class M, class V>
void parevolve (const G& grid, const M& mapper, V& c, double t, double& dt)
{
  // check data partitioning
9
  assert(grid.overlapSize(0)>0 || (grid.ghostSize(0)>0)); /*@\label{peh:assert}@*/
Peter Bastian's avatar
Peter Bastian committed
10 11 12 13 14 15 16 17

  // first we extract the dimensions of the grid
  const int dim = G::dimension;
  const int dimworld = G::dimensionworld;

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

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

Peter Bastian's avatar
Peter Bastian committed
21
  // iterator type
22 23
  typedef typename GridView::template Codim<0>::
  template Partition<Dune::All_Partition>::Iterator LeafIterator;       /*@\label{peh:pit}@*/
Peter Bastian's avatar
Peter Bastian committed
24 25

  // intersection iterator type
26
  typedef typename GridView::IntersectionIterator IntersectionIterator;
Peter Bastian's avatar
Peter Bastian committed
27

28 29 30
  // type of intersection
  typedef typename IntersectionIterator::Intersection Intersection;

Peter Bastian's avatar
Peter Bastian committed
31 32 33 34
  // entity pointer type
  typedef typename G::template Codim<0>::EntityPointer EntityPointer;

  // allocate a temporary vector for the update
Peter Bastian's avatar
Peter Bastian committed
35
  V update(c.size());
Christian Engwer's avatar
Christian Engwer committed
36
  for (typename V::size_type i=0; i<c.size(); i++) update[i] = 0;
Peter Bastian's avatar
Peter Bastian committed
37 38 39 40

  // initialize dt very large
  dt = 1E100;

41 42 43
  // get grid view instance on leaf grid
  GridView gridView = grid.leafView();

Peter Bastian's avatar
Peter Bastian committed
44 45
  // compute update vector and optimum dt in one grid traversal
  // iterate over all entities, but update is only used on interior entities
46 47
  LeafIterator endit = gridView.template end<0,Dune::All_Partition>(); /*@\label{peh:end}@*/
  for (LeafIterator it = gridView.template begin<0,Dune::All_Partition>(); it!=endit; ++it) /*@\label{peh:begin}@*/
Peter Bastian's avatar
Peter Bastian committed
48 49
  {
    // cell geometry type
50
    Dune::GeometryType gt = it->type();
Peter Bastian's avatar
Peter Bastian committed
51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70

    // cell center in reference element
    const Dune::FieldVector<ct,dim>&
    local = Dune::ReferenceElements<ct,dim>::general(gt).position(0,0);

    // cell center in global coordinates
    Dune::FieldVector<ct,dimworld>
    global = it->geometry().global(local);

    // cell volume, assume linear map here
    double volume = it->geometry().integrationElement(local)
                    *Dune::ReferenceElements<ct,dim>::general(gt).volume();

    // cell index
    int indexi = mapper.map(*it);

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

    // run through all intersections with neighbors and boundary
71 72
    const IntersectionIterator isend = gridView.iend(*it);
    for( IntersectionIterator is = gridView.ibegin(*it); is != isend; ++is )
Peter Bastian's avatar
Peter Bastian committed
73
    {
74 75
      const Intersection &intersection = *is;

Peter Bastian's avatar
Peter Bastian committed
76
      // get geometry type of face
Martin Nolte's avatar
Martin Nolte committed
77
      Dune::GeometryType gtf = intersection.type();
78 79 80

      const Dune::ReferenceElement< ct, dim-1 > &refElement
        = Dune::ReferenceElements< ct, dim-1 >::general( gtf );
Peter Bastian's avatar
Peter Bastian committed
81 82

      // center in face's reference element
83
      const Dune::FieldVector< ct, dim-1 > &facelocal = refElement.position( 0, 0 );
Peter Bastian's avatar
Peter Bastian committed
84

85
      // get normal vector scaled with volume
86 87 88
      Dune::FieldVector< ct, dimworld > integrationOuterNormal
        = intersection.integrationOuterNormal( facelocal );
      integrationOuterNormal *= refElement.volume();
Peter Bastian's avatar
Peter Bastian committed
89 90

      // center of face in global coordinates
91
      Dune::FieldVector< ct, dimworld > faceglobal
Martin Nolte's avatar
Martin Nolte committed
92
        = intersection.geometry().global( facelocal );
Peter Bastian's avatar
Peter Bastian committed
93 94 95 96 97 98 99 100 101 102 103

      // evaluate velocity at face center
      Dune::FieldVector<double,dim> velocity = u(faceglobal,t);

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

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

      // handle interior face
104
      if( intersection.neighbor() )
Peter Bastian's avatar
Peter Bastian committed
105 106
      {
        // access neighbor
107
        EntityPointer outside = intersection.outside();
108
        int indexj = mapper.map(*outside);
Peter Bastian's avatar
Peter Bastian committed
109

110 111 112
        const int insideLevel = it->level();
        const int outsideLevel = outside->level();

113
        // handle face from one side
114 115
        if( (insideLevel > outsideLevel)
            || ((insideLevel == outsideLevel) && (indexi < indexj)) )
Peter Bastian's avatar
Peter Bastian committed
116
        {
117
          // compute factor in neighbor
118
          Dune::GeometryType nbgt = outside->type();
119 120 121 122 123 124
          const Dune::FieldVector<ct,dim>&
          nblocal = Dune::ReferenceElements<ct,dim>::general(nbgt).position(0,0);
          double nbvolume = outside->geometry().integrationElement(nblocal)
                            *Dune::ReferenceElements<ct,dim>::general(nbgt).volume();
          double nbfactor = velocity*integrationOuterNormal/nbvolume;

125
          if( factor < 0 )       // inflow
Peter Bastian's avatar
Peter Bastian committed
126
          {
127 128 129
            update[indexi] -= c[indexj]*factor;
            update[indexj] += c[indexj]*nbfactor;
          }
130
          else       // outflow
131 132 133
          {
            update[indexi] -= c[indexi]*factor;
            update[indexj] += c[indexi]*nbfactor;
Peter Bastian's avatar
Peter Bastian committed
134 135 136 137 138
          }
        }
      }

      // handle boundary face
139 140 141
      if( intersection.boundary() )
      {
        if( factor < 0 )       // inflow, apply boundary condition
Peter Bastian's avatar
Peter Bastian committed
142
          update[indexi] -= b(faceglobal,t)*factor;
143
        else       // outflow
Peter Bastian's avatar
Peter Bastian committed
144
          update[indexi] -= c[indexi]*factor;
145 146
      }
    }       // end all intersections
Peter Bastian's avatar
Peter Bastian committed
147 148

    // compute dt restriction
149
    if (it->partitionType()==Dune::InteriorEntity)       /*@\label{peh:inter}@*/
Peter Bastian's avatar
Peter Bastian committed
150 151 152 153
      dt = std::min(dt,1.0/sumfactor);

  }       // end grid traversal

154 155
  // global min over all partitions
  dt = grid.comm().min(dt);                            /*@\label{peh:min}@*/
Peter Bastian's avatar
Peter Bastian committed
156 157 158 159
  // scale dt with safety factor
  dt *= 0.99;

  // exchange update
160
  VectorExchange<M,V> dh(mapper,update);               /*@\label{peh:dist0}@*/
161 162
  grid.template
  communicate<VectorExchange<M,V> >(dh,Dune::InteriorBorder_All_Interface,
163
                                    Dune::ForwardCommunication);                                       /*@\label{peh:dist1}@*/
Peter Bastian's avatar
Peter Bastian committed
164 165 166 167 168 169 170

  // update the concentration vector
  for (unsigned int i=0; i<c.size(); ++i)
    c[i] += dt*update[i];

  return;
}