heat.cc

 
#include <config.h>
 
#include <iostream>
 
#include <amdis/AMDiS.hpp>
#include <amdis/AdaptInstationary.hpp>
#include <amdis/LocalOperators.hpp>
#include <amdis/ProblemInstat.hpp>
#include <amdis/ProblemStat.hpp>
#include <amdis/GridFunctions.hpp>
 
using namespace AMDiS;
 
// 1 component with polynomial degree 1
//using Grid = Dune::AlbertaGrid<GRIDDIM, WORLDDIM>;
using HeatParam   = YaspGridBasis<GRIDDIM, 2>;
using HeatProblem = ProblemStat<HeatParam>;
using HeatProblemInstat = ProblemInstat<HeatParam>;
 
int main(int argc, char** argv)
{
  Environment env(argc, argv);
 
  HeatProblem prob("heat");
  prob.initialize(INIT_ALL);
 
  HeatProblemInstat probInstat("heat", prob);
  probInstat.initialize(INIT_UH_OLD);
 
  AdaptInfo adaptInfo("adapt");
 
  auto invTau = std::ref(probInstat.invTau());
 
  // 1/τ*⟨u,v⟩
  auto opTimeLhs = makeOperator(tag::test_trial{}, invTau);
  prob.addMatrixOperator(opTimeLhs);
 
  // ⟨∇u,∇v⟩
  auto opL = makeOperator(tag::gradtest_gradtrial{}, 1.0);
  prob.addMatrixOperator(opL);
 
  // 1/τ*⟨u_old,v⟩
  auto opTimeRhs = makeOperator(tag::test{},
    invokeAtQP([invTau](double u) { return u * invTau.get(); }, prob.solution()), 2);
  prob.addVectorOperator(opTimeRhs);
 
  // (-1.0,v)
  auto opForce = makeOperator(tag::test{}, [](auto const& x) { return -1.0; }, 0);
  prob.addVectorOperator(opForce);
 
 
  // 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, dbcValues);
 
  AdaptInstationary adapt("adapt", prob, adaptInfo, probInstat, adaptInfo);
  adapt.adapt();
 
  return 0;
}