Element.cc 12.9 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
#include "Element.h"
#include "DOFAdmin.h"
#include "Mesh.h"
#include "CoarseningManager.h"
#include "FixVec.h"
#include "ElementRegion_ED.h"

namespace AMDiS {


  int Element::getRegion() const 
  {
    ElementRegion_ED* red_;

15
16
17
18
19
    if (!elementData) 
      return -1;

    red_ = dynamic_cast<ElementRegion_ED*>(elementData->getElementData(ELEMENT_REGION));

20
21
22
23
24
25
26
27
    if (red_) 
      return red_->getRegion();
    
    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
39
40
    dof = mesh->createDOFPtrs();  
  }

  Element::Element(Mesh *aMesh)
  {
    mesh = aMesh;

    index = mesh ? mesh->getNextElementIndex() : -1; 
  
41
42
    child[0] = NULL;
    child[1] = NULL;
43
44
45

    newCoord = NULL;
    elementData = NULL;
46
47
48
49
50
51

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

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

64
65
66
    if (newCoord) {
      DELETE newCoord;
    }
67
68
  }

69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
  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++) {
	  if (dof[j] != NULL) {
	    if (Mesh::serializedDOFs.count(dof[j][0]) == 0) {
	      //	      std::cout << "FREE INNER: " << ndof << std::endl;
	      FREE_MEMORY(dof[j], DegreeOfFreedom, ndof);
	      Mesh::serializedDOFs[dof[j][0]] = NULL;
	    }
	  }
	  
	  j++;
	}
      }
    }

    //    std::cout << "FREE OUTER: " << mesh->getNumberOfNodes() << std::endl;
    FREE_MEMORY(dof, DegreeOfFreedom*, mesh->getNumberOfNodes());

    if (child[0]) {
      child[0]->deleteElementDOFs();
    }
    if (child[1]) {
      child[1]->deleteElementDOFs();
    }
  }

108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
  Element* Element::cloneWithDOFs()
  {
    Element *el;

    if (isLine()) {
      el = NEW Line(mesh);
    } else if (isTriangle()) {
      el = NEW Triangle(mesh);
    } else {
      el = NEW Tetrahedron(mesh);
    }
    el->index = index;
    el->mark = mark;
    if (newCoord) {
      WorldVector<double> *nc = NEW WorldVector<double>();
      *nc = *newCoord;
      el->newCoord = nc;
    }

    /* =========== And here we clone the DOFs =========== */

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

    /* =========== And clone the children ============= */

    if (child[0]) {
165
      el->child[0] = child[0]->cloneWithDOFs();
166
167
    }
    if (child[1]) {
168
      el->child[1] = child[1]->cloneWithDOFs();
169
170
171
172
173
    }

    return el;
  }

174
175
176
177
178
179
180
181
182
  /****************************************************************************/
  /*  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)						\
183
  ldof = el->dof[n0 + i] + nd0;						\
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
    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)
  {
202
203
    int j, k, n0, nd, nd0;
    DegreeOfFreedom *ldof;  
204
205
206
207
208
209
210
    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);
211
      for (int i = 0; i < vertices; i++) {
212
213
214
215
	CHANGE_DOFS_1(this);
      }
    }

216
    if (mesh->getDim() > 1) {
217
218
219
      if ((nd = admin->getNumberOfDOFs(EDGE)))  {
	nd0 = admin->getNumberOfPreDOFs(EDGE);
	n0 = admin->getMesh()->getNode(EDGE);
220
	for (int i = 0; i < edges; i++) {
221
222
223
224
225
	  CHANGE_DOFS_1(this);
	}
      }
    }

226
    if (mesh->getDim() == 3) {
227
228
229
      if ((nd = admin->getNumberOfDOFs(FACE)))  {
	nd0 = admin->getNumberOfPreDOFs(FACE);
	n0 = admin->getMesh()->getNode(FACE);
230
	for (int i = 0; i < faces; i++) {
231
232
233
234
235
236
237
238
	  CHANGE_DOFS_1(this);
	}
      }
    }

    if ((nd = admin->getNumberOfDOFs(CENTER)))  {
      nd0 = admin->getNumberOfPreDOFs(CENTER);
      n0 = admin->getMesh()->getNode(CENTER);
239
      int i = 0;          /* only one center */
240
241
242
243
244
245
246
      CHANGE_DOFS_1(this);
    }
  }


  void Element::newDOFFct2(const DOFAdmin* admin)
  {
247
248
    int i, j, k, n0, nd, nd0;
    DegreeOfFreedom  *ldof;
249
250
251
252
253
254
255
256
257
258
259
260
    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);
      }
    }

261
    if (mesh->getDim() > 1) {
262
263
264
265
266
267
268
269
270
      if ((nd = admin->getNumberOfDOFs(EDGE)))  {
	nd0 = admin->getNumberOfPreDOFs(EDGE);
	n0 = admin->getMesh()->getNode(EDGE);
	for (i = 0; i < edges; i++) {
	  CHANGE_DOFS_2(this);
	}
      }
    }

