neumann.cc 2.12 KB
Newer Older
1
2
3
4
5
6
7
8
#include "AMDiS.h"

using namespace AMDiS;
using namespace std;

class N : public AbstractFunction<double, WorldVector<double> >
{
public:
9
10
  double operator()(const WorldVector<double>& x) const 
  {
11
12
    return 1.0;
  }
13
14
};

15
/// Dirichlet boundary function
16
17
18
class G : public AbstractFunction<double, WorldVector<double> >
{
public:
19
20
21
  /// Implementation of AbstractFunction::operator().
  double operator()(const WorldVector<double>& x) const 
  {
22
23
    return exp(-10.0 * (x * x));
  }
24
25
};

26
/// RHS function
27
28
29
class F : public AbstractFunction<double, WorldVector<double> >
{
public:
30
  F(int degree) : AbstractFunction<double, WorldVector<double> >(degree) {}
31

32
33
34
  /// Implementation of AbstractFunction::operator().
  double operator()(const WorldVector<double>& x) const 
  {
35
36
37
    int dow = x.getSize();
    double r2 = (x*x);
    double ux = exp(-10.0*r2);
38
39
    return -(400.0*r2 - 20.0*dow)*ux;
  }
40
41
42
43
44
45
46
47
48
49
50
};

// // ===== main program // 
int main(int argc, char* argv[])
{
  FUNCNAME("main");

  // ===== check for init file =====
  TEST_EXIT(argc == 2)("usage: neumann initfile\n");

  // ===== init parameters =====
51
  Parameters::init(argv[1]);
52
53

  // ===== create and init the scalar problem ===== 
54
  ProblemStat neumann("neumann");
55
56
57
  neumann.initialize(INIT_ALL);

  // === create adapt info ===
58
  AdaptInfo *adaptInfo = new AdaptInfo("neumann->adapt", 1);
59
60

  // === create adapt ===
61
  AdaptStationary *adapt = new AdaptStationary("neumann->adapt",
62
63
64
65
					       &neumann,
					       adaptInfo);
  
  // ===== create matrix operator =====
66
  Operator matrixOperator(neumann.getFeSpace());
67
68
  matrixOperator.addSecondOrderTerm(new Simple_SOT);
  neumann.addMatrixOperator(&matrixOperator, 0, 0);
69
70

  // ===== create rhs operator =====
71
  int degree = neumann.getFeSpace()->getBasisFcts()->getDegree();
72
  Operator rhsOperator(neumann.getFeSpace());
73
  rhsOperator.addZeroOrderTerm(new CoordsAtQP_ZOT(new F(degree)));
74
  neumann.addVectorOperator(&rhsOperator, 0);
75

76
  // ===== add boundary conditions =====
77
78
  neumann.addNeumannBC(1, 0, 0, new N);
  neumann.addDirichletBC(2, 0, 0, new G);
79

80
81
82
83
84
85
86
  // ===== start adaption loop =====
  adapt->adapt();

  neumann.writeFiles(adaptInfo, true);
}