Element.cc 13.3 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
44
45
46
47
48

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

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

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

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

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

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

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

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

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

    return el;
  }

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

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

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

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


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

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

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

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

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

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

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

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
404
405
406
407
      if (ndof > 0) {
	for (int i = 0; i < mesh->getGeo(position); i++) {
	  if (dof[j] != NULL) {
	    if (Mesh::serializedDOFs[dof[j][0]] == NULL) {
408
409
	      Mesh::serializedDOFs[dof[j][0]] = dof[j];
	      out.write(reinterpret_cast<const char*>(&ndof), sizeof(int));
410
	      out.write(reinterpret_cast<const char*>(dof[j]), ndof * sizeof(DegreeOfFreedom));
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
	    } 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
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
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
450
    std::string typeName;
451
452
453
454
455

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

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

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

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

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
548
549
550
551
552
553
554
  int Element::calcMemoryUsage()
  {
    int result = 0;

    result += sizeof(Element);
    if (child[0]) {
      result += child[0]->calcMemoryUsage() + child[1]->calcMemoryUsage();
    }

    return result;
  }

555
}