HL_SignedDist.cc 9.99 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
//
// Software License for AMDiS
//
// Copyright (c) 2010 Dresden University of Technology 
// All rights reserved.
// Authors: Simon Vey, Thomas Witkowski et al.
//
// This file is part of AMDiS
//
// See also license.opensource.txt in the distribution.


13
14
15
16
17
#include "HL_SignedDist.h"

const Flag HL_SignedDist::VEL_EXT                = 0X01L;
const Flag HL_SignedDist::VEL_EXT_FROM_VEL_FIELD = 0X02L;

18
19
20
21
void HL_SignedDist::calcSignedDistFct(AdaptInfo *adaptInfo_,
				      const DOFVector<double> *lS_DOF_, 
				      DOFVector<double> *sD_DOF_,
				      ElementFunction<double> *elFct) 
22
23
24
25
26
{
  adaptInfo = adaptInfo_;

  if (elFct == NULL) {
    TEST_EXIT(lS_DOF_)("illegal level set function !\n");
27
    TEST_EXIT(lS_DOF_->getFeSpace() == sD_DOF_->getFeSpace())("DOF vectors do not have the same FE space!\n");
28
29
  }
  TEST_EXIT(sD_DOF_)("illegal DOF vector for signed distance function !\n");
30
  TEST_EXIT(sD_DOF_->getFeSpace()->getBasisFcts()->getDegree() == 1)("does it work for non-linear finite elements ?\n");
31
32
33
  
  lS_DOF = lS_DOF_;
  sD_DOF = sD_DOF_;
34
  feSpace = sD_DOF->getFeSpace();
35
36
37
38
39
40
41
42
43
44

  // ===== Initialization. =====
  initialize(elFct);
  
  // ===== Print level set function. =====
  if (lS_DOF) 
    printLevelSetFct();
  
  // ===== Create DOF vector to mark boundary vertices. =====
  if (bound_DOF) 
45
46
    delete bound_DOF;
  bound_DOF = new DOFVector<double>(feSpace, "bound_DOF");
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
  bound_DOF->set(0.0);
  
  // ===== Boundary vertex and boundary value initialization. =====
  initializeBoundary();
  
  // ===== Print boundary initialization. =====
  printBoundInitFct();
  
    // ===== Calculate distance function. =====
  HL_updateIteration();
  
  // ===== Transformation to signed distance function. =====
  if (elFct == NULL)
    setSign();
  
  // ===== Print calculated signed distance function. =====
  printSignedDistFct();
}

66
67

void HL_SignedDist::calcSignedDistFct(AdaptInfo *adaptInfo_, DOFVector<double> *lS_DOF_)
68
69
70
{
  TEST_EXIT(lS_DOF_)("illegal level set function lS_DOF_ !\n");

71
  sD_DOF = new DOFVector<double>(lS_DOF_->getFeSpace(), "sD_DOF");
72
73
74
75
76

  calcSignedDistFct(adaptInfo_, lS_DOF_, sD_DOF, NULL);

  *lS_DOF_ = *sD_DOF;

77
  delete sD_DOF;
78
79
}

80
81
82
83
84
85
86
87

void HL_SignedDist::calcVelocityExt(AdaptInfo *adaptInfo_,
				    DOFVector<double> *origVel_DOF_, 
				    DOFVector<double> *vel_DOF_,
				    const DOFVector<double> *lS_DOF_,
				    DOFVector<double> *sD_DOF_,
				    bool calcSDFct,
				    ElementFunction<double> *elFct)
