ellipt.cc 3.45 KB
Newer Older
1
2
3
4
5
#include "AMDiS.h"

using namespace AMDiS;
using namespace std;

6
7
8
9
10
// ===========================================================================
// ===== function definitions ================================================
// ===========================================================================

/// Dirichlet boundary function
11
12
13
14
class G : public AbstractFunction<double, WorldVector<double> >
{
public:

15
16
17
  /// Implementation of AbstractFunction::operator().
  double operator()(const WorldVector<double>& x) const 
  {
18
19
    return exp(-10.0 * (x * x));
  }
20
21
};

22
/// RHS function
23
24
25
26
class F : public AbstractFunction<double, WorldVector<double> >
{
public:

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

29
30
31
  /// Implementation of AbstractFunction::operator().
  double operator()(const WorldVector<double>& x) const 
  {
32
33
34
    int dow = x.getSize();
    double r2 = (x * x);
    double ux = exp(-10.0 * r2);
35
36
    return -(400.0 * r2 - 20.0 * dow) * ux;
  }
37
38
};

39
40
41
42
// ===========================================================================
// ===== main program ========================================================
// ===========================================================================

43
44
45
46
47
48
49
int main(int argc, char* argv[])
{
  FUNCNAME("main");

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

50

51
52
53
  // ===== init parameters =====
  Parameters::init(true, argv[1]);

54

55
56
57
58
  // ===== create and init the scalar problem ===== 
  ProblemScal ellipt("ellipt");
  ellipt.initialize(INIT_ALL);

59

60
  // === create adapt info ===
61
62
  AdaptInfo adaptInfo("ellipt->adapt");

63
64

  // === create adapt ===
65
66
  AdaptStationary adapt("ellipt->adapt", ellipt, adaptInfo);

67
68
  
  // ===== create matrix operator =====
69
  Operator matrixOperator(ellipt.getFeSpace());
70
  matrixOperator.addSecondOrderTerm(new Laplace_SOT);
71
72
  ellipt.addMatrixOperator(matrixOperator);

73
74

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

80

81
82
83
84
  // ===== add boundary conditions =====
  ellipt.addDirichletBC(1, new G);


85
  // ===== start adaption loop =====
Thomas Witkowski's avatar
Thomas Witkowski committed
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
115
116
117
118
119
  //  adapt.adapt();



  // Speichert hier bei jedem 2D-Element die DOFs auf dem Element.
  std::vector<DegreeOfFreedom> localDofs(3);
  // Damit wir jeden DOF nur einmal aufrufen, speichern wir hier jeden
  // besuchten DOF.
  std::set<DegreeOfFreedom> dofsVisited;

  // Okay, einfach das Gitter traversieren.
  TraverseStack stack;
  ElInfo *elInfo = stack.traverseFirst(ellipt.getMesh(), -1, Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS);
  while (elInfo) {
    // Zuerst lassen wir uns die DOFs des aktuellen Elements in localDofs abspeichern.
    ellipt.getFeSpace()->getBasisFcts()->getLocalIndices(elInfo->getElement(), 
							 ellipt.getFeSpace()->getAdmin(),
							 localDofs);

    // Dann durchlaufen wir die drei DOFs (ist also fix für lineare Elemente in 2D!!!!)
    for (int i = 0; i < 3; i++) {
      // Falls der DOF noch nicht ausgegeben wurde.
      if (dofsVisited.count(localDofs[i]) == 0) {
	// Gib globalen DOF Index und seine Welt-Koordinaten aus.
	std::cout << "DOF " << localDofs[i] << ": " 
		  << elInfo->getCoord(i)[0] << " " << elInfo->getCoord(i)[1] << "\n";
	
	// Merke dir den DOF, dass er schon ausgegebe  wurde.
	dofsVisited.insert(localDofs[i]);
      }
    }

    elInfo = stack.traverseNext(elInfo);
  }
120
121
122
123

}