parametric.cc 4.64 KB
Newer Older
1
2
3
4
5
6
7
8
9
#include "AMDiS.h"

using namespace std;
using namespace AMDiS;

// ===========================================================================
// ===== function definitions ================================================
// ===========================================================================

10
/// RHS function
11
12
13
class F : public AbstractFunction<double, WorldVector<double> >
{
public:
14
  F(int degree) : AbstractFunction<double, WorldVector<double> >(degree) {}
15

16
17
18
  /// Implementation of AbstractFunction::operator().
  double operator()(const WorldVector<double>& x) const 
  {
19
20
    return -2.0 * x[0];
  }
21
22
23
24
25
26
27
};

/** \brief
 * Implementation of a scalar problem. In \ref buildBeforeCoarsen() parametric
 * coordinates for the vertices are filled in a DOFVector. This DOFVector is
 * used in \ref parametric to parametrize elements while mesh traversal. 
 */
28
class ParametricSphere : public ProblemStat
29
30
{
public:
31
  /// constructor
32
  ParametricSphere(const char *name) 
33
    : ProblemStat(name),
34
      parametric(NULL)
35
  {}
36

37
38
39
40
41
42
43
  /// destructor
  ~ParametricSphere() 
  {
    for (int i = 0; i < Global::getGeo(WORLD); i++)
      delete parametricCoords[i];

    delete parametric;
44
45
  }

46
  /// initialization of the base class and creation of \ref parametric.
47
  void initialize(Flag initFlag,
48
		  ProblemStat *adoptProblem = NULL,
49
50
		  Flag adoptFlag = INIT_NOTHING)
  {
51
    ProblemStat::initialize(initFlag, adoptProblem, adoptFlag);
52
53
54
55

    // ===== create projection =====
    WorldVector<double> ballCenter;
    ballCenter.set(0.0);
56
    new BallProject(1, 
57
58
59
60
61
		    VOLUME_PROJECTION, 
		    ballCenter, 
		    1.0);

    // ===== create coords vector =====
Thomas Witkowski's avatar
Thomas Witkowski committed
62
    for (int i = 0; i < Global::getGeo(WORLD); i++)
63
      parametricCoords[i] = new DOFVector<double>(this->getFeSpace(), 
64
65
66
						  "parametric coords");

    // ===== create parametric object =====
67
    parametric = new ParametricFirstOrder(&parametricCoords);
68
69
70

    // ===== enable parametric traverse =====
    this->getMesh()->setParametric(parametric);
71
72
  }

73
74
75
76

  /** \brief
   * Implementation of ProblemStatBase::buildBeforeCoarsen().
   */
77
78
  void buildBeforeCoarsen(AdaptInfo *adaptInfo, Flag flag) 
  {
79
80
    FUNCNAME("ParametricSphere::buildAfterCoarsen()");
    MSG("calculation of parametric coordinates\n");
Marcel Schiffel's avatar
Marcel Schiffel committed
81
    int preDOFs = this->getFeSpace()->getAdmin()->getNumberOfPreDofs(VERTEX);
82
83
84
85
86
87
88
89
90
91
92
93
94
    int dim = this->getMesh()->getDim();
    int dow = Global::getGeo(WORLD);
    WorldVector<double> vertexCoords;
    const DegreeOfFreedom **dof;
    DegreeOfFreedom vertexIndex;

    // ===== disable parametric traverse =====
    this->getMesh()->setParametric(NULL);

    TraverseStack stack;
    ElInfo *elInfo = NULL;
    elInfo = stack.traverseFirst(this->getMesh(), -1, 
				 Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS);
95
    while (elInfo) {
96
      dof = elInfo->getElement()->getDof();
97
      for (int i = 0; i < dim + 1; i++) {
98
99
	vertexCoords = elInfo->getCoord(i);
	Projection::getProjection(1)->project(vertexCoords);
100
	for (int j = 0; j < dow; j++)
101
102
103
104
105
106
107
	  (*(parametricCoords[j]))[dof[i][preDOFs]] = vertexCoords[j];
      }
      elInfo = stack.traverseNext(elInfo);
    }
    
    // ===== enable parametric traverse =====
    this->getMesh()->setParametric(parametric);
108
  }
109

110

111
protected:
112
  /// DOFVector storing parametric coordinates.
113
114
  WorldVector<DOFVector<double>*> parametricCoords;

115
  /// Parametrizes vertex coordinates while mesh traversal.
116
117
118
119
120
121
122
123
124
125
126
  ParametricFirstOrder *parametric;
};

// ===========================================================================
// ===== main program ========================================================
// ===========================================================================

int main(int argc, char* argv[])
{
  FUNCNAME("parametric main");

127
  AMDiS::init(argc, argv);
128
129
130
131
132
133

  // ===== create and init the scalar problem ===== 
  ParametricSphere parametric("parametric");
  parametric.initialize(INIT_ALL);

  // === create adapt info ===
134
  AdaptInfo *adaptInfo = new AdaptInfo("parametric->adapt", parametric.getNumComponents());
135
136

  // === create adapt ===
137
  AdaptStationary *adapt = new AdaptStationary("parametric->adapt",
138
139
140
141
					       &parametric,
					       adaptInfo);
  
  // ===== create matrix operator =====
142
  Operator matrixOperator(parametric.getFeSpace());
143
  matrixOperator.addTerm(new Simple_SOT);
144
  parametric.addMatrixOperator(&matrixOperator, 0, 0);
145
146

  // ===== create rhs operator =====
147
  Operator rhsOperator(parametric.getFeSpace());
148

149
  int degree = parametric.getFeSpace()->getBasisFcts()->getDegree();
150

151
  rhsOperator.addTerm(new CoordsAtQP_ZOT(new F(degree)));
152
  parametric.addVectorOperator(&rhsOperator, 0);
153
154
155
156
157
158

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

  parametric.writeFiles(adaptInfo, true);

159
  AMDiS::finalize();
160
161
162
}