heat.cc 2.18 KB
Newer Older
1
2
3
4
5
6
// -*- 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>
7
#include <dune/amdis/AdaptInstationary.hpp>
8
9
#include <dune/amdis/ProblemInstat.hpp>
#include <dune/amdis/ProblemStat.hpp>
10
#include <dune/amdis/Terms.hpp>
11
#include <dune/amdis/common/Literals.hpp>
12
13
14
15

using namespace AMDiS;


Praetorius, Simon's avatar
Praetorius, Simon committed
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
// class HeatParam
// {
//   template <class... Bs>
//   using FeSpaceTuple = std::tuple<Bs...>;
//   
//   template <class Mesh, int deg>
//   using Lagrange = Dune::Functions::PQkNodalBasis<typename Mesh::LeafGridView, deg>;
//   
// public:
//   static constexpr int dim = AMDIS_DIM;
//   static constexpr int dimworld = AMDIS_DOW;
//   static constexpr int nComponents = 1;
//   
//   // default values
//   using Mesh = Dune::UGGrid<dim>;
//   using FeSpaces = FeSpaceTuple<Lagrange<Mesh, 2>>;
// };
33
34
35


// 1 component with polynomial degree 1
Praetorius, Simon's avatar
Praetorius, Simon committed
36
using HeatParam   = DefaultTraitsMesh<Dune::AlbertaGrid<AMDIS_DIM, AMDIS_DOW>, AMDIS_DIM, AMDIS_DOW, 2>;
37
38
39
40
41
using HeatProblem = ProblemStat<HeatParam>;
using HeatProblemInstat = ProblemInstat<HeatParam>;

int main(int argc, char** argv)
{
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
  AMDiS::init(argc, argv);
  
  HeatProblem prob("heat");
  prob.initialize(INIT_ALL);
  
  HeatProblemInstat probInstat("heat", prob);
  probInstat.initialize(INIT_UH_OLD);
  
  AdaptInfo adaptInfo("adapt");
  
  using Op = HeatProblem::OperatorType;
  Op opTimeLhs, opL, opTimeRhs, opForce; 
  
  opTimeLhs.addZOT( constant(1.0) );
  prob.addMatrixOperator(opTimeLhs, 0, 0, probInstat.getInvTau());
  
  opL.addSOT( constant(1.0) );
  prob.addMatrixOperator(opL, 0, 0);
  
  opTimeRhs.addZOT( valueOf(prob.getSolution(0_c)) );
  prob.addVectorOperator(opTimeRhs, 0, probInstat.getInvTau());
  
  opForce.addZOT( eval([](auto const& x) { return -1.0; }) );
  prob.addVectorOperator(opForce, 0);
  
  
  // set boundary condition
  auto predicate = [](auto const& p){ return p[0] < 1.e-8 || p[1] < 1.e-8; };
  auto dbcValues = [](auto const& p){ return 0.0; };
  prob.addDirichletBC(predicate, 0, 0, dbcValues);
      
  *prob.getSolution() = 0.0;
  
  AdaptInstationary adapt("adapt", prob, adaptInfo, probInstat, adaptInfo);
  adapt.adapt();
  
  AMDiS::finalize();
  return 0;
80
}