periodic.cc 2.73 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif
#include <iostream>

#if HAVE_DUNE_SPGRID
#include <dune/grid/spgrid.hh>
#endif

#include <amdis/AMDiS.hpp>
#include <amdis/Integrate.hpp>
#include <amdis/ProblemStat.hpp>
13
#include <amdis/LocalOperators.hpp>
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
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
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
#include <amdis/common/Literals.hpp>

using namespace AMDiS;


template <class Grid>
void print(Grid const& grid)
{
  auto const& indexSet = grid.leafIndexSet();
  std::cout << "vertices:\n";
  for (auto const& v : vertices(grid.leafGridView())) {
    auto coord = v.geometry().corner(0);
    std::cout << indexSet.index(v) << ": [" << coord[0] << ", " << coord[1] << "]\n";
  }
  std::cout << "\n";

  std::cout << "elements:\n";
  for (auto const& e : elements(grid.leafGridView())) {
    std::cout << indexSet.index(e) << ": {" << indexSet.subIndex(e,0,2);
    for (unsigned int i = 1; i < e.subEntities(2); ++i) {
      std::cout << ", " << indexSet.subIndex(e,i,2);
     }
     std::cout << "}\n";
  }
  std::cout << "\n";

  std::cout << "boundary-intersections:\n";
  for (auto const& e : elements(grid.leafGridView())) {
    for (auto const& inter : intersections(grid.leafGridView(), e)) {
      if (!inter.boundary())
        continue;

      std::cout << inter.boundarySegmentIndex() << ": [" << inter.geometry().center() << "] {" << indexSet.index(inter.inside());
      if (inter.neighbor())
        std::cout << ", " << indexSet.index(inter.outside());
      std::cout << "}\n";
    }
  }
}

template <class Grid>
void run(Grid& grid)
{
  using Traits = LagrangeBasis<typename Grid::LeafGridView, 1>;

  ProblemStat<Traits> prob("ellipt", grid);
  prob.initialize(INIT_ALL);
  prob.boundaryManager().setBoxBoundary({-1,-1,1,1});

  print(grid);

  using BC = PeriodicBC<FieldVector<double,2>, typename Traits::GlobalBasis::MultiIndex>;
  BC periodicBC(prob.boundaryManagerPtr(),-1,{{{1.0,0.0}, {0.0,1.0}}, {1.0, 0.0}});

  periodicBC.init(prob.globalBasis(), prob.globalBasis());
  std::cout << "periodicNodes:\n";
  for (std::size_t i = 0; i < periodicBC.periodicNodes().size(); ++i)
    std::cout << i << ": " << periodicBC.periodicNodes()[i] << "\n";
  std::cout << "\n";

  std::cout << "associations:\n";
  for (auto const& a : periodicBC.associations())
    std::cout << a.first << " -> " << a.second << "\n";
  std::cout << "\n";
}

int main(int argc, char** argv)
{
  AMDiS::init(argc, argv);

#if HAVE_DUNE_SPGRID
  Dune::SPCube<double,2> cube({0.0,0.0},{1.0,1.0});
  Dune::SPDomain<double,2> domain({cube}, Dune::SPTopology<2>(0b01));
  Dune::SPGrid<double,2> grid1(domain, Dune::SPMultiIndex<2>({2,2}));
  run(grid1);
#endif

  // NOTE: 'periodic' YaspGrid not supported yet, but a simple YaspGrid can be made periodic
  Dune::YaspGrid<2> grid2({1.0,1.0}, {2,2});
  run(grid2);

  AMDiS::finalize();
  return 0;
}