HL_SignedDist.cc 9.74 KB
Newer Older
1
2
3
4
5
#include "HL_SignedDist.h"

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

6
7
8
9
void HL_SignedDist::calcSignedDistFct(AdaptInfo *adaptInfo_,
				      const DOFVector<double> *lS_DOF_, 
				      DOFVector<double> *sD_DOF_,
				      ElementFunction<double> *elFct) 
10
11
12
13
14
{
  adaptInfo = adaptInfo_;

  if (elFct == NULL) {
    TEST_EXIT(lS_DOF_)("illegal level set function !\n");
15
    TEST_EXIT(lS_DOF_->getFeSpace() == sD_DOF_->getFeSpace())("DOF vectors do not have the same FE space!\n");
16
17
  }
  TEST_EXIT(sD_DOF_)("illegal DOF vector for signed distance function !\n");
18
  TEST_EXIT(sD_DOF_->getFeSpace()->getBasisFcts()->getDegree() == 1)("does it work for non-linear finite elements ?\n");
19
20
21
  
  lS_DOF = lS_DOF_;
  sD_DOF = sD_DOF_;
22
  feSpace = sD_DOF->getFeSpace();
23
24
25
26
27
28
29
30
31
32

  // ===== Initialization. =====
  initialize(elFct);
  
  // ===== Print level set function. =====
  if (lS_DOF) 
    printLevelSetFct();
  
  // ===== Create DOF vector to mark boundary vertices. =====
  if (bound_DOF) 
33
34
    delete bound_DOF;
  bound_DOF = new DOFVector<double>(feSpace, "bound_DOF");
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
  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();
}

54
55

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

59
  sD_DOF = new DOFVector<double>(lS_DOF_->getFeSpace(), "sD_DOF");
60
61
62
63
64

  calcSignedDistFct(adaptInfo_, lS_DOF_, sD_DOF, NULL);

  *lS_DOF_ = *sD_DOF;

65
  delete sD_DOF;
66
67
}

68
69
70
71
72
73
74
75

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)
76
77
78
79
80
81
82
83
84
85
{
  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);
86
  } else {
87
88
89

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

90
    sD_DOF = new DOFVector<double>(vel_DOF_->getFeSpace(), "sD_DOF");
91
92
93

    calcSignedDistFct(adaptInfo_, lS_DOF_, sD_DOF, elFct);

94
    delete sD_DOF;
95
96
97
98
99
  }

  velExt->printVelDOF(adaptInfo_);
}

100
101
102
103
104
105

void HL_SignedDist::calcVelocityExt(AdaptInfo *adaptInfo_,
				    DOFVector<double> *vel_DOF_,
				    DOFVector<double> *lS_DOF_,
				    bool calcSDFct,
				    ElementFunction<double> *elFct)
106
107
108
109
110
111
112
113
114
{
  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_ = 
115
    new DOFVector<double>(vel_DOF_->getFeSpace(), "vel_DOF_");
116
117
118
119
120

  velExt->setVelocity(vel_DOF_, newVel_DOF_);

  velExt->printOrigVelDOF(adaptInfo_);

121
  if (calcSDFct && elFct == NULL) {
122
    calcSignedDistFct(adaptInfo_, lS_DOF_);
123
  } else {
124
    sD_DOF = new DOFVector<double>(vel_DOF_->getFeSpace(), "sD_DOF");
125
126
127
128
129
130

    calcSignedDistFct(adaptInfo_, lS_DOF_, sD_DOF, elFct);

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

131
    delete sD_DOF;
132
133
134
135
136
137
  }

  vel_DOF_->copy(*newVel_DOF_);

  velExt->printVelDOF(adaptInfo_);

138
  delete newVel_DOF_;
139
140
}

141
142
143
144
145
146
147
148

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)
149
150
151
152
153
154
155
156
157
158
{
  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);
159
  } else {
160
161
162

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

163
    sD_DOF = new DOFVector<double>(vel_DOF_[0]->getFeSpace(), "sD_DOF");
164
165
166

    calcSignedDistFct(adaptInfo_, lS_DOF_, sD_DOF, elFct);

167
    delete sD_DOF;
168
169
170
171
172
  }

  velExt->printVelDOF(adaptInfo_);
}

173
174
175
176
177
178

void HL_SignedDist::calcVelocityExt(AdaptInfo *adaptInfo_,
				    std::vector<DOFVector<double> *> vel_DOF_,
				    DOFVector<double> *lS_DOF_,
				    bool calcSDFct,
				    ElementFunction<double> *elFct)
