Element.cc 14.8 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
17
18
  Element::Element(Mesh *aMesh)
  {
    mesh = aMesh;
    index = mesh ? mesh->getNextElementIndex() : -1; 
19
20
    child[0] = NULL;
    child[1] = NULL;
21
22
    newCoord = NULL;
    elementData = NULL;
23
    mark = 0;
24

25
26
27
28
    if (mesh)
      createNewDofPtrs();
    else
      mesh = NULL;    
29
30
  }

31

32
33
  // call destructor through Mesh::freeElement !!!
  Element::~Element()
34
  {  
35
    if (child[0])
Thomas Witkowski's avatar
Thomas Witkowski committed
36
      delete child[0];
37
38
    
    if (child[1])
Thomas Witkowski's avatar
Thomas Witkowski committed
39
      delete child[1];   
40

41
    if (newCoord)
Thomas Witkowski's avatar
Thomas Witkowski committed
42
      delete newCoord;    
43
44
45

    if (elementData) {
      elementData->deleteDecorated();
Thomas Witkowski's avatar
Thomas Witkowski committed
46
      delete elementData;
47
48
49
    }
  }

50

Thomas Witkowski's avatar
Thomas Witkowski committed
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
  int Element::getRegion() const 
  {
    if (!elementData) 
      return -1;

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

    if (red)
      return red->getRegion();    
    
    return -1;
  }


66
  void Element::createNewDofPtrs()
Thomas Witkowski's avatar
Thomas Witkowski committed
67
68
69
70
71
72
73
74
75
  {
    FUNCNAME("Element::setDofPtrs()");

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

    dof = mesh->createDofPtrs();  
  }


76
77
78
79
80
81
82
83
  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
84
	delete tmp;
85
86
87
88
89
90
91
	tmp = NULL;
	return true;
      } else {
	return elementData->deleteDecorated(typeID);
      }
    }
    return false;
92
93
  }

94

95
96
97
98
99
100
101
102
103
  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
104
      for (int i = 0; i < mesh->getNumberOfDOFAdmin(); i++)
105
	ndof += mesh->getDofAdmin(i).getNumberOfDofs(position);
106
107
108

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

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

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

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

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

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

160
      for (int i = 0; i < mesh->getNumberOfDOFAdmin(); i++)
161
	ndof += mesh->getDofAdmin(i).getNumberOfDofs(position);
162
163
164
165

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

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

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

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

    return el;
  }

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

201
  void Element::newDofFct1(const DOFAdmin* admin)
202
  {
203
    int n0, nd, nd0;
204

205
    if ((nd = admin->getNumberOfDofs(VERTEX)))  {
206
      int vertices = mesh->getGeo(VERTEX);
207
208
      nd0 = admin->getNumberOfPreDOFs(VERTEX);
      n0 = admin->getMesh()->getNode(VERTEX);
209
210
      for (int i = 0; i < vertices; i++)
	changeDofs1(admin, n0, nd0, nd, i);
211
212
    }

213
    if (mesh->getDim() > 1) {
214
      if ((nd = admin->getNumberOfDofs(EDGE)))  {
215
	int edges = mesh->getGeo(EDGE); 
216
217
	nd0 = admin->getNumberOfPreDOFs(EDGE);
	n0 = admin->getMesh()->getNode(EDGE);
218
219
	for (int i = 0; i < edges; i++)
	  changeDofs1(admin, n0, nd0, nd, i);
220
221
222
      }
    }

223
    if (mesh->getDim() == 3) {
224
      if ((nd = admin->getNumberOfDofs(FACE)))  {
225
	int faces = mesh->getGeo(FACE);
226
227
	nd0 = admin->getNumberOfPreDOFs(FACE);
	n0 = admin->getMesh()->getNode(FACE);
228
229
	for (int i = 0; i < faces; i++)
	  changeDofs1(admin, n0, nd0, nd, i);
230
231
232
      }
    }

233
    if ((nd = admin->getNumberOfDofs(CENTER)))  {
234
      nd0 = admin->getNumberOfPreDOFs(CENTER);
235
236
      n0 = admin->getMesh()->getNode(CENTER);      
      changeDofs1(admin, n0, nd0, nd, 0);
237
238
239
240
    }
  }


241
  void Element::newDofFct2(const DOFAdmin* admin)
