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

namespace AMDiS {

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

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

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

19
20
21
    if (red) {
      return red->getRegion();
    }
22
23
24
25
26
27
    
    return -1;
  }

  void Element::setDOFPtrs()
  {
28
    FUNCNAME("Element::setDOFPtrs()");
29
30

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

32
33
34
35
36
37
38
    dof = mesh->createDOFPtrs();  
  }

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

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

  // call destructor through Mesh::freeElement !!!
  Element::~Element()
55
  {  
56
    if (child[0]) {
57
      DELETE child[0];
58
59
    }
    if (child[1]) {
60
      DELETE child[1];
61
    }
62

63
64
65
    if (newCoord) {
      DELETE newCoord;
    }
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88

    if (elementData) {
      elementData->deleteDecorated();
      DELETE elementData;
    }
  }

  bool Element::deleteElementData(int typeID)
  {
    FUNCNAME("Element::deleteElementData()");

    if (elementData) {
      if (elementData->isOfType(typeID)) {
	ElementData *tmp = elementData;
	elementData = elementData->getDecorated();
	DELETE tmp;
	tmp = NULL;
	return true;
      } else {
	return elementData->deleteDecorated(typeID);
      }
    }
    return false;
89
90
  }

91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
  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;
     
      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++) {
106
107
108
	  if (dof[j]) {
	    if (deletedDOFs.count(dof[j]) == 0) {
	      deletedDOFs[dof[j]] = true;
109
110
	      FREE_MEMORY(dof[j], DegreeOfFreedom, ndof);
	    }
111
	  }  
112
113
114
115
116
117
	  j++;
	}
      }
    }

    FREE_MEMORY(dof, DegreeOfFreedom*, mesh->getNumberOfNodes());
118
    
119
120
121
122
123
124
125
126
    if (child[0]) {
      child[0]->deleteElementDOFs();
    }
    if (child[1]) {
      child[1]->deleteElementDOFs();
    }
  }

127
128
129
  Element* Element::cloneWithDOFs()
  {
    Element *el;
130
    
131
    if (isLine()) {
132
      el = NEW Line(NULL);
133
    } else if (isTriangle()) {
134
      el = NEW Triangle(NULL);
135
    } else {
136
      el = NEW Tetrahedron(NULL);
137
    }
138
139

    el->mesh = mesh;
140
141
    el->index = index;
    el->mark = mark;
142
    el->hidden = hidden;
143
144
145
146
147
    if (newCoord) {
      WorldVector<double> *nc = NEW WorldVector<double>();
      *nc = *newCoord;
      el->newCoord = nc;
    }
148
    
149
    /* =========== And here we clone the DOFs =========== */
150
   
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
    el->dof = GET_MEMORY(DegreeOfFreedom*, mesh->getNumberOfNodes());

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

    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) {
	    if (Mesh::serializedDOFs[dof[j][0]] == NULL) {
	      el->dof[j] = GET_MEMORY(DegreeOfFreedom, ndof);
	      for (int k = 0; k < ndof; k++) {
		el->dof[j][k] = dof[j][k];
	      }
	      Mesh::serializedDOFs[dof[j][0]] = el->dof[j];
	    } else {
	      el->dof[j] = Mesh::serializedDOFs[dof[j][0]];
	    }
	  } else {
	    el->dof[j] = NULL;
	  }
	  j++;
	}
      }
    }
183
    
184
185
186
    /* =========== And clone the children ============= */

    if (child[0]) {
187
      el->child[0] = child[0]->cloneWithDOFs();
188
189
    }
    if (child[1]) {
190
      el->child[1] = child[1]->cloneWithDOFs();
191
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)
  {
269
270
    int i, j, k, n0, nd, nd0;
    DegreeOfFreedom  *ldof;
271
272
273
274
275
276
277
278
279
280
281
282
    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);
      for (i = 0; i < vertices; i++) {
	CHANGE_DOFS_2(this);
      }
    }

283
    if (mesh->getDim() > 1) {
284
285
286
287
288
289
290
291
292
      if ((nd = admin->getNumberOfDOFs(EDGE)))  {
	nd0 = admin->getNumberOfPreDOFs(EDGE);
	n0 = admin->getMesh()->getNode(EDGE);
	for (i = 0; i < edges; i++) {
	  CHANGE_DOFS_2(this);
	}
      }
    }

