Element.cc 13.2 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
18
    ElementRegion_ED* red = 
      dynamic_cast<ElementRegion_ED*>(elementData->getElementData(ELEMENT_REGION));
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
45
46
47
48
49

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

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

61
    if (newCoord)
Thomas Witkowski's avatar
Thomas Witkowski committed
62
      delete newCoord;    
63
64
65

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

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

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

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

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

    el->mesh = mesh;
137
138
139
    el->index = index;
    el->mark = mark;
    if (newCoord) {
Thomas Witkowski's avatar
Thomas Witkowski committed
140
      WorldVector<double> *nc = new WorldVector<double>();
141
142
143
      *nc = *newCoord;
      el->newCoord = nc;
    }
144
    
145
    /* =========== And here we clone the DOFs =========== */
146
   
147
148
149
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
    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++;
	}
      }
    }
179
    
180
181
182
    /* =========== And clone the children ============= */

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

    return el;
  }

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

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

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

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


  void Element::newDOFFct2(const DOFAdmin* admin)
  {
265
266
    int i, j, k, n0, nd, nd0;
    DegreeOfFreedom  *ldof;
267
268
269
270
271
272
273
274
275
276
277
278
    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);
      }
    }

279
    if (mesh->getDim() > 1) {
280
281
282
283
284
285
286
287
288
      if ((nd = admin->getNumberOfDOFs(EDGE)))  {
	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
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
      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
  {
317
318
    int nv = 0;
    int ov = 0;
319
320
321
    int vertices = mesh->getGeo(VERTEX);
    int dim = mesh->getDim();

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

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

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
390
      for (int i = 0; i < mesh->getNumberOfDOFAdmin(); i++) {
391
	ndof += mesh->getDOFAdmin(i).getNumberOfDOFs(position);	
392
393
      }

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

    // write newCoord
Thomas Witkowski's avatar
Thomas Witkowski committed
422
423
    if (newCoord) {
      out << "WorldVector\n";
424
425
      newCoord->serialize(out);
    } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
426
      out << "NULL\n";
427
428
429
    }

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

Thomas Witkowski's avatar
Thomas Witkowski committed
438
  void Element::deserialize(std::istream &in)
439
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
440
    std::string typeName;
441
442
443
444
445

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

Thomas Witkowski's avatar
Thomas Witkowski committed
446
447
    if (typeName != "NULL") {
      if (typeName == "Line") {
Thomas Witkowski's avatar
Thomas Witkowski committed
448
449
	child[0] = new Line(NULL);
	child[1] = new Line(NULL);      
450
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
451
      if (typeName == "Triangle") {
Thomas Witkowski's avatar
Thomas Witkowski committed
452
453
	child[0] = new Triangle(NULL);
	child[1] = new Triangle(NULL);      
454
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
455
      if (typeName == "Tetrahedron") {
Thomas Witkowski's avatar
Thomas Witkowski committed
456
457
	child[0] = new Tetrahedron(NULL);
	child[1] = new Tetrahedron(NULL);      
458
      }
459
460
461
462
463
464
465
466
467
468
469
470
      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); 

471
    for (int i = 0; i < nodes; i++) {
472
473
474
      int dofs;
      in.read(reinterpret_cast<char*>(&dofs), sizeof(int));

475
476
      if (dofs) {
	if (dofs != -1) {
477
478
	  dof[i] = GET_MEMORY(DegreeOfFreedom, dofs);
	  in.read(reinterpret_cast<char*>(dof[i]), dofs * sizeof(DegreeOfFreedom));
479
	  if (Mesh::serializedDOFs[dof[i][0]] != NULL) {
480
481
482
	    DegreeOfFreedom *dofPtr = Mesh::serializedDOFs[dof[i][0]];
	    FREE_MEMORY(dof[i], DegreeOfFreedom, dofs);
	    dof[i] = dofPtr;
483
484
	  } else {
	    Mesh::serializedDOFs[dof[i][0]] = dof[i];
485
486
487
488
489
490
491
492
493
494
495
496
497
498
	  }
	} 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));

499
    // read mark
500
501
502
503
504
505
    in.read(reinterpret_cast<char*>(&mark), sizeof(signed char));

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

506
507
    if (typeName != "NULL") {
      if (typeName == "WorldVector") {
Thomas Witkowski's avatar
Thomas Witkowski committed
508
	newCoord = new WorldVector<double>;
509
510
511
512
513
514
515
516
517
518
519
520
	newCoord->deserialize(in);
      } else {
	ERROR_EXIT("unexpected type name\n");
      }
    } else {
      newCoord = NULL;
    }

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

521
    if (typeName != "NULL") {
522
      elementData = CreatorMap<ElementData>::getCreator(typeName)->create();
523
      if (elementData) {
524
525
526
527
528
529
530
531
532
	elementData->deserialize(in);
      } else {
	ERROR_EXIT("unexpected type name\n");
      }
    } else {
      elementData = NULL;
    }
  }

533
534
535
536
537
  int Element::calcMemoryUsage()
  {
    int result = 0;

    result += sizeof(Element);
538
539
    result += mesh->getNumberOfNodes() * sizeof(DegreeOfFreedom*);

540
541
542
543
544
545
546
    if (child[0]) {
      result += child[0]->calcMemoryUsage() + child[1]->calcMemoryUsage();
    }

    return result;
  }

547
}