DataCollector.cc 13.7 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
  DataCollector::DataCollector(const FiniteElemSpace *feSpace,
			       DOFVector<double> *values,
			       int level,
			       Flag traverseFlag,
			       bool (*writeElem)(ElInfo*))
    : writeElem_(writeElem)
  {
17
    FUNCNAME("DataCollector::DataCollector()");
18
19
20

    TEST_EXIT(feSpace)("no feSpace\n");
    TEST_EXIT(values)("no value Vector\n");
21
   
22
23
24
25
26
27
    // === get mesh  ===    
    mesh_ = feSpace->getMesh();

    // === get admin ===    
    localAdmin_ = const_cast<DOFAdmin*>(feSpace->getAdmin());
    // === create vertex info vector ===
Thomas Witkowski's avatar
Thomas Witkowski committed
28
    vertexInfos_ = new DOFVector< std::list<VertexInfo> >(feSpace, "vertex infos");
29

Thomas Witkowski's avatar
Thomas Witkowski committed
30
31
32
    interpPointInd_ = new DOFVector<int>(feSpace, "interpolation point indices");
    interpPointCoords_ = new DOFVector< std::list<WorldVector<double> > >(feSpace, "interpolation point coordinates");
    dofCoords_ = new DOFVector< std::list<WorldVector<double> > >(feSpace, "dof coords");
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49

    dim_ = mesh_->getDim();
    
    nPreDofs_ = localAdmin_->getNumberOfPreDOFs(VERTEX);
    nVertices_ = 0;
    nElements_ = 0;
    nInterpPoints_ = 0;
    nConnections_ = 0;

    feSpace_ = feSpace;
    values_ = values;
    level_ = level;
    traverseFlag_ = traverseFlag;

    elementDataCollected_ = false;
    valueDataCollected_ = false;
    periodicDataCollected_ = false;
Thomas Witkowski's avatar
Thomas Witkowski committed
50

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

  DataCollector::~DataCollector()
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
56
57
58
59
60
    delete vertexInfos_;
    delete interpPointInd_;
    delete interpPointCoords_;
    delete dofCoords_;
    delete vertexCoords;
61
62
  }

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

68
    if (!periodicDataCollected_)
Thomas Witkowski's avatar
Thomas Witkowski committed
69
      startCollectingPeriodicData();
70

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

75
  void DataCollector::startCollectingElementData()
76
  {
77
    FUNCNAME("DataCollector::startCollectingElementData()");
78
79
80
81
82
83
84
85
86
87
88
89
90
91

    Flag flag = traverseFlag_;
    flag |= 
      Mesh::FILL_NEIGH      | 
      Mesh::FILL_COORDS     |
      Mesh::FILL_OPP_COORDS |
      Mesh::FILL_BOUND;

    elements_.resize(0, dim_);

    TraverseStack stack;

    // Traverse elements to create continuous element indices
    ElInfo *elInfo = stack.traverseFirst(mesh_, level_, flag);
92
93
    while (elInfo) {
      if (!writeElem_ || writeElem_(elInfo))
94
95
96
97
98
99
	outputIndices_[elInfo->getElement()->getIndex()] = nElements_++;
      elInfo = stack.traverseNext(elInfo);
    }

    // Traverse elements to create element information
    elInfo = stack.traverseFirst(mesh_, level_, flag);
100

101
    while (elInfo) {
102
      if (!writeElem_ || writeElem_(elInfo))
103
	addElementData(elInfo);
104

105
106
107
108
109
110
      elInfo = stack.traverseNext(elInfo);
    }

    elementDataCollected_ = true;
  }

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

    DOFVector<int>::Iterator intPointIt(interpPointInd_, USED_DOFS);
116
    for (intPointIt.reset(); !intPointIt.end(); ++intPointIt)
117
118
119
120
      (*intPointIt) = -1;

    interpPoints_.clear();

Thomas Witkowski's avatar
Thomas Witkowski committed
121
122
    basisFcts_ = const_cast<BasisFunction*>(feSpace_->getBasisFcts());
    nBasisFcts_ = basisFcts_->getNumber();
123
    localDOFs_ = new DegreeOfFreedom[nBasisFcts_];
124
  
125
126
127
128
    TraverseStack stack;

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

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

160
  void DataCollector::startCollectingPeriodicData()