179
180
181
182
183
184
185
186
187
188
{
  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);
189
190
191
  for (int nV = 0; nV < nVelDOFs; ++nV) {    
    TEST_EXIT(vel_DOF_[nV])("velocity vector vel_DOF_ not defined !\n");
    
192
    newVel_DOF_[nV] = new DOFVector<double>((vel_DOF_[nV])->getFeSpace(), "vel_DOF_");
193
194
195
196
197
198
  }

  velExt->setVelocity(vel_DOF_, newVel_DOF_);

  velExt->printOrigVelDOF(adaptInfo_);

199
  if (calcSDFct && elFct == NULL) {
200
    calcSignedDistFct(adaptInfo_, lS_DOF_);
201
  } else {
202
    sD_DOF = new DOFVector<double>(vel_DOF_[0]->getFeSpace(), "sD_DOF");
203
204
205
206
207
208

    calcSignedDistFct(adaptInfo_, lS_DOF_, sD_DOF, elFct);

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

209
    delete sD_DOF;
210
211
212
213
214
215
216
  }

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

  velExt->printVelDOF(adaptInfo_);

217
218
  for (int nV = 0; nV < nVelDOFs; ++nV)
    delete newVel_DOF_[nV];
219
220
}

221
222
223
224
225
226
227
228

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)
229
230
231
232
233
234
235
236
237
238
{
  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);
239
  } else {
240
241
242

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

243
    sD_DOF = new DOFVector<double>(vel_DOF_->getFeSpace(), "sD_DOF");
244
245
246

    calcSignedDistFct(adaptInfo_, lS_DOF_, sD_DOF, elFct);

247
    delete sD_DOF;
248
249
250
251
252
  }

  velExt->printVelDOF(adaptInfo_);
}

253
254
255
256
257
258
259

void HL_SignedDist::calcVelocityExtFromVelocityField(AdaptInfo *adaptInfo_,
						     std::vector<DOFVector<double> *> &velField_, 
						     DOFVector<double> *vel_DOF_,
						     DOFVector<double> *lS_DOF_,
						     bool calcSDFct,
						     ElementFunction<double> *elFct)
260
261
262
263
264
265
266
267
268
269
270
271
272
{
  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_);
273
  } else {
274
275
276

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

277
    sD_DOF = new DOFVector<double>(vel_DOF_->getFeSpace(), "sD_DOF");
278
279
280
281
282
283

    calcSignedDistFct(adaptInfo_, lS_DOF_, sD_DOF, elFct);

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

284
    delete sD_DOF;
285
286
287
288
289
  }

  velExt->printVelDOF(adaptInfo_);
}

290
291

void HL_SignedDist::initialize(ElementFunction<double> *elFct)
292
293
{
  if (elUpdate)
294
    delete elUpdate;
295
  if (bndElDist)
296
    delete bndElDist;
297
298
  
  if (elLS)
299
    delete elLS;
300
  if (phi)
301
    delete phi;
302
303
304
  
  // ===== Create ElementLevelSet. =====
  if (elFct == NULL) {
305
306
307
308
    phi = new ElementFunctionDOFVec<double>(lS_DOF);
    elLS = new ElementLevelSet("ElementLevelSet", phi, feSpace->getMesh());
  } else {
    elLS = new ElementLevelSet("ElementLevelSet", elFct, feSpace->getMesh());
309
310
311
312
313
314
315
316
  }

  
  // ===== Define boundary initialization strategy. =====
  int boundInitFlag;
  GET_PARAMETER(0, name + "->boundary initialization", "%d", 
		&boundInitFlag);
  switch(boundInitFlag) {
317
  case 0: bndElDist = new BoundaryElementLevelSetDist(elLS, dim);
318
    break;
319
  case 1: bndElDist = new BoundaryElementTopDist(elLS, dim, velExt);
320
    break;
321
  case 2: bndElDist = new BoundaryElementEdgeDist(elLS, dim);
322
    break;
323
  case 3: bndElDist = new BoundaryElementNormalDist(elLS, dim);
324
325
326
327
328
329
330
331
332
333
    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) {
334
  case 2: elUpdate = new ElementUpdate_2d(velExt);
335
    break;
336
  case 3: elUpdate = new ElementUpdate_3d(velExt);
337
338
339
340
341
342
    break;
  default: ERROR_EXIT("illegal dimension !\n");
    break;
  }
}

343
344

void HL_SignedDist::setSign() 
345
{
346
347
348
349
350
351
352
  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;      
353
}