HL_SignedDist.cc 10 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
{
  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");

123
  if (calcSDFct) {
124
    TEST_EXIT(lS_DOF_)("illegal level set function lS_DOF_ !\n");
125
  }
126
127

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

  velExt->setVelocity(vel_DOF_, newVel_DOF_);

  velExt->printOrigVelDOF(adaptInfo_);

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

    calcSignedDistFct(adaptInfo_, lS_DOF_, sD_DOF, elFct);

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

144
    delete sD_DOF;
145
146
147
148
149
150
  }

  vel_DOF_->copy(*newVel_DOF_);

  velExt->printVelDOF(adaptInfo_);

151
  delete newVel_DOF_;
152
153
}

154
155
156
157
158
159
160
161

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)
162
163
164
165
166
167
168
169
170
171
{
  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);
172
  } else {
173
174
175

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

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

    calcSignedDistFct(adaptInfo_, lS_DOF_, sD_DOF, elFct);

180
    delete sD_DOF;
181
182
183
184
185
  }

  velExt->printVelDOF(adaptInfo_);
}

186
187
188
189
190
191

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

  int nVelDOFs = vel_DOF_.size();

198
  if (calcSDFct) {
199
    TEST_EXIT(lS_DOF_)("illegal level set function lS_DOF_ !\n");
200
  }
201
202

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

  velExt->setVelocity(vel_DOF_, newVel_DOF_);

  velExt->printOrigVelDOF(adaptInfo_);

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

    calcSignedDistFct(adaptInfo_, lS_DOF_, sD_DOF, elFct);

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

223
    delete sD_DOF;
224
225
226
227
228
229
230
  }

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

  velExt->printVelDOF(adaptInfo_);

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

235
236
237
238
239
240
241
242

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)
243
244
245
246
247
248
249
250
251
252
{
  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);
253
  } else {
254
255
256

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

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

    calcSignedDistFct(adaptInfo_, lS_DOF_, sD_DOF, elFct);

261
    delete sD_DOF;
262
263
264
265
266
  }

  velExt->printVelDOF(adaptInfo_);
}

267
268
269
270
271
272
273

void HL_SignedDist::calcVelocityExtFromVelocityField(AdaptInfo *adaptInfo_,
						     std::vector<DOFVector<double> *> &velField_, 
						     DOFVector<double> *vel_DOF_,
						     DOFVector<double> *lS_DOF_,
						     bool calcSDFct,
						     ElementFunction<double> *elFct)
274
275
276
277
278
279
{
  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");

280
  if (calcSDFct) {
281
    TEST_EXIT(lS_DOF_)("illegal level set function lS_DOF_ !\n");
282
  }
283
284
285
286
287

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

  if (calcSDFct && elFct == NULL) {
    calcSignedDistFct(adaptInfo_, lS_DOF_);
288
  } else {
289
290
291

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

292
    sD_DOF = new DOFVector<double>(vel_DOF_->getFeSpace(), "sD_DOF");
293
294
295
296
297
298

    calcSignedDistFct(adaptInfo_, lS_DOF_, sD_DOF, elFct);

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

299
    delete sD_DOF;
300
301
302
303
304
  }

  velExt->printVelDOF(adaptInfo_);
}

305
306

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

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

352
353
354
  switch (dim) {
  case 2: 
    elUpdate = new ElementUpdate_2d(velExt);
355
    break;
356
357
  case 3: 
    elUpdate = new ElementUpdate_3d(velExt);
358
    break;
359
360
  default: 
    ERROR_EXIT("illegal dimension !\n");
361
362
363
364
    break;
  }
}

365
366

void HL_SignedDist::setSign() 
367
{
368
369
370
371
372
373
374
  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;      
375
}