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
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]) {
56
      DELETE child[0];
57
58
    }
    if (child[1]) {
59
      DELETE child[1];
60
    }
61

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

    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;
88
89
  }

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

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

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

    el->mesh = mesh;
139
140
141
142
143
144
145
    el->index = index;
    el->mark = mark;
    if (newCoord) {
      WorldVector<double> *nc = NEW WorldVector<double>();
      *nc = *newCoord;
      el->newCoord = nc;
    }
146
    
147
    /* =========== And here we clone the DOFs =========== */
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
179
180
    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++;
	}
      }
    }
181
    
182
183
184
    /* =========== And clone the children ============= */

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

    return el;
  }

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

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

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

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


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

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

291
    if (mesh->getDim() == 3) {
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
317
318
      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
  {
319
320
    int nv = 0;
    int ov = 0;
321
322
323
    int vertices = mesh->getGeo(VERTEX);
    int dim = mesh->getDim();

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

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

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

  double Element::getNewCoord(int j) const {
365
    if (j >= 0) {
366
      TEST_EXIT_DBG(newCoord)("newCoord = NULL\n");
367
368
      return (*newCoord)[j];
    } else { 
369
      return (newCoord != NULL);
370
371
372
373
    }
  }

  void Element::eraseNewCoord() {
374
    if (newCoord != NULL) {
375
      DELETE newCoord;
376
377
      newCoord = NULL;
    }
378
  }
379
 
Thomas Witkowski's avatar
Thomas Witkowski committed
380
  void Element::serialize(std::ostream &out) 
381
382
  {
    // write children
383
    if (child[0]) {
Thomas Witkowski's avatar
Thomas Witkowski committed
384
      out << child[0]->getTypeName() << "\n";
385
386
387
      child[0]->serialize(out);
      child[1]->serialize(out);
    } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
388
      out << "NULL\n";
389
390
391
392
393
    }

    // write dofs
    int dim = mesh->getDim();
    int nodes = mesh->getNumberOfNodes();
Thomas Witkowski's avatar
Thomas Witkowski committed
394
    int j = 0;
395
    out.write(reinterpret_cast<const char*>(&nodes), sizeof(int));
Thomas Witkowski's avatar
Thomas Witkowski committed
396
397
   
    for (int pos = 0; pos <= dim; pos++) {
398
      GeoIndex position = INDEX_OF_DIM(pos, dim);
Thomas Witkowski's avatar
Thomas Witkowski committed
399
      int ndof = 0;
400

Thomas Witkowski's avatar
Thomas Witkowski committed
401
      for (int i = 0; i < mesh->getNumberOfDOFAdmin(); i++) {
402
	ndof += mesh->getDOFAdmin(i).getNumberOfDOFs(position);	
403
404
      }

Thomas Witkowski's avatar
Thomas Witkowski committed
405
406
407
408
      if (ndof > 0) {
	for (int i = 0; i < mesh->getGeo(position); i++) {
	  if (dof[j] != NULL) {
	    if (Mesh::serializedDOFs[dof[j][0]] == NULL) {
409
410
	      Mesh::serializedDOFs[dof[j][0]] = dof[j];
	      out.write(reinterpret_cast<const char*>(&ndof), sizeof(int));
411
	      out.write(reinterpret_cast<const char*>(dof[j]), ndof * sizeof(DegreeOfFreedom));
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
	    } 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
433
434
    if (newCoord) {
      out << "WorldVector\n";
435
436
      newCoord->serialize(out);
    } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
437
      out << "NULL\n";
438
439
440
    }

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

Thomas Witkowski's avatar
Thomas Witkowski committed
449
  void Element::deserialize(std::istream &in)
450
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
451
    std::string typeName;
452
453
454
455
456

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

Thomas Witkowski's avatar
Thomas Witkowski committed
457
458
    if (typeName != "NULL") {
      if (typeName == "Line") {
459
460
461
	child[0] = NEW Line(NULL);
	child[1] = NEW Line(NULL);      
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
462
      if (typeName == "Triangle") {
463
464
465
	child[0] = NEW Triangle(NULL);
	child[1] = NEW Triangle(NULL);      
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
466
      if (typeName == "Tetrahedron") {
467
468
469
	child[0] = NEW Tetrahedron(NULL);
	child[1] = NEW Tetrahedron(NULL);      
      }
470
471
472
473
474
475
476
477
478
479
480
481
      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); 

482
    for (int i = 0; i < nodes; i++) {
483
484
485
      int dofs;
      in.read(reinterpret_cast<char*>(&dofs), sizeof(int));

486
487
      if (dofs) {
	if (dofs != -1) {
488
489
	  dof[i] = GET_MEMORY(DegreeOfFreedom, dofs);
	  in.read(reinterpret_cast<char*>(dof[i]), dofs * sizeof(DegreeOfFreedom));
490
	  if (Mesh::serializedDOFs[dof[i][0]] != NULL) {
491
492
493
	    DegreeOfFreedom *dofPtr = Mesh::serializedDOFs[dof[i][0]];
	    FREE_MEMORY(dof[i], DegreeOfFreedom, dofs);
	    dof[i] = dofPtr;
494
495
	  } else {
	    Mesh::serializedDOFs[dof[i][0]] = dof[i];
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
	  }
	} 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));

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

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

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

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

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

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

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

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

    return result;
  }

558
}