mpcci1.cc 2.59 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
#include "AMDiS.h"
#include "MpCCIAdapter.h"
#include "MultiGridSolver.h"
#include "GSSmoother.h"

using namespace std;
using namespace AMDiS;

// ===== function definitions // 
/** \brief
 * Dirichlet boundary function
 */
class G : public AbstractFunction<double, WorldVector<double> >
{
public:
  MEMORY_MANAGED(G);

  /** \brief
   * Implementation of AbstractFunction::operator().
   */
21
22
23
  double operator()(const WorldVector<double>& x) const
    return exp(-10.0 * (x * x));
  }
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
};

/** \brief
 * RHS function
 */
class F : public AbstractFunction<double, WorldVector<double> >
{
public:
  MEMORY_MANAGED(F);

  F(int degree) : AbstractFunction<double, WorldVector<double> >(degree) {};

  /** \brief
   * Implementation of AbstractFunction::operator().
   */
39
  double operator()(const WorldVector<double>& x) const {
40
    int dim = x.getSize();
41
42
43
44
    double r2 = x * x;
    double ux = exp(-10.0 * r2);
    return -(400.0 * r2 - 20.0 * dim) * ux;
  }
45
46
47
48
49
};

class WaitingProblem : public ProblemScal
{
public:
50
51
52
  WaitingProblem(const char *name) 
    : ProblemScal(name) 
  {}
53
54
55
56
57
58

  void solve(AdaptInfo *adaptInfo) {
    FUNCNAME("WaitingProblem::solve()");

    static MultiGridSolver *mgSolver = NULL;

59
60
61
62
63
64
    if (!mgSolver) 
      mgSolver = NEW MultiGridSolver(feSpace_, 
				     NEW GSSmoother(1.0),
				     systemMatrix_,
				     solution_,
				     rhs_);
65
66

    mgSolver->solve();
67
  }
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
};

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

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

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

  // ===== create and init the scalar problem ===== 
  WaitingProblem mpcci1("mpcci1");

  mpcci1.initialize(INIT_ALL);

  // === create adapt info ===
  AdaptInfo *adaptInfo = NEW AdaptInfo("mpcci1->adapt", 1);

  // === create adapt ===
  AdaptStationary *adapt = NEW AdaptStationary("mpcci1->adapt",
					       &mpcci1,
					       adaptInfo);
  
  //mpcci1.setAdapt(adapt);

  // ===== add boundary conditions =====
  mpcci1.addDirichletBC(1, NEW G);

  // ===== create matrix operator =====
  Operator matrixOperator(Operator::MATRIX_OPERATOR, mpcci1.getFESpace());
  matrixOperator.addSecondOrderTerm(NEW Laplace_SOT);
  mpcci1.addMatrixOperator(&matrixOperator);

  // ===== create rhs operator =====
  int degree = mpcci1.getFESpace()->getBasisFcts()->getDegree();
  Operator rhsOperator(Operator::VECTOR_OPERATOR, mpcci1.getFESpace());
  rhsOperator.addZeroOrderTerm(NEW CoordsAtQP_ZOT(NEW F(degree)));
  mpcci1.addVectorOperator(&rhsOperator);

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