161
  {    
162
    FUNCNAME("DataCollector::startCollectingPeriodicData()");
163
164
165
166
167

    periodicConnections_.clear();
    
    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(mesh_, level_, traverseFlag_);
168
169
    while (elInfo) {
      if (!writeElem_ || writeElem_(elInfo)) {
170
171
172
173
174
	LeafDataPeriodic *ldp = dynamic_cast<LeafDataPeriodic*>
	  (elInfo->getElement()->
	   getElementData()->
	   getElementData(PERIODIC));
	
175
176
	if (ldp) {
	  nConnections_ += dynamic_cast<LeafDataPeriodic*>(ldp)->getInfoList().size();
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
	}

	periodicConnections_.push_back(DimVec<bool>(dim_, DEFAULT_VALUE, false));
      }
      elInfo = stack.traverseNext(elInfo);
    }   

    nConnections_ /= 2;

    periodicInfos_.clear();

    Flag flag = traverseFlag_;
    flag |= 
      Mesh::FILL_COORDS    |
      Mesh::FILL_OPP_COORDS|
      Mesh::FILL_NEIGH     | 
      Mesh::FILL_BOUND;

    elInfo = stack.traverseFirst(mesh_, level_, flag);
196
197
    while (elInfo) {
      if (!writeElem_ || writeElem_(elInfo))
198
199
200
201
202
203
204
	addPeriodicData(elInfo);
      elInfo = stack.traverseNext(elInfo);
    }   

    periodicDataCollected_ = true;
  }

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

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

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

    // for all vertices
235
    for (int i = 0; i < dim_ + 1; i++) {
236
237
238
239
240
      // get coords of this vertex
      vertexCoords = elInfo->getCoord(i);
      
      // get dof index of this vertex
      vertexDOF = dof[i][nPreDofs_];
241
     
242
      // search for coords at this dof
Thomas Witkowski's avatar
Thomas Witkowski committed
243
      std::list<VertexInfo>::iterator it =
244
245
246
247
248
	  find((*vertexInfos_)[vertexDOF].begin(),
	       (*vertexInfos_)[vertexDOF].end(),
	       vertexCoords);
      
      // coords not yet in list?
249
      if (it == (*vertexInfos_)[vertexDOF].end()) {
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
	// create new vertex info and increment number of vertices
	VertexInfo newVertexInfo = {vertexCoords, nVertices_};

	// add new vertex info to list
	(*vertexInfos_)[vertexDOF].push_front(newVertexInfo);
	
	// set iterator to new vertex info
	it = (*vertexInfos_)[vertexDOF].begin();	

	nVertices_++;
      }
      
      // 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) ?
	outputIndices_[elInfo->getNeighbour(i)->getIndex()] :
	-1;
    }
    
272
    if (dim_ == 3)
273
274
275
276
277
278
      elementInfo.type = (dynamic_cast<ElInfo3d*>(elInfo)->getType());

    // remember element info
    elements_.push_back(elementInfo);
  }

279
  void DataCollector::addValueData(ElInfo *elInfo)
280
  {
281
    FUNCNAME("DataCollector::addValueData()");
282
    
Thomas Witkowski's avatar
Thomas Witkowski committed
283
    basisFcts_->getLocalIndices(elInfo->getElement(), localAdmin_, localDOFs_);
284
285

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

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

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

314
      elInfo->coordToWorld(*basisFcts_->getCoords(i), *vertexCoords);
315
      
Thomas Witkowski's avatar
Thomas Witkowski committed
316
      if ((*interpPointInd_)[dofi] == -1) {
317
	// mark as interpolation point
Thomas Witkowski's avatar
Thomas Witkowski committed
318
	(*interpPointInd_)[dofi] = -3; 
319
320
321
	
	// 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
322
	std::list<WorldVector<double> >::iterator it =
Thomas Witkowski's avatar
Thomas Witkowski committed
323
324
		find((*interpPointCoords_)[dofi].begin(),
		     (*interpPointCoords_)[dofi].end(),
Thomas Witkowski's avatar
Thomas Witkowski committed
325
		     *vertexCoords);
326
	
Thomas Witkowski's avatar
Thomas Witkowski committed
327
	if (it == (*interpPointCoords_)[dofi].end()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
328
	  (*interpPointCoords_)[dofi].push_back(*vertexCoords); 
329
	  nInterpPoints_++;
330
	}
331
      }      
332
    }  
333
334
  }

335
  void DataCollector::addInterpData(ElInfo *elInfo)
336
  {
337
    FUNCNAME("DataCollector::addInterpData()");
338
    
Thomas Witkowski's avatar
Thomas Witkowski committed
339
    basisFcts_->getLocalIndices(elInfo->getElement(), localAdmin_, localDOFs_);
340

341
    std::vector<int> elemInterpPoints(0);
342
    for (int i = mesh_->getGeo(VERTEX); i < nBasisFcts_; i++)
343
      elemInterpPoints.push_back((*interpPointInd_)[localDOFs_[i]]);
Thomas Witkowski's avatar
Thomas Witkowski committed
344
    
345
    interpPoints_.push_back(elemInterpPoints); 
346
347
  }

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

    LeafDataPeriodic *ldp = dynamic_cast<LeafDataPeriodic*>
      (elInfo->getElement()->
       getElementData()->
       getElementData(PERIODIC));
    
