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

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

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

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

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

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

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

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

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

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

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

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

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

508
509
    if (typeName != "NULL") {
      if (typeName == "WorldVector") {
510
511
512
513
514
515
516
517
518
519
520
521
522
	newCoord = NEW WorldVector<double>;
	newCoord->deserialize(in);
      } else {
	ERROR_EXIT("unexpected type name\n");
      }
    } else {
      newCoord = NULL;
    }

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

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

535
536
537
538
539
  int Element::calcMemoryUsage()
  {
    int result = 0;

    result += sizeof(Element);
540
541
    result += mesh->getNumberOfNodes() * sizeof(DegreeOfFreedom*);

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

    return result;
  }

549
}