Element.cc 14.9 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
#include "Element.h"
#include "DOFAdmin.h"
#include "Mesh.h"
#include "CoarseningManager.h"
#include "FixVec.h"
#include "ElementRegion_ED.h"
19
#include "Serializer.h"
20
#include "MeshStructure.h"
21
#include "BasisFunction.h"
22
23
24

namespace AMDiS {

25
  std::map<DegreeOfFreedom*, bool> Element::deletedDOFs;
26

27

28
29
30
31
  Element::Element(Mesh *aMesh)
  {
    mesh = aMesh;
    index = mesh ? mesh->getNextElementIndex() : -1; 
32
33
    child[0] = NULL;
    child[1] = NULL;
34
35
    newCoord = NULL;
    elementData = NULL;
36
    mark = 0;
37
    dofValid = true;
38

39
40
41
42
    if (mesh)
      createNewDofPtrs();
    else
      mesh = NULL;    
43
44
  }

45

46
47
  // call destructor through Mesh::freeElement !!!
  Element::~Element()
48
  {  
49
    if (child[0])
Thomas Witkowski's avatar
Thomas Witkowski committed
50
      delete child[0];
51
52
    
    if (child[1])
Thomas Witkowski's avatar
Thomas Witkowski committed
53
      delete child[1];   
54

55
    if (newCoord)
Thomas Witkowski's avatar
Thomas Witkowski committed
56
      delete newCoord;    
57
58
59

    if (elementData) {
      elementData->deleteDecorated();
Thomas Witkowski's avatar
Thomas Witkowski committed
60
      delete elementData;
61
62
63
    }
  }

64

Thomas Witkowski's avatar
Thomas Witkowski committed
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
  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;
  }


80
  void Element::createNewDofPtrs()
Thomas Witkowski's avatar
Thomas Witkowski committed
81
82
83
84
85
86
87
88
89
  {
    FUNCNAME("Element::setDofPtrs()");

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

    dof = mesh->createDofPtrs();  
  }


90
91
92
93
94
95
96
97
  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
98
	delete tmp;
99
100
101
102
103
104
105
	tmp = NULL;
	return true;
      } else {
	return elementData->deleteDecorated(typeID);
      }
    }
    return false;
106
107
  }

108

109
110
111
112
113
114
115
116
117
  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
118
      for (int i = 0; i < mesh->getNumberOfDOFAdmin(); i++)
119
	ndof += mesh->getDofAdmin(i).getNumberOfDofs(position);
120
121
122

      if (ndof > 0) {
	for (int i = 0; i < mesh->getGeo(position); i++) {
123
124
125
	  if (dof[j]) {
	    if (deletedDOFs.count(dof[j]) == 0) {
	      deletedDOFs[dof[j]] = true;
126
	      delete [] dof[j];
127
	    }
128
	  }  
129
130
131
132
133
	  j++;
	}
      }
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
134
    delete [] dof;
135
    
136
    if (child[0])
137
      child[0]->deleteElementDOFs();
138
    if (child[1])
139
140
141
      child[1]->deleteElementDOFs();
  }

142

143
144
145
  Element* Element::cloneWithDOFs()
  {
    Element *el;
146
    
147
    if (isLine()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
148
      el = new Line(NULL);
149
    } else if (isTriangle()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
150
      el = new Triangle(NULL);
151
    } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
152
      el = new Tetrahedron(NULL);
153
    }
154
155

    el->mesh = mesh;
156
157
158
    el->index = index;
    el->mark = mark;
    if (newCoord) {
Thomas Witkowski's avatar
Thomas Witkowski committed
159
      WorldVector<double> *nc = new WorldVector<double>();
160
161
162
      *nc = *newCoord;
      el->newCoord = nc;
    }
163
    
164
    /* =========== And here we clone the DOFs =========== */
165
   
166
    el->dof = new DegreeOfFreedom*[mesh->getNumberOfNodes()];
