adaptivefinitevolume.cc 3.99 KB
Newer Older
Peter Bastian's avatar
Peter Bastian committed
1 2 3 4 5 6
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#include "config.h"               // know what grids are present
#include <iostream>               // for input/output to shell
#include <fstream>                // for input/output to files
#include <vector>                 // STL vector class
Christian Engwer's avatar
Christian Engwer committed
7
#include <dune/grid/common/mcmgmapper.hh> // mapper class
Peter Bastian's avatar
Peter Bastian committed
8 9

#include "unitcube.hh"
Peter Bastian's avatar
Peter Bastian committed
10
#include "transportproblem2.hh"
Peter Bastian's avatar
Peter Bastian committed
11 12 13
#include "initialize.hh"
#include "evolve.hh"
#include "finitevolumeadapt.hh"
Peter Bastian's avatar
Peter Bastian committed
14
#include "vtkout.hh"
Peter Bastian's avatar
Peter Bastian committed
15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30

//===============================================================
// the time loop function working for all types of grids
//===============================================================

//! Parameter for mapper class
template<int dim>
struct P0Layout
{
  bool contains (int codim, Dune::GeometryType gt)
  {
    if (codim==0) return true;
    return false;
  }
};

Peter Bastian's avatar
Peter Bastian committed
31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67
template<class G>
void gnuplot (G& grid, std::vector<double>& c)
{
  // 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
  // such a type is exported by every grid implementation
  typedef typename G::ctype ct;

  // the grid has an iterator providing the access to
  // all elements (better codim 0 entities) which are leafs
  // of the refinement tree.
  // Note the use of the typename keyword and the traits class
  typedef typename G::template Codim<0>::LeafIterator ElementLeafIterator;

  // make a mapper for codim 0 entities in the leaf grid
  Dune::LeafMultipleCodimMultipleGeomTypeMapper<G,P0Layout>
  mapper(grid);

  // iterate through all entities of codim 0 at the leafs
  int count = 0;
  for (ElementLeafIterator it = grid.template leafbegin<0>();
       it!=grid.template leafend<0>(); ++it)
  {
    Dune::GeometryType gt = it->geometry().type();
    const Dune::FieldVector<ct,dim>&
    local = Dune::ReferenceElements<ct,dim>::general(gt).position(0,0);
    Dune::FieldVector<ct,dimworld>
    global = it->geometry().global(local);
    std::cout << global[0] << " " << c[mapper.map(*it)] << std::endl;
    count++;
  }
}


Peter Bastian's avatar
Peter Bastian committed
68 69 70 71
template<class G>
void timeloop (G& grid, double tend, int lmin, int lmax)
{
  // make a mapper for codim 0 entities in the leaf grid
Peter Bastian's avatar
Peter Bastian committed
72 73
  Dune::LeafMultipleCodimMultipleGeomTypeMapper<G,P0Layout>
  mapper(grid);
Peter Bastian's avatar
Peter Bastian committed
74 75 76 77 78 79 80 81

  // allocate a vector for the concentration
  std::vector<double> c(mapper.size());

  // initialize concentration with initial values
  initialize(grid,mapper,c);
  for (int i=grid.maxLevel(); i<lmax; i++)
  {
Peter Bastian's avatar
Peter Bastian committed
82
    if (grid.maxLevel()>=lmax) break;
Peter Bastian's avatar
Peter Bastian committed
83
    finitevolumeadapt(grid,mapper,c,lmin,lmax,0);
Peter Bastian's avatar
Peter Bastian committed
84 85
    initialize(grid,mapper,c);
  }
Peter Bastian's avatar
Peter Bastian committed
86
  vtkout(grid,c,"concentration",0);
Peter Bastian's avatar
Peter Bastian committed
87 88 89

  double dt, t=0;
  int k=0;
Peter Bastian's avatar
Peter Bastian committed
90
  int modulo=5;
Peter Bastian's avatar
Peter Bastian committed
91 92
  std::cout << "s=" << grid.size(0) << " k=" << k
            << " t=" << t << std::endl;
Peter Bastian's avatar
Peter Bastian committed
93 94 95 96
  while (t<tend)
  {
    k++;
    evolve(grid,mapper,c,t,dt);
Peter Bastian's avatar
Peter Bastian committed
97
    return;
Peter Bastian's avatar
Peter Bastian committed
98
    t += dt;
Peter Bastian's avatar
Peter Bastian committed
99
    if (k%modulo==0) vtkout(grid,c,"concentration",k/modulo);
Peter Bastian's avatar
Peter Bastian committed
100 101
    std::cout << "s=" << grid.size(0) << " k=" << k
              << " t=" << t << " dt=" << dt << std::endl;
Peter Bastian's avatar
Peter Bastian committed
102
    //            finitevolumeadapt(grid,mapper,c,lmin,lmax,k);
Peter Bastian's avatar
Peter Bastian committed
103
  }
Peter Bastian's avatar
Peter Bastian committed
104
  vtkout(grid,c,"concentration",k/modulo);
Peter Bastian's avatar
Peter Bastian committed
105 106 107 108 109 110 111 112 113 114 115 116
}

//===============================================================
// The main function creates objects and does the time loop
//===============================================================

int main (int argc , char ** argv)
{
#if HAVE_MPI
  MPI_Init(&argc,&argv);
#endif

Peter Bastian's avatar
Peter Bastian committed
117 118
  UnitCube<Dune::OneDGrid<1,1>,1> uc0;
  UnitCube<Dune::SGrid<1,1>,1> uc1;
Peter Bastian's avatar
Peter Bastian committed
119 120
  UnitCube<Dune::YaspGrid<2,2>,1> uc;
#if HAVE_UG
121
  UnitCube<Dune::UGGrid<3,3>,1> uc2;
Peter Bastian's avatar
Peter Bastian committed
122 123
#endif
#if HAVE_ALBERTA
Peter Bastian's avatar
update  
Peter Bastian committed
124
#if DUNE_PROBLEM_DIM==2
Peter Bastian's avatar
Peter Bastian committed
125
  UnitCube<Dune::AlbertaGrid<2,2>,1> uc3;
Peter Bastian's avatar
update  
Peter Bastian committed
126
#endif
Peter Bastian's avatar
Peter Bastian committed
127
#endif
Peter Bastian's avatar
Peter Bastian committed
128 129
  //    uc3.grid().globalRefine(8);
  //    timeloop(uc3.grid(),0.5,8,18);
130 131 132
  //  uc2.grid().globalRefine(3);
  //  timeloop(uc2.grid(),0.5,3,7);
  uc0.grid().globalRefine(4);
Peter Bastian's avatar
Peter Bastian committed
133
  timeloop(uc0.grid(),0.5,4,8);
Peter Bastian's avatar
Peter Bastian committed
134 135 136 137 138 139 140 141

#if HAVE_MPI
  MPI_Finalize();
#endif

  // done
  return 0;
}