DataCollector.cc 13.5 KB
Newer Older
1
2
3
4
5
6
7
8
#include "DataCollector.h"
#include "Traverse.h"
#include "DOFVector.h"
#include "SurfaceRegion_ED.h"
#include "ElementRegion_ED.h"
#include "Projection.h"

namespace AMDiS {
9

10
11
12
13
14
15
16
17
18
19
  DataCollector::DataCollector(const FiniteElemSpace *fe,
			       DOFVector<double> *val,
			       int l,
			       Flag flag,
			       bool (*writeElemFct)(ElInfo*))
    : values(val),
      level(l),
      traverseFlag(flag),
      feSpace(fe),
      writeElem(writeElemFct)
20
  {
21
    FUNCNAME("DataCollector::DataCollector()");
22
23
24

    TEST_EXIT(feSpace)("no feSpace\n");
    TEST_EXIT(values)("no value Vector\n");
25
   
26
    // === get mesh  ===    
27
    mesh = feSpace->getMesh();
28
29

    // === get admin ===    
30
    localAdmin = const_cast<DOFAdmin*>(feSpace->getAdmin());
31
    // === create vertex info vector ===
32
    vertexInfos = new DOFVector< std::list<VertexInfo> >(feSpace, "vertex infos");
33

34
35
36
    interpPointInd = new DOFVector<int>(feSpace, "interpolation point indices");
    interpPointCoords = new DOFVector< std::list<WorldVector<double> > >(feSpace, "interpolation point coordinates");
    dofCoord = new DOFVector< std::list<WorldVector<double> > >(feSpace, "dof coords");
37

38
    dim = mesh->getDim();
39
    
40
41
42
43
44
    nPreDofs = localAdmin->getNumberOfPreDOFs(VERTEX);
    nVertices = 0;
    nElements = 0;
    nInterpPoints = 0;
    nConnection = 0;
45

46
47
48
    elementDataCollected = false;
    valueDataCollected = false;
    periodicDataCollected = false;
Thomas Witkowski's avatar
Thomas Witkowski committed
49

Thomas Witkowski's avatar
Thomas Witkowski committed
50
    vertexCoords = new WorldVector<double>;
Thomas Witkowski's avatar
Thomas Witkowski committed
51
  } 
52
53
54

