1 #include <dune/common/version.hh>
2 #include<dune/geometry/referenceelements.hh>
4 template<
class G,
class M,
class V>
5 void evolve(
const G& grid,
const M& mapper, V& c,
double t,
double& dt)
8 const int dim = G::dimension;
9 const int dimworld = G::dimensionworld;
12 typedef typename G::ctype ct;
15 typedef typename G::LeafGridView GridView;
18 typedef typename GridView::template Codim<0>::Iterator LeafIterator;
21 typedef typename GridView::IntersectionIterator IntersectionIterator;
24 typedef typename G::template Codim<0>::EntityPointer EntityPointer;
27 GridView gridView = grid.leafGridView();
31 for (
typename V::size_type i=0; i<c.size(); i++) update[i] = 0;
37 LeafIterator endit = gridView.template end<0>();
38 for (LeafIterator it = gridView.template begin<0>(); it!=endit; ++it) {
40 Dune::GeometryType gt = it->type();
43 const Dune::FieldVector<ct,dim>&
44 #if DUNE_VERSION_NEWER(DUNE_GRID, 2, 3)
45 local = Dune::ReferenceElements<ct,dim>::general(gt).position(0,0);
47 local = Dune::GenericReferenceElements<ct,dim>::general(gt).position(0,0);
55 #if DUNE_VERSION_NEWER(DUNE_GRID, 2, 3)
56 double volume = it->geometry().integrationElement(local)
57 *Dune::ReferenceElements<ct,dim>::general(gt).volume();
59 double volume = it->geometry().integrationElement(local)
60 *Dune::GenericReferenceElements<ct,dim>::general(gt).volume();
64 #if DUNE_VERSION_NEWER(DUNE_GRID, 2, 4)
66 int indexi = mapper.index(*it);
69 int indexi = mapper.map(*it);
73 double sumfactor = 0.0;
75 IntersectionIterator isend = gridView.iend(*it);
76 for (IntersectionIterator is = gridView.ibegin(*it); is!=isend; ++is) {
80 Dune::GeometryType gtf = is->type();
83 const Dune::FieldVector<ct,dim-1>&
84 #if DUNE_VERSION_NEWER(DUNE_GRID, 2, 3)
85 facelocal = Dune::ReferenceElements<ct,dim-1>::general(gtf).position(0,0);
87 facelocal = Dune::GenericReferenceElements<ct,dim-1>::general(gtf).position(0,0);
91 Dune::FieldVector<ct,dimworld> integrationOuterNormal
92 = is->integrationOuterNormal(facelocal);
93 integrationOuterNormal
94 #if DUNE_VERSION_NEWER(DUNE_GRID, 2, 3)
95 *= Dune::ReferenceElements<ct,dim-1>::general(gtf).volume();
97 *= Dune::GenericReferenceElements<ct,dim-1>::general(gtf).volume();
101 Dune::FieldVector<ct,dimworld>
102 faceglobal = is->geometry().center();
106 Dune::FieldVector<double,dim> velocity =
u(faceglobal,t);
109 double factor = velocity*integrationOuterNormal/volume;
112 if (factor>=0) sumfactor += factor;
118 EntityPointer outside = is->outside();
119 #if DUNE_VERSION_NEWER(DUNE_GRID, 2, 4)
120 int indexj = mapper.index(*outside);
122 int indexj = mapper.map(*outside);
127 if (it->level()>outside->level() ||
128 (it->level()==outside->level() && indexi<indexj))
131 Dune::GeometryType nbgt = outside->type();
132 const Dune::FieldVector<ct,dim>&
133 #if DUNE_VERSION_NEWER(DUNE_GRID, 2, 3)
134 nblocal = Dune::ReferenceElements<ct,dim>::general(nbgt).position(0,0);
135 double nbvolume = outside->geometry().integrationElement(nblocal)
136 *Dune::ReferenceElements<ct,dim>::general(nbgt).volume();
138 nblocal = Dune::GenericReferenceElements<ct,dim>::general(nbgt).position(0,0);
139 double nbvolume = outside->geometry().integrationElement(nblocal)
140 *Dune::GenericReferenceElements<ct,dim>::general(nbgt).volume();
142 double nbfactor = velocity*integrationOuterNormal/nbvolume;
146 update[indexi] -= c[indexj]*factor;
147 update[indexj] += c[indexj]*nbfactor;
150 update[indexi] -= c[indexi]*factor;
151 update[indexj] += c[indexi]*nbfactor;
160 update[indexi] -=
b(faceglobal,t)*factor;
162 update[indexi] -= c[indexi]*factor;
166 dt = std::min(dt,1.0/sumfactor);
174 for (
unsigned int i=0; i<c.size(); ++i)
175 c[i] += dt*update[i];
Dune::FieldVector< double, dimworld > u(const Dune::FieldVector< ct, dimworld > &, double)
Definition: transportproblem2.hh:36
void evolve(const G &grid, const M &mapper, V &c, double t, double &dt)
Definition: evolve.hh:5
double b(const Dune::FieldVector< ct, dimworld > &x, double t)
Definition: transportproblem2.hh:25