242
  {
243
    int n0, nd0;
244

245
    int nd = admin->getNumberOfDofs(VERTEX);
246
    if (nd) {
247
      int vertices = mesh->getGeo(VERTEX);
248
249
      nd0 = admin->getNumberOfPreDOFs(VERTEX);
      n0 = admin->getMesh()->getNode(VERTEX);
250
251
      for (int i = 0; i < vertices; i++)
	changeDofs2(n0, nd0, nd, i);
252
253
    }

254
    if (mesh->getDim() > 1) {
255
      nd = admin->getNumberOfDofs(EDGE);
256
      if (nd) {
257
	int edges = mesh->getGeo(EDGE); 
258
259
	nd0 = admin->getNumberOfPreDOFs(EDGE);
	n0 = admin->getMesh()->getNode(EDGE);
260
261
	for (int i = 0; i < edges; i++)
	  changeDofs2(n0, nd0, nd, i);
262
263
264
      }
    }

265
    if (mesh->getDim() == 3) {
266
      nd = admin->getNumberOfDofs(FACE);
267
      if (nd) {
268
	int faces = mesh->getGeo(FACE);
269
270
	nd0 = admin->getNumberOfPreDOFs(FACE);
	n0 = admin->getMesh()->getNode(FACE);
271
272
	for (int i = 0; i < faces; i++)
	  changeDofs2(n0, nd0, nd, i);
273
274
275
      }
    }

276
    nd = admin->getNumberOfDofs(CENTER);
277
    if (nd) {
278
279
      nd0 = admin->getNumberOfPreDOFs(CENTER);
      n0 = admin->getMesh()->getNode(CENTER);
280
      // only one center
281
      changeDofs2(n0, nd0, nd, 0);	  
282
283
284
    }
  }

285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305

  void Element::changeDofs1(const DOFAdmin* admin, int n0, int nd0, int nd, int pos)
  {
    DegreeOfFreedom *ldof = dof[n0 + pos] + nd0;
    for (int j = 0; j < nd; j++) {
      int k = ldof[j];
      if (k >= 0)
	ldof[j] = -admin->getMesh()->newDOF[k] - 1;
    }
  }
  
  
  void Element::changeDofs2(int n0, int nd0, int nd, int pos)
  {
    DegreeOfFreedom *ldof = dof[n0 + pos] + nd0;
    for (int j = 0; j < nd; j++) {
      int k = ldof[j];
      if (k < 0)
	ldof[j] = -k - 1;      
    }
  }
306

307

308
309
310
311
312
313
314
  /****************************************************************************/
  /* 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
  {
315
316
    int nv = 0;
    int ov = 0;
317
318
319
    int vertices = mesh->getGeo(VERTEX);
    int dim = mesh->getDim();

320
    for (int i = 0; i < vertices; i++) {
321
      if (nv < i - 1)  
322
323
324
325
326
327
328
329
330
331
332
	return(-1);

      for (int j = 0; j < dim; j++) {
	if (dof[i] == pdof[j]) {
	  /****************************************************************************/
	  /* i is a common vertex                                                     */
	  /****************************************************************************/
	  ov += i;
	  nv++;
	  break;
	}
333
      }
334
335
336
337
338

    }
    
    if (nv != mesh->getDim()) 
      return(-1);
