ellipt.cc 2.04 KB
Newer Older
Praetorius, Simon's avatar
Praetorius, Simon committed
1
2
3
4
5
6
7
8
9
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:

#include <iostream>

#include <dune/amdis/AMDiS.hpp>
#include <dune/amdis/AdaptInstationary.hpp>
#include <dune/amdis/ProblemInstat.hpp>
#include <dune/amdis/ProblemStat.hpp>
10
#include <dune/amdis/Terms.hpp>
11
#include <dune/amdis/common/Literals.hpp>
Praetorius, Simon's avatar
Praetorius, Simon committed
12
#include <dune/amdis/BoundaryElementIterator.hpp>
Praetorius, Simon's avatar
Praetorius, Simon committed
13

14
15
using namespace AMDiS;

Praetorius, Simon's avatar
Praetorius, Simon committed
16
// 1 component with polynomial degree 1
Praetorius, Simon's avatar
Praetorius, Simon committed
17
using ElliptParam   = DefaultTraitsMesh<Dune::AlbertaGrid<AMDIS_DIM, AMDIS_DOW>, AMDIS_DIM, AMDIS_DOW, 1>;
Praetorius, Simon's avatar
Praetorius, Simon committed
18
19
20
21
22
23
using ElliptProblem = ProblemStat<ElliptParam>;
using ElliptProblemInstat = ProblemInstat<ElliptParam>;

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

Praetorius, Simon's avatar
Praetorius, Simon committed
25
26
  ElliptProblem prob("ellipt");
  prob.initialize(INIT_ALL);
27

Praetorius, Simon's avatar
Praetorius, Simon committed
28
29
  ElliptProblemInstat probInstat("ellipt", prob);
  probInstat.initialize(INIT_UH_OLD);
30

Praetorius, Simon's avatar
Praetorius, Simon committed
31
  AdaptInfo adaptInfo("adapt");
32

Praetorius, Simon's avatar
Praetorius, Simon committed
33
  using Op = ElliptProblem::OperatorType;
34
35
  Op opL, opForce;

Praetorius, Simon's avatar
Praetorius, Simon committed
36
37
  opL.addSOT(1.0);
  prob.addMatrixOperator(opL, 0, 0);
38

Praetorius, Simon's avatar
Praetorius, Simon committed
39
40
  opForce.addZOT( eval([](auto const& x) { return -1.0; }) );
  prob.addVectorOperator(opForce, 0);
41
42


Praetorius, Simon's avatar
Praetorius, Simon committed
43
44
45
46
  // set boundary condition
  auto predicate = [](auto const& x){ return x[0] < 1.e-8 || x[1] < 1.e-8; }; // define boundary
  auto dbcValues = [](auto const& x){ return 0.0; }; // set value
  prob.addDirichletBC(predicate, 0, 0, dbcValues);
47

Praetorius, Simon's avatar
Praetorius, Simon committed
48
  *prob.getSolution() = 0.0; // maybe not necessary
49

Praetorius, Simon's avatar
Praetorius, Simon committed
50
  prob.buildAfterCoarsen(adaptInfo, Flag(0));
51

Praetorius, Simon's avatar
Praetorius, Simon committed
52
53
54
55
56
  // write matrix to file
  {
    mtl::io::matrix_market_ostream out("matrix.mtx");
    out << prob.getSystemMatrix()->getMatrix<0,0>();
  }
57

Praetorius, Simon's avatar
Praetorius, Simon committed
58
  prob.solve(adaptInfo);
59
  prob.writeFiles(adaptInfo, true);
60
61
62


  msg("----------------------------");
63
  for (auto const& b : boundary_elements(prob.leafGridView()))
64
65
66
    msg("b = ", b.geometry().center());

  msg("----------------------------");
Praetorius, Simon's avatar
Praetorius, Simon committed
67
  size_t i = 0;
68
  for (auto const& f : boundary_facets(prob.leafGridView())) {
69
    msg("f = ", f.geometry().center()/*, ", id = ", f.boundaryId()*/);
Praetorius, Simon's avatar
Praetorius, Simon committed
70
71
72
    if (i++ >= 10)
      break;
  }
73

Praetorius, Simon's avatar
Praetorius, Simon committed
74
75
  AMDiS::finalize();
  return 0;
76
}