Element.cc 15.1 KB
Newer Older
1
2
3
4
5
6
#include "Element.h"
#include "DOFAdmin.h"
#include "Mesh.h"
#include "CoarseningManager.h"
#include "FixVec.h"
#include "ElementRegion_ED.h"
7
#include "Serializer.h"
8
#include "MeshStructure.h"
9
10
11

namespace AMDiS {

12
  std::map<DegreeOfFreedom*, bool> Element::deletedDOFs;
13

14

15
16
  int Element::getRegion() const 
  {
17
18
19
    if (!elementData) 
      return -1;

20
21
    ElementRegion_ED* red = 
      dynamic_cast<ElementRegion_ED*>(elementData->getElementData(ELEMENT_REGION));
22

23
24
    if (red)
      return red->getRegion();    
25
26
27
28
    
    return -1;
  }

29

30
31
  void Element::setDOFPtrs()
  {
32
    FUNCNAME("Element::setDOFPtrs()");
33
34

    TEST_EXIT_DBG(mesh)("no mesh!\n");
35

36
37
38
    dof = mesh->createDOFPtrs();  
  }

39

40
41
42
43
  Element::Element(Mesh *aMesh)
  {
    mesh = aMesh;
    index = mesh ? mesh->getNextElementIndex() : -1; 
44
45
    child[0] = NULL;
    child[1] = NULL;
46
47
    newCoord = NULL;
    elementData = NULL;
48
    mark = 0;
49
50
51
52
53
54

    if (mesh) {
      setDOFPtrs();
    } else {
      mesh = NULL;
    }
55
56
  }

57

58
59
  // call destructor through Mesh::freeElement !!!
  Element::~Element()
60
  {  
61
    if (child[0])
Thomas Witkowski's avatar
Thomas Witkowski committed
62
      delete child[0];
63
64
    
    if (child[1])
Thomas Witkowski's avatar
Thomas Witkowski committed
65
      delete child[1];   
66

67
    if (newCoord)
Thomas Witkowski's avatar
Thomas Witkowski committed
68
      delete newCoord;    
69
70
71

    if (elementData) {
      elementData->deleteDecorated();
Thomas Witkowski's avatar
Thomas Witkowski committed
72
      delete elementData;
73
74
75
    }
  }

76

77
78
79
80
81
82
83
84
  bool Element::deleteElementData(int typeID)
  {
    FUNCNAME("Element::deleteElementData()");

    if (elementData) {
      if (elementData->isOfType(typeID)) {
	ElementData *tmp = elementData;
	elementData = elementData->getDecorated();
Thomas Witkowski's avatar
Thomas Witkowski committed
85
	delete tmp;
86
87
88
89
90
91
92
	tmp = NULL;
	return true;
      } else {
	return elementData->deleteDecorated(typeID);
      }
    }
    return false;
93
94
  }

95

96
97
98
99
100
101
102
103
104
  void Element::deleteElementDOFs()
  {
    int dim = mesh->getDim();
    int j = 0;

    for (int pos = 0; pos <= dim; pos++) {
      GeoIndex position = INDEX_OF_DIM(pos, dim);
      int ndof = 0;
     
Thomas Witkowski's avatar
Thomas Witkowski committed
105
      for (int i = 0; i < mesh->getNumberOfDOFAdmin(); i++)
106
107
108
109
	ndof += mesh->getDOFAdmin(i).getNumberOfDOFs(position);

      if (ndof > 0) {
	for (int i = 0; i < mesh->getGeo(position); i++) {
110
111
112
	  if (dof[j]) {
	    if (deletedDOFs.count(dof[j]) == 0) {
	      deletedDOFs[dof[j]] = true;
113
	      delete [] dof[j];
114
	    }
115
	  }  
116
117
118
119
120
	  j++;
	}
      }
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
121
    delete [] dof;
122
    
123
    if (child[0])
124
      child[0]->deleteElementDOFs();
125
    if (child[1])
126
127
128
      child[1]->deleteElementDOFs();
  }

129
130
131
  Element* Element::cloneWithDOFs()
  {
    Element *el;
132
    
133
    if (isLine()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
134
      el = new Line(NULL);
135
    } else if (isTriangle()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
136
      el = new Triangle(NULL);
137
    } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
138
      el = new Tetrahedron(NULL);
139
    }
140
141

    el->mesh = mesh;
142
143
144
    el->index = index;
    el->mark = mark;
    if (newCoord) {
Thomas Witkowski's avatar
Thomas Witkowski committed
145
      WorldVector<double> *nc = new WorldVector<double>();
146
147
148
      *nc = *newCoord;
      el->newCoord = nc;
    }
149
    
150
    /* =========== And here we clone the DOFs =========== */
151
   
152
    el->dof = new DegreeOfFreedom*[mesh->getNumberOfNodes()];
153
154
155
156
157
158
159
160

    int dim = mesh->getDim();
    int j = 0;

    for (int pos = 0; pos <= dim; pos++) {
      GeoIndex position = INDEX_OF_DIM(pos, dim);
      int ndof = 0;

161
      for (int i = 0; i < mesh->getNumberOfDOFAdmin(); i++)
162
163
164
165
166
	ndof += mesh->getDOFAdmin(i).getNumberOfDOFs(position);

      if (ndof > 0) {
	for (int i = 0; i < mesh->getGeo(position); i++) {
	  if (dof[j] != NULL) {
167
168
169
	    std::pair<DegreeOfFreedom, int> idx = std::make_pair(dof[j][0], pos);

	    if (Mesh::serializedDOFs[idx] == NULL) {
170
171
	      el->dof[j] = new DegreeOfFreedom[ndof];
	      for (int k = 0; k < ndof; k++)
172
		el->dof[j][k] = dof[j][k];
173

174
	      Mesh::serializedDOFs[idx] = el->dof[j];
175
	    } else {
176
	      el->dof[j] = Mesh::serializedDOFs[idx];
177
178
179
180
181
182
183
184
	    }
	  } else {
	    el->dof[j] = NULL;
	  }
	  j++;
	}
      }
    }
185
    
186
187
    /* =========== And clone the children ============= */

188
    if (child[0]) 
189
      el->child[0] = child[0]->cloneWithDOFs();
190
    if (child[1])
191
      el->child[1] = child[1]->cloneWithDOFs();
192
193
194
195

    return el;
  }

196
197
198
199
200
201
202
203
204
  /****************************************************************************/
  /*  ATTENTION:                                                              */
  /*  new_dof_fct() destroys new_dof !!!!!!!!!!                               */
  /*  should be used only at the end of dof_compress()!!!!!                   */
  /****************************************************************************/