339
340
341
342
343
344
345
346
347
348
    /****************************************************************************/
    /*  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:
349
      return 3 - ov;
350
351
      break;
    case 3:
352
      return 6 - ov;
353
354
355
356
357
358
359
      break;
    default:
      ERROR_EXIT("invalid dim\n");
      return 0;
    }
  }

360
361
362

  void Element::eraseNewCoord() 
  {
363
    if (newCoord != NULL) {
Thomas Witkowski's avatar
Thomas Witkowski committed
364
      delete newCoord;
365
366
      newCoord = NULL;
    }
367
  }
368
 
369

Thomas Witkowski's avatar
Thomas Witkowski committed
370
  void Element::serialize(std::ostream &out) 
371
372
  {
    // write children
373
    if (child[0]) {
Thomas Witkowski's avatar
Thomas Witkowski committed
374
      out << child[0]->getTypeName() << "\n";
375
376
377
      child[0]->serialize(out);
      child[1]->serialize(out);
    } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
378
      out << "NULL\n";
379
380
381
382
383
    }

    // write dofs
    int dim = mesh->getDim();
    int nodes = mesh->getNumberOfNodes();
Thomas Witkowski's avatar
Thomas Witkowski committed
384
    int j = 0;
385
    SerUtil::serialize(out, nodes);
Thomas Witkowski's avatar
Thomas Witkowski committed
386
387
   
    for (int pos = 0; pos <= dim; pos++) {
388
      GeoIndex position = INDEX_OF_DIM(pos, dim);
Thomas Witkowski's avatar
Thomas Witkowski committed
389
      int ndof = 0;
390

391
      for (int i = 0; i < mesh->getNumberOfDOFAdmin(); i++)
392
	ndof += mesh->getDofAdmin(i).getNumberOfDofs(position);	      
393

Thomas Witkowski's avatar
Thomas Witkowski committed
394
395
396
      if (ndof > 0) {
	for (int i = 0; i < mesh->getGeo(position); i++) {
	  if (dof[j] != NULL) {
397
398
399
400
401
	    // 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];
402
403
	      SerUtil::serialize(out, ndof);
	      SerUtil::serialize(out, pos);
404
405
	      out.write(reinterpret_cast<const char*>(dof[j]), 
			ndof * sizeof(DegreeOfFreedom));
406
407
	    } else {
	      int minusOne = -1;
408
409
	      SerUtil::serialize(out, minusOne);
	      SerUtil::serialize(out, pos);
410
411
	      out.write(reinterpret_cast<const char*>(&(dof[j][0])), 
			sizeof(DegreeOfFreedom));
412
413
414
	    }
	  } else {
	    int zero = 0;
415
416
	    SerUtil::serialize(out, zero);
	    SerUtil::serialize(out, pos);
417
418
419
420
421
422
423
	  }
	  j++;
	}
      }
    }

    // write index
424
    SerUtil::serialize(out, index);
425
426

    // write mark
427
    SerUtil::serialize(out, mark);
428
429

    // write newCoord
Thomas Witkowski's avatar
Thomas Witkowski committed
430
431
    if (newCoord) {
      out << "WorldVector\n";
432
433
      newCoord->serialize(out);
    } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
434
      out << "NULL\n";
435
436
437
    }

    // write element data
Thomas Witkowski's avatar
Thomas Witkowski committed
438
439
    if (elementData) {
      out << elementData->getTypeName() << "\n";
440
441
      elementData->serialize(out);
    } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
442
      out << "NULL\n";
443
444
445
    }
  }

446

Thomas Witkowski's avatar
Thomas Witkowski committed
447
  void Element::deserialize(std::istream &in)
448
  {
449
450
    FUNCNAME("Element::deserialize()");

451
    std::string typeName = "";
452
453
454
455
456

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

Thomas Witkowski's avatar
Thomas Witkowski committed
457
458
    if (typeName != "NULL") {
      if (typeName == "Line") {
Thomas Witkowski's avatar
Thomas Witkowski committed
459
460
	child[0] = new Line(NULL);
	child[1] = new Line(NULL);      
461
      }else  if (typeName == "Triangle") {
Thomas Witkowski's avatar
Thomas Witkowski committed
462
463
	child[0] = new Triangle(NULL);
	child[1] = new Triangle(NULL);      
464
      } else  if (typeName == "Tetrahedron") {
Thomas Witkowski's avatar
Thomas Witkowski committed
465
466
	child[0] = new Tetrahedron(NULL);
	child[1] = new Tetrahedron(NULL);      
467
468
      } else {
	ERROR_EXIT("Wrong element type!\n");
469
      }
470

471
472
473
474
475
476
477
478
      child[0]->deserialize(in);
      child[1]->deserialize(in);
    } else {
      child[0] = child[1] = NULL;
    }

    // read dofs
    int nodes;
479
    SerUtil::deserialize(in, nodes);
480

481
482
    TEST_EXIT_DBG(nodes > 0)("Should not happen!\n");

483
    dof = new DegreeOfFreedom*[nodes]; 
484

485
    for (int i = 0; i < nodes; i++) {
486
      int nDofs, pos;
487
488
      SerUtil::deserialize(in, nDofs);
      SerUtil::deserialize(in, pos);
489

490
491
492
493
      if (nDofs) {
	if (nDofs != -1) {
	  dof[i] = new DegreeOfFreedom[nDofs];
	  in.read(reinterpret_cast<char*>(dof[i]), nDofs * sizeof(DegreeOfFreedom));
494
495
496
497
498
499

	  // 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];
500
	    delete [] dof[i];
501
	    dof[i] = dofPtr;
502
	  } else {
503
	    Mesh::serializedDOFs[idx] = dof[i];
504
	  }
505

506
507
	} else {
	  DegreeOfFreedom index;
508
509
	  SerUtil::deserialize(in, index);
      
510
511
	  std::pair<DegreeOfFreedom, int> idx = std::make_pair(index, pos);
	  TEST_EXIT(Mesh::serializedDOFs.find(idx) !=  Mesh::serializedDOFs.end())
512
	    ("This should never happen!\n");
513
	  dof[i] = Mesh::serializedDOFs[idx];
514
515
516
517
518
	}
      } else {
	dof[i] = NULL;
      }
    }
519
   
520
    // read index
521
    SerUtil::deserialize(in, index);
522

523
    // read mark
524
    SerUtil::deserialize(in, mark);
525
526
527
528
529

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

530
531
    if (typeName != "NULL") {
      if (typeName == "WorldVector") {
Thomas Witkowski's avatar
Thomas Witkowski committed
532
	newCoord = new WorldVector<double>;
533
534
535
536
537
538
539
540
541
542
543
544
	newCoord->deserialize(in);
      } else {
	ERROR_EXIT("unexpected type name\n");
      }
    } else {
      newCoord = NULL;
    }

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

545
    if (typeName != "NULL") {
546
      elementData = CreatorMap<ElementData>::getCreator(typeName)->create();
547
548

      if (elementData)
549
	elementData->deserialize(in);
550
551
      else
	ERROR_EXIT("unexpected type name\n");      
552
553
554
555
556
    } else {
      elementData = NULL;
    }
  }

557

558
559
560
561
562
  int Element::calcMemoryUsage()
  {
    int result = 0;

    result += sizeof(Element);
563
564
    result += mesh->getNumberOfNodes() * sizeof(DegreeOfFreedom*);

565
566
    if (child[0])
      result += child[0]->calcMemoryUsage() + child[1]->calcMemoryUsage();    
567
568
569
570

    return result;
  }

571
572
573
574
575
576
577
578
579
580
581
582
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

  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();
  }

622
}