Element.cc 11.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
56
  }

  // call destructor through Mesh::freeElement !!!
  Element::~Element()
  {
57
58
59
60
    if (child[0]) 
      DELETE child[0];
    if (child[1]) 
      DELETE child[1];
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
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
  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]) {
124
      el->child[0] = child[0]->cloneWithDOFs();
125
126
    }
    if (child[1]) {
127
      el->child[1] = child[1]->cloneWithDOFs();
128
129
130
131
132
    }

    return el;
  }

133
134
135
136
137
138
139
140
141
  /****************************************************************************/
  /*  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)						\
142
  ldof = el->dof[n0 + i] + nd0;						\
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
    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)
  {
161
162
    int j, k, n0, nd, nd0;
    DegreeOfFreedom *ldof;  
163
164
165
166
167
168
169
    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);
170
      for (int i = 0; i < vertices; i++) {
171
172
173
174
	CHANGE_DOFS_1(this);
      }
    }

175
    if (mesh->getDim() > 1) {
176
177
178
      if ((nd = admin->getNumberOfDOFs(EDGE)))  {
	nd0 = admin->getNumberOfPreDOFs(EDGE);
	n0 = admin->getMesh()->getNode(EDGE);
179
	for (int i = 0; i < edges; i++) {
180
181
182
183
184
	  CHANGE_DOFS_1(this);
	}
      }
    }

185
    if (mesh->getDim() == 3) {
186
187
188
      if ((nd = admin->getNumberOfDOFs(FACE)))  {
	nd0 = admin->getNumberOfPreDOFs(FACE);
	n0 = admin->getMesh()->getNode(FACE);
189
	for (int i = 0; i < faces; i++) {
190
191
192
193
194
195
196
197
	  CHANGE_DOFS_1(this);
	}
      }
    }

    if ((nd = admin->getNumberOfDOFs(CENTER)))  {
      nd0 = admin->getNumberOfPreDOFs(CENTER);
      n0 = admin->getMesh()->getNode(CENTER);
198
      int i = 0;          /* only one center */
199
200
201
202
203
204
205
      CHANGE_DOFS_1(this);
    }
  }


  void Element::newDOFFct2(const DOFAdmin* admin)
  {
206
207
    int i, j, k, n0, nd, nd0;
    DegreeOfFreedom  *ldof;
208
209
210
211
212
213
214
215
216
217
218
219
    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);
      }
    }

220
    if (mesh->getDim() > 1) {
221
222
223
224
225
226
227
228
229
      if ((nd = admin->getNumberOfDOFs(EDGE)))  {
	nd0 = admin->getNumberOfPreDOFs(EDGE);
	n0 = admin->getMesh()->getNode(EDGE);
	for (i = 0; i < edges; i++) {
	  CHANGE_DOFS_2(this);
	}
      }
    }

