vecellipt.cc 2.98 KB
Newer Older
1
2
3
4
5
6
7
8
9
#include "AMDiS.h"

using namespace std;
using namespace AMDiS;

// ===========================================================================
// ===== function definitions ================================================
// ===========================================================================

10
/// Dirichlet boundary function
11
12
13
class G : public AbstractFunction<double, WorldVector<double> >
{
public:
14
15
16
  /// Implementation of AbstractFunction::operator().
  double operator()(const WorldVector<double>& x) const 
  {
17
18
    return exp(-10.0 * (x * x));
  }
19
20
};

21
/// RHS function
22
23
24
class F : public AbstractFunction<double, WorldVector<double> >
{
public:
25
  F(int degree) : AbstractFunction<double, WorldVector<double> >(degree) {}
26

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

// ===========================================================================
// ===== main program ========================================================
// ===========================================================================

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

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

  // ===== init parameters =====
  Parameters::init(false, argv[1]);

  // ===== create and init the scalar problem ===== 
  ProblemVec vecellipt("vecellipt");
  vecellipt.initialize(INIT_ALL);

  // === create adapt info ===
56
  AdaptInfo *adaptInfo = new AdaptInfo("vecellipt->adapt",
57
58
59
				       vecellipt.getNumComponents());

  // === create adapt ===
60
  AdaptStationary *adapt = new AdaptStationary("vecellipt->adapt",
61
62
63
64
65
					       &vecellipt,
					       adaptInfo);
  

  // ===== add boundary conditions =====
66
  vecellipt.addDirichletBC(1, 0, 0, new G);
67
68
69
70
71

  // ===== create operators =====
  Operator matrixOperator00(Operator::MATRIX_OPERATOR,
			    vecellipt.getFESpace(0),
			    vecellipt.getFESpace(0));
72
  matrixOperator00.addSecondOrderTerm(new Laplace_SOT);
73
74
75
76
77
78
79
  vecellipt.addMatrixOperator(&matrixOperator00, 0, 0);

  Operator rhsOperator0(Operator::VECTOR_OPERATOR,
			vecellipt.getFESpace(0));

  int degree = vecellipt.getFESpace(0)->getBasisFcts()->getDegree();

80
  rhsOperator0.addZeroOrderTerm(new CoordsAtQP_ZOT(new F(degree)));
81
82
83
84
85
86
87
88
89
90
91

  vecellipt.addVectorOperator(&rhsOperator0, 0);

  Operator matrixOperator10(Operator::MATRIX_OPERATOR,
			    vecellipt.getFESpace(1),
			    vecellipt.getFESpace(0));

  Operator matrixOperator11(Operator::MATRIX_OPERATOR,
			    vecellipt.getFESpace(1),
			    vecellipt.getFESpace(1));

92
  matrixOperator10.addZeroOrderTerm(new Simple_ZOT);
93
94
95

  vecellipt.addMatrixOperator(&matrixOperator10, 1, 0);

96
  matrixOperator11.addZeroOrderTerm(new Simple_ZOT(-1.0));
97
98
99
100
101
102
103
104
105
106
107
108

  vecellipt.addMatrixOperator(&matrixOperator11, 1, 1);

  // ===== start adaption loop =====
  adapt->adapt();

  vecellipt.writeFiles(adaptInfo, true);

  return 0;
}