SHARE
TWEET

Untitled

a guest Dec 8th, 2019 71 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. #ifdef HAVE_CONFIG_H
  2. # include "config.h"
  3. #endif
  4. #include <iostream>
  5. #include <vector>
  6. #include <array>
  7.  
  8. #include <dune/common/exceptions.hh>
  9. #include <dune/common/parallel/mpihelper.hh> // An initializer of MPI
  10. #include <dune/grid/uggrid.hh>             // Koristimo UGGrid
  11. #include <dune/grid/io/file/gmshreader.hh> // GmshReader klasa
  12. #include <dune/grid/io/file/vtk.hh>
  13. #include <dune/geometry/quadraturerules.hh>
  14.  
  15. template <int dim>
  16. Dune::FieldVector<double, dim> q(Dune::FieldVector<double, dim>const & x)
  17. {
  18.     auto y = x;
  19.     y /= x.dot(x);
  20.     return y;
  21. }
  22.  
  23. template<typename GV>
  24. double integration(GV gv, int p)
  25. {
  26.     const int dim = 2;
  27.     double integral = 0.0;
  28.     for (auto const & element : elements(gv))
  29.     {
  30.         for (auto const & side : intersections(gv, element))
  31.         {
  32.             if (side.boundary())
  33.             {
  34.                 const auto side_geom = side.geometry();
  35.                 auto outer_normal = side.centerUnitOuterNormal();
  36.                 const auto &rule = Dune::QuadratureRules<double, dim - 1>::rule(side_geom.type(), p);
  37.  
  38.                 for (auto const &qpoint : rule)
  39.                 {
  40.                     auto fval = q(side_geom.global(qpoint.position()));;
  41.                     double weight = qpoint.weight();
  42.                     double detjac = side_geom.integrationElement(qpoint.position());
  43.                     integral += (fval * outer_normal) * weight * detjac;
  44.                 }
  45.             }
  46.         }
  47.     }
  48.  
  49.     return integral;
  50. }
  51.  
  52. int main(int argc, char** argv)
  53. {
  54.     const int dim = 2;
  55.    
  56.     Dune::MPIHelper::instance(argc, argv);
  57.     using GridType = Dune::UGGrid<dim>;
  58.     using GridView = GridType::LeafGridView;
  59.  
  60.     if (argc < 2)
  61.         std::cout << "Need msh file name!" << std::endl, std::exit(-1);
  62.  
  63.     GridType *grid = Dune::GmshReader<GridType>::read(argv[1]);
  64.     grid->loadBalance();
  65.    
  66.     GridView gv = grid->leafGridView();
  67.    
  68.     auto &set = gv.indexSet();
  69.     std::vector<double> q_data(dim * gv.size(dim));
  70.    
  71.     for (auto const & vertex : vertices(gv))
  72.         for (unsigned int i = 0; i < dim; ++i)
  73.             q_data[dim * set.index(vertex) + i] = q(vertex.geometry().center())[i];
  74.    
  75.     Dune::VTKWriter<GridView> writer(gv);
  76.     writer.addVertexData(q_data, "q", dim);
  77.     writer.write("flux");
  78.  
  79.     std::cout << "Flux = " << integration(gv, 5) << std::endl;
  80.     return 0;
  81. }
RAW Paste Data
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy. OK, I Understand
 
Top