88
89
90
91
92
93
94
95
96
97
{
  TEST_EXIT(velExtType.isSet(VEL_EXT))("illegal velocity extension type !\n");
  TEST_EXIT(velExt)("velExt not defined !\n");

  velExt->setVelocity(origVel_DOF_, vel_DOF_);

  velExt->printOrigVelDOF(adaptInfo_);

  if (calcSDFct || sD_DOF_ != NULL) {
    calcSignedDistFct(adaptInfo_, lS_DOF_, sD_DOF_, elFct);
98
  } else {
99
100
101

    TEST_EXIT(vel_DOF_)("illegal level set function lS_DOF_ !\n");

102
    sD_DOF = new DOFVector<double>(vel_DOF_->getFeSpace(), "sD_DOF");
103
104
105

    calcSignedDistFct(adaptInfo_, lS_DOF_, sD_DOF, elFct);

106
    delete sD_DOF;
107
108
109
110
111
  }

  velExt->printVelDOF(adaptInfo_);
}

112
113
114
115
116
117

void HL_SignedDist::calcVelocityExt(AdaptInfo *adaptInfo_,
				    DOFVector<double> *vel_DOF_,
				    DOFVector<double> *lS_DOF_,
				    bool calcSDFct,
				    ElementFunction<double> *elFct)
118
119
120
121
122
123
124
125
126
{
  TEST_EXIT(velExtType.isSet(VEL_EXT))("illegal velocity extension type !\n");
  TEST_EXIT(velExt)("velExt not defined !\n");
  TEST_EXIT(vel_DOF_)("velocity vector vel_DOF_ not defined !\n");

  if (calcSDFct)
    TEST_EXIT(lS_DOF_)("illegal level set function lS_DOF_ !\n");

  DOFVector<double> *newVel_DOF_ = 
127
    new DOFVector<double>(vel_DOF_->getFeSpace(), "vel_DOF_");
128
129
130
131
132

  velExt->setVelocity(vel_DOF_, newVel_DOF_);

  velExt->printOrigVelDOF(adaptInfo_);

133
  if (calcSDFct && elFct == NULL) {
134
    calcSignedDistFct(adaptInfo_, lS_DOF_);
135
  } else {
136
    sD_DOF = new DOFVector<double>(vel_DOF_->getFeSpace(), "sD_DOF");
137
138
139
140
141
142

    calcSignedDistFct(adaptInfo_, lS_DOF_, sD_DOF, elFct);

    if (calcSDFct) 
      lS_DOF_->copy(*sD_DOF);

143
    delete sD_DOF;
144
145
146
147
148
149
  }

  vel_DOF_->copy(*newVel_DOF_);

  velExt->printVelDOF(adaptInfo_);

150
  delete newVel_DOF_;
151
152
}

153
154
155
156
157
158
159
160

void HL_SignedDist::calcVelocityExt(AdaptInfo *adaptInfo_,
				    std::vector<DOFVector<double> *> origVel_DOF_, 
				    std::vector<DOFVector<double> *> vel_DOF_,
				    const DOFVector<double> *lS_DOF_,
				    DOFVector<double> *sD_DOF_,
				    bool calcSDFct,
				    ElementFunction<double> *elFct)
161
162
163
164
165
166
167
168
169
170
{
  TEST_EXIT(velExtType.isSet(VEL_EXT))("illegal velocity extension type !\n");
  TEST_EXIT(velExt)("velExt not defined !\n");

  velExt->setVelocity(origVel_DOF_, vel_DOF_);

  velExt->printOrigVelDOF(adaptInfo_);

  if (calcSDFct || sD_DOF_ != NULL) {
    calcSignedDistFct(adaptInfo_, lS_DOF_, sD_DOF_, elFct);
171
  } else {
172
173
174

    TEST_EXIT(vel_DOF_[0])("illegal level set function lS_DOF_ !\n");

175
    sD_DOF = new DOFVector<double>(vel_DOF_[0]->getFeSpace(), "sD_DOF");
176
177
178

    calcSignedDistFct(adaptInfo_, lS_DOF_, sD_DOF, elFct);

179
    delete sD_DOF;
180
181
182
183
184
  }

  velExt->printVelDOF(adaptInfo_);
}

185
186
187
188
189
190

void HL_SignedDist::calcVelocityExt(AdaptInfo *adaptInfo_,
				    std::vector<DOFVector<double> *> vel_DOF_,
				    DOFVector<double> *lS_DOF_,
				    bool calcSDFct,
				    ElementFunction<double> *elFct)