  /* CHANGE_DOFS_1 changes old dofs to NEGATIVE new dofs */

#define CHANGE_DOFS_1(el)						\
205
  ldof = el->dof[n0 + i] + nd0;						\
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
    for (j = 0; j < nd; j++) {						\
      if ((k = ldof[j]) >= 0) {						\
	/* do it only once! (dofs are visited more than once) */	\
	ldof[j] = - admin->getMesh()->newDOF[k] - 1;			\
      } }

  /* CHANGE_DOFS_2 changes NEGATIVE new dofs to POSITIVE */

#define CHANGE_DOFS_2(el)						\
  ldof = el->dof[n0+i] + nd0;						\
    for (j = 0; j < nd; j++) {						\
      if ((k = ldof[j]) < 0) {						\
	/* do it only once! (dofs are visited more than once) */	\
	ldof[j] = - k - 1;						\
      } }

  void Element::newDOFFct1(const DOFAdmin* admin)
  {
224
225
    int j, k, n0, nd, nd0;
    DegreeOfFreedom *ldof;  
226
227
228
229
230
231
232
    int vertices = mesh->getGeo(VERTEX);
    int edges = mesh->getGeo(EDGE); 
    int faces = mesh->getGeo(FACE);

    if ((nd = admin->getNumberOfDOFs(VERTEX)))  {
      nd0 = admin->getNumberOfPreDOFs(VERTEX);
      n0 = admin->getMesh()->getNode(VERTEX);
233
      for (int i = 0; i < vertices; i++) {
234
235
236
237
	CHANGE_DOFS_1(this);
      }
    }

238
    if (mesh->getDim() > 1) {
239
240
241
      if ((nd = admin->getNumberOfDOFs(EDGE)))  {
	nd0 = admin->getNumberOfPreDOFs(EDGE);
	n0 = admin->getMesh()->getNode(EDGE);
242
	for (int i = 0; i < edges; i++) {
243
244
245
246
247
	  CHANGE_DOFS_1(this);
	}
      }
    }

248
    if (mesh->getDim() == 3) {
249
250
251
      if ((nd = admin->getNumberOfDOFs(FACE)))  {
	nd0 = admin->getNumberOfPreDOFs(FACE);
	n0 = admin->getMesh()->getNode(FACE);
252
	for (int i = 0; i < faces; i++) {
253
254
255
256
257
258
259
260
	  CHANGE_DOFS_1(this);
	}
      }
    }

    if ((nd = admin->getNumberOfDOFs(CENTER)))  {
      nd0 = admin->getNumberOfPreDOFs(CENTER);
      n0 = admin->getMesh()->getNode(CENTER);
261
      int i = 0;          /* only one center */
262
263
264
265
266
267
268
      CHANGE_DOFS_1(this);
    }
  }


