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

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

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

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

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

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

34
35
36
37
38
39
40
    dof = mesh->createDOFPtrs();  
  }

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

    if (mesh) {
      setDOFPtrs();
    } else {
      mesh = NULL;
    }
52
53
54
55
  }

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

63
    if (newCoord)
Thomas Witkowski's avatar
Thomas Witkowski committed
64
      delete newCoord;    
65
66
67

    if (elementData) {
      elementData->deleteDecorated();
Thomas Witkowski's avatar
Thomas Witkowski committed
68
      delete elementData;
69
70
71
72
73
74
75
76
77
78
79
    }
  }

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

90
91
92
93
94
95
96
97
98
  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
99
      for (int i = 0; i < mesh->getNumberOfDOFAdmin(); i++)
100
101
102
103
	ndof += mesh->getDOFAdmin(i).getNumberOfDOFs(position);

      if (ndof > 0) {
	for (int i = 0; i < mesh->getGeo(position); i++) {
104
105
106
	  if (dof[j]) {
	    if (deletedDOFs.count(dof[j]) == 0) {
	      deletedDOFs[dof[j]] = true;
107
	      delete [] dof[j];
108
	    }
109
	  }  
110
111
112
113
114
	  j++;
	}
      }
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
115
    delete [] dof;
116
    
117
    if (child[0])
118
      child[0]->deleteElementDOFs();
119
    if (child[1])
120
121
122
      child[1]->deleteElementDOFs();
  }

123
124
125
  Element* Element::cloneWithDOFs()
  {
    Element *el;
126
    
127
    if (isLine()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
128
      el = new Line(NULL);
129
    } else if (isTriangle()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
130
      el = new Triangle(NULL);
131
    } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
132
      el = new Tetrahedron(NULL);
133
    }
134
135

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

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

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

155
      for (int i = 0; i < mesh->getNumberOfDOFAdmin(); i++)
156
157
158
159
160
	ndof += mesh->getDOFAdmin(i).getNumberOfDOFs(position);

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

	    if (Mesh::serializedDOFs[idx] == NULL) {
164
165
	      el->dof[j] = new DegreeOfFreedom[ndof];
	      for (int k = 0; k < ndof; k++)
166
		el->dof[j][k] = dof[j][k];
167

168
	      Mesh::serializedDOFs[idx] = el->dof[j];
169
	    } else {
170
	      el->dof[j] = Mesh::serializedDOFs[idx];
171
172
173
174
175
176
177
178
	    }
	  } else {
	    el->dof[j] = NULL;
	  }
	  j++;
	}
      }
    }
179
    
180
181
    /* =========== And clone the children ============= */

182
    if (child[0]) 
183
      el->child[0] = child[0]->cloneWithDOFs();
184
    if (child[1])
185
      el->child[1] = child[1]->cloneWithDOFs();
186
187
188
189

    return el;
  }

190
191
192
193
194
195
196
197
198
  /****************************************************************************/
  /*  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)						\
199
  ldof = el->dof[n0 + i] + nd0;						\
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
    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)
  {
218
219
    int j, k, n0, nd, nd0;
    DegreeOfFreedom *ldof;  
220
221
222
223
224
225
226
    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);
227
      for (int i = 0; i < vertices; i++) {
228
229
230
231
	CHANGE_DOFS_1(this);
      }
    }

232
    if (mesh->getDim() > 1) {
233
234
235
      if ((nd = admin->getNumberOfDOFs(EDGE)))  {
	nd0 = admin->getNumberOfPreDOFs(EDGE);
	n0 = admin->getMesh()->getNode(EDGE);
236
	for (int i = 0; i < edges; i++) {
237
238
239
240
241
	  CHANGE_DOFS_1(this);
	}
      }
    }

242
    if (mesh->getDim() == 3) {
243
244
245
      if ((nd = admin->getNumberOfDOFs(FACE)))  {
	nd0 = admin->getNumberOfPreDOFs(FACE);
	n0 = admin->getMesh()->getNode(FACE);
246
	for (int i = 0; i < faces; i++) {
247
248
249
250
251
252
253
254
	  CHANGE_DOFS_1(this);
	}
      }
    }

    if ((nd = admin->getNumberOfDOFs(CENTER)))  {
      nd0 = admin->getNumberOfPreDOFs(CENTER);
      n0 = admin->getMesh()->getNode(CENTER);
255
      int i = 0;          /* only one center */