191
192
193
194
195
196
197
198
199
200
{
  TEST_EXIT(velExtType.isSet(VEL_EXT))("illegal velocity extension type !\n");
  TEST_EXIT(velExt)("velExt not defined !\n");

  int nVelDOFs = vel_DOF_.size();

  if (calcSDFct)
    TEST_EXIT(lS_DOF_)("illegal level set function lS_DOF_ !\n");

  std::vector<DOFVector<double> *> newVel_DOF_(nVelDOFs);
201
202
203
  for (int nV = 0; nV < nVelDOFs; ++nV) {    
    TEST_EXIT(vel_DOF_[nV])("velocity vector vel_DOF_ not defined !\n");
    
204
    newVel_DOF_[nV] = new DOFVector<double>((vel_DOF_[nV])->getFeSpace(), "vel_DOF_");
205
206
207
208
209
210
  }

  velExt->setVelocity(vel_DOF_, newVel_DOF_);

  velExt->printOrigVelDOF(adaptInfo_);

211
  if (calcSDFct && elFct == NULL) {
212
    calcSignedDistFct(adaptInfo_, lS_DOF_);
213
  } else {
214
    sD_DOF = new DOFVector<double>(vel_DOF_[0]->getFeSpace(), "sD_DOF");
215
216
217
218
219
220

    calcSignedDistFct(adaptInfo_, lS_DOF_, sD_DOF, elFct);

    if (calcSDFct) 
      lS_DOF_->copy(*sD_DOF);

221
    delete sD_DOF;
222
223
224
225
226
227
228
  }

  for (int nV=0; nV<nVelDOFs; ++nV)
    (vel_DOF_[nV])->copy(*(newVel_DOF_[nV]));

  velExt->printVelDOF(adaptInfo_);

229
230
  for (int nV = 0; nV < nVelDOFs; ++nV)
    delete newVel_DOF_[nV];
231
232
}

233
234
235
236
237
238
239
240

void HL_SignedDist::calcVelocityExtFromVelocityField(AdaptInfo *adaptInfo_,
						     std::vector<DOFVector<double> *> &velField_, 
						     DOFVector<double> *vel_DOF_,
						     const DOFVector<double> *lS_DOF_,
						     DOFVector<double> *sD_DOF_,
						     bool calcSDFct,
						     ElementFunction<double> *elFct)
241
242
243
244
245
246
247
248
249
250
{
  TEST_EXIT(velExtType.isSet(VEL_EXT_FROM_VEL_FIELD))
    ("illegal velocity extension type !\n");
  TEST_EXIT(velExt)("velExt not defined !\n");
  TEST_EXIT(elFct == NULL)("not implemented yet for elFct != NULL !\n");

  ((VelocityExtFromVelocityField *)(velExt))->setVelocityField(velField_, lS_DOF_, vel_DOF_);

  if (calcSDFct) {
    calcSignedDistFct(adaptInfo_, lS_DOF_, sD_DOF_, elFct);
251
  } else {
252
253
254

    TEST_EXIT(vel_DOF_)("illegal velocity vector vel_DOF_ !\n");

255
    sD_DOF = new DOFVector<double>(vel_DOF_->getFeSpace(), "sD_DOF");
256
257
258

    calcSignedDistFct(adaptInfo_, lS_DOF_, sD_DOF, elFct);

259
    delete sD_DOF;
260
261
262
263
264
  }

  velExt->printVelDOF(adaptInfo_);
}

265
266
267
268
269
270
271

void HL_SignedDist::calcVelocityExtFromVelocityField(AdaptInfo *adaptInfo_,
						     std::vector<DOFVector<double> *> &velField_, 
						     DOFVector<double> *vel_DOF_,
						     DOFVector<double> *lS_DOF_,
						     bool calcSDFct,
						     ElementFunction<double> *elFct)
