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
124
125
126
127
128
129
130
131
132
  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]) {
      el->child[0] = child[0]->clone();
    }
    if (child[1]) {
      el->child[1] = child[1]->clone();
    }

    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
}