DataCollector.cc 13.7 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
18
19
20
#include "DataCollector.h"
#include "Traverse.h"
#include "DOFVector.h"
#include "SurfaceRegion_ED.h"
#include "ElementRegion_ED.h"
#include "Projection.h"

namespace AMDiS {
21

22
23
24
25
26
27
28
29
30
  DataCollector::DataCollector(const FiniteElemSpace *fe,
			       DOFVector<double> *val,
			       int l,
			       Flag flag,
			       bool (*writeElemFct)(ElInfo*))
    : values(val),
      level(l),
      traverseFlag(flag),
      feSpace(fe),
31
32
33
34
35
36
37
      nVertices(0),
      nElements(0),
      nInterpPoints(0),
      nConnection(0),
      elementDataCollected(false),
      valueDataCollected(false),
      periodicDataCollected(false),
38
      writeElem(writeElemFct)
39
  {
40
    FUNCNAME("DataCollector::DataCollector()");
41

42
    TEST_EXIT(feSpace)("No finite elem space defined!\n");
43
   
44
    // === get mesh  ===    
45
    mesh = feSpace->getMesh();
46
47

    // === get admin ===    
48
    localAdmin = const_cast<DOFAdmin*>(feSpace->getAdmin());
49
    // === create vertex info vector ===
50
    vertexInfos = new DOFVector<std::list<VertexInfo> >(feSpace, "vertex infos");
51

52
53
54
    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");
55

56
    dim = mesh->getDim();    
57
    nPreDofs = localAdmin->getNumberOfPreDofs(VERTEX);
Thomas Witkowski's avatar
Thomas Witkowski committed
58
  } 
59

60

61
62
  DataCollector::~DataCollector()
  {
63
64
65
66
    delete vertexInfos;
    delete interpPointInd;
    delete interpPointCoords;
    delete dofCoord;
67
68
  }

69

Thomas Witkowski's avatar
Thomas Witkowski committed
70
71
  void DataCollector::fillAllData()
  {
72
    if (!elementDataCollected)
Thomas Witkowski's avatar
Thomas Witkowski committed
73
      startCollectingElementData();
74

75
    if (!periodicDataCollected)
Thomas Witkowski's avatar
Thomas Witkowski committed
76
      startCollectingPeriodicData();
77

78
    if (!valueDataCollected)
Thomas Witkowski's avatar
Thomas Witkowski committed
79
80
81
      startCollectingValueData();
  }

82

83
  void DataCollector::startCollectingElementData()
84
  {
85
    FUNCNAME("DataCollector::startCollectingElementData()");
86

87
    Flag flag = traverseFlag;
88
89
90
91
92
93
    flag |= 
      Mesh::FILL_NEIGH      | 
      Mesh::FILL_COORDS     |
      Mesh::FILL_OPP_COORDS |
      Mesh::FILL_BOUND;

94
    elements.resize(0, dim);
95
96
97
98

    TraverseStack stack;

    // Traverse elements to create continuous element indices
99
    ElInfo *elInfo = stack.traverseFirst(mesh, level, flag);
100
    while (elInfo) {
101
102
      if (!writeElem || writeElem(elInfo))
	outputIndices[elInfo->getElement()->getIndex()] = nElements++;
103
104
105
106
      elInfo = stack.traverseNext(elInfo);
    }

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

109
    while (elInfo) {
110
      if (!writeElem || writeElem(elInfo))
111
	addElementData(elInfo);
112

113
114
115
      elInfo = stack.traverseNext(elInfo);
    }

116
    elementDataCollected = true;
117
118
  }

119

120
  void DataCollector::startCollectingValueData()
121
  {
122
    FUNCNAME("DataCollector::startCollectingValueData()");
123

124
    DOFVector<int>::Iterator intPointIt(interpPointInd, USED_DOFS);
125
    for (intPointIt.reset(); !intPointIt.end(); ++intPointIt)
126
127
      (*intPointIt) = -1;

128
    interpPoints.clear();
129

130
131
132
    basisFcts = const_cast<BasisFunction*>(feSpace->getBasisFcts());
    nBasisFcts = basisFcts->getNumber();
    localDOFs = new DegreeOfFreedom[nBasisFcts];
133
  
134
135
136
137
    TraverseStack stack;

    // Traverse elements to add value information and to mark
    // interpolation points.
138
139
    ElInfo *elInfo = stack.traverseFirst(mesh, level, 
					 traverseFlag | Mesh::FILL_COORDS);
140
    while (elInfo) {
141
      if (!writeElem || writeElem(elInfo))
Thomas Witkowski's avatar
Thomas Witkowski committed
142
	addValueData(elInfo);
143
144
145
      elInfo = stack.traverseNext(elInfo);
    }

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

169

170
  void DataCollector::startCollectingPeriodicData()
171
  {    
172
    FUNCNAME("DataCollector::startCollectingPeriodicData()");
173

174
    periodicConnections.clear();
175
176
    
    TraverseStack stack;
177
    ElInfo *elInfo = stack.traverseFirst(mesh, level, traverseFlag);
178
    while (elInfo) {
179
      if (!writeElem || writeElem(elInfo)) {
180
181
182
183
184
	LeafDataPeriodic *ldp = dynamic_cast<LeafDataPeriodic*>
	  (elInfo->getElement()->
	   getElementData()->
	   getElementData(PERIODIC));
	
185
	if (ldp)
186
	  nConnection += dynamic_cast<LeafDataPeriodic*>(ldp)->getInfoList().size();
187

188
	periodicConnections.push_back(DimVec<bool>(dim, DEFAULT_VALUE, false));
189
190
191
192
      }
      elInfo = stack.traverseNext(elInfo);
    }   

193
    nConnection /= 2;
194

195
    periodicInfos.clear();
196

197
    Flag flag = traverseFlag;
198
199
200
201
202
203
    flag |= 
      Mesh::FILL_COORDS    |
      Mesh::FILL_OPP_COORDS|
      Mesh::FILL_NEIGH     | 
      Mesh::FILL_BOUND;

204
    elInfo = stack.traverseFirst(mesh, level, flag);
205
    while (elInfo) {
206
      if (!writeElem || writeElem(elInfo))
207
208
209
210
	addPeriodicData(elInfo);
      elInfo = stack.traverseNext(elInfo);
    }   

211
    periodicDataCollected = true;
212
213
  }

214

215
  void DataCollector::addElementData(ElInfo* elInfo)
216
  {
217
    FUNCNAME("DataCollector::addElementData()");
218

219
    const DegreeOfFreedom **dof = elInfo->getElement()->getDof();
220

221
222
223
224
//     MSG("ELEMENT: %d %d\n", 
// 	elInfo->getElement()->getIndex(),
// 	elInfo->getMacroElement()->getIndex());    
//     MSG("DOFs: %d %d %d %d\n", dof[0][0], dof[1][0], dof[2][0], dof[3][0]);
225

226
    // create ElementInfo
227
    ElementInfo elementInfo(dim);
228
229
    
    // read element region
230
    ElementData *ed = elInfo->getElement()->getElementData(ELEMENT_REGION);
231
    
Thomas Witkowski's avatar
Thomas Witkowski committed
232
    if (ed)
233
      elementInfo.elementRegion = dynamic_cast<ElementRegion_ED*>(ed)->getRegion();
Thomas Witkowski's avatar
Thomas Witkowski committed
234
235
    else 
      elementInfo.elementRegion = -1;    
236
   
237
238
239
    // read surface regions to element info
    ed = elInfo->getElement()->getElementData(SURFACE_REGION);
    elementInfo.surfaceRegions.set(-1);
240
    while (ed) {
241
242
243
244
245
246
      SurfaceRegion_ED *sr = dynamic_cast<SurfaceRegion_ED*>(ed);
      elementInfo.surfaceRegions[sr->getSide()] = sr->getRegion();
      ed = ed->getDecorated(SURFACE_REGION);
    }

    // for all vertices
247
    for (int i = 0; i < mesh->getGeo(VERTEX); i++) {
248
      // get coords of this vertex
249
      WorldVector<double> &vertexCoords = elInfo->getCoord(i);
250
251
      
      // get dof index of this vertex
252
      DegreeOfFreedom vertexDof = dof[i][nPreDofs];
253

254
      // search for coords at this dof
Thomas Witkowski's avatar
Thomas Witkowski committed
255
      std::list<VertexInfo>::iterator it =
256
	  find((*vertexInfos)[vertexDof].begin(), (*vertexInfos)[vertexDof].end(),
257
258
259
	       vertexCoords);
      
      // coords not yet in list?
260
      if (it == (*vertexInfos)[vertexDof].end()) {
261
	// create new vertex info and increment number of vertices
262
	VertexInfo newVertexInfo = {vertexCoords, nVertices};
263
264

	// add new vertex info to list
265
	(*vertexInfos)[vertexDof].push_front(newVertexInfo);
266

267
	// set iterator to new vertex info
268
	it = (*vertexInfos)[vertexDof].begin();	
269

270
	nVertices++;
271
272
273
274
275
276
277
278
      }
      
      // 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) ?
279
	outputIndices[elInfo->getNeighbour(i)->getIndex()] :
280
281
282
	-1;
    }
    
283
    if (dim == 3)
284
      elementInfo.type = elInfo->getType();
285
286

    // remember element info
287
    elements.push_back(elementInfo);
288
289
  }

290

291
  void DataCollector::addValueData(ElInfo *elInfo)
292
  {
293
    FUNCNAME("DataCollector::addValueData()");
294
    
295
    WorldVector<double> vertexCoords;
296
    basisFcts->getLocalIndices(elInfo->getElement(), localAdmin, localDOFs);
297
298

    // First, traverse all DOFs at the vertices of the element, determine
299
300
301
    // 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
302

303
      (*interpPointInd)[dofi] = -2; // mark as vertex
304
305
      
      // get coords of this vertex
306
      vertexCoords = elInfo->getCoord(i);
307

308
      // search for coords at this dof
Thomas Witkowski's avatar
Thomas Witkowski committed
309
      std::list<WorldVector<double> >::iterator it =
310
311
	  find((*dofCoord)[dofi].begin(),
	       (*dofCoord)[dofi].end(),
312
	       vertexCoords);
313
314
      
      // coords not yet in list?
315
      if (it == (*dofCoord)[dofi].end()) {
316
	// add new coords to list
317
	(*dofCoord)[dofi].push_back(vertexCoords);
318
319
      }
    }
320
   
321
322
    // Then, traverse all interpolation DOFs of the element, determine
    // their coordinates and add them to the corresponding entry in 
323
324
325
    // interpPointCoords.
    for (int i = mesh->getGeo(VERTEX); i < nBasisFcts; i++) {
      DegreeOfFreedom dofi = localDOFs[i];
326

327
      elInfo->coordToWorld(*basisFcts->getCoords(i), vertexCoords);
328
      
329
      if ((*interpPointInd)[dofi] == -1) {
330
	// mark as interpolation point
331
	(*interpPointInd)[dofi] = -3; 
332
333
334
	
	// 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
335
	std::list<WorldVector<double> >::iterator it =
336
		find((*interpPointCoords)[dofi].begin(), 
337
		     (*interpPointCoords)[dofi].end(),
338
		     vertexCoords);
339
	
340
	if (it == (*interpPointCoords)[dofi].end()) {
341
	  (*interpPointCoords)[dofi].push_back(vertexCoords); 
342
	  nInterpPoints++;
343
	}
344
      }      
345
    }  
346
347
  }

348

349
  void DataCollector::addInterpData(ElInfo *elInfo)
350
  {
351
    FUNCNAME("DataCollector::addInterpData()");
352
    
353
    basisFcts->getLocalIndices(elInfo->getElement(), localAdmin, localDOFs);
354

355
    std::vector<int> elemInterpPoints(0);
356
357
    for (int i = mesh->getGeo(VERTEX); i < nBasisFcts; i++)
      elemInterpPoints.push_back((*interpPointInd)[localDOFs[i]]);
Thomas Witkowski's avatar
Thomas Witkowski committed
358
    
359
    interpPoints.push_back(elemInterpPoints); 
360
361
  }

362

363
364
  void DataCollector::addPeriodicData(ElInfo *elInfo) 
  {
365
366
367
368
369
370
371
    FUNCNAME("DataCollector::addPeriodicData");

    LeafDataPeriodic *ldp = dynamic_cast<LeafDataPeriodic*>
      (elInfo->getElement()->
       getElementData()->
       getElementData(PERIODIC));
    
372
    if (ldp) {
Thomas Witkowski's avatar
Thomas Witkowski committed
373
      std::list<LeafDataPeriodic::PeriodicInfo>::iterator it;
374
      
375
376
377
378
      for (it = dynamic_cast<LeafDataPeriodic*>(ldp)->getInfoList().begin();
	   it != dynamic_cast<LeafDataPeriodic*>(ldp)->getInfoList().end();
	   ++it) {

379
380
	int outputIndex = outputIndices[elInfo->getElement()->getIndex()];
	int neighIndex  = outputIndices[elInfo->
381
382
383
					 getNeighbour(it->elementSide)->
					 getIndex()];
	
384
	if (!periodicConnections[outputIndex][it->elementSide]) {
385
	  PeriodicInfo periodicInfo;
386
	  
387
388
389
390
391
392
	  periodicInfo.mode = it->periodicMode;
	  periodicInfo.type = it->type;
	  periodicInfo.outputIndex = outputIndex;
	  periodicInfo.neighIndex = neighIndex;
	  periodicInfo.vertexMap.clear();
	  
393
394
	  periodicConnections[outputIndex][it->elementSide] = true;
	  periodicConnections
395
396
397
	    [neighIndex][elInfo->getOppVertex(it->elementSide)] = true;
	    
  
398
399
	  for (int i = 0; i < dim; i++) {
	    int index1 = elInfo->getElement()->getVertexOfPosition(INDEX_OF_DIM(dim - 1, dim),
400
401
								   it->elementSide,
								   i);
402
	    int dof1 = elInfo->getElement()->getDof(index1, nPreDofs);
403
	    
404
405
	    for (int j = 0; j < dim; j++) {
	      int index2 = elInfo->getElement()->getVertexOfPosition(INDEX_OF_DIM(dim - 1, dim),
406
407
								     elInfo->getOppVertex(it->elementSide),
								     j);
408
	      int dof2 = elInfo->getNeighbour(it->elementSide)->getDof(index2, nPreDofs);
409
	      
410
	      if ((dof1 == dof2) || (mesh->associated(dof1, dof2))) {
411
412
413
		periodicInfo.vertexMap[index1] = index2;
		break;
	      } 
414
415
	    }
	  }
416
	  
417
	  periodicInfos.push_back(periodicInfo);
418
	}
419
      }
420
421
422
    }
  }

423

Thomas Witkowski's avatar
Thomas Witkowski committed
424
  std::list<ElementInfo>* DataCollector::getElementInfos()
425
  {        
426
    if (!elementDataCollected)
427
428
      startCollectingElementData();
      
429
    return &elements;
430
431
  }

432

Thomas Witkowski's avatar
Thomas Witkowski committed
433
  DOFVector< std::list<VertexInfo> >* DataCollector::getVertexInfos()
434
  {
435
    if (!elementDataCollected)
436
437
      startCollectingElementData();

438
    return vertexInfos;
439
440
  }

441

Thomas Witkowski's avatar
Thomas Witkowski committed
442
  std::list<PeriodicInfo>* DataCollector::getPeriodicInfos()
443
  {
444
    if (!periodicDataCollected)
445
446
      startCollectingPeriodicData();

447
    return &periodicInfos;
448
449
  }

450

451
452
  int DataCollector::getNumberVertices()
  {
453
    if (!elementDataCollected)
454
455
      startCollectingElementData();

456
    return nVertices;
457
458
  }

459

460
461
  int DataCollector::getNumberElements()
  {
462
    if (!elementDataCollected)
463
464
      startCollectingElementData();

465
    return nElements;
466
467
  }

468

469
470
  int DataCollector::getNumberInterpPoints()
  {
471
    if (!valueDataCollected)
472
473
      startCollectingValueData();

474
    return nInterpPoints;
475
  }
476

477
478
479
  
  int DataCollector::getNumberConnections()
  {
480
    if (!periodicDataCollected)
481
482
      startCollectingPeriodicData();

483
    return nConnection;
484
485
486
487
488
  }


  DOFVector<double>* DataCollector::getValues()
  {
489
    if (!valueDataCollected)
490
491
      startCollectingValueData();

492
    return values;
493
494
  }

495

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

501
    return dofCoord;
502
503
  }

504

505
506
  DOFVector<int>* DataCollector::getInterpPointInd()
  {
507
    if (!valueDataCollected)
508
509
      startCollectingValueData();

510
    return interpPointInd;
511
512
  }

513

Thomas Witkowski's avatar
Thomas Witkowski committed
514
  DOFVector< std::list<WorldVector<double> > >* DataCollector::getInterpPointCoords()
515
  {
516
    if (!valueDataCollected)
517
518
      startCollectingValueData();

519
    return interpPointCoords;
520
521
  }

522

Thomas Witkowski's avatar
Thomas Witkowski committed
523
  std::vector< std::vector<int> >* DataCollector::getInterpPoints()
524
  {
525
    if (!valueDataCollected)
526
527
      startCollectingValueData();

528
    return &interpPoints;
529
530
  }
}