ElliptProject.cpp 2.15 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
#include "ElliptProject.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 dow = x.getSize();
    double r2 = (x * x);
    double ux = exp(-10.0 * r2);
    return -(400.0 * r2 - 20.0 * dow) * ux;
  }
};

Elliptdemo::Elliptdemo():
  ellipt("ellipt"),
  adaptInfo(NULL),
  adapt(NULL),
  matrixOperator(NULL),
  rhsOperator(NULL)
{}

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

void Elliptdemo::create(std::string filename) {
    // ===== init parameters =====
  Parameters::init(true, filename);


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


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


  // === create adapt ===
  adapt = new AdaptStationary("ellipt->adapt", &ellipt, adaptInfo);

  
  // ===== create matrix operator =====
  matrixOperator = new Operator(ellipt.getFeSpace());
  matrixOperator->addSecondOrderTerm(new Laplace_SOT);
  ellipt.addMatrixOperator(matrixOperator);


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


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


}

int Elliptdemo::solve(SolutionInformation& info) {
  int retCode = adapt->adapt();
  info.dofVec = ellipt.getSolution();
  return retCode;
}

void createProjectList(ProjectList& list) {
}