272
273
274
275
276
277
278
279
280
281
282
283
284
{
  TEST_EXIT(velExtType.isSet(VEL_EXT_FROM_VEL_FIELD))
    ("illegal velocity extension type !\n");
  TEST_EXIT(velExt)("velExt not defined !\n");
  TEST_EXIT(elFct == NULL)("not implemented yet for elFct != NULL !\n");

  if (calcSDFct)
    TEST_EXIT(lS_DOF_)("illegal level set function lS_DOF_ !\n");

  ((VelocityExtFromVelocityField *)(velExt))->setVelocityField(velField_, lS_DOF_, vel_DOF_);

  if (calcSDFct && elFct == NULL) {
    calcSignedDistFct(adaptInfo_, lS_DOF_);
285
  } else {
286
287
288

    TEST_EXIT(vel_DOF_)("illegal velocity vector vel_DOF_ !\n");

289
    sD_DOF = new DOFVector<double>(vel_DOF_->getFeSpace(), "sD_DOF");
290
291
292
293
294
295

    calcSignedDistFct(adaptInfo_, lS_DOF_, sD_DOF, elFct);

    if (calcSDFct)
      lS_DOF_->copy(*sD_DOF);

296
    delete sD_DOF;
297
298
299
300
301
  }

  velExt->printVelDOF(adaptInfo_);
}

302
303

void HL_SignedDist::initialize(ElementFunction<double> *elFct)
304
305
{
  if (elUpdate)
306
    delete elUpdate;
307
  if (bndElDist)
308
    delete bndElDist;
309
310
  
  if (elLS)
311
    delete elLS;
312
  if (phi)
313
    delete phi;
314
315
316
  
  // ===== Create ElementLevelSet. =====
  if (elFct == NULL) {
317
318
319
320
    phi = new ElementFunctionDOFVec<double>(lS_DOF);
    elLS = new ElementLevelSet("ElementLevelSet", phi, feSpace->getMesh());
  } else {
    elLS = new ElementLevelSet("ElementLevelSet", elFct, feSpace->getMesh());
321
322
323
324
325
326
327
328
  }

  
  // ===== Define boundary initialization strategy. =====
  int boundInitFlag;
  GET_PARAMETER(0, name + "->boundary initialization", "%d", 
		&boundInitFlag);
  switch(boundInitFlag) {
329
  case 0: bndElDist = new BoundaryElementLevelSetDist(elLS, dim);
330
    break;
331
  case 1: bndElDist = new BoundaryElementTopDist(elLS, dim, velExt);
332
    break;
333
  case 2: bndElDist = new BoundaryElementEdgeDist(elLS, dim);
334
    break;
335
  case 3: bndElDist = new BoundaryElementNormalDist(elLS, dim);
336
337
338
339
340
341
342
343
344
345
    break;
  default: ERROR_EXIT("illegal boundary initialization !\n");
    break;
  }
  
  // ===== Create elUpdate. =====
  TEST_EXIT(velExt == NULL || boundInitFlag == 1)
    ("velocity extension only works with topological boundary element initialization method !\n");

  switch(dim) {
346
  case 2: elUpdate = new ElementUpdate_2d(velExt);
347
    break;
348
  case 3: elUpdate = new ElementUpdate_3d(velExt);
349
350
351
352
353
354
    break;
  default: ERROR_EXIT("illegal dimension !\n");
    break;
  }
}

355
356

void HL_SignedDist::setSign() 
357
{
358
359
360
361
362
363
364
  DOFVector<double>::Iterator it_sD(const_cast<DOFVector<double> *>(sD_DOF), USED_DOFS);
  DOFVector<double>::Iterator it_lS(const_cast<DOFVector<double> *>(lS_DOF), USED_DOFS);

  // Vertex lies inside the zero level set ?  
  for (it_sD.reset(), it_lS.reset(); !it_sD.end(); ++it_sD, ++it_lS)    
    if ((*it_lS) < 0)
      (*it_sD) *= -1;      
365
}