NeumannProject.cpp 2.33 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
#include "NeumannProject.h"
#include "AdaptInfo.h"
#include "AdaptStationary.h"

class N : public AbstractFunction<double, WorldVector<double> >
{
public:
  double operator()(const WorldVector<double>& x) const 
  {
    return 1.0;
  }
};

/// 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;
  }
};

Neumanndemo::Neumanndemo():
  neumann("neumann"),
  adaptInfo(NULL),
  adapt(NULL),
  matrixOperator(NULL),
  rhsOperator(NULL)
{}

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

60
void Neumanndemo::create(string& filename) {
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
   // ===== init parameters =====
  Parameters::init(true, filename);

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

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

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

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

  // ===== add boundary conditions =====
  neumann.addNeumannBC(1, new N);
  neumann.addDirichletBC(2, new G);

}

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