271
    if (mesh->getDim() == 3) {
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
      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
  {
299
300
    int nv = 0;
    int ov = 0;
301
302
303
    int vertices = mesh->getGeo(VERTEX);
    int dim = mesh->getDim();

304
    for (int i = 0; i < vertices; i++) {
305
      if (nv < i - 1)  
306
307
308
309
310
311
312
313
314
315
316
	return(-1);

      for (int j = 0; j < dim; j++) {
	if (dof[i] == pdof[j]) {
	  /****************************************************************************/
	  /* i is a common vertex                                                     */
	  /****************************************************************************/
	  ov += i;
	  nv++;
	  break;
	}
317
      }
318
319
320
321
322

    }
    
    if (nv != mesh->getDim()) 
      return(-1);
323
324
325
326
327
328
329
330
331
332
    /****************************************************************************/
    /*  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:
333
      return 3 - ov;
334
335
      break;
    case 3:
336
      return 6 - ov;
337
338
339
340
341
342
343
344
      break;
    default:
      ERROR_EXIT("invalid dim\n");
      return 0;
    }
  }

  double Element::getNewCoord(int j) const {
345
    if (j >= 0) {
346
      TEST_EXIT_DBG(newCoord)("newCoord = NULL\n");
347
348
      return (*newCoord)[j];
    } else { 
349
      return (newCoord != NULL);
350
351
352
353
    }
  }

  void Element::eraseNewCoord() {
354
    if (newCoord != NULL) {
355
      DELETE newCoord;
356
357
      newCoord = NULL;
    }
358
  }
359
 
Thomas Witkowski's avatar
Thomas Witkowski committed
360
  void Element::serialize(std::ostream &out) 
361
362
  {
    // write children
363
    if (child[0]) {
Thomas Witkowski's avatar
Thomas Witkowski committed
364
      out << child[0]->getTypeName() << "\n";
365
366
367
      child[0]->serialize(out);
      child[1]->serialize(out);
    } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
368
      out << "NULL\n";
369
370
371
372
373
    }

    // write dofs
    int dim = mesh->getDim();
    int nodes = mesh->getNumberOfNodes();
Thomas Witkowski's avatar
Thomas Witkowski committed
374
    int j = 0;
375
    out.write(reinterpret_cast<const char*>(&nodes), sizeof(int));
Thomas Witkowski's avatar
Thomas Witkowski committed
376
377
   
    for (int pos = 0; pos <= dim; pos++) {
378
      GeoIndex position = INDEX_OF_DIM(pos, dim);
Thomas Witkowski's avatar
Thomas Witkowski committed
379
      int ndof = 0;
380

Thomas Witkowski's avatar
Thomas Witkowski committed
381
      for (int i = 0; i < mesh->getNumberOfDOFAdmin(); i++) {
382
	ndof += mesh->getDOFAdmin(i).getNumberOfDOFs(position);	
383
384
      }

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

    // write element data
Thomas Witkowski's avatar
Thomas Witkowski committed
421
422
    if (elementData) {
      out << elementData->getTypeName() << "\n";
423
424
      elementData->serialize(out);
    } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
425
      out << "NULL\n";
426
427
428
    }
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
429
  void Element::deserialize(std::istream &in)
430
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
431
    std::string typeName;
432
433
434
435
436

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

Thomas Witkowski's avatar
Thomas Witkowski committed
437
438
    if (typeName != "NULL") {
      if (typeName == "Line") {
439
440
441
	child[0] = NEW Line(NULL);
	child[1] = NEW Line(NULL);      
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
442
      if (typeName == "Triangle") {
443
444
445
	child[0] = NEW Triangle(NULL);
	child[1] = NEW Triangle(NULL);      
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
446
      if (typeName == "Tetrahedron") {
447
448
449
	child[0] = NEW Tetrahedron(NULL);
	child[1] = NEW Tetrahedron(NULL);      
      }
450
451
452
453
454
455
456
457
458
459
460
461
      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); 

462
    for (int i = 0; i < nodes; i++) {
463
464
465
      int dofs;
      in.read(reinterpret_cast<char*>(&dofs), sizeof(int));

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

497
498
    if (typeName != "NULL") {
      if (typeName == "WorldVector") {
499
500
501
502
503
504
505
506
507
508
509
510
511
	newCoord = NEW WorldVector<double>;
	newCoord->deserialize(in);
      } else {
	ERROR_EXIT("unexpected type name\n");
      }
    } else {
      newCoord = NULL;
    }

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

512
    if (typeName != "NULL") {
513
      elementData = CreatorMap<ElementData>::getCreator(typeName)->create();
514
      if (elementData) {
515
516
517
518
519
520
521
522
523
	elementData->deserialize(in);
      } else {
	ERROR_EXIT("unexpected type name\n");
      }
    } else {
      elementData = NULL;
    }
  }

524
525
526
527
528
529
530
531
532
533
534
535
  int Element::calcMemoryUsage()
  {
    int result = 0;

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

    return result;
  }

536
}