256
257
258
259
260
261
262
      CHANGE_DOFS_1(this);
    }
  }


  void Element::newDOFFct2(const DOFAdmin* admin)
  {
263
    int i, j, k, n0, nd0;
264
    DegreeOfFreedom  *ldof;
265
266
267
268
    int vertices = mesh->getGeo(VERTEX);
    int edges = mesh->getGeo(EDGE); 
    int faces = mesh->getGeo(FACE);

269
270
    int nd = admin->getNumberOfDOFs(VERTEX);
    if (nd) {
271
272
273
274
275
276
277
      nd0 = admin->getNumberOfPreDOFs(VERTEX);
      n0 = admin->getMesh()->getNode(VERTEX);
      for (i = 0; i < vertices; i++) {
	CHANGE_DOFS_2(this);
      }
    }

278
    if (mesh->getDim() > 1) {
279
280
      nd = admin->getNumberOfDOFs(EDGE);
      if (nd) {
281
282
283
284
285
286
287
288
	nd0 = admin->getNumberOfPreDOFs(EDGE);
	n0 = admin->getMesh()->getNode(EDGE);
	for (i = 0; i < edges; i++) {
	  CHANGE_DOFS_2(this);
	}
      }
    }

289
    if (mesh->getDim() == 3) {
290
291
      nd = admin->getNumberOfDOFs(FACE);
      if (nd) {
292
293
294
295
296
297
298
299
	nd0 = admin->getNumberOfPreDOFs(FACE);
	n0 = admin->getMesh()->getNode(FACE);
	for (i = 0; i < faces; i++) {
	  CHANGE_DOFS_2(this);
	}
      }
    }

300
301
    nd = admin->getNumberOfDOFs(CENTER);
    if (nd) {
302
303
      nd0 = admin->getNumberOfPreDOFs(CENTER);
      n0 = admin->getMesh()->getNode(CENTER);
304
305
      // only one center
      i = 0;   
306
307
308
309
310
311
312
313
314
315
316
317
318
319
      CHANGE_DOFS_2(this);
    }
  }

#undef CHANGE_DOFS_1
#undef CHANGE_DOFS_2

  /****************************************************************************/
  /* 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
  {
320
321
    int nv = 0;
    int ov = 0;
322
323
324
    int vertices = mesh->getGeo(VERTEX);
    int dim = mesh->getDim();

325
    for (int i = 0; i < vertices; i++) {
326
      if (nv < i - 1)  
327
328
329
330
331
332
333
334
335
336
337
	return(-1);

      for (int j = 0; j < dim; j++) {
	if (dof[i] == pdof[j]) {
	  /****************************************************************************/
	  /* i is a common vertex                                                     */
	  /****************************************************************************/
	  ov += i;
	  nv++;
	  break;
	}
338
      }
339
340
341
342
343

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

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

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

393
394
      for (int i = 0; i < mesh->getNumberOfDOFAdmin(); i++)
	ndof += mesh->getDOFAdmin(i).getNumberOfDOFs(position);	      
395

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

    // write index
426
    SerUtil::serialize(out, index);
427
428

    // write mark
429
    SerUtil::serialize(out, mark);
430
431

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

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

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

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

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

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

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

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

482
    dof = new DegreeOfFreedom*[nodes]; 
483

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

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

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

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

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

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

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

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

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

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

556
557
558
559
560
  int Element::calcMemoryUsage()
  {
    int result = 0;

    result += sizeof(Element);
561
562
    result += mesh->getNumberOfNodes() * sizeof(DegreeOfFreedom*);

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

    return result;
  }

569
570
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
  void fitElementToMeshCode(RefinementManager *refineManager, MeshStructure &code, 
			    Element *el, int ithSide, int elType)
  {
    if (code.isLeafElement())
      return;

    if (el->isLeaf()) {
      el->setMark(1);
      refineManager->refineMesh(el->getMesh());
    }

    int s1 = el->getSideOfChild(0, ithSide, elType);
    int s2 = el->getSideOfChild(1, ithSide, elType);

    code.nextElement();
    
    if (s1 != -1) 
      fitElementToMeshCode(refineManager, code, 
			   el->getFirstChild(), s1, el->getChildType(elType));

    code.nextElement();
    
    if (s2 != -1)
      fitElementToMeshCode(refineManager, code, 
			   el->getSecondChild(), s2, el->getChildType(elType));
  }

596
}