  DataCollector::~DataCollector()
  {
55
56
57
58
    delete vertexInfos;
    delete interpPointInd;
    delete interpPointCoords;
    delete dofCoord;
Thomas Witkowski's avatar
Thomas Witkowski committed
59
    delete vertexCoords;
60
61
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
62
63
  void DataCollector::fillAllData()
  {
64
    if (!elementDataCollected)
Thomas Witkowski's avatar
Thomas Witkowski committed
65
      startCollectingElementData();
66

67
    if (!periodicDataCollected)
Thomas Witkowski's avatar
Thomas Witkowski committed
68
      startCollectingPeriodicData();
69

70
    if (!valueDataCollected)
Thomas Witkowski's avatar
Thomas Witkowski committed
71
72
73
      startCollectingValueData();
  }

74
  void DataCollector::startCollectingElementData()
75
  {
76
    FUNCNAME("DataCollector::startCollectingElementData()");
77

78
    Flag flag = traverseFlag;
79
80
81
82
83
84
    flag |= 
      Mesh::FILL_NEIGH      | 
      Mesh::FILL_COORDS     |
      Mesh::FILL_OPP_COORDS |
      Mesh::FILL_BOUND;

85
    elements.resize(0, dim);
86
87
88
89

    TraverseStack stack;

    // Traverse elements to create continuous element indices
90
    ElInfo *elInfo = stack.traverseFirst(mesh, level, flag);
91
    while (elInfo) {
92
93
      if (!writeElem || writeElem(elInfo))
	outputIndices[elInfo->getElement()->getIndex()] = nElements++;
94
95
96
97
      elInfo = stack.traverseNext(elInfo);
    }

    // Traverse elements to create element information
98
    elInfo = stack.traverseFirst(mesh, level, flag);
99

100
    while (elInfo) {
101
      if (!writeElem || writeElem(elInfo))
102
	addElementData(elInfo);
103

104
105
106
      elInfo = stack.traverseNext(elInfo);
    }

107
    elementDataCollected = true;
108
109
  }

110
  void DataCollector::startCollectingValueData()
111
  {
112
    FUNCNAME("DataCollector::startCollectingValueData()");
113

114
    DOFVector<int>::Iterator intPointIt(interpPointInd, USED_DOFS);
115
    for (intPointIt.reset(); !intPointIt.end(); ++intPointIt)
116
117
      (*intPointIt) = -1;

118
    interpPoints.clear();
119

120
121
122
    basisFcts = const_cast<BasisFunction*>(feSpace->getBasisFcts());
    nBasisFcts = basisFcts->getNumber();
    localDOFs = new DegreeOfFreedom[nBasisFcts];
123
  
124
125
126
127
    TraverseStack stack;

    // Traverse elements to add value information and to mark
    // interpolation points.
128
129
    ElInfo *elInfo = stack.traverseFirst(mesh, level, 
					 traverseFlag | Mesh::FILL_COORDS);
130
    while (elInfo) {
131
      if (!writeElem || writeElem(elInfo))
Thomas Witkowski's avatar
Thomas Witkowski committed
132
	addValueData(elInfo);
133
134
135
      elInfo = stack.traverseNext(elInfo);
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
136
137
    // If there are interpolation points, add them to the corresponding 
    // data array.
138
    if (nInterpPoints > 0) {
Thomas Witkowski's avatar
Thomas Witkowski committed
139
140
141
      // Remove all interpolation marks and, instead, set to each
      // interpolation point its continous index starting from 0.
      int i = 0;
142
143
      for (intPointIt.reset(); !intPointIt.end(); ++intPointIt)
	if (*intPointIt == -3)
Thomas Witkowski's avatar
Thomas Witkowski committed
144
145
146
	  *intPointIt = i++;
      
      // Traverse elements to create interpolation values.
147
      elInfo = stack.traverseFirst(mesh, level, traverseFlag | Mesh::FILL_COORDS);
Thomas Witkowski's avatar
Thomas Witkowski committed
148
      while (elInfo) {
149
	if (!writeElem || writeElem(elInfo))
Thomas Witkowski's avatar
Thomas Witkowski committed
150
151
152
	  addInterpData(elInfo);
	elInfo = stack.traverseNext(elInfo);
      }
153
    }
154
   
155
156
    delete [] localDOFs;
    valueDataCollected = true;
157
158
  }

159
  void DataCollector::startCollectingPeriodicData()
160
  {    
161
    FUNCNAME("DataCollector::startCollectingPeriodicData()");
162

163
    periodicConnections.clear();
164
165
    
    TraverseStack stack;
166
    ElInfo *elInfo = stack.traverseFirst(mesh, level, traverseFlag);
167
    while (elInfo) {
168
      if (!writeElem || writeElem(elInfo)) {
169
170
171
172
173
	LeafDataPeriodic *ldp = dynamic_cast<LeafDataPeriodic*>
	  (elInfo->getElement()->
	   getElementData()->
	   getElementData(PERIODIC));
	
174
	if (ldp) {
175
	  nConnection += dynamic_cast<LeafDataPeriodic*>(ldp)->getInfoList().size();
176
177
	}

178
	periodicConnections.push_back(DimVec<bool>(dim, DEFAULT_VALUE, false));
179
180
181
182
      }
      elInfo = stack.traverseNext(elInfo);
    }   

183
    nConnection /= 2;
184

185
    periodicInfos.clear();
186

187
    Flag flag = traverseFlag;
188
189
190
191
192
193
    flag |= 
      Mesh::FILL_COORDS    |
      Mesh::FILL_OPP_COORDS|
      Mesh::FILL_NEIGH     | 
      Mesh::FILL_BOUND;

194
    elInfo = stack.traverseFirst(mesh, level, flag);
195
    while (elInfo) {
196
      if (!writeElem || writeElem(elInfo))
197
198
199
200
	addPeriodicData(elInfo);
      elInfo = stack.traverseNext(elInfo);
    }   

201
    periodicDataCollected = true;
202
203
  }

204
  void DataCollector::addElementData(ElInfo* elInfo)
205
  {
206
    FUNCNAME("DataCollector::addElementData()");
207
208
209
210

    const DegreeOfFreedom **dof = elInfo->getElement()->getDOF();
    DegreeOfFreedom vertexDOF;
    WorldVector<double> vertexCoords;
211

212
    // create ElementInfo
213
    ElementInfo elementInfo(dim);
214
215
    
    // read element region
216
    ElementData *ed = elInfo->getElement()->getElementData(ELEMENT_REGION);
217
    
Thomas Witkowski's avatar
Thomas Witkowski committed
218
    if (ed)
219
      elementInfo.elementRegion = dynamic_cast<ElementRegion_ED*>(ed)->getRegion();
Thomas Witkowski's avatar
Thomas Witkowski committed
220
221
    else 
      elementInfo.elementRegion = -1;    
222
   
223
224
225
    // read surface regions to element info
    ed = elInfo->getElement()->getElementData(SURFACE_REGION);
    elementInfo.surfaceRegions.set(-1);
226
    while (ed) {
227
228
229
230
231
232
      SurfaceRegion_ED *sr = dynamic_cast<SurfaceRegion_ED*>(ed);
      elementInfo.surfaceRegions[sr->getSide()] = sr->getRegion();
      ed = ed->getDecorated(SURFACE_REGION);
    }

    // for all vertices
233
    for (int i = 0; i < dim + 1; i++) {
234
235
236
237
      // get coords of this vertex
      vertexCoords = elInfo->getCoord(i);
      
      // get dof index of this vertex
238
      vertexDOF = dof[i][nPreDofs];
239
     
240
      // search for coords at this dof
Thomas Witkowski's avatar
Thomas Witkowski committed
241
      std::list<VertexInfo>::iterator it =
242
243
	  find((*vertexInfos)[vertexDOF].begin(),
	       (*vertexInfos)[vertexDOF].end(),
244
245
246
	       vertexCoords);
      
      // coords not yet in list?
247
      if (it == (*vertexInfos)[vertexDOF].end()) {
248
	// create new vertex info and increment number of vertices
249
	VertexInfo newVertexInfo = {vertexCoords, nVertices};
250
251

	// add new vertex info to list
252
	(*vertexInfos)[vertexDOF].push_front(newVertexInfo);
253
254
	
	// set iterator to new vertex info
255
	it = (*vertexInfos)[vertexDOF].begin();	
256

257
	nVertices++;
258
259
260
261
262
263
264
265
      }
      
      // fill element info
      elementInfo.vertexInfo[i] = it;
      elementInfo.boundary[i] = elInfo->getBoundary(i);
      elementInfo.projection[i] = elInfo->getProjection(i);
      elementInfo.neighbour[i] = 
	elInfo->getNeighbour(i) ?
266
	outputIndices[elInfo->getNeighbour(i)->getIndex()] :
267
268
269
	-1;
    }
    
270
    if (dim == 3)
271
272
273
      elementInfo.type = (dynamic_cast<ElInfo3d*>(elInfo)->getType());

    // remember element info
274
    elements.push_back(elementInfo);
275
276
  }

277
  void DataCollector::addValueData(ElInfo *elInfo)
278
  {
279
    FUNCNAME("DataCollector::addValueData()");
280
    
281
    basisFcts->getLocalIndices(elInfo->getElement(), localAdmin, localDOFs);
282
283

    // First, traverse all DOFs at the vertices of the element, determine
284
285
286
    // their coordinates and add them to the corresponding entry in dofCoords.
    for (int i = 0; i < mesh->getGeo(VERTEX); i++) {
      DegreeOfFreedom dofi = localDOFs[i];
Thomas Witkowski's avatar
Thomas Witkowski committed
287

288
      (*interpPointInd)[dofi] = -2; // mark as vertex
289
290
      
      // get coords of this vertex
Thomas Witkowski's avatar
Thomas Witkowski committed
291
      *vertexCoords = elInfo->getCoord(i);
292

293
      // search for coords at this dof
Thomas Witkowski's avatar
Thomas Witkowski committed
294
      std::list<WorldVector<double> >::iterator it =
295
296
	  find((*dofCoord)[dofi].begin(),
	       (*dofCoord)[dofi].end(),
Thomas Witkowski's avatar
Thomas Witkowski committed
297
	       *vertexCoords);
298
299
      
      // coords not yet in list?
300
      if (it == (*dofCoord)[dofi].end()) {
301
	// add new coords to list
302
	(*dofCoord)[dofi].push_back(*vertexCoords);
303
304
      }
    }
305
   
306
307
    // Then, traverse all interpolation DOFs of the element, determine
    // their coordinates and add them to the corresponding entry in 
308
309
310
    // interpPointCoords.
    for (int i = mesh->getGeo(VERTEX); i < nBasisFcts; i++) {
      DegreeOfFreedom dofi = localDOFs[i];
311

312
      elInfo->coordToWorld(*basisFcts->getCoords(i), *vertexCoords);
313
      
314
      if ((*interpPointInd)[dofi] == -1) {
315
	// mark as interpolation point
316
	(*interpPointInd)[dofi] = -3; 
317
318
319
	
	// search for interpolation point coordinates, and insert them to the 
	// dof-entry, if not contained in the list
Thomas Witkowski's avatar
Thomas Witkowski committed
320
	std::list<WorldVector<double> >::iterator it =
321
322
		find((*interpPointCoords)[dofi].begin(),
		     (*interpPointCoords)[dofi].end(),
Thomas Witkowski's avatar
Thomas Witkowski committed
323
		     *vertexCoords);
324
	
325
326
327
	if (it == (*interpPointCoords)[dofi].end()) {
	  (*interpPointCoords)[dofi].push_back(*vertexCoords); 
	  nInterpPoints++;
328
	}
329
      }      
330
    }  
331
332
  }

333
  void DataCollector::addInterpData(ElInfo *elInfo)
334
  {
335
    FUNCNAME("DataCollector::addInterpData()");
336
    
337
    basisFcts->getLocalIndices(elInfo->getElement(), localAdmin, localDOFs);
338

339
    std::vector<int> elemInterpPoints(0);
340
341
    for (int i = mesh->getGeo(VERTEX); i < nBasisFcts; i++)
      elemInterpPoints.push_back((*interpPointInd)[localDOFs[i]]);
Thomas Witkowski's avatar
Thomas Witkowski committed
342
    
343
    interpPoints.push_back(elemInterpPoints); 
344
345
  }

346
347
  void DataCollector::addPeriodicData(ElInfo *elInfo) 
  {
348
349
350
351
352
353
354
    FUNCNAME("DataCollector::addPeriodicData");

    LeafDataPeriodic *ldp = dynamic_cast<LeafDataPeriodic*>
      (elInfo->getElement()->
       getElementData()->
       getElementData(PERIODIC));
    
355
    if (ldp) {
Thomas Witkowski's avatar
Thomas Witkowski committed
356
      std::list<LeafDataPeriodic::PeriodicInfo>::iterator it;
357
      
358
359
360
361
      for (it = dynamic_cast<LeafDataPeriodic*>(ldp)->getInfoList().begin();
	   it != dynamic_cast<LeafDataPeriodic*>(ldp)->getInfoList().end();
	   ++it) {

362
363
	int outputIndex = outputIndices[elInfo->getElement()->getIndex()];
	int neighIndex  = outputIndices[elInfo->
364
365
366
					 getNeighbour(it->elementSide)->
					 getIndex()];
	
367
	if (!periodicConnections[outputIndex][it->elementSide]) {
368
	  PeriodicInfo periodicInfo;
369
	  
370
371
372
373
374
375
	  periodicInfo.mode = it->periodicMode;
	  periodicInfo.type = it->type;
	  periodicInfo.outputIndex = outputIndex;
	  periodicInfo.neighIndex = neighIndex;
	  periodicInfo.vertexMap.clear();
	  
376
377
	  periodicConnections[outputIndex][it->elementSide] = true;
	  periodicConnections
378
379
380
	    [neighIndex][elInfo->getOppVertex(it->elementSide)] = true;
	    
  
381
382
	  for (int i = 0; i < dim; i++) {
	    int index1 = elInfo->getElement()->getVertexOfPosition(INDEX_OF_DIM(dim - 1, dim),
383
384
								   it->elementSide,
								   i);
385
	    int dof1 = elInfo->getElement()->getDOF(index1, nPreDofs);
386
	    
387
388
	    for (int j = 0; j < dim; j++) {
	      int index2 = elInfo->getElement()->getVertexOfPosition(INDEX_OF_DIM(dim - 1, dim),
389
390
								     elInfo->getOppVertex(it->elementSide),
								     j);
391
	      int dof2 = elInfo->getNeighbour(it->elementSide)->getDOF(index2, nPreDofs);
392
	      
393
	      if ((dof1 == dof2) || (mesh->associated(dof1, dof2))) {
394
395
396
		periodicInfo.vertexMap[index1] = index2;
		break;
	      } 
397
398
	    }
	  }
399
	  
400
	  periodicInfos.push_back(periodicInfo);
401
	}
402
      }
403
404
405
    }
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
406
  std::list<ElementInfo>* DataCollector::getElementInfos()
407
  {        
408
    if (!elementDataCollected)
409
410
      startCollectingElementData();
      
411
    return &elements;
412
413
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
414
  DOFVector< std::list<VertexInfo> >* DataCollector::getVertexInfos()
415
  {
416
    if (!elementDataCollected)
417
418
      startCollectingElementData();

419
    return vertexInfos;
420
421
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
422
  std::list<PeriodicInfo>* DataCollector::getPeriodicInfos()
423
  {
424
    if (!periodicDataCollected)
425
426
      startCollectingPeriodicData();

427
    return &periodicInfos;
428
429
430
431
  }

  int DataCollector::getNumberVertices()
  {
432
    if (!elementDataCollected)
433
434
      startCollectingElementData();

435
    return nVertices;
436
437
438
439
  }

  int DataCollector::getNumberElements()
  {
440
    if (!elementDataCollected)
441
442
      startCollectingElementData();

443
    return nElements;
444
445
446
447
  }

  int DataCollector::getNumberInterpPoints()
  {
448
    if (!valueDataCollected)
449
450
      startCollectingValueData();

451
    return nInterpPoints;
452
453
454
455
  }
  
  int DataCollector::getNumberConnections()
  {
456
    if (!periodicDataCollected)
457
458
      startCollectingPeriodicData();

459
    return nConnection;
460
461
462
463
  }

  const FiniteElemSpace* DataCollector::getFeSpace()
  {
464
    return feSpace;
465
466
467
468
  }

  Mesh* DataCollector::getMesh()
  {
469
    return mesh;
470
471
472
473
  }

  DOFVector<double>* DataCollector::getValues()
  {
474
    if (!valueDataCollected)
475
476
      startCollectingValueData();

477
    return values;
478
479
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
480
  DOFVector< std::list<WorldVector<double> > >* DataCollector::getDofCoords()
481
  {
482
    if (!valueDataCollected)
483
484
      startCollectingValueData();

485
    return dofCoord;
486
487
488
489
  }

  DOFVector<int>* DataCollector::getInterpPointInd()
  {
490
    if (!valueDataCollected)
491
492
      startCollectingValueData();

493
    return interpPointInd;
494
495
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
496
  DOFVector< std::list<WorldVector<double> > >* DataCollector::getInterpPointCoords()
497
  {
498
    if (!valueDataCollected)
499
500
      startCollectingValueData();

501
    return interpPointCoords;
502
503
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
504
  std::vector< std::vector<int> >* DataCollector::getInterpPoints()
505
  {
506
    if (!valueDataCollected)
507
508
      startCollectingValueData();

509
    return &interpPoints;
510
511
  }
}