ElementObjectDatabase.cc 31.8 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
#include "VertexVector.h"
14
#include "parallel/ElementObjectDatabase.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
15
16
17

namespace AMDiS {

18
  void ElementObjectDatabase::create(map<int, int>& rankMap, MeshLevelData& ld)
19
  {
20
21
22
23
24
    FUNCNAME("ElementObjectDatabase::create()");

    macroElementRankMap = &rankMap;
    levelData = &ld;

25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
    // === Fills macro element data structures. ===
    
    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, 
					 Mesh::CALL_LEAF_EL | 
					 Mesh::FILL_NEIGH | 
					 Mesh::FILL_BOUND);
    while (elInfo) {
      TEST_EXIT_DBG(elInfo->getLevel() == 0)("Should not happen!\n");

      Element *el = elInfo->getElement();
      macroElIndexMap.insert(make_pair(el->getIndex(), el));
      macroElIndexTypeMap.insert(make_pair(el->getIndex(), elInfo->getType()));
     
      // Add all sub object of the element to the variable elObjDb.
      addElement(elInfo);

      elInfo = stack.traverseNext(elInfo);
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
44
45
46
47
48
49
    
    // Handle periodic boudaries
    if (removePeriodicBoundary == false) 
      createPeriodicData();
    else
      removePeriodicData();
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68

    // Create data about the reverse modes of neighbouring elements.
    createReverseModeData();
  }


