torus.cc 4.92 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
#include "AMDiS.h"

using namespace AMDiS;

class YRotation
{
public:
  static WorldVector<double>& rotate(WorldVector<double> &x, double angle) 
  {
    double x0 = x[0] * cos(angle) + x[2] * sin(angle);
    x[2] = -x[0] * sin(angle) + x[2] * cos(angle);
    x[0] = x0;
    return x;
Thomas Witkowski's avatar
Thomas Witkowski committed
14
  }
15
16
};

17
/// RHS function
18
19
20
21
22
23
class F : public AbstractFunction<double, WorldVector<double> >
{
public:
  F(int degree) 
    : AbstractFunction<double, WorldVector<double> >(degree),
      rotation(0.0)
24
  {}
25

26
27
28
  /// Implementation of AbstractFunction::operator().
  double operator()(const WorldVector<double>& x) const 
  {
29
30
    WorldVector<double> myX = x;
    YRotation::rotate(myX, -rotation);
31
32
    return -2.0 * myX[0];
  }
33

34
35
  void rotate(double r) 
  { 
36
    rotation += r; 
37
  }
38
39
40
41
42
43
44
45

private:
  double rotation;
};

class TorusProject : public Projection
{
public:
46
  /// Constructor.
47
48
49
50
51
52
53
  TorusProject(int id, 
	       ProjectionType type,
	       double radius1,
	       double radius2) 
    : Projection(id, type),
      radius1_(radius1),
      radius2_(radius2)
54
  {}
55

56
57
  /// Destructor.
  virtual ~TorusProject() {}
58

59
60
61
  /// Implementation of Projection::project();
  void project(WorldVector<double> &x) 
  {
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78

    WorldVector<double> xPlane = x;
    xPlane[2] = 0.0;

    double norm = sqrt(xPlane*xPlane);
    TEST_EXIT(norm != 0.0)("can't project vector x\n");
    
    WorldVector<double> center = xPlane;
    center *= radius1_ / norm;

    x -= center;

    norm = sqrt(x*x);
    TEST_EXIT(norm != 0.0)("can't project vector x\n");
    x *= radius2_/norm;

    x += center;
Thomas Witkowski's avatar
Thomas Witkowski committed
79
  }
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97

protected:
  double radius1_;

  double radius2_;
};

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

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

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

  // ===== create projection =====
98
99
  double r2 = (1.5 - 1.0 / sqrt(2.0)) / 2.0;
  double r1 = 1.0 / sqrt(2.0) + r2;
100

101
  new TorusProject(1, 
102
103
104
105
106
107
108
109
110
		   VOLUME_PROJECTION, 
		   r1,
		   r2);

  // ===== create and init the scalar problem ===== 
  ProblemScal torus("torus");
  torus.initialize(INIT_ALL);

  // === create adapt info ===
111
  AdaptInfo *adaptInfo = new AdaptInfo("torus->adapt");
112
113

  // === create adapt ===
114
  AdaptStationary *adapt = new AdaptStationary("torus->adapt",
115
116
117
118
119
					       &torus,
					       adaptInfo);
  
  // ===== create matrix operator =====
  Operator matrixOperator(Operator::MATRIX_OPERATOR,
120
			  torus.getFeSpace());
121
  matrixOperator.addSecondOrderTerm(new Laplace_SOT);
122
123
124
125
  torus.addMatrixOperator(&matrixOperator);

  // ===== create rhs operator =====
  Operator rhsOperator(Operator::VECTOR_OPERATOR,
126
		       torus.getFeSpace());
127

128
  int degree = torus.getFeSpace()->getBasisFcts()->getDegree();
129
130

  F f(degree);
131
  rhsOperator.addZeroOrderTerm(new CoordsAtQP_ZOT(&f));
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
  torus.addVectorOperator(&rhsOperator);

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

  torus.writeFiles(adaptInfo, true);


  double rotation = M_PI/3.0;
  int dim = torus.getMesh()->getDim();
  int dow = Global::getGeo(WORLD);

  DegreeOfFreedom dof;
  WorldVector<double> x;
 
147
  const FiniteElemSpace *feSpace = torus.getFeSpace();
148
149
  const BasisFunction *basFcts = feSpace->getBasisFcts();
  int numBasFcts = basFcts->getNumber();
150
  DegreeOfFreedom *localIndices = new DegreeOfFreedom[numBasFcts];
151
152
153
  DOFAdmin *admin = feSpace->getAdmin();
 
  WorldVector<DOFVector<double>*> parametricCoords;
154
155
  for (int i = 0; i < dow; i++)
    parametricCoords[i] = new DOFVector<double>(feSpace, "parametric coords");
156
157
158
159
160
161

  std::map<DegreeOfFreedom, bool> visited;

  TraverseStack stack;
  ElInfo *elInfo = stack.traverseFirst(torus.getMesh(), -1, 
				       Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS);
162
  while (elInfo) {
163
    basFcts->getLocalIndices(elInfo->getElement(), admin, localIndices);
164
    for (int i = 0; i < dim + 1; i++) {
165
166
167
      dof = localIndices[i];
      x = elInfo->getCoord(i);
      YRotation::rotate(x, rotation);
168
169
      if (!visited[dof]) {
	for (int j = 0; j < dow; j++)
170
	  (*(parametricCoords[j]))[dof] = x[j];
171

172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
	visited[dof] = true;
      }
    }
    elInfo = stack.traverseNext(elInfo);
  }

  ParametricFirstOrder parametric(&parametricCoords);
  torus.getMesh()->setParametric(&parametric);

  f.rotate(rotation);
  adaptInfo->reset();
  adapt->adapt();

  visited.clear();
  elInfo = stack.traverseFirst(torus.getMesh(), -1, 
			       Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS);
188
  while (elInfo) {
189
    basFcts->getLocalIndices(elInfo->getElement(), admin, localIndices);
190
    for (int i = 0; i < dim + 1; i++) {
191
192
193
      dof = localIndices[i];
      x = elInfo->getCoord(i);
      YRotation::rotate(x, rotation);
194
195
      if (!visited[dof]) {
	for (int j = 0; j < dow; j++)
196
	  (*(parametricCoords[j]))[dof] = x[j];
197

198
199
200
201
202
203
204
205
206
207
	visited[dof] = true;
      }
    }
    elInfo = stack.traverseNext(elInfo);
  }

  f.rotate(rotation);
  adaptInfo->reset();
  adapt->adapt();

208
209
  torus.writeFiles(adaptInfo, true);

Thomas Witkowski's avatar
Thomas Witkowski committed
210
  for (int i = 0; i < dow; i++)
211
    delete parametricCoords[i];
Thomas Witkowski's avatar
Thomas Witkowski committed
212

213
  delete [] localIndices;
214
215
216
}