vecellipt.cc 1.47 KB
Newer Older
Praetorius, Simon's avatar
Praetorius, Simon committed
1
#include <config.h>
2

3
4
#include <iostream>

5
#include <amdis/AMDiS.hpp>
6
#include <amdis/LocalOperators.hpp>
7
#include <amdis/ProblemStat.hpp>
8
9

using namespace AMDiS;
10
using namespace Dune::Functions::BasisFactory;
11
12
13

int main(int argc, char** argv)
{
14
  Environment env(argc, argv);
15

16
17
18
  using Grid = Dune::YaspGrid<GRIDDIM>;
  auto grid = MeshCreator<Grid>{"elliptMesh"}.create();
  ProblemStat prob("ellipt", *grid, power<2>(lagrange<2>()));
19
20
21
22
  prob.initialize(INIT_ALL);

  AdaptInfo adaptInfo("adapt");

23
  auto opL = makeOperator(tag::gradtest_gradtrial{}, 1.0);
24
25
  prob.addMatrixOperator(opL, 1, 1);

26
27
28
  auto opM = makeOperator(tag::test_trial{}, 1.0);
  prob.addMatrixOperator(std::ref(opM), 0, 0);
  prob.addMatrixOperator(std::ref(opM), 1, 0);
29

30
  auto opForce = makeOperator(tag::test{}, [](auto const& x) { return -1.0; }, 0);
31
32
33
34
35
36
37
38
  prob.addVectorOperator(opForce, 0);


  // 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, 1, 1, dbcValues);

39
  prob.buildAfterAdapt(adaptInfo, Flag(0));
40

41
#ifdef DEBUG_MTL
42
43
44
  // write matrix to file
  if (Parameters::get<int>("elliptMesh->global refinements").value_or(0) < 4) {
    mtl::io::matrix_market_ostream out("matrix.mtx");
Praetorius, Simon's avatar
Praetorius, Simon committed
45
    out << prob.systemMatrix().matrix();
46

Praetorius, Simon's avatar
Praetorius, Simon committed
47
    std::cout << prob.systemMatrix().matrix() << '\n';
48
  }
49
#endif
50
51

  prob.solve(adaptInfo);
Praetorius, Simon's avatar
Praetorius, Simon committed
52
  prob.writeFiles(adaptInfo);
53
54
55

  return 0;
}