  void ElementObjectDatabase::createMacroElementInfo(vector<MacroElement*> &mel)
  {
    macroElIndexMap.clear();
    macroElIndexTypeMap.clear();

    for (vector<MacroElement*>::iterator it = mel.begin(); 
	 it != mel.end(); ++it) {
      macroElIndexMap.insert(make_pair((*it)->getIndex(), (*it)->getElement()));
      macroElIndexTypeMap.insert(make_pair((*it)->getIndex(), (*it)->getElType()));
    }
  }

      
69
  void ElementObjectDatabase::addElement(ElInfo *elInfo)
70
  {
71
    FUNCNAME("ElementObjectDatabase::addElement()");
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97

    TEST_EXIT_DBG(mesh)("Mesh not set!\n");

    Element *el = elInfo->getElement();


    // === First, add all element objects to the database. ===

    for (int i = 0; i < el->getGeo(VERTEX); i++)
      addVertex(el, i);
    
    for (int i = 0; i < el->getGeo(EDGE); i++)
      addEdge(el, i);      
    
    for (int i = 0; i < el->getGeo(FACE); i++)
      addFace(el, i);


    // === Get periodic boundary information. === 
    
    switch (mesh->getDim()) {
    case 2:
      for (int i = 0; i < el->getGeo(EDGE); i++) {
	if (mesh->isPeriodicAssociation(elInfo->getBoundary(EDGE, i))) {
	  // The current element's i-th edge is periodic.
	  Element *neigh = elInfo->getNeighbour(i);
98
99

	  TEST_EXIT_DBG(neigh)("Should not happen!\n");
100
101
102
103
	  
	  DofEdge edge0 = el->getEdge(i);
	  DofEdge edge1 = neigh->getEdge(elInfo->getOppVertex(i));
	  BoundaryType boundaryType = elInfo->getBoundary(EDGE, i);
104

Thomas Witkowski's avatar
Thomas Witkowski committed
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
	  // Add the periodic edge.
	  periodicEdges[make_pair(edge0, edge1)] = boundaryType;
	  periodicEdgeAssoc[edge0].insert(edge1);
	  
	  // Add both vertices of the edge to be periodic.
	  periodicVertices[make_pair(edge0.first, edge1.first)] = boundaryType;
	  periodicVertices[make_pair(edge0.second, edge1.second)] = boundaryType;
	  periodicDofAssoc[edge0.first].insert(boundaryType);
	  periodicDofAssoc[edge0.second].insert(boundaryType);
	  
	  TEST_EXIT_DBG(edge0.first == 
			mesh->getPeriodicAssociations(boundaryType)[edge1.first] &&
			edge0.second == 
			mesh->getPeriodicAssociations(boundaryType)[edge1.second])
	    ("Should not happen!\n");
120
121
122
123
124
125
126
127
	}
      }
      break;
    case 3:
      for (int i = 0; i < el->getGeo(FACE); i++) {
	if (mesh->isPeriodicAssociation(elInfo->getBoundary(FACE, i))) {
	  // The current element's i-th face is periodic.
	  Element *neigh = elInfo->getNeighbour(i);
128
129

	  TEST_EXIT_DBG(neigh)("Should not happen!\n");
130
131
132
133
	  
	  DofFace face0 = el->getFace(i);
	  DofFace face1 = neigh->getFace(elInfo->getOppVertex(i));
	  BoundaryType boundaryType = elInfo->getBoundary(FACE, i);
134

135
136
137
138
	  // Add the periodic face.
	  periodicFaces[make_pair(face0, face1)] = elInfo->getBoundary(i);
	  
	  /// Add all three vertices of the face to be periodic.
139
140
141
142
143
144
	  periodicVertices[make_pair(face0.get<0>(), face1.get<0>())] = 
	    boundaryType;
	  periodicVertices[make_pair(face0.get<1>(), face1.get<1>())] = 
	    boundaryType;
	  periodicVertices[make_pair(face0.get<2>(), face1.get<2>())] = 
	    boundaryType;
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
	  
	  periodicDofAssoc[face0.get<0>()].insert(boundaryType);
	  periodicDofAssoc[face0.get<1>()].insert(boundaryType);
	  periodicDofAssoc[face0.get<2>()].insert(boundaryType);
	  
	  TEST_EXIT_DBG(face0.get<0>() == 
			mesh->getPeriodicAssociations(boundaryType)[face1.get<0>()] &&
			face0.get<1>() == 
			mesh->getPeriodicAssociations(boundaryType)[face1.get<1>()] &&
			face0.get<2>() == 
			mesh->getPeriodicAssociations(boundaryType)[face1.get<2>()])
	    ("Should not happen!\n");
	  
	  // Create all three edges of the element and add them to be periodic.
	  DofEdge elEdge0 = make_pair(face0.get<0>(), face0.get<1>());
	  DofEdge elEdge1 = make_pair(face0.get<0>(), face0.get<2>());
	  DofEdge elEdge2 = make_pair(face0.get<1>(), face0.get<2>());
	  DofEdge neighEdge0 = make_pair(face1.get<0>(), face1.get<1>());
	  DofEdge neighEdge1 = make_pair(face1.get<0>(), face1.get<2>());
	  DofEdge neighEdge2 = make_pair(face1.get<1>(), face1.get<2>());
	  	  
	  periodicEdges[make_pair(elEdge0, neighEdge0)] = boundaryType;
	  periodicEdges[make_pair(elEdge1, neighEdge1)] = boundaryType;
	  periodicEdges[make_pair(elEdge2, neighEdge2)] = boundaryType;
	  
	  periodicEdgeAssoc[elEdge0].insert(neighEdge0);
	  periodicEdgeAssoc[elEdge1].insert(neighEdge1);
	  periodicEdgeAssoc[elEdge2].insert(neighEdge2);
	}
      }
      break;
    default:
      ERROR_EXIT("Should not happen!\n");
    }       
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
  void ElementObjectDatabase::addVertex(Element *el, int ith)
  {
    DegreeOfFreedom vertex = el->getDof(ith, 0);
    int elIndex = el->getIndex();
    ElementObjectData elObj(elIndex, ith);

    if (elIndex == 53 && ith == 0)
      MSG("A: 53/0 ON DOF %d\n", vertex);

    if (elIndex == 229 && ith == 0)
      MSG("A: 229/0 ON DOF %d\n", vertex);
    
    vertexElements[vertex].push_back(elObj);
    vertexLocalMap[elObj] = vertex;
  }


  void ElementObjectDatabase::addEdge(Element *el, int ith)
  {
    FUNCNAME("ElementObjectDatabase::addEdge()");

    DofEdge edge = el->getEdge(ith);
    int elIndex = el->getIndex();
    ElementObjectData elObj(elIndex, ith);
    
    edgeElements[edge].push_back(elObj);
    edgeLocalMap[elObj] = edge;
  }

  
  void ElementObjectDatabase::addFace(Element *el, int ith)
  {
    DofFace face = el->getFace(ith);
    int elIndex = el->getIndex();
    ElementObjectData elObj(elIndex, ith);
    
    faceElements[face].push_back(elObj);
    faceLocalMap[elObj] = face;
  }

221

222
  void ElementObjectDatabase::createPeriodicData()
223
  {
224
    FUNCNAME("ElementObjectDatabase::createPeriodicData()");
225
226
227

    TEST_EXIT_DBG(mesh)("Mesh not set!\n");

228
229
    // === Return, if there are no periodic vertices, i.e., there are no  ===
    // === periodic boundaries in the mesh.                               ===
230
231
232
    
    if (periodicVertices.size() == 0)
      return;
233
 
234
235
236
237
238
239
    // === Calculate smallest periodic boundary ID in mesh. ===

    smallestPeriodicBcType = 0;
    for (map<BoundaryType, VertexVector*>::iterator it = mesh->getPeriodicAssociations().begin();
	 it != mesh->getPeriodicAssociations().end(); ++it)
      smallestPeriodicBcType = std::min(smallestPeriodicBcType, it->first);
240
241
242
243

    // === Get all vertex DOFs that have multiple periodic associations. ===
    
    // We group all vertices together, that have either two or three periodic
244
245
246
    // associations. For rectangular domains in 2D, the four corner vertices have
    // all two periodic associations. For box domains in 3D, the eight corner
    // vertices have all three periodic associations.
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266

    vector<DegreeOfFreedom> multPeriodicDof2, multPeriodicDof3;
    for (map<DegreeOfFreedom, std::set<BoundaryType> >::iterator it = periodicDofAssoc.begin();
	 it != periodicDofAssoc.end(); ++it) {
      TEST_EXIT_DBG((mesh->getDim() == 2 && it->second.size() <= 2) ||
		    (mesh->getDim() == 3 && it->second.size() <= 3))
	("Should not happen!\n");
      
      if (it->second.size() == 2)
	multPeriodicDof2.push_back(it->first);
      if (it->second.size() == 3)
	multPeriodicDof3.push_back(it->first);
    }

    if (mesh->getDim() == 2) {
      TEST_EXIT_DBG(multPeriodicDof2.size() == 0 || 
		    multPeriodicDof2.size() == 4)
	("Should not happen (%d)!\n", multPeriodicDof2.size());
      TEST_EXIT_DBG(multPeriodicDof3.size() == 0)("Should not happen!\n");
    }
267

268
269
270
271
272
    if (mesh->getDim() == 3) {
      TEST_EXIT_DBG(multPeriodicDof3.size() == 0 || 
		    multPeriodicDof3.size() == 8)
	("Should not happen (%d)!\n", multPeriodicDof3.size());
    }
273

274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
    if (multPeriodicDof2.size() > 0) {
      for (unsigned int i = 0; i < multPeriodicDof2.size(); i++) {
	DegreeOfFreedom dof0 = multPeriodicDof2[i];
	if (dof0 == -1)
	  continue;
	
	DegreeOfFreedom dof1 = -1;
	DegreeOfFreedom dof2 = -1;
	BoundaryType type0 = *(periodicDofAssoc[dof0].begin());
	BoundaryType type1 = *(++(periodicDofAssoc[dof0].begin()));
	
	for (PerBoundMap<DegreeOfFreedom>::iterator it = periodicVertices.begin();
	     it != periodicVertices.end(); ++it) {
	  if (it->first.first == dof0 && it->second == type0)
	    dof1 = it->first.second;
	  if (it->first.first == dof0 && it->second == type1)
	    dof2 = it->first.second;
	  
	  if (dof1 != -1 && dof2 != -1)
	    break;
	}
	
	TEST_EXIT_DBG(dof1 != -1 && dof2 != -1)("Should not happen!\n");
	
	DegreeOfFreedom dof3 = -1;
	for (PerBoundMap<DegreeOfFreedom>::iterator it = periodicVertices.begin();
	     it != periodicVertices.end(); ++it) {
	  if (it->first.first == dof1 && it->second == type1) {
	    dof3 = it->first.second;
	    
	    TEST_EXIT_DBG(periodicVertices[make_pair(dof2, dof3)] == type0)
	      ("Should not happen!\n");
	    
	    break;
	  }
	}
	  
	TEST_EXIT_DBG(dof3 != -1)("Should not happen!\n");
	TEST_EXIT_DBG(periodicVertices.count(make_pair(dof0, dof3)) == 0)
	  ("Should not happen!\n");
	TEST_EXIT_DBG(periodicVertices.count(make_pair(dof3, dof0)) == 0)
	  ("Should not happen!\n");
316
317

 	periodicVertices[make_pair(dof0, dof3)] = 
318
	  provideConnectedPeriodicBoundary(type0, type1);
319
 	periodicVertices[make_pair(dof3, dof0)] = 
320
 	  provideConnectedPeriodicBoundary(type0, type1);
321

322
323
324
325
326
327
328
329
	for (unsigned int j = i + 1; j < multPeriodicDof2.size(); j++)
	  if (multPeriodicDof2[j] == dof3)
	    multPeriodicDof2[j] = -1;
      }
    }

    if (multPeriodicDof3.size() > 0) {
      int nMultPeriodicDofs = multPeriodicDof3.size();
330

331
332
333
334
335
336
337
338
      for (int i = 0; i < nMultPeriodicDofs; i++) {
	for (int j = i + 1; j < nMultPeriodicDofs; j++) {
	  pair<DegreeOfFreedom, DegreeOfFreedom> perDofs0 = 
	    make_pair(multPeriodicDof3[i], multPeriodicDof3[j]);
	  pair<DegreeOfFreedom, DegreeOfFreedom> perDofs1 = 
	    make_pair(multPeriodicDof3[j], multPeriodicDof3[i]);
	  
	  if (periodicVertices.count(perDofs0) == 0) {
339
	    BoundaryType b = getNewBoundaryType();
340
341
	    periodicVertices[perDofs0] = b;
	    periodicVertices[perDofs1] = b;
342
343
344
345
	  }
	}
      }
    }
346
347
  
   
348
349
350
351
352
353
354
355
    // === Get all edges that have multiple periodic associations (3D only!). ===
    
    for (map<DofEdge, std::set<DofEdge> >::iterator it = periodicEdgeAssoc.begin();
	 it != periodicEdgeAssoc.end(); ++it) {
      if (it->second.size() > 1) {
	TEST_EXIT_DBG(mesh->getDim() == 3)("Should not happen!\n");
	TEST_EXIT_DBG(it->second.size() == 2)("Should not happen!\n");
	
356
357
358
	DofEdge edge0 = it->first;
	DofEdge edge1 = *(it->second.begin());
	DofEdge edge2 = *(++(it->second.begin()));
359
	
360
361
362
363
364
365
366
367
368
369
	pair<DofEdge, DofEdge> perEdge0 = make_pair(edge1, edge2);
	pair<DofEdge, DofEdge> perEdge1 = make_pair(edge2, edge1);

	TEST_EXIT_DBG(periodicEdges.count(make_pair(edge0, edge1)) == 1)
	  ("Should not happen!\n");
	TEST_EXIT_DBG(periodicEdges.count(make_pair(edge1, edge0)) == 1)
	  ("Should not happen!\n");

	BoundaryType type0 = periodicEdges[make_pair(edge0, edge1)];
	BoundaryType type1 = periodicEdges[make_pair(edge0, edge2)];
370
	BoundaryType type2 = provideConnectedPeriodicBoundary(type0, type1);
371
372
373

 	periodicEdges[perEdge0] = type2;
 	periodicEdges[perEdge1] = type2;
374
375
376
377
378
379
380
381
382
383
384
385
386
387
      }
    }


    // === In debug mode we make some tests, if the periodic structures are set ===
    // === in a symmetric way, i.e., if A -> B for a specific boundary type,    ===
    // === there must be a mapping B -> A with the same boundary type.          ===
    
#if (DEBUG != 0)
    for (PerBoundMap<DegreeOfFreedom>::iterator it = periodicVertices.begin();
	 it != periodicVertices.end(); ++it) {
      pair<DegreeOfFreedom, DegreeOfFreedom> testVertex = 
	make_pair(it->first.second, it->first.first);
      
388
389
390
391
      TEST_EXIT_DBG(periodicVertices.count(testVertex) == 1)
	("Should not happen!\n");
      TEST_EXIT_DBG(periodicVertices[testVertex] == it->second)
	("Should not happen!\n");
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
    }

    for (PerBoundMap<DofEdge>::iterator it = periodicEdges.begin();
	 it != periodicEdges.end(); ++it) {
      pair<DofEdge, DofEdge> testEdge = 
	make_pair(it->first.second, it->first.first);
      
      TEST_EXIT_DBG(periodicEdges.count(testEdge) == 1)("Should not happen!\n");
      TEST_EXIT_DBG(periodicEdges[testEdge] == it->second)("Should not happen!\n");
    }
    
    for (PerBoundMap<DofFace>::iterator it = periodicFaces.begin();
	 it != periodicFaces.end(); ++it) {
      pair<DofFace, DofFace> testFace = 
	make_pair(it->first.second, it->first.first);
      
      TEST_EXIT_DBG(periodicFaces.count(testFace) == 1)("Should not happen!\n");
      TEST_EXIT_DBG(periodicFaces[testFace] == it->second)("Should not happen!\n");
    }
#endif       
  }


Thomas Witkowski's avatar
Thomas Witkowski committed
415
416
417
418
void ElementObjectDatabase::removePeriodicData()
{
}

419
  BoundaryType ElementObjectDatabase::getNewBoundaryType()
420
  {
421
    FUNCNAME("ElementObjectDatabase::getNewBoundaryType()");
422
423
424
425
426
427
428
429
430
431

    BoundaryType newPeriodicBoundaryType = 0;
    for (map<BoundaryType, VertexVector*>::iterator it = mesh->getPeriodicAssociations().begin();
	 it != mesh->getPeriodicAssociations().end(); ++it)
      newPeriodicBoundaryType = std::min(newPeriodicBoundaryType, it->first);
    
    TEST_EXIT_DBG(newPeriodicBoundaryType < 0)("Should not happen!\n");
    newPeriodicBoundaryType--;
    
    mesh->getPeriodicAssociations()[newPeriodicBoundaryType] = 
432
      new VertexVector(feSpace->getAdmin(), "");
433
434
435
436
437
    
    return newPeriodicBoundaryType;
  }


438
  BoundaryType 
439
  ElementObjectDatabase::provideConnectedPeriodicBoundary(BoundaryType b0, 
440
							  BoundaryType b1)
441
  {
442
    FUNCNAME("ElementObjectDatabase::provideConnectedPeriodicBoundary()");
443
444

    std::pair<BoundaryType, BoundaryType> bConn = 
445
      (b0 <= b1 ? make_pair(b0, b1) : make_pair(b1, b0));
446
447

    if (bConnMap.count(bConn) == 0) {
448
      BoundaryType newPeriodicBoundaryType = getNewBoundaryType();
449
450
451
452
453

      VertexVector &vecB0 = mesh->getPeriodicAssociations(b0);
      VertexVector &vecB1 = mesh->getPeriodicAssociations(b1);
      VertexVector &vecC = mesh->getPeriodicAssociations(newPeriodicBoundaryType);

454
      DOFIteratorBase it(const_cast<DOFAdmin*>(feSpace->getAdmin()), USED_DOFS);
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472

      for (it.reset(); !it.end(); ++it) {
	if (!it.isDofFree()) {
	  TEST_EXIT_DBG(vecB1[vecB0[it.getDOFIndex()]] == vecB0[vecB1[it.getDOFIndex()]])
	    ("Should not happen!\n");

	  vecC[it.getDOFIndex()] = vecB1[vecB0[it.getDOFIndex()]];
	}
      }
      
      bConnMap[bConn] = newPeriodicBoundaryType;
    }


    return bConnMap[bConn];
  }


473
  void ElementObjectDatabase::updateRankData()
Thomas Witkowski's avatar
Thomas Witkowski committed
474
  {
475
    FUNCNAME("ElementObjectDatabase::updateRankData()");
476
    
Thomas Witkowski's avatar
Thomas Witkowski committed
477
    vertexInRank.clear();
478
    for (map<DegreeOfFreedom, vector<ElementObjectData> >::iterator it = vertexElements.begin();
Thomas Witkowski's avatar
Thomas Witkowski committed
479
	 it != vertexElements.end(); ++it) {
480
      for (vector<ElementObjectData>::iterator it2 = it->second.begin(); 
Thomas Witkowski's avatar
Thomas Witkowski committed
481
	   it2 != it->second.end(); ++it2) {
482
	int elementInRank =(* macroElementRankMap)[it2->elIndex];	
Thomas Witkowski's avatar
Thomas Witkowski committed
483
484
	if (it2->elIndex > vertexInRank[it->first][elementInRank].elIndex)
	  vertexInRank[it->first][elementInRank] = *it2;
Thomas Witkowski's avatar
Thomas Witkowski committed
485
486
487
      }
    }
    
Thomas Witkowski's avatar
Thomas Witkowski committed
488
    edgeInRank.clear();
489
    for (map<DofEdge, vector<ElementObjectData> >::iterator it = edgeElements.begin();
Thomas Witkowski's avatar
Thomas Witkowski committed
490
	 it != edgeElements.end(); ++it) {
491
      for (vector<ElementObjectData>::iterator it2 = it->second.begin(); 
Thomas Witkowski's avatar
Thomas Witkowski committed
492
	   it2 != it->second.end(); ++it2) {
493
	int elementInRank = (*macroElementRankMap)[it2->elIndex];
Thomas Witkowski's avatar
Thomas Witkowski committed
494
495
	if (it2->elIndex > edgeInRank[it->first][elementInRank].elIndex)
	  edgeInRank[it->first][elementInRank] = *it2;
Thomas Witkowski's avatar
Thomas Witkowski committed
496
497
498
      }
    }
    
Thomas Witkowski's avatar
Thomas Witkowski committed
499
    faceInRank.clear();
500
    for (map<DofFace, vector<ElementObjectData> >::iterator it = faceElements.begin();
Thomas Witkowski's avatar
Thomas Witkowski committed
501
	 it != faceElements.end(); ++it) {
502
      for (vector<ElementObjectData>::iterator it2 = it->second.begin(); 
Thomas Witkowski's avatar
Thomas Witkowski committed
503
	   it2 != it->second.end(); ++it2) {
504
	int elementInRank = (*macroElementRankMap)[it2->elIndex];
Thomas Witkowski's avatar
Thomas Witkowski committed
505
506
	if (it2->elIndex > faceInRank[it->first][elementInRank].elIndex)
	  faceInRank[it->first][elementInRank] = *it2;
Thomas Witkowski's avatar
Thomas Witkowski committed
507
508
509
510
      }
    }
  }

511

512
  void ElementObjectDatabase::createReverseModeData()
513
  {
514
    FUNCNAME("ElementObjectDatabase::createReverseModeData()");
515

516
517
    // === In 2D, all reverse modes are always true! ===

518
519
520
    if (mesh->getDim() == 2)
      return;

521
522
523

    // === First, create reverse modes for all "directly" neighbouring elements. ===

524
525
526
527
528
529
    for (map<DofEdge, vector<ElementObjectData> >::iterator edgeIt = edgeElements.begin();
	 edgeIt != edgeElements.end(); ++edgeIt) {

      vector<ElementObjectData>& els = edgeIt->second;

      for (unsigned int i = 0; i < els.size(); i++) {
530
531
	BoundaryObject obj0(macroElIndexMap[els[i].elIndex], 
			    macroElIndexTypeMap[els[i].elIndex], 
532
533
534
			    EDGE, els[i].ithObject);

	for (unsigned int j = i + 1; j < els.size(); j++) {
535
536
	  BoundaryObject obj1(macroElIndexMap[els[j].elIndex], 
			      macroElIndexTypeMap[els[j].elIndex], 
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
			      EDGE, els[j].ithObject);

	  bool reverseMode = 
	    BoundaryObject::computeReverseMode(obj0, obj1, feSpace, INTERIOR);

	  edgeReverseMode[make_pair(els[i], els[j])] = reverseMode;
	  edgeReverseMode[make_pair(els[j], els[i])] = reverseMode;
	}
      }
    }

    for (map<DofFace, vector<ElementObjectData> >::iterator faceIt = faceElements.begin();
	 faceIt != faceElements.end(); ++faceIt) {

      vector<ElementObjectData>& els = faceIt->second;

      for (unsigned int i = 0; i < els.size(); i++) {
554
555
	BoundaryObject obj0(macroElIndexMap[els[i].elIndex], 
			    macroElIndexTypeMap[els[i].elIndex], 
556
557
558
			    FACE, els[i].ithObject);

	for (unsigned int j = i + 1; j < els.size(); j++) {
559
560
	  BoundaryObject obj1(macroElIndexMap[els[j].elIndex], 
			      macroElIndexTypeMap[els[j].elIndex], 
561
562
563
564
565
566
567
568
569
570
			      FACE, els[j].ithObject);

	  bool reverseMode = 
	    BoundaryObject::computeReverseMode(obj0, obj1, feSpace, INTERIOR);

	  faceReverseMode[make_pair(els[i], els[j])] = reverseMode;
	  faceReverseMode[make_pair(els[j], els[i])] = reverseMode;
	}
      }
    }
571
572
573
574
575
576
577
578
579
580


    // === And create reverse modes for periodic neighbouring elements. ===

    for (PerBoundMap<DofEdge>::iterator edgeIt = periodicEdges.begin();
	 edgeIt != periodicEdges.end(); ++edgeIt) {
      vector<ElementObjectData> &edges0 = edgeElements[edgeIt->first.first];
      vector<ElementObjectData> &edges1 = edgeElements[edgeIt->first.second];

      for (unsigned int i = 0; i < edges0.size(); i++) {
581
582
	BoundaryObject obj0(macroElIndexMap[edges0[i].elIndex], 
			    macroElIndexTypeMap[edges0[i].elIndex], 
583
584
585
			    EDGE, edges0[i].ithObject);

	for (unsigned int j = 0; j < edges1.size(); j++) {
586
587
	  BoundaryObject obj1(macroElIndexMap[edges1[j].elIndex], 
			      macroElIndexTypeMap[edges1[j].elIndex], 
588
589
590
			      EDGE, edges1[j].ithObject);

	  bool reverseMode = 
591
592
	    BoundaryObject::computeReverseMode(obj0, obj1, feSpace, 
					       edgeIt->second);
593
594
595
596
597
598
599
600
601
602
603
604
605
606

	  edgeReverseMode[make_pair(edges0[i], edges1[j])] = reverseMode;
	  edgeReverseMode[make_pair(edges1[j], edges0[i])] = reverseMode;
	}
      }
    }

    for (PerBoundMap<DofFace>::iterator faceIt = periodicFaces.begin();
	 faceIt != periodicFaces.end(); ++faceIt) {
      vector<ElementObjectData> &faces0 = faceElements[faceIt->first.first];
      vector<ElementObjectData> &faces1 = faceElements[faceIt->first.second];

      TEST_EXIT_DBG(faces0.size() == faces1.size() == 1)("Should not happen!\n");

607
608
      BoundaryObject obj0(macroElIndexMap[faces0[0].elIndex], 
			  macroElIndexTypeMap[faces0[0].elIndex], 
609
			  FACE, faces0[0].ithObject);      
610
611
      BoundaryObject obj1(macroElIndexMap[faces1[0].elIndex], 
			  macroElIndexTypeMap[faces1[0].elIndex], 
612
613
614
615
616
617
618
619
			  FACE, faces1[0].ithObject);
      
      bool reverseMode = 
	BoundaryObject::computeReverseMode(obj0, obj1, feSpace, faceIt->second);
      
      faceReverseMode[make_pair(faces0[0], faces1[0])] = reverseMode;
      faceReverseMode[make_pair(faces1[0], faces0[0])] = reverseMode;
    }  
620
621
622
  }


623
  int ElementObjectDatabase::getIterateOwner(int level)
624
625
626
627
628
629
  {
    FUNCNAME("ElementObjectDatabase::getIterateOwner()");

    TEST_EXIT_DBG(macroElementRankMap)("Should not happen!\n");

    int owner = -1;
630
631
    vector<ElementObjectData> *objData;
    
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
    switch (iterGeoPos) {
    case VERTEX:
      objData = &(vertexElements[vertexIter->first]);
      break;
    case EDGE:
      objData = &(edgeElements[edgeIter->first]);
      break;
    case FACE:
      objData = &(faceElements[faceIter->first]);
      break;
    }
    
    std::set<int> &levelRanks = levelData->getLevelRanks(level);
    bool allRanks = (levelRanks.size() == 1 && *(levelRanks.begin()) == -1);
    
    for (vector<ElementObjectData>::iterator it = objData->begin();
	 it != objData->end(); ++it) {
      int elRank = (*macroElementRankMap)[it->elIndex];
      if (allRanks || levelRanks.count(elRank))
	owner = std::max(owner, elRank);
    }
    
    TEST_EXIT_DBG(owner >= 0)("Cannot find owner on level %d\n", level);
    
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
    return owner;
  }


  int ElementObjectDatabase::getIterateMaxLevel()
  {
    FUNCNAME("ElementObjectDatabase::getIterateMaxLevel()");

    int nLevel = levelData->getLevelNumber();
    TEST_EXIT_DBG(nLevel > 0)("Should not happen!\n");

    vector<std::set<int> > ranksInLevel;
    ranksInLevel.resize(nLevel);

    switch (iterGeoPos) {
    case VERTEX:
      {
	vector<ElementObjectData>& vertexData = vertexElements[vertexIter->first];
	for (vector<ElementObjectData>::iterator it = vertexData.begin();
	     it != vertexData.end(); ++it)
	  ranksInLevel[0].insert((*macroElementRankMap)[it->elIndex]);
      }
      break;
    case EDGE:
      {
	vector<ElementObjectData>& edgeData = edgeElements[edgeIter->first];
	for (vector<ElementObjectData>::iterator it = edgeData.begin();
	     it != edgeData.end(); ++it)
	  ranksInLevel[0].insert((*macroElementRankMap)[it->elIndex]);
      }
      break;
    case FACE:
688
689
690
691
692
693
      {
	vector<ElementObjectData>& faceData = faceElements[faceIter->first];
	for (vector<ElementObjectData>::iterator it = faceData.begin();
	     it != faceData.end(); ++it)
	  ranksInLevel[0].insert((*macroElementRankMap)[it->elIndex]);
      }
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
      break;
    default:
      ERROR_EXIT("Should not happen!\n");
    }

    for (std::set<int>::iterator it = ranksInLevel[0].begin(); 
	 it != ranksInLevel[0].end(); ++it) 
      for (int level = 1; level < nLevel; level++)
	ranksInLevel[level].insert(levelData->getLevelId(level, *it));

    int maxLevel = 0;
    for (int level = 1; level < nLevel; level++)
      if (ranksInLevel[level].size() > 1)
	maxLevel = level;
    
    return maxLevel;
  }


713
  void ElementObjectDatabase::serialize(ostream &out)
714
  {
715
    FUNCNAME("ElementObjectDatabase::serialize()");
716
717
718

    int nSize = vertexElements.size();
    SerUtil::serialize(out, nSize);
719
    for (map<DegreeOfFreedom, vector<ElementObjectData> >::iterator it = vertexElements.begin();
720
721
722
723
724
725
726
	 it != vertexElements.end(); ++it) {
      SerUtil::serialize(out, it->first);
      serialize(out, it->second);
    }

    nSize = edgeElements.size();
    SerUtil::serialize(out, nSize);
727
    for (map<DofEdge, vector<ElementObjectData> >::iterator it = edgeElements.begin();
728
729
730
731
732
733
734
	 it != edgeElements.end(); ++it) {
      SerUtil::serialize(out, it->first);
      serialize(out, it->second);
    }

    nSize = faceElements.size();
    SerUtil::serialize(out, nSize);
735
    for (map<DofFace, vector<ElementObjectData> >::iterator it = faceElements.begin();
736
737
738
739
740
741
	 it != faceElements.end(); ++it) {
      SerUtil::serialize(out, it->first);
      serialize(out, it->second);
    }


742

743
744
    nSize = vertexLocalMap.size();
    SerUtil::serialize(out, nSize);
745
    for (map<ElementObjectData, DegreeOfFreedom>::iterator it = vertexLocalMap.begin();
746
747
748
749
750
751
752
	 it != vertexLocalMap.end(); ++it) {
      it->first.serialize(out);
      SerUtil::serialize(out, it->second);
    }

    nSize = edgeLocalMap.size();
    SerUtil::serialize(out, nSize);
753
    for (map<ElementObjectData, DofEdge>::iterator it = edgeLocalMap.begin();
754
755
756
757
758
759
760
	 it != edgeLocalMap.end(); ++it) {
      it->first.serialize(out);
      SerUtil::serialize(out, it->second);
    }

    nSize = faceLocalMap.size();
    SerUtil::serialize(out, nSize);
761
    for (map<ElementObjectData, DofFace>::iterator it = faceLocalMap.begin();
762
763
764
765
766
767
	 it != faceLocalMap.end(); ++it) {
      it->first.serialize(out);
      SerUtil::serialize(out, it->second);
    }


768
    nSize = vertexInRank.size();
769
    SerUtil::serialize(out, nSize);
770
    for (map<DegreeOfFreedom, map<int, ElementObjectData> >::iterator it = vertexInRank.begin();
771
772
773
774
775
776
	 it != vertexInRank.end(); ++it) {
      SerUtil::serialize(out, it->first);
      serialize(out, it->second);
    }

    nSize = edgeInRank.size();
777
    SerUtil::serialize(out, nSize);
778
    for (map<DofEdge, map<int, ElementObjectData> >::iterator it = edgeInRank.begin();
779
780
781
782
783
784
	 it != edgeInRank.end(); ++it) {
      SerUtil::serialize(out, it->first);
      serialize(out, it->second);
    }

    nSize = faceInRank.size();
785
    SerUtil::serialize(out, nSize);
786
    for (map<DofFace, map<int, ElementObjectData> >::iterator it = faceInRank.begin();
787
788
789
790
791
	 it != faceInRank.end(); ++it) {
      SerUtil::serialize(out, it->first);
      serialize(out, it->second);
    }

792
793
794
    SerUtil::serialize(out, periodicVertices);
    SerUtil::serialize(out, periodicEdges);
    SerUtil::serialize(out, periodicFaces);
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830


    nSize = periodicDofAssoc.size();
    SerUtil::serialize(out, nSize);
    for (map<DegreeOfFreedom, std::set<BoundaryType> >::iterator it = periodicDofAssoc.begin();
	 it != periodicDofAssoc.end(); ++it) {
      SerUtil::serialize(out, it->first);
      SerUtil::serialize(out, it->second);
    }

    nSize = periodicEdgeAssoc.size();
    SerUtil::serialize(out, nSize);
    for (map<DofEdge, std::set<DofEdge> >::iterator it = periodicEdgeAssoc.begin();
	 it != periodicEdgeAssoc.end(); ++it) {
      SerUtil::serialize(out, it->first);
      SerUtil::serialize(out, it->second);
    }


    nSize = edgeReverseMode.size();
    SerUtil::serialize(out, nSize);
    for (map<pair<ElementObjectData, ElementObjectData>, bool>::iterator it = edgeReverseMode.begin();
	 it != edgeReverseMode.end(); ++it) {
      it->first.first.serialize(out);
      it->first.second.serialize(out);
      SerUtil::serialize(out, it->second);
    }

    nSize = faceReverseMode.size();
    SerUtil::serialize(out, nSize);
    for (map<pair<ElementObjectData, ElementObjectData>, bool>::iterator it = faceReverseMode.begin();
	 it != faceReverseMode.end(); ++it) {
      it->first.first.serialize(out);
      it->first.second.serialize(out);
      SerUtil::serialize(out, it->second);
    }
831
832
833
  }


834
  void ElementObjectDatabase::deserialize(istream &in)
835
  {
836
    FUNCNAME("ElementObjectDatabase::deserialize()");
837
838
839
840
841
842

    int nSize;
    SerUtil::deserialize(in, nSize);
    vertexElements.clear();
    for (int i = 0; i < nSize; i++) {
      DegreeOfFreedom dof;
843
      vector<ElementObjectData> data;
844
845
846
847
848
849
850
851
852
      SerUtil::deserialize(in, dof);
      deserialize(in, data);
      vertexElements[dof] = data;
    }

    SerUtil::deserialize(in, nSize);
    edgeElements.clear();
    for (int i = 0; i < nSize; i++) {
      DofEdge edge;
853
      vector<ElementObjectData> data;
854
855
856
857
858
859
860
861
862
      SerUtil::deserialize(in, edge);
      deserialize(in, data);
      edgeElements[edge] = data;
    }

    SerUtil::deserialize(in, nSize);
    faceElements.clear();
    for (int i = 0; i < nSize; i++) {
      DofFace face;
863
      vector<ElementObjectData> data;
864
865
866
867
      SerUtil::deserialize(in, face);
      deserialize(in, data);
      faceElements[face] = data;
    }
868
   
869

870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900

    SerUtil::deserialize(in, nSize);
    vertexLocalMap.clear();
    for (int i = 0; i < nSize; i++) {
      ElementObjectData data;
      DegreeOfFreedom dof;
      data.deserialize(in);
      SerUtil::deserialize(in, dof);
      vertexLocalMap[data] = dof;
    }

    SerUtil::deserialize(in, nSize);
    edgeLocalMap.clear();
    for (int i = 0; i < nSize; i++) {
      ElementObjectData data;
      DofEdge edge;
      data.deserialize(in);
      SerUtil::deserialize(in, edge);
      edgeLocalMap[data] = edge;
    }

    SerUtil::deserialize(in, nSize);
    faceLocalMap.clear();
    for (int i = 0; i < nSize; i++) {
      ElementObjectData data;
      DofFace face;
      data.deserialize(in);
      SerUtil::deserialize(in, face);
      faceLocalMap[data] = face;
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
901

902
903
904
905
906

    SerUtil::deserialize(in, nSize);
    vertexInRank.clear();
    for (int i = 0; i < nSize; i++) {
      DegreeOfFreedom dof;
907
      map<int, ElementObjectData> data;
908
909
910
911
912
913
914
915
916
      SerUtil::deserialize(in, dof);
      deserialize(in, data);
      vertexInRank[dof] = data;
    }

    SerUtil::deserialize(in, nSize);
    edgeInRank.clear();
    for (int i = 0; i < nSize; i++) {
      DofEdge edge;
917
      map<int, ElementObjectData> data;
918
919
920
921
922
923
924
925
926
      SerUtil::deserialize(in, edge);
      deserialize(in, data);
      edgeInRank[edge] = data;
    }

    SerUtil::deserialize(in, nSize);
    faceInRank.clear();
    for (int i = 0; i < nSize; i++) {
      DofFace face;
927
      map<int, ElementObjectData> data;
928
929
930
931
      SerUtil::deserialize(in, face);
      deserialize(in, data);
      faceInRank[face] = data;
    }
932
933
934

    SerUtil::deserialize(in, periodicVertices);
    SerUtil::deserialize(in, periodicEdges);
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
    SerUtil::deserialize(in, periodicFaces);

    

    SerUtil::deserialize(in, nSize);
    periodicDofAssoc.clear();
    for (int i = 0; i < nSize; i++) {
      DegreeOfFreedom dof;
      std::set<DegreeOfFreedom> dofs;
      SerUtil::deserialize(in, dof);
      SerUtil::deserialize(in, dofs);
      periodicDofAssoc[dof] = dofs;
    }

    SerUtil::deserialize(in, nSize);
    periodicEdgeAssoc.clear();
    for (int i = 0; i < nSize; i++) {
      DofEdge edge;
      std::set<DofEdge> edges;
      SerUtil::deserialize(in, edge);
      SerUtil::deserialize(in, edges);
      periodicEdgeAssoc[edge] = edges;
    }



    SerUtil::deserialize(in, nSize);
    edgeReverseMode.clear();
    for (int i = 0; i < nSize; i++) {
      ElementObjectData obj0, obj1;
      bool reverseMode;
      obj0.deserialize(in);
      obj1.deserialize(in);
      SerUtil::deserialize(in, reverseMode);

      edgeReverseMode[make_pair(obj0, obj1)] = reverseMode;
    }

    SerUtil::deserialize(in, nSize);
    faceReverseMode.clear();
    for (int i = 0; i < nSize; i++) {
      ElementObjectData obj0, obj1;
      bool reverseMode;
      obj0.deserialize(in);
      obj1.deserialize(in);
      SerUtil::deserialize(in, reverseMode);

      faceReverseMode[make_pair(obj0, obj1)] = reverseMode;
    }
984
985
986
  }


987
988
  void ElementObjectDatabase::serialize(ostream &out, 
					vector<ElementObjectData>& elVec)
989
990
991
992
993
994
995
996
  {
    int nSize = elVec.size();
    SerUtil::serialize(out, nSize);
    for (int i = 0; i < nSize; i++)
      elVec[i].serialize(out);
  }


997
998
  void ElementObjectDatabase::deserialize(istream &in, 
					  vector<ElementObjectData>& elVec)
999
1000
1001
1002
1003
1004
1005
1006
1007
  {
    int nSize;
    SerUtil::deserialize(in, nSize);
    elVec.resize(nSize);
    for (int i = 0; i < nSize; i++)
      elVec[i].deserialize(in);
  }

  
1008
1009
  void ElementObjectDatabase::serialize(ostream &out, 
					map<int, ElementObjectData>& data)
1010
1011
1012
  {
    int nSize = data.size();
    SerUtil::serialize(out, nSize);
1013
    for (map<int, ElementObjectData>::iterator it = data.begin();
1014
1015
1016
1017
1018
1019
1020
	 it != data.end(); ++it) {
      SerUtil::serialize(out, it->first);
      it->second.serialize(out);
    }
  }

  
1021
1022
  void ElementObjectDatabase::deserialize(istream &in, 
					  map<int, ElementObjectData>& data)
1023
1024
1025
1026
1027
1028
1029
1030
1031
1032
1033
1034
1035
  {
    int nSize;
    SerUtil::deserialize(in, nSize);
    for (int i = 0; i < nSize; i++) {
      int index;
      ElementObjectData elObj;
      SerUtil::deserialize(in, index);
      elObj.deserialize(in);

      data[index] = elObj;
    }
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
1036
}