230
    if (mesh->getDim() == 3) {
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
      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
  {
258
259
    int nv = 0;
    int ov = 0;
260
261
262
    int vertices = mesh->getGeo(VERTEX);
    int dim = mesh->getDim();

263
    for (int i = 0; i < vertices; i++) {
264
      if (nv < i - 1)  
265
266
267
268
269
270
271
272
273
274
275
	return(-1);

      for (int j = 0; j < dim; j++) {
	if (dof[i] == pdof[j]) {
	  /****************************************************************************/
	  /* i is a common vertex                                                     */
	  /****************************************************************************/
	  ov += i;
	  nv++;
	  break;
	}
276
      }
277
278
279
280
281

    }
    
    if (nv != mesh->getDim()) 
      return(-1);
282
283
284
285
286
287
288
289
290
291
    /****************************************************************************/
    /*  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:
292
      return 3 - ov;
293
294
      break;
    case 3:
295
      return 6 - ov;
296
297
298
299
300
301
302
303
      break;
    default:
      ERROR_EXIT("invalid dim\n");
      return 0;
    }
  }

  double Element::getNewCoord(int j) const {
304
    if (j >= 0) {
305
      TEST_EXIT_DBG(newCoord)("newCoord = NULL\n");
306
307
      return (*newCoord)[j];
    } else { 
308
      return (newCoord != NULL);
309
310
311
312
    }
  }

  void Element::eraseNewCoord() {
313
    if (newCoord != NULL) {
314
      DELETE newCoord;
315
316
      newCoord = NULL;
    }
317
  }
318
 
Thomas Witkowski's avatar
Thomas Witkowski committed
319
  void Element::serialize(std::ostream &out) 
320
321
  {
    // write children
322
    if (child[0]) {
Thomas Witkowski's avatar
Thomas Witkowski committed
323
      out << child[0]->getTypeName() << "\n";
324
325
326
      child[0]->serialize(out);
      child[1]->serialize(out);
    } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
327
      out << "NULL\n";
328
329
330
331
332
    }

    // write dofs
    int dim = mesh->getDim();
    int nodes = mesh->getNumberOfNodes();
Thomas Witkowski's avatar
Thomas Witkowski committed
333
    int j = 0;
334
    out.write(reinterpret_cast<const char*>(&nodes), sizeof(int));
Thomas Witkowski's avatar
Thomas Witkowski committed
335
336
   
    for (int pos = 0; pos <= dim; pos++) {
337
      GeoIndex position = INDEX_OF_DIM(pos, dim);
Thomas Witkowski's avatar
Thomas Witkowski committed
338
      int ndof = 0;
339

Thomas Witkowski's avatar
Thomas Witkowski committed
340
      for (int i = 0; i < mesh->getNumberOfDOFAdmin(); i++) {
341
	ndof += mesh->getDOFAdmin(i).getNumberOfDOFs(position);	
342
343
      }

Thomas Witkowski's avatar
Thomas Witkowski committed
344
345
346
347
      if (ndof > 0) {
	for (int i = 0; i < mesh->getGeo(position); i++) {
	  if (dof[j] != NULL) {
	    if (Mesh::serializedDOFs[dof[j][0]] == NULL) {
348
349
	      Mesh::serializedDOFs[dof[j][0]] = dof[j];
	      out.write(reinterpret_cast<const char*>(&ndof), sizeof(int));
350
	      out.write(reinterpret_cast<const char*>(dof[j]), ndof * sizeof(DegreeOfFreedom));
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
	    } 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
372
373
    if (newCoord) {
      out << "WorldVector\n";
374
375
      newCoord->serialize(out);
    } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
376
      out << "NULL\n";
377
378
379
    }

    // write element data
Thomas Witkowski's avatar
Thomas Witkowski committed
380
381
    if (elementData) {
      out << elementData->getTypeName() << "\n";
382
383
      elementData->serialize(out);
    } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
384
      out << "NULL\n";
385
386
387
    }
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
388
  void Element::deserialize(std::istream &in)
389
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
390
    std::string typeName;
391
392
393
394
395

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

Thomas Witkowski's avatar
Thomas Witkowski committed
396
397
    if (typeName != "NULL") {
      if (typeName == "Line") {
398
399
400
	child[0] = NEW Line(NULL);
	child[1] = NEW Line(NULL);      
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
401
      if (typeName == "Triangle") {
402
403
404
	child[0] = NEW Triangle(NULL);
	child[1] = NEW Triangle(NULL);      
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
405
      if (typeName == "Tetrahedron") {
406
407
408
	child[0] = NEW Tetrahedron(NULL);
	child[1] = NEW Tetrahedron(NULL);      
      }
409
410
411
412
413
414
415
416
417
418
419
420
      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); 

421
    for (int i = 0; i < nodes; i++) {
422
423
424
      int dofs;
      in.read(reinterpret_cast<char*>(&dofs), sizeof(int));

425
426
      if (dofs) {
	if (dofs != -1) {
427
428
	  dof[i] = GET_MEMORY(DegreeOfFreedom, dofs);
	  in.read(reinterpret_cast<char*>(dof[i]), dofs * sizeof(DegreeOfFreedom));
429
	  if (Mesh::serializedDOFs[dof[i][0]] != NULL) {
430
431
432
	    DegreeOfFreedom *dofPtr = Mesh::serializedDOFs[dof[i][0]];
	    FREE_MEMORY(dof[i], DegreeOfFreedom, dofs);
	    dof[i] = dofPtr;
433
434
	  } else {
	    Mesh::serializedDOFs[dof[i][0]] = dof[i];
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
	  }
	} 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();

456
457
    if (typeName != "NULL") {
      if (typeName == "WorldVector") {
458
459
460
461
462
463
464
465
466
467
468
469
470
	newCoord = NEW WorldVector<double>;
	newCoord->deserialize(in);
      } else {
	ERROR_EXIT("unexpected type name\n");
      }
    } else {
      newCoord = NULL;
    }

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

471
    if (typeName != "NULL") {
472
      elementData = CreatorMap<ElementData>::getCreator(typeName)->create();
473
      if (elementData) {
474
475
476
477
478
479
480
481
482
	elementData->deserialize(in);
      } else {
	ERROR_EXIT("unexpected type name\n");
      }
    } else {
      elementData = NULL;
    }
  }

483
484
485
486
487
488
489
490
491
492
493
494
  int Element::calcMemoryUsage()
  {
    int result = 0;

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

    return result;
  }

495
}