167
168
169
170
171
172
173
174

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

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

175
      for (int i = 0; i < mesh->getNumberOfDOFAdmin(); i++)
176
	ndof += mesh->getDofAdmin(i).getNumberOfDofs(position);
177
178
179
180

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

	    if (Mesh::serializedDOFs[idx] == NULL) {
184
185
	      el->dof[j] = new DegreeOfFreedom[ndof];
	      for (int k = 0; k < ndof; k++)
186
		el->dof[j][k] = dof[j][k];
187

188
	      Mesh::serializedDOFs[idx] = el->dof[j];
189
	    } else {
190
	      el->dof[j] = Mesh::serializedDOFs[idx];
191
192
193
194
195
196
197
198
	    }
	  } else {
	    el->dof[j] = NULL;
	  }
	  j++;
	}
      }
    }
199
    
200
201
    /* =========== And clone the children ============= */

202
    if (child[0]) 
203
      el->child[0] = child[0]->cloneWithDOFs();
204
    if (child[1])
205
      el->child[1] = child[1]->cloneWithDOFs();
206
207
208
209

    return el;
  }

210

211
212
213
214
215
216
  /****************************************************************************/
  /*  ATTENTION:                                                              */
  /*  new_dof_fct() destroys new_dof !!!!!!!!!!                               */
  /*  should be used only at the end of dof_compress()!!!!!                   */
  /****************************************************************************/

217
  void Element::newDofFct1(const DOFAdmin* admin, std::vector<int>& newDofIndex)
218
  {
219
    int n0, nd, nd0;
220

221
    if ((nd = admin->getNumberOfDofs(VERTEX)))  {
222
      int vertices = mesh->getGeo(VERTEX);
223
      nd0 = admin->getNumberOfPreDofs(VERTEX);
224
      n0 = admin->getMesh()->getNode(VERTEX);
225
      for (int i = 0; i < vertices; i++)
226
	changeDofs1(admin, newDofIndex, n0, nd0, nd, i);
227
228
    }

229
    if (mesh->getDim() > 1) {
230
      if ((nd = admin->getNumberOfDofs(EDGE)))  {
231
	int edges = mesh->getGeo(EDGE); 
232
	nd0 = admin->getNumberOfPreDofs(EDGE);
233
	n0 = admin->getMesh()->getNode(EDGE);
234
	for (int i = 0; i < edges; i++)
235
	  changeDofs1(admin, newDofIndex, n0, nd0, nd, i);
236
237
238
      }
    }

239
    if (mesh->getDim() == 3) {
240
      if ((nd = admin->getNumberOfDofs(FACE)))  {
241
	int faces = mesh->getGeo(FACE);
242
	nd0 = admin->getNumberOfPreDofs(FACE);
243
	n0 = admin->getMesh()->getNode(FACE);
244
	for (int i = 0; i < faces; i++)
245
	  changeDofs1(admin, newDofIndex, n0, nd0, nd, i);
246
247
248
      }
    }

249
    if ((nd = admin->getNumberOfDofs(CENTER)))  {
250
      nd0 = admin->getNumberOfPreDofs(CENTER);
251
      n0 = admin->getMesh()->getNode(CENTER);      
252
      changeDofs1(admin, newDofIndex, n0, nd0, nd, 0);
253
254
255
256
    }
  }


257
  void Element::newDofFct2(const DOFAdmin* admin)
