BallProject.cpp 2.72 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
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
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
#include "BallProject.h"
#include "AdaptStationary.h"
/// Dirichlet boundary function
class G : public AbstractFunction<double, WorldVector<double> >
{
public:
  
  /// Implementation of AbstractFunction::operator().
  double operator()(const WorldVector<double>& x) const 
  {
    return exp(-10.0 * (x * x));
  }
};

/// RHS function
class F : public AbstractFunction<double, WorldVector<double> >
{
public:
  F(int degree) : AbstractFunction<double, WorldVector<double> >(degree) {}

  /// Implementation of AbstractFunction::operator().
  double operator()(const WorldVector<double>& x) const 
  {
    int dim = x.getSize();
    double r2 = x * x;
    double ux = exp(-10.0*r2);
    return -(400.0 * r2 - 20.0 * dim) * ux;
  }
};

BallDemo::BallDemo():
  ball("ball"),
  ballCenter(NULL),
  adaptInfo(NULL),
  adapt(NULL),
  matrixOperator(NULL),
  rhsOperator(NULL)
{
}

BallDemo::~BallDemo() {
  if(matrixOperator != NULL)
    delete matrixOperator;
  if(rhsOperator != NULL)
    delete rhsOperator;
  if(ballCenter != NULL)
    delete ballCenter;
  if(adapt != NULL)
    delete adapt;
  if(adaptInfo != NULL)
    delete adaptInfo;
}

void BallDemo::create(std::string& filename)  {
  // ===== init parameters =====
  Parameters::init(false, filename);
  ballCenter = new WorldVector< double >();
  ballCenter->set(0.0);

  // ===== create projection =====
  new AMDiS::BallProject(1, 
		  BOUNDARY_PROJECTION, 
		  *ballCenter, 
		  1.0);

  // ===== create and init the scalar problem ===== 
  ball.initialize(INIT_ALL);

  // === create adapt info ===
  adaptInfo = new AdaptInfo("ball->adapt");

  // === create adapt ===
  adapt = new AdaptStationary("ball->adapt", &ball, adaptInfo);
  
  // ===== create matrix operator =====
  matrixOperator=new Operator(ball.getFeSpace());
  matrixOperator->addSecondOrderTerm(new Laplace_SOT);
  ball.addMatrixOperator(matrixOperator);

  // ===== create rhs operator =====
  int degree = ball.getFeSpace()->getBasisFcts()->getDegree();
  rhsOperator = new Operator(ball.getFeSpace());
  rhsOperator->addZeroOrderTerm(new CoordsAtQP_ZOT(new F(degree)));
  ball.addVectorOperator(rhsOperator);

  // ===== add boundary conditions =====
  ball.addDirichletBC(1, new G);


}

int BallDemo::solve(SolutionInformation& solinfo) {
  assert(adaptInfo!=NULL);
  assert(adapt!=NULL);
  assert(matrixOperator!=NULL);
  assert(rhsOperator!=NULL);

  int retCode = adapt->adapt();
  solinfo.dofVec = ball.getSolution();
  return retCode;
}

void createProjectList(ProjectList& list) {
  list.clear();

  Project* demo = new BallDemo();
  ProjectInfo ballInfo(demo, "init/ball.dat.2d", "../testdata/balldata_2d");
  list.push_back(ballInfo);
  ballInfo = ProjectInfo(demo, "init/ball.dat.3d", "../testdata/balldata_3d");
  list.push_back(ballInfo);
}