  void Element::newDOFFct2(const DOFAdmin* admin)
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
269
    int j, k, n0, nd0;
270
    DegreeOfFreedom  *ldof;
271
272
273
274
    int vertices = mesh->getGeo(VERTEX);
    int edges = mesh->getGeo(EDGE); 
    int faces = mesh->getGeo(FACE);

275
276
    int nd = admin->getNumberOfDOFs(VERTEX);
    if (nd) {
277
278
      nd0 = admin->getNumberOfPreDOFs(VERTEX);
      n0 = admin->getMesh()->getNode(VERTEX);
Thomas Witkowski's avatar
Thomas Witkowski committed
279
      for (int i = 0; i < vertices; i++) {
280
281
282
283
	CHANGE_DOFS_2(this);
      }
    }

284
    if (mesh->getDim() > 1) {
285
286
      nd = admin->getNumberOfDOFs(EDGE);
      if (nd) {
287
288
	nd0 = admin->getNumberOfPreDOFs(EDGE);
	n0 = admin->getMesh()->getNode(EDGE);
Thomas Witkowski's avatar
Thomas Witkowski committed
289
	for (int i = 0; i < edges; i++) {
290
291
292
293
294
	  CHANGE_DOFS_2(this);
	}
      }
    }

295
    if (mesh->getDim() == 3) {
296
297
      nd = admin->getNumberOfDOFs(FACE);
      if (nd) {
298
299
	nd0 = admin->getNumberOfPreDOFs(FACE);
	n0 = admin->getMesh()->getNode(FACE);
Thomas Witkowski's avatar
Thomas Witkowski committed
300
	for (int i = 0; i < faces; i++) {
301
302
303
304
305
	  CHANGE_DOFS_2(this);
	}
      }
    }

306
307
    nd = admin->getNumberOfDOFs(CENTER);
    if (nd) {
308
309
      nd0 = admin->getNumberOfPreDOFs(CENTER);
      n0 = admin->getMesh()->getNode(CENTER);
310
      // only one center
Thomas Witkowski's avatar
Thomas Witkowski committed
311
      int i = 0;   
312
313
314
315
316
317
318
      CHANGE_DOFS_2(this);
    }
  }

#undef CHANGE_DOFS_1
#undef CHANGE_DOFS_2

319

320
321
322
323
324
325
326
  /****************************************************************************/
  /* opp_vertex checks whether the face with vertices dof[0],..,dof[DIM-1] is */
  /* part of mel's boundary. returns the opposite vertex if true, -1 else     */
  /****************************************************************************/

  int Element::oppVertex(FixVec<DegreeOfFreedom*, DIMEN> pdof) const
  {
327
328
    int nv = 0;
    int ov = 0;
329
330
331
    int vertices = mesh->getGeo(VERTEX);
    int dim = mesh->getDim();

332
    for (int i = 0; i < vertices; i++) {
333
      if (nv < i - 1)  
334
335
336
337
338
339
340
341
342
343
344
	return(-1);

      for (int j = 0; j < dim; j++) {
	if (dof[i] == pdof[j]) {
	  /****************************************************************************/
	  /* i is a common vertex                                                     */
	  /****************************************************************************/
	  ov += i;
	  nv++;
	  break;
	}
345
      }
346
347
348
349
350

    }
    
    if (nv != mesh->getDim()) 
      return(-1);
351
352
353
354
355
356
357
358
359
360
    /****************************************************************************/
    /*  the opposite vertex is 3(6) - (sum of indices of common vertices) in    */
    /*  2d(3d)                                                                  */
    /****************************************************************************/

    switch(mesh->getDim()) {
    case 1:
      return ov;
      break;
    case 2:
361
      return 3 - ov;
362
363
      break;
    case 3:
364
      return 6 - ov;
365
366
367
368
369
370
371
      break;
    default:
      ERROR_EXIT("invalid dim\n");
      return 0;
    }
  }

372
373
374

  void Element::eraseNewCoord() 
  {
375
    if (newCoord != NULL) {
Thomas Witkowski's avatar
Thomas Witkowski committed
376
      delete newCoord;
377
378
      newCoord = NULL;
    }
379
  }
380
 
381

Thomas Witkowski's avatar
Thomas Witkowski committed
382
  void Element::serialize(std::ostream &out) 