258
  {
259
    int n0, nd0;
260

261
    int nd = admin->getNumberOfDofs(VERTEX);
262
    if (nd) {
263
      int vertices = mesh->getGeo(VERTEX);
264
      nd0 = admin->getNumberOfPreDofs(VERTEX);
265
      n0 = admin->getMesh()->getNode(VERTEX);
266
267
      for (int i = 0; i < vertices; i++)
	changeDofs2(n0, nd0, nd, i);
268
269
    }

270
    if (mesh->getDim() > 1) {
271
      nd = admin->getNumberOfDofs(EDGE);
272
      if (nd) {
273
	int edges = mesh->getGeo(EDGE); 
274
	nd0 = admin->getNumberOfPreDofs(EDGE);
275
	n0 = admin->getMesh()->getNode(EDGE);
276
277
	for (int i = 0; i < edges; i++)
	  changeDofs2(n0, nd0, nd, i);
278
279
280
      }
    }

281
    if (mesh->getDim() == 3) {
282
      nd = admin->getNumberOfDofs(FACE);
283
      if (nd) {
284
	int faces = mesh->getGeo(FACE);
285
	nd0 = admin->getNumberOfPreDofs(FACE);
286
	n0 = admin->getMesh()->getNode(FACE);
287
288
	for (int i = 0; i < faces; i++)
	  changeDofs2(n0, nd0, nd, i);
289
290
291
      }
    }

292
    nd = admin->getNumberOfDofs(CENTER);
293
    if (nd) {
294
      nd0 = admin->getNumberOfPreDofs(CENTER);
295
      n0 = admin->getMesh()->getNode(CENTER);
296
      // only one center
297
      changeDofs2(n0, nd0, nd, 0);	  
298
299
300
    }
  }

301

302
303
  void Element::changeDofs1(const DOFAdmin* admin, std::vector<int>& newDofIndex,
			    int n0, int nd0, int nd, int pos)
304
305
306
307
308
  {
    DegreeOfFreedom *ldof = dof[n0 + pos] + nd0;
    for (int j = 0; j < nd; j++) {
      int k = ldof[j];
      if (k >= 0)
309
	ldof[j] = -newDofIndex[k] - 1;
310
311
312
313
314
315
316
317
318
319
320
321
322
    }
  }
  
  
  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;      
    }
  }
323

324

