HL_SignedDistTraverse.cc 9.31 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
/******************************************************************************
 *
 * AMDiS - Adaptive multidimensional simulations
 *
 * Copyright (C) 2013 Dresden University of Technology. All Rights Reserved.
 * Web: https://fusionforge.zih.tu-dresden.de/projects/amdis
 *
 * Authors: 
 * Simon Vey, Thomas Witkowski, Andreas Naumann, Simon Praetorius, et al.
 *
 * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
 * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
 *
 *
 * This file is part of AMDiS
 *
 * See also license.opensource.txt in the distribution.
 * 
 ******************************************************************************/
20
21


22
23
24
#include "HL_SignedDistTraverse.h"
#include "VelocityExtFromVelocityField.h"

25
26
namespace reinit {

27
void HL_SignedDistTraverse::initializeBoundary() 
28
29
{
  FUNCNAME("HL_SignedDistTraverse::initializeBoundary()");
Thomas Witkowski's avatar
Thomas Witkowski committed
30
 
31
32
33
34
35
36
37
38
39
  // ===== All non-boundary vertices are initialized with "infinity". =====
  sD_DOF->set(inftyValue);
  
  // ===== Traverse mesh and initialize boundary elements. =====
  TraverseStack stack;
  FixVec<double, VERTEX> distVec(dim, NO_INIT);
  int elStatus;
  
  const int nBasFcts = feSpace->getBasisFcts()->getNumber();
40
  if (static_cast<int>(locInd.size()) < nBasFcts)
41
    locInd.resize(nBasFcts);
42
43
44
45
46
47
48
49
50
  
  ElInfo *elInfo;
  if (velExt && velExtType.isSet(VEL_EXT_FROM_VEL_FIELD)) {
    elInfo = stack.traverseFirst(feSpace->getMesh(),
				 -1, 
				 Mesh::CALL_LEAF_EL | 
				 Mesh::FILL_BOUND |
				 Mesh::FILL_COORDS |
				 Mesh::FILL_GRD_LAMBDA);
51
  } else {
52
53
54
55
56
57
    elInfo = stack.traverseFirst(feSpace->getMesh(),
				 -1, 
				 Mesh::CALL_LEAF_EL | 
				 Mesh::FILL_BOUND |
				 Mesh::FILL_COORDS);
  }
58

59
60
61
  while(elInfo) {

    // Set elInfo in case velocity extension from velocity field is used.
62
63
    if (velExt && velExtType.isSet(VEL_EXT_FROM_VEL_FIELD))
      ((VelocityExtFromVelocityField *)velExt)->setElInfo(elInfo);    
64
65
    
    // Get local indices of vertices.
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
66
67
68
    feSpace->getBasisFcts()->getLocalIndices(const_cast<Element*>(elInfo->getElement()),
					     const_cast<DOFAdmin*>(feSpace->getAdmin()),
					     locInd); 
69
70
71
72
73
74
75
76
    
    // Get element status.
    elStatus = elLS->createElementLevelSet(elInfo);
    
    // Is element cut by the interface ?
    if (elStatus == ElementLevelSet::LEVEL_SET_BOUNDARY) {
      
      // Reset element distance vector.
77
78
      for (int i=0; i <= dim; i++)
	distVec[i] = inftyValue;      
79
80
      
      // Mark all vertices as boundary vertices.
81
82
      for (int i = 0; i <= dim; i++)
	(*bound_DOF)[locInd[i]] = 1.0;      
83
84
85
86
87
      
      // Calculate distance for all vertices.
      if (bndElDist->calcDistOnBoundaryElement(elInfo, distVec) !=
	  ElementLevelSet::LEVEL_SET_BOUNDARY) {
	ERROR_EXIT("error in distance calculation !\n");
88
      } else {
89
90
	
	// If distance is smaller, correct to new distance.
91
	for (int i = 0; i <= dim; i++) {
92
93
	  
	  // --> for test purposes:
Thomas Witkowski's avatar
blub    
Thomas Witkowski committed
94
95
96
97
98
99
	  if (distVec[i] > 1000) {
	    MSG("Element %d  Knoten %d  keine Randwertinitialisierung!\n",
		elInfo->getElement()->getIndex(), i);
	  }
		  
	  // --> end: for test purposes	  
100
101
102
	  if ((*sD_DOF)[locInd[i]] > distVec[i]) {
	    (*sD_DOF)[locInd[i]] = distVec[i];
	    //If Distance is corrected, calculate new velocity.
103
	    if (velExt != NULL)	      
104
	      velExt->calcVelocityBoundary(&(locInd[0]), i);	      
105
106
107
108
109
110
	  }
	}
      }
    }  // end of: elStatus == ElementLevelSet::LEVEL_SET_BOUNDARY
    
    elInfo = stack.traverseNext(elInfo);
111
  }  // end of: mesh traverse 
112
113
114
115
116
117
  
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
  // In parallel AMDiS synchronize the bound_DOF DOFVector with the max-assigner and the sD_DOF DOFVector with the min-assigner
  AMDiS::Parallel::MeshDistributor::globalMeshDistributor->synchVector(*bound_DOF, max_assigner());
  AMDiS::Parallel::MeshDistributor::globalMeshDistributor->synchVector(*sD_DOF, min_to_zero_assigner());
#endif
118
119
}

120

121
void HL_SignedDistTraverse::HL_updateIteration()
122
{
Thomas Witkowski's avatar
blub    
Thomas Witkowski committed
123
124
  FUNCNAME("HL_SignedDistTraverse::HL_updateIteration()");

125
126
  // ===== Create DOF vector for the last iteration step. =====
  if (sDOld_DOF)
127
128
129
    delete sDOld_DOF;

  sDOld_DOF = new DOFVector<double>(feSpace, "sDOld_DOF");
130
131
  sDOld_DOF->copy(const_cast<DOFVector<double> &>(*sD_DOF));
  
132
133
134
135
136
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
  // Update sDOld_DOF on interior domain boundaries
  AMDiS::Parallel::MeshDistributor::globalMeshDistributor->synchVector(*sDOld_DOF, min_to_zero_assigner());
#endif  
  
137
  // ===== Gauss-Seidel or Jacobi iteration ? =====
138
139
140
141
  if (GaussSeidelFlag)
    update_DOF = sD_DOF;  
  else
    update_DOF = sDOld_DOF;  
142
143
144
145
146
147
  
  // ===== Iteration loop: proceed until tolerance is reached. =====
  TraverseStack stack;
  ElInfo *elInfo;
  tol_reached = false;
  int itCntr = 0;
148
  while (!tol_reached && itCntr != maxIt) {   
149
150
151
    ++itCntr;
    tol_reached = true;
    
152
153
154
155
156
    // ===== Traverse mesh: perform Hopf-Lax element update on each element. =====
    elInfo = 
      stack.traverseFirst(feSpace->getMesh(), -1, 
			  Mesh::CALL_LEAF_EL | Mesh::FILL_BOUND | Mesh::FILL_COORDS);
    while (elInfo) {
157
158
159
160
      HL_elementUpdate(elInfo);
      elInfo = stack.traverseNext(elInfo);
    }
    
161
162
163
164
165
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    // Update sD_DOF on interior domain boundaries
    AMDiS::Parallel::MeshDistributor::globalMeshDistributor->synchVector(*sD_DOF, min_to_zero_assigner());
#endif    
    
166
167
    // ===== Is tolerance reached ? =====
    tol_reached = checkTol();
168
169
170
171
172
173
174
175
176
177
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    // Convert bool to int that the datatype MPI_INT can be used in the MPI_Allreduce() command below.
    int int_tol_reached = (tol_reached ? 1 : 0);
    int int_tol_reached_overall;
    // Here MPI_MIN operation in the MPI_Allreduce() command is used that each processor runs the cycle again until each
    // processor has reached the tolerance. Otherwise the synchVector() method from above would fail.
    MPI_Allreduce(&int_tol_reached, &int_tol_reached_overall, 1, MPI_INT, MPI_MIN, MPI_COMM_WORLD);
    int_tol_reached = int_tol_reached_overall;
    tol_reached = (int_tol_reached == 0 ? false : true);
#endif
178
179
180
181
    
    sDOld_DOF->copy(const_cast<DOFVector<double> &>(*sD_DOF));
  }
  
182

Thomas Witkowski's avatar
blub    
Thomas Witkowski committed
183
  MSG("Calculation of signed distance function via mesh traverse iteration:\n");
184
  if (GaussSeidelFlag)
Thomas Witkowski's avatar
blub    
Thomas Witkowski committed
185
    MSG("\tGauss-Seidel iteration\n");
186
  else
Thomas Witkowski's avatar
blub    
Thomas Witkowski committed
187
    MSG("\tJacobi iteration\n");
188
  
Thomas Witkowski's avatar
blub    
Thomas Witkowski committed
189
  MSG("\tnumber of iterations needed: %d\n", itCntr);
190
191
}

192

193
void HL_SignedDistTraverse::HL_elementUpdate(ElInfo *elInfo) 
194
195
196
{
  // ===== Get global indices of vertices of element. =====
  const int nBasFcts = feSpace->getBasisFcts()->getNumber();
197
  if (static_cast<int>(locInd.size()) < nBasFcts)
198
    locInd.resize(nBasFcts);
199
  
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
200
201
202
  feSpace->getBasisFcts()->getLocalIndices(const_cast<Element *>(elInfo->getElement()),
					   const_cast<DOFAdmin *>(feSpace->getAdmin()),
					   locInd);
203
204
  
  // ===== Hopf-Lax element update for each vertex of element. =====
205
  for (int i = 0; i <= dim; i++) {
206
207
208
209
    
    // ===== Calculate update for non-boundary vertex. =====
    if ((*bound_DOF)[locInd[i]] < 1.e-15) {
      //save permutation of vertexes for calculation of the velocity
210
      if (velExt != NULL)
211
212
213
	velExt->setPermutation(i, 1);
	
      double update = calcElementUpdate(elInfo, i, &(locInd[0]));
214
215
216
217
218
219
220
221
      // ---> for test purposes: count number of calculated updates
      ++calcUpdate_Cntr;
      // ---> end: for test purposes
      
      // Calculates minimum of all element updates on elements
      // containing vertex i.
      if (update < (*sD_DOF)[locInd[i]]) {
	(*sD_DOF)[locInd[i]] = update;
222
223

	// If Distance is corrected, calculate new velocity.
224
	if(velExt != NULL)
225
226
	  velExt->calcVelocity(&(locInd[0]), i);
	
227
228
229
230
231
	// ---> for test purposes: count number of calculated updates
	++setUpdate_Cntr;
	// ---> end: for test purposes
      }
    }
232
  } 
233
234
}

235

236
237
238
double HL_SignedDistTraverse::calcElementUpdate(ElInfo *elInfo, 
						int nXh,
						const DegreeOfFreedom *locInd) 
239
240
{
  // ===== Get local indices of element vertices (except xh). =====
241
242
243
244
245
  int nWh = 0;  
  int nYh = (nXh + 1) % (dim+1);
  int nZh = (nXh + 2) % (dim+1);
  if (dim == 3) 
    nWh = (nXh + 3) % (dim+1);
246
247
248
249
250
251
  
  // ===== Get world coordinates of vertices of element and their values
  //       of uh. 
  //       The coordinates of the vertex the update is calculated for
  //       are stored at the end of the vector elVert. =====
  switch (dim) {
252
253
  case 2: 
    elVert[0] = &(elInfo->getCoord(nYh));
254
255
256
257
258
259
260
261
    elVert[1] = &(elInfo->getCoord(nZh));
    elVert[2] = &(elInfo->getCoord(nXh));
    
    uhVal[0] = (*update_DOF)[locInd[nYh]];
    uhVal[1] = (*update_DOF)[locInd[nZh]];
    uhVal[2] = (*update_DOF)[locInd[nXh]];
    
    break;
262
263
  case 3: 
    elVert[0] = &(elInfo->getCoord(nYh));
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
    elVert[1] = &(elInfo->getCoord(nZh));
    elVert[2] = &(elInfo->getCoord(nWh));
    elVert[3] = &(elInfo->getCoord(nXh));

    uhVal[0] = (*update_DOF)[locInd[nYh]];
    uhVal[1] = (*update_DOF)[locInd[nZh]];
    uhVal[2] = (*update_DOF)[locInd[nWh]];
    uhVal[3] = (*update_DOF)[locInd[nXh]];
    
    break;
  default: ERROR_EXIT("illegal dimension !\n");
    break;
  }
  
  // ===== Calculate Hopf-Lax element update for vertex. =====
279
  return elUpdate->calcElementUpdate(elVert, uhVal);
280
281
}

282

283
bool HL_SignedDistTraverse::checkTol()
284
285
286
287
{
  DOFVector<double>::Iterator it_sD(sD_DOF, USED_DOFS);
  DOFVector<double>::Iterator it_sDOld(sDOld_DOF, USED_DOFS);
  
288
289
  for (it_sD.reset(), it_sDOld.reset(); !it_sD.end(); ++it_sD, ++it_sDOld)
    if ((*it_sDOld) - (*it_sD) > tol || (*it_sDOld) - (*it_sD) < 0)
290
291
292
293
      return false;
  
  return true;
}
294
295

}