383
384
  {
    // write children
385
    if (child[0]) {
Thomas Witkowski's avatar
Thomas Witkowski committed
386
      out << child[0]->getTypeName() << "\n";
387
388
389
      child[0]->serialize(out);
      child[1]->serialize(out);
    } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
390
      out << "NULL\n";
391
392
393
394
395
    }

    // write dofs
    int dim = mesh->getDim();
    int nodes = mesh->getNumberOfNodes();
Thomas Witkowski's avatar
Thomas Witkowski committed
396
    int j = 0;
397
    SerUtil::serialize(out, nodes);
Thomas Witkowski's avatar
Thomas Witkowski committed
398
399
   
    for (int pos = 0; pos <= dim; pos++) {
400
      GeoIndex position = INDEX_OF_DIM(pos, dim);
Thomas Witkowski's avatar
Thomas Witkowski committed
401
      int ndof = 0;
402

403
404
      for (int i = 0; i < mesh->getNumberOfDOFAdmin(); i++)
	ndof += mesh->getDOFAdmin(i).getNumberOfDOFs(position);	      
405

Thomas Witkowski's avatar
Thomas Witkowski committed
406
407
408
      if (ndof > 0) {
	for (int i = 0; i < mesh->getGeo(position); i++) {
	  if (dof[j] != NULL) {
409
410
411
412
413
	    // Create index to check if the dofs were already written.
	    std::pair<DegreeOfFreedom, int> idx = std::make_pair(dof[j][0], pos);

	    if (Mesh::serializedDOFs[idx] == NULL) {
	      Mesh::serializedDOFs[idx] = dof[j];
414
415
	      SerUtil::serialize(out, ndof);
	      SerUtil::serialize(out, pos);
416
417
	      out.write(reinterpret_cast<const char*>(dof[j]), 
			ndof * sizeof(DegreeOfFreedom));
418
419
	    } else {
	      int minusOne = -1;
420
421
	      SerUtil::serialize(out, minusOne);
	      SerUtil::serialize(out, pos);
422
423
	      out.write(reinterpret_cast<const char*>(&(dof[j][0])), 
			sizeof(DegreeOfFreedom));
424
425
426
	    }
	  } else {
	    int zero = 0;
427
428
	    SerUtil::serialize(out, zero);
	    SerUtil::serialize(out, pos);
429
430
431
432
433
434
435
	  }
	  j++;
	}
      }
    }

    // write index
436
    SerUtil::serialize(out, index);
437
438

    // write mark
439
    SerUtil::serialize(out, mark);
440
441

    // write newCoord
Thomas Witkowski's avatar
Thomas Witkowski committed
442
443
    if (newCoord) {
      out << "WorldVector\n";
444
445
      newCoord->serialize(out);
    } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
446
      out << "NULL\n";
447
448
449
    }

    // write element data
Thomas Witkowski's avatar
Thomas Witkowski committed
450
451
    if (elementData) {
      out << elementData->getTypeName() << "\n";
452
453
      elementData->serialize(out);
    } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
454
      out << "NULL\n";
455
456
457
    }
  }

458

Thomas Witkowski's avatar
Thomas Witkowski committed
459
  void Element::deserialize(std::istream &in)
