vecellipt.cc 1.55 KB
Newer Older
1 2 3 4
#include <iostream>

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

using namespace AMDiS;

// 1 component with polynomial degree 1
10 11
//using Grid = Dune::AlbertaGrid<AMDIS_DIM, AMDIS_DOW>;
using ElliptParam   = YaspGridBasis<AMDIS_DIM, 2,2>;
12 13 14 15 16 17 18 19 20 21 22
using ElliptProblem = ProblemStat<ElliptParam>;

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

  ElliptProblem prob("ellipt");
  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 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56
  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);

  *prob.getSolution() = 0.0; // maybe not necessary

  prob.buildAfterCoarsen(adaptInfo, Flag(0));

  // write matrix to file
  if (Parameters::get<int>("elliptMesh->global refinements").value_or(0) < 4) {
    mtl::io::matrix_market_ostream out("matrix.mtx");
    out << prob.getSystemMatrix()->getMatrix();

    std::cout << prob.getSystemMatrix()->getMatrix() << '\n';
  }

  prob.solve(adaptInfo);
  prob.writeFiles(adaptInfo, true);

  AMDiS::finalize();
  return 0;
}