ellipt.cc 1.63 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

13
14
using namespace AMDiS;

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

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

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

Praetorius, Simon's avatar
Praetorius, Simon committed
27
  AdaptInfo adaptInfo("adapt");
28

29
  using Op = ElliptProblem::ElementOperator;
30
31
  Op opL, opForce;

Praetorius, Simon's avatar
Praetorius, Simon committed
32
  opL.addSOT(1.0);
33
  prob.addMatrixOperator(opL, _0, _0);
34

Praetorius, Simon's avatar
Praetorius, Simon committed
35
  opForce.addZOT([](auto const& x) { return -1.0; });
36
  prob.addVectorOperator(opForce, _0);
37
38


Praetorius, Simon's avatar
Praetorius, Simon committed
39
40
41
  // 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
42
  prob.addDirichletBC(predicate, _0, _0, dbcValues);
43

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

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

Praetorius, Simon's avatar
Praetorius, Simon committed
48
  // write matrix to file
49
  if (Parameters::get<int>("elliptMesh->global refinements").value_or(0) < 5) {
Praetorius, Simon's avatar
Praetorius, Simon committed
50
    mtl::io::matrix_market_ostream out("matrix.mtx");
51
    out << prob.getSystemMatrix()->getMatrix();
52

53
    std::cout << prob.getSystemMatrix()->getMatrix() << '\n';
Praetorius, Simon's avatar
Praetorius, Simon committed
54
  }
55

Praetorius, Simon's avatar
Praetorius, Simon committed
56
  prob.solve(adaptInfo);
57
  prob.writeFiles(adaptInfo, true);
58

Praetorius, Simon's avatar
Praetorius, Simon committed
59
60
  AMDiS::finalize();
  return 0;
61
}