460
  {
461
462
    FUNCNAME("Element::deserialize()");

463
    std::string typeName = "";
464
465
466
467
468

    // read children
    in >> typeName;
    in.get();

Thomas Witkowski's avatar
Thomas Witkowski committed
469
470
    if (typeName != "NULL") {
      if (typeName == "Line") {
Thomas Witkowski's avatar
Thomas Witkowski committed
471
472
	child[0] = new Line(NULL);
	child[1] = new Line(NULL);      
473
      }else  if (typeName == "Triangle") {
Thomas Witkowski's avatar
Thomas Witkowski committed
474
475
	child[0] = new Triangle(NULL);
	child[1] = new Triangle(NULL);      
476
      } else  if (typeName == "Tetrahedron") {
Thomas Witkowski's avatar
Thomas Witkowski committed
477
478
	child[0] = new Tetrahedron(NULL);
	child[1] = new Tetrahedron(NULL);      
479
480
      } else {
	ERROR_EXIT("Wrong element type!\n");
481
      }
482

483
484
485
486
487
488
489
490
      child[0]->deserialize(in);
      child[1]->deserialize(in);
    } else {
      child[0] = child[1] = NULL;
    }

    // read dofs
    int nodes;
491
    SerUtil::deserialize(in, nodes);
492

493
494
    TEST_EXIT_DBG(nodes > 0)("Should not happen!\n");

495
    dof = new DegreeOfFreedom*[nodes]; 
496

497
    for (int i = 0; i < nodes; i++) {
498
      int nDofs, pos;
499
500
      SerUtil::deserialize(in, nDofs);
      SerUtil::deserialize(in, pos);
501

502
503
504
505
      if (nDofs) {
	if (nDofs != -1) {
	  dof[i] = new DegreeOfFreedom[nDofs];
	  in.read(reinterpret_cast<char*>(dof[i]), nDofs * sizeof(DegreeOfFreedom));
506
507
508
509
510
511

	  // Create index to check if the dofs were alread read from file.
	  std::pair<DegreeOfFreedom, int> idx = std::make_pair(dof[i][0], pos);

	  if (Mesh::serializedDOFs[idx] != NULL) {
	    DegreeOfFreedom *dofPtr = Mesh::serializedDOFs[idx];
512
	    delete [] dof[i];
513
	    dof[i] = dofPtr;
514
	  } else {
515
	    Mesh::serializedDOFs[idx] = dof[i];
516
	  }
517

518
519
	} else {
	  DegreeOfFreedom index;
520
521
	  SerUtil::deserialize(in, index);
      
522
523
	  std::pair<DegreeOfFreedom, int> idx = std::make_pair(index, pos);
	  TEST_EXIT(Mesh::serializedDOFs.find(idx) !=  Mesh::serializedDOFs.end())
524
	    ("This should never happen!\n");
525
	  dof[i] = Mesh::serializedDOFs[idx];
526
527
528
529
530
	}
      } else {
	dof[i] = NULL;
      }
    }
531
   
532
    // read index
533
    SerUtil::deserialize(in, index);
534

535
    // read mark
536
    SerUtil::deserialize(in, mark);
537
538
539
540
541

    // read newCoord
    in >> typeName;
    in.get();

542
543
    if (typeName != "NULL") {
      if (typeName == "WorldVector") {
Thomas Witkowski's avatar
Thomas Witkowski committed
544
	newCoord = new WorldVector<double>;
545
546
547
548
549
550
551
552
553
554
555
556
	newCoord->deserialize(in);
      } else {
	ERROR_EXIT("unexpected type name\n");
      }
    } else {
      newCoord = NULL;
    }

    // read element data
    in >> typeName;
    in.get();

557
    if (typeName != "NULL") {
558
      elementData = CreatorMap<ElementData>::getCreator(typeName)->create();
559
560

      if (elementData)
561
	elementData->deserialize(in);
562
563
      else
	ERROR_EXIT("unexpected type name\n");      
564
565
566
567
568
    } else {
      elementData = NULL;
    }
  }

569

570
571
572
573
574
  int Element::calcMemoryUsage()
  {
    int result = 0;

    result += sizeof(Element);
575
576
    result += mesh->getNumberOfNodes() * sizeof(DegreeOfFreedom*);

577
578
    if (child[0])
      result += child[0]->calcMemoryUsage() + child[1]->calcMemoryUsage();    
579
580
581
582

    return result;
  }

583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633

  void writeDotFile(Element *el, std::string filename, int maxLevels)
  {
    std::vector<int> graphNodes;
    std::vector<std::pair<int, int> > graphEdges;

    std::vector<Element*> traverseElements, nextTraverseElements;
    traverseElements.push_back(el);

    int nLevels = 0;

    while (!traverseElements.empty() && (maxLevels == -1 || nLevels < maxLevels)) {
      nextTraverseElements.clear();

      for (unsigned int i = 0; i < traverseElements.size(); i++) {
	graphNodes.push_back(traverseElements[i]->getIndex());

	if (!traverseElements[i]->isLeaf() && 
	    (maxLevels == -1 || nLevels + 1 < maxLevels)) {
	  graphEdges.push_back(std::make_pair(traverseElements[i]->getIndex(),
					      traverseElements[i]->getChild(0)->getIndex()));
	  graphEdges.push_back(std::make_pair(traverseElements[i]->getIndex(),
					      traverseElements[i]->getChild(1)->getIndex()));
	  nextTraverseElements.push_back(traverseElements[i]->getChild(0));
	  nextTraverseElements.push_back(traverseElements[i]->getChild(1));
	}
      }

      traverseElements = nextTraverseElements;
      nLevels++;
    }

    std::ofstream file;
    file.open(filename.c_str());

    file << "digraph G\n";
    file << "{\n";

    for (unsigned int i = 0; i < graphNodes.size(); i++)
      file << "   node" << graphNodes[i] 
	   << " [ label = \"" << graphNodes[i] << "\"];\n";    

    for (unsigned int i = 0; i < graphEdges.size(); i++) 
      file << "   \"node" << graphEdges[i].first << "\" -> \"node"
	   << graphEdges[i].second << "\";\n";

    file << "}\n";

    file.close();
  }

634
}