293
    if (mesh->getDim() == 3) {
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
      if ((nd = admin->getNumberOfDOFs(FACE)))  {
	nd0 = admin->getNumberOfPreDOFs(FACE);
	n0 = admin->getMesh()->getNode(FACE);
	for (i = 0; i < faces; i++) {
	  CHANGE_DOFS_2(this);
	}
      }
    }

    if ((nd = admin->getNumberOfDOFs(CENTER)))  {
      nd0 = admin->getNumberOfPreDOFs(CENTER);
      n0 = admin->getMesh()->getNode(CENTER);
      i = 0;          /* only one center */
      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
  {
321
322
    int nv = 0;
    int ov = 0;
323
324
325
    int vertices = mesh->getGeo(VERTEX);
    int dim = mesh->getDim();

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

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

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

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

    // write dofs
    int dim = mesh->getDim();
    int nodes = mesh->getNumberOfNodes();
Thomas Witkowski's avatar
Thomas Witkowski committed
387
    int j = 0;
388
    out.write(reinterpret_cast<const char*>(&nodes), sizeof(int));
Thomas Witkowski's avatar
Thomas Witkowski committed
389
390
   
    for (int pos = 0; pos <= dim; pos++) {
391
      GeoIndex position = INDEX_OF_DIM(pos, dim);
Thomas Witkowski's avatar
Thomas Witkowski committed
392
      int ndof = 0;
393

Thomas Witkowski's avatar
Thomas Witkowski committed
394
      for (int i = 0; i < mesh->getNumberOfDOFAdmin(); i++) {
395
	ndof += mesh->getDOFAdmin(i).getNumberOfDOFs(position);	
396
397
      }

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

    // write index
    out.write(reinterpret_cast<const char*>(&index), sizeof(int));

    // write mark
    out.write(reinterpret_cast<const char*>(&mark), sizeof(signed char));

425
426
427
    // write hidden flag
    out.write(reinterpret_cast<const char*>(&hidden), sizeof(hidden));

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
445
  void Element::deserialize(std::istream &in)
446
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
447
    std::string typeName;
448
449
450
451
452

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

Thomas Witkowski's avatar
Thomas Witkowski committed
453
454
    if (typeName != "NULL") {
      if (typeName == "Line") {
455
456
457
	child[0] = NEW Line(NULL);
	child[1] = NEW Line(NULL);      
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
458
      if (typeName == "Triangle") {
459
460
461
	child[0] = NEW Triangle(NULL);
	child[1] = NEW Triangle(NULL);      
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
462
      if (typeName == "Tetrahedron") {
463
464
465
	child[0] = NEW Tetrahedron(NULL);
	child[1] = NEW Tetrahedron(NULL);      
      }
466
467
468
469
470
471
472
473
474
475
476
477
      child[0]->deserialize(in);
      child[1]->deserialize(in);
    } else {
      child[0] = child[1] = NULL;
    }

    // read dofs
    int nodes;
    in.read(reinterpret_cast<char*>(&nodes), sizeof(int));

    dof = GET_MEMORY(DegreeOfFreedom*, nodes); 

478
    for (int i = 0; i < nodes; i++) {
479
480
481
      int dofs;
      in.read(reinterpret_cast<char*>(&dofs), sizeof(int));

482
483
      if (dofs) {
	if (dofs != -1) {
484
485
	  dof[i] = GET_MEMORY(DegreeOfFreedom, dofs);
	  in.read(reinterpret_cast<char*>(dof[i]), dofs * sizeof(DegreeOfFreedom));
486
	  if (Mesh::serializedDOFs[dof[i][0]] != NULL) {
487
488
489
	    DegreeOfFreedom *dofPtr = Mesh::serializedDOFs[dof[i][0]];
	    FREE_MEMORY(dof[i], DegreeOfFreedom, dofs);
	    dof[i] = dofPtr;
490
491
	  } else {
	    Mesh::serializedDOFs[dof[i][0]] = dof[i];
492
493
494
495
496
497
498
499
500
501
502
503
504
505
	  }
	} else {
	  DegreeOfFreedom index;
	  in.read(reinterpret_cast<char*>(&index), sizeof(DegreeOfFreedom));
	  dof[i] = Mesh::serializedDOFs[index];
	}
      } else {
	dof[i] = NULL;
      }
    }
    
    // read index
    in.read(reinterpret_cast<char*>(&index), sizeof(int));

506
    // read mark
507
508
    in.read(reinterpret_cast<char*>(&mark), sizeof(signed char));

509
510
511
    // read hidden
    in.read(reinterpret_cast<char*>(&hidden), sizeof(bool));

512
513
514
515
    // read newCoord
    in >> typeName;
    in.get();

516
517
    if (typeName != "NULL") {
      if (typeName == "WorldVector") {
518
519
520
521
522
523
524
525
526
527
528
529
530
	newCoord = NEW WorldVector<double>;
	newCoord->deserialize(in);
      } else {
	ERROR_EXIT("unexpected type name\n");
      }
    } else {
      newCoord = NULL;
    }

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

531
    if (typeName != "NULL") {
532
      elementData = CreatorMap<ElementData>::getCreator(typeName)->create();
533
      if (elementData) {
534
535
536
537
538
539
540
541
542
	elementData->deserialize(in);
      } else {
	ERROR_EXIT("unexpected type name\n");
      }
    } else {
      elementData = NULL;
    }
  }

543
544
545
546
547
  int Element::calcMemoryUsage()
  {
    int result = 0;

    result += sizeof(Element);
548
549
    result += mesh->getNumberOfNodes() * sizeof(DegreeOfFreedom*);

550
551
552
553
554
555
556
    if (child[0]) {
      result += child[0]->calcMemoryUsage() + child[1]->calcMemoryUsage();
    }

    return result;
  }

557
}