parametric.cc 4.98 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
#include "AMDiS.h"

using namespace std;
using namespace AMDiS;

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

/** \brief
 * RHS function
 */
class F : public AbstractFunction<double, WorldVector<double> >
{
public:
  MEMORY_MANAGED(F);

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

  /** \brief
   * Implementation of AbstractFunction::operator().
   */
23
24
25
  double operator()(const WorldVector<double>& x) const {
    return -2.0 * x[0];
  }
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
};

/** \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. 
 */
class ParametricSphere : public ProblemScal 
{
public:
  /** \brief
   * constructor
   */
  ParametricSphere(const char *name) 
    : ProblemScal(name),
      parametric(NULL)
  {};

  /** \brief
   * destructor
   */
  ~ParametricSphere() {
48
    for (int i = 0; i < Global::getGeo(WORLD); i++) {
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
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
      DELETE parametricCoords[i];
    }
    DELETE parametric;
  }

  /** \brief
   * initialization of the base class and creation of \ref parametric.
   */
  void initialize(Flag initFlag,
		  ProblemScal *adoptProblem = NULL,
		  Flag adoptFlag = INIT_NOTHING)
  {
    ProblemScal::initialize(initFlag, adoptProblem, adoptFlag);

    // ===== create projection =====
    WorldVector<double> ballCenter;
    ballCenter.set(0.0);
    NEW BallProject(1, 
		    VOLUME_PROJECTION, 
		    ballCenter, 
		    1.0);

    // ===== create coords vector =====
    int i;
    for(i = 0; i < Global::getGeo(WORLD); i++) {
      parametricCoords[i] = NEW DOFVector<double>(this->getFESpace(), 
						  "parametric coords");
    }

    // ===== create parametric object =====
    parametric = NEW ParametricFirstOrder(&parametricCoords);

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

  /** \brief
   * Implementation of ProblemStatBase::buildBeforeCoarsen().
   */
  void buildBeforeCoarsen(AdaptInfo *adaptInfo, Flag flag) {
    FUNCNAME("ParametricSphere::buildAfterCoarsen()");
    MSG("calculation of parametric coordinates\n");
    int i, j;
    int preDOFs = this->getFESpace()->getAdmin()->getNumberOfPreDOFs(VERTEX);
    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);
    while(elInfo) {
      dof = elInfo->getElement()->getDOF();
      for(i = 0; i < dim + 1; i++) {
	vertexCoords = elInfo->getCoord(i);
	Projection::getProjection(1)->project(vertexCoords);
	for(j = 0; j < dow; j++) {
	  (*(parametricCoords[j]))[dof[i][preDOFs]] = vertexCoords[j];
	}
      }
      elInfo = stack.traverseNext(elInfo);
    }
    

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

protected:
  /** \brief
   * DOFVector storing parametric coordinates.
   */
  WorldVector<DOFVector<double>*> parametricCoords;

  /** \brief
   * Parametrizes vertex coordinates while mesh traversal.
   */
  ParametricFirstOrder *parametric;
};

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

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

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

  // ===== init parameters =====
  Parameters::init(false, argv[1]);

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

  // === create adapt info ===
  AdaptInfo *adaptInfo = NEW AdaptInfo("parametric->adapt", 1);

  // === create adapt ===
  AdaptStationary *adapt = NEW AdaptStationary("parametric->adapt",
					       &parametric,
					       adaptInfo);
  
  // ===== create matrix operator =====
  Operator matrixOperator(Operator::MATRIX_OPERATOR,
			  parametric.getFESpace());
  matrixOperator.addSecondOrderTerm(NEW Laplace_SOT);
  parametric.addMatrixOperator(&matrixOperator);

  // ===== create rhs operator =====
  Operator rhsOperator(Operator::VECTOR_OPERATOR,
		       parametric.getFESpace());

  int degree = parametric.getFESpace()->getBasisFcts()->getDegree();

  rhsOperator.addZeroOrderTerm(NEW CoordsAtQP_ZOT(NEW F(degree)));
  parametric.addVectorOperator(&rhsOperator);

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

  parametric.writeFiles(adaptInfo, true);

  return 0;
}