357
    if (ldp) {
Thomas Witkowski's avatar
Thomas Witkowski committed
358
      std::list<LeafDataPeriodic::PeriodicInfo>::iterator it;
359
      
360
361
362
363
364
365
366
367
368
369
370
      for (it = dynamic_cast<LeafDataPeriodic*>(ldp)->getInfoList().begin();
	   it != dynamic_cast<LeafDataPeriodic*>(ldp)->getInfoList().end();
	   ++it) {

	int outputIndex = outputIndices_[elInfo->getElement()->getIndex()];
	int neighIndex  = outputIndices_[elInfo->
					 getNeighbour(it->elementSide)->
					 getIndex()];
	
	if (!periodicConnections_[outputIndex][it->elementSide]) {
	  PeriodicInfo periodicInfo;
371
	  
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
	  periodicInfo.mode = it->periodicMode;
	  periodicInfo.type = it->type;
	  periodicInfo.outputIndex = outputIndex;
	  periodicInfo.neighIndex = neighIndex;
	  periodicInfo.vertexMap.clear();
	  
	  periodicConnections_[outputIndex][it->elementSide] = true;
	  periodicConnections_
	    [neighIndex][elInfo->getOppVertex(it->elementSide)] = true;
	    
  
	  for (int i = 0; i < dim_; i++) {
	    int index1 = elInfo->getElement()->getVertexOfPosition(INDEX_OF_DIM(dim_ - 1, dim_),
								   it->elementSide,
								   i);
	    int dof1 = elInfo->getElement()->getDOF(index1, nPreDofs_);
	    
	    for (int j = 0; j < dim_; j++) {
	      int index2 = elInfo->getElement()->getVertexOfPosition(INDEX_OF_DIM(dim_ - 1, dim_),
								     elInfo->getOppVertex(it->elementSide),
								     j);
	      int dof2 = elInfo->getNeighbour(it->elementSide)->getDOF(index2, nPreDofs_);
	      
	      if ((dof1 == dof2) || (mesh_->associated(dof1, dof2))) {
		periodicInfo.vertexMap[index1] = index2;
		break;
	      } 
399
400
	    }
	  }
401
402
	  
	  periodicInfos_.push_back(periodicInfo);
403
	}
404
      }
405
406
407
    }
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
408
  std::list<ElementInfo>* DataCollector::getElementInfos()
409
  {        
410
    if (!elementDataCollected_)
411
412
413
414
415
      startCollectingElementData();
      
    return &elements_;
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
416
  DOFVector< std::list<VertexInfo> >* DataCollector::getVertexInfos()
417
  {
418
    if (!elementDataCollected_)
419
420
421
422
423
      startCollectingElementData();

    return vertexInfos_;
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
424
  std::list<PeriodicInfo>* DataCollector::getPeriodicInfos()
425
  {
426
    if (!periodicDataCollected_)
427
428
429
430
431
432
433
      startCollectingPeriodicData();

    return &periodicInfos_;
  }

  int DataCollector::getNumberVertices()
  {
434
    if (!elementDataCollected_)
435
436
437
438
439
440
441
      startCollectingElementData();

    return nVertices_;
  }

  int DataCollector::getNumberElements()
  {
442
    if (!elementDataCollected_)
443
444
445
446
447
448
449
      startCollectingElementData();

    return nElements_;
  }

  int DataCollector::getNumberInterpPoints()
  {
450
    if (!valueDataCollected_)
451
452
453
454
455
456
457
      startCollectingValueData();

    return nInterpPoints_;
  }
  
  int DataCollector::getNumberConnections()
  {
458
    if (!periodicDataCollected_)
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
      startCollectingPeriodicData();

    return nConnections_;
  }

  const FiniteElemSpace* DataCollector::getFeSpace()
  {
    return feSpace_;
  }

  Mesh* DataCollector::getMesh()
  {
    return mesh_;
  }

  DOFVector<double>* DataCollector::getValues()
  {
476
    if (!valueDataCollected_)
477
478
479
480
481
      startCollectingValueData();

    return values_;
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
482
  DOFVector< std::list<WorldVector<double> > >* DataCollector::getDofCoords()
483
  {
484
    if (!valueDataCollected_)
485
486
487
488
489
490
491
      startCollectingValueData();

    return dofCoords_;
  }

  DOFVector<int>* DataCollector::getInterpPointInd()
  {
492
    if (!valueDataCollected_)
493
494
495
496
497
      startCollectingValueData();

    return interpPointInd_;
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
498
  DOFVector< std::list<WorldVector<double> > >* DataCollector::getInterpPointCoords()
499
  {
500
    if (!valueDataCollected_)
501
502
503
504
505
      startCollectingValueData();

    return interpPointCoords_;
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
506
  std::vector< std::vector<int> >* DataCollector::getInterpPoints()
507
  {
508
    if (!valueDataCollected_)
509
510
511
512
513
      startCollectingValueData();

    return &interpPoints_;
  }
}