325
326
327
328
329
330
331
  /****************************************************************************/
  /* 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
  {
332
333
    int nv = 0;
    int ov = 0;
334
335
336
    int vertices = mesh->getGeo(VERTEX);
    int dim = mesh->getDim();

337
    for (int i = 0; i < vertices; i++) {
338
      if (nv < i - 1)  
339
340
341
342
	return(-1);

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

    }
    
    if (nv != mesh->getDim()) 
      return(-1);
354
355

    // The opposite vertex is 3(6) - (sum of indices of common vertices) in 2D(3D)
356
357
358
359
360
361

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

373
374
375

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

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

    // write dofs
    int dim = mesh->getDim();
    int nodes = mesh->getNumberOfNodes();
Thomas Witkowski's avatar
Thomas Witkowski committed
397
    int j = 0;
398
    SerUtil::serialize(out, nodes);
399
    SerUtil::serialize(out, dofValid);
Thomas Witkowski's avatar
Thomas Witkowski committed
400
   
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
    if (dofValid) {
      for (int pos = 0; pos <= dim; pos++) {
	GeoIndex position = INDEX_OF_DIM(pos, dim);
	int ndof = 0;
	
	for (int i = 0; i < mesh->getNumberOfDOFAdmin(); i++)
	  ndof += mesh->getDofAdmin(i).getNumberOfDofs(position);	      
	
	if (ndof > 0) {
	  for (int i = 0; i < mesh->getGeo(position); i++) {
	    if (dof[j] != NULL) {
	      // 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];
		SerUtil::serialize(out, ndof);
		SerUtil::serialize(out, pos);
		out.write(reinterpret_cast<const char*>(dof[j]), 
			  ndof * sizeof(DegreeOfFreedom));
	      } else {
		int minusOne = -1;
		SerUtil::serialize(out, minusOne);
		SerUtil::serialize(out, pos);
		out.write(reinterpret_cast<const char*>(&(dof[j][0])), 
			  sizeof(DegreeOfFreedom));
	      }
428
	    } else {
429
430
	      int zero = 0;
	      SerUtil::serialize(out, zero);
431
	      SerUtil::serialize(out, pos);
432
	    }
433
	    j++;
434
435
436
437
438
439
	  }
	}
      }
    }

    // write index
440
    SerUtil::serialize(out, index);
441
442

    // write mark
443
    SerUtil::serialize(out, mark);
444
445

    // write newCoord
Thomas Witkowski's avatar
Thomas Witkowski committed
446
447
    if (newCoord) {
      out << "WorldVector\n";
448
449
      newCoord->serialize(out);
    } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
450
      out << "NULL\n";
451
452
453
    }

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

462

Thomas Witkowski's avatar
Thomas Witkowski committed
463
  void Element::deserialize(std::istream &in)
464
  {
465
466
    FUNCNAME("Element::deserialize()");

467
    std::string typeName = "";
468
469
470
471
472

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

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

487
488
489
490
491
492
493
494
      child[0]->deserialize(in);
      child[1]->deserialize(in);
    } else {
      child[0] = child[1] = NULL;
    }

    // read dofs
    int nodes;
495
    SerUtil::deserialize(in, nodes);
496
    SerUtil::deserialize(in, dofValid);
497

498
    TEST_EXIT_DBG(nodes > 0)("Should not happen!\n");      
499
    dof = new DegreeOfFreedom*[nodes]; 
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
     
    if (dofValid) {
      for (int i = 0; i < nodes; i++) {
	int nDofs, pos;
	SerUtil::deserialize(in, nDofs);
	SerUtil::deserialize(in, pos);
	
	if (nDofs) {
	  if (nDofs != -1) {
	    dof[i] = new DegreeOfFreedom[nDofs];
	    in.read(reinterpret_cast<char*>(dof[i]), nDofs * sizeof(DegreeOfFreedom));
	    
	    // 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];
	      delete [] dof[i];
	      dof[i] = dofPtr;
	    } else {
	      Mesh::serializedDOFs[idx] = dof[i];
	    }
	    
523
	  } else {
524
525
526
527
528
529
530
	    DegreeOfFreedom index;
	    SerUtil::deserialize(in, index);
	    
	    std::pair<DegreeOfFreedom, int> idx = std::make_pair(index, pos);
	    TEST_EXIT(Mesh::serializedDOFs.find(idx) !=  Mesh::serializedDOFs.end())
	      ("This should never happen!\n");
	    dof[i] = Mesh::serializedDOFs[idx];
531
532
	  }
	} else {
533
	  dof[i] = NULL;
534
535
536
	}
      }
    }
537
      
538
    // read index
539
    SerUtil::deserialize(in, index);
540

541
    // read mark
542
    SerUtil::deserialize(in, mark);
543
544
545
546
547

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

548
549
    if (typeName != "NULL") {
      if (typeName == "WorldVector") {
Thomas Witkowski's avatar
Thomas Witkowski committed
550
	newCoord = new WorldVector<double>;
551
552
553
554
555
556
557
558
559
560
561
562
	newCoord->deserialize(in);
      } else {
	ERROR_EXIT("unexpected type name\n");
      }
    } else {
      newCoord = NULL;
    }

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

563
    if (typeName != "NULL") {
564
      elementData = CreatorMap<ElementData>::getCreator(typeName)->create();
565
566

      if (elementData)
567
	elementData->deserialize(in);
568
569
      else
	ERROR_EXIT("unexpected type name\n");      
570
571
572
573
574
    } else {
      elementData = NULL;
    }
  }

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
622
623
624
625

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

626
627
628
629
630
631
632
633
634
635
636

  void Element::getAllDofs(FiniteElemSpace* feSpace, 
			   BoundaryObject bound, 
			   DofContainer& dofs)
  {
    getNodeDofs(feSpace, bound, dofs);

    if (feSpace->getBasisFcts()->getDegree() > 1)
      getHigherOrderDofs(feSpace, bound, dofs);
  }

637
}