Element.cc 15.1 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

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

44

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

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

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

63

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


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

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

    dof = mesh->createDofPtrs();  
  }


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

107

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

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

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

141

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

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

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

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

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

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

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

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

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

    return el;
  }

209

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

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

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

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

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

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


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

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

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

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

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

300

301
302
  void Element::changeDofs1(const DOFAdmin* admin, std::vector<int>& newDofIndex,
			    int n0, int nd0, int nd, int pos)
303
  {
304
305
    FUNCNAME("Element::changeDofs1()");

306
307
308
309
    DegreeOfFreedom *ldof = dof[n0 + pos] + nd0;
    for (int j = 0; j < nd; j++) {
      int k = ldof[j];
      if (k >= 0)
310
	ldof[j] = -newDofIndex[k] - 1;
311
312
313
314
315
316
317
318
319
320
321
322
323
    }
  }
  
  
  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;      
    }
  }
324

325

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

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

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

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

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

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

374
375
376

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

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

    // write dofs
    int dim = mesh->getDim();
    int nodes = mesh->getNumberOfNodes();
Thomas Witkowski's avatar
Thomas Witkowski committed
398
    int j = 0;
399
    SerUtil::serialize(out, nodes);
400
    int dofValid = (dof != NULL ? 1 : 0);
401
    SerUtil::serialize(out, dofValid);
Thomas Witkowski's avatar
Thomas Witkowski committed
402
   
403
    if (dof != NULL) {
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
      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));
	      }
430
	    } else {
431
432
	      int zero = 0;
	      SerUtil::serialize(out, zero);
433
	      SerUtil::serialize(out, pos);
434
	    }
435
	    j++;
436
437
438
439
440
441
	  }
	}
      }
    }

    // write index
442
    SerUtil::serialize(out, index);
443
444

    // write mark
445
    SerUtil::serialize(out, mark);
446
447

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

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

464

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

469
    std::string typeName = "";
470
471
472
473
474

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

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

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

    // read dofs
    int nodes;
497
    SerUtil::deserialize(in, nodes);
498
    int dofValid = 0;
499
    SerUtil::deserialize(in, dofValid);
500

501
502
    TEST_EXIT_DBG(nodes > 0)("Should not happen!\n");    
    dof = (dofValid ? new DegreeOfFreedom*[nodes] : NULL);
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
     
    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];
	    }
	    
526
	  } else {
527
528
529
530
531
532
533
	    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];
534
535
	  }
	} else {
536
	  dof[i] = NULL;
537
538
539
	}
      }
    }
540
      
541
    // read index
542
    SerUtil::deserialize(in, index);
543

544
    // read mark
545
    SerUtil::deserialize(in, mark);
546
547
548
549
550

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

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

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

566
    if (typeName != "NULL") {
567
568
      elementData = 
	CreatorMap<ElementData>::getCreator(typeName, "deserializetion from file")->create();
569
570

      if (elementData)
571
	elementData->deserialize(in);
572
573
      else
	ERROR_EXIT("unexpected type name\n");      
574
575
576
577
578
    } else {
      elementData = NULL;
    }
  }

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
626
627
628
629

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

630

631
  void Element::getAllDofs(const FiniteElemSpace* feSpace, 
632
			   BoundaryObject bound, 
633
634
			   DofContainer& dofs,
			   bool baseDofPtr)
635
  {
636
    getNodeDofs(feSpace, bound, dofs, baseDofPtr);
637
638

    if (feSpace->getBasisFcts()->getDegree() > 1)
639
      getHigherOrderDofs(feSpace, bound, dofs, baseDofPtr);
640
641
  }

642
}