Mesh.cc 34.6 KB
Newer Older
Thomas Witkowski's avatar
Thomas Witkowski committed
1
2
3
4
5
6
#include <algorithm>
#include <set>
#include <map>

#include "time.h"

7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
#include "AdaptStationary.h"
#include "AdaptInstationary.h"
#include "FiniteElemSpace.h"
#include "ElementData.h"
#include "MacroElement.h"
#include "MacroReader.h"
#include "Mesh.h"
#include "Traverse.h"
#include "Parameters.h"
#include "FixVec.h"
#include "DOFVector.h"
#include "CoarseningManager.h"
#include "DOFIterator.h"
#include "VertexVector.h"
#include "MacroWriter.h"
#include "PeriodicMap.h"
#include "Projection.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
24
#include "ElInfoStack.h"
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44

namespace AMDiS {

#define TIME_USED(f,s) ((double)((s)-(f))/(double)CLOCKS_PER_SEC)

  //**************************************************************************
  //  flags, which information should be present in the elInfo structure     
  //**************************************************************************

  const Flag Mesh::FILL_NOTHING    = 0X00L;
  const Flag Mesh::FILL_COORDS     = 0X01L;
  const Flag Mesh::FILL_BOUND      = 0X02L;
  const Flag Mesh::FILL_NEIGH      = 0X04L;
  const Flag Mesh::FILL_OPP_COORDS = 0X08L;
  const Flag Mesh::FILL_ORIENTATION= 0X10L;
  const Flag Mesh::FILL_DET        = 0X20L;
  const Flag Mesh::FILL_GRD_LAMBDA = 0X40L;
  const Flag Mesh::FILL_ADD_ALL    = 0X80L;


45
46
47
  const Flag Mesh::FILL_ANY_1D = (0X01L|0X02L|0X04L|0X08L|0x20L|0X40L|0X80L);
  const Flag Mesh::FILL_ANY_2D = (0X01L|0X02L|0X04L|0X08L|0x20L|0X40L|0X80L);
  const Flag Mesh::FILL_ANY_3D = (0X01L|0X02L|0X04L|0X08L|0X10L|0x20L|0X40L|0X80L);
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65

  //**************************************************************************
  //  flags for Mesh traversal                                                
  //**************************************************************************

  const Flag Mesh::CALL_EVERY_EL_PREORDER  = 0X0100L;
  const Flag Mesh::CALL_EVERY_EL_INORDER   = 0X0200L;
  const Flag Mesh::CALL_EVERY_EL_POSTORDER = 0X0400L;
  const Flag Mesh::CALL_LEAF_EL            = 0X0800L;
  const Flag Mesh::CALL_LEAF_EL_LEVEL      = 0X1000L;
  const Flag Mesh::CALL_EL_LEVEL           = 0X2000L;
  const Flag Mesh::CALL_MG_LEVEL           = 0X4000L ; // used in mg methods 


  // const Flag Mesh::USE_PARAMETRIC          = 0X8000L ; // used in mg methods 

  DOFAdmin* Mesh::compressAdmin = NULL;
  Mesh* Mesh::traversePtr = NULL;
Thomas Witkowski's avatar
Thomas Witkowski committed
66
  std::vector<DegreeOfFreedom> Mesh::dof_used;
67
  const int Mesh::MAX_DOF = 100;
Thomas Witkowski's avatar
Thomas Witkowski committed
68
  std::map<DegreeOfFreedom, DegreeOfFreedom*> Mesh::serializedDOFs;
69
70
71

  struct delmem { 
    DegreeOfFreedom* ptr;
72
    int len;
73
74
75
  };


Thomas Witkowski's avatar
Thomas Witkowski committed
76
  Mesh::Mesh(const std::string& aName, int dimension) 
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
    : name(aName), 
      dim(dimension), 
      nVertices(0),
      nEdges(0),
      nLeaves(0), 
      nElements(0),
      parametric(NULL), 
      preserveCoarseDOFs(false),
      nDOFEl(0),
      nDOF(dimension, DEFAULT_VALUE, 0),
      nNodeEl(0),
      node(dimension, DEFAULT_VALUE, 0),
      elementPrototype(NULL),
      elementDataPrototype(NULL),
      elementIndex(-1),
      initialized(false),
93
      macroFileInfo(NULL),
94
95
96
      final_lambda(dimension, DEFAULT_VALUE, 0.0)
  {

97
    FUNCNAME("Mesh::Mesh()");
98
99
100
101

    // set default element prototype
    switch(dim) {
    case 1:
Thomas Witkowski's avatar
Thomas Witkowski committed
102
      elementPrototype = new Line(this);
103
104
      break;
    case 2:
Thomas Witkowski's avatar
Thomas Witkowski committed
105
      elementPrototype = new Triangle(this);
106
107
      break;
    case 3:
Thomas Witkowski's avatar
Thomas Witkowski committed
108
      elementPrototype = new Tetrahedron(this);
109
110
111
112
113
114
115
      break;
    default:
      ERROR_EXIT("invalid dimension\n");
    }

    elementPrototype->setIndex(-1);

116
117
    elementIndex = 0;
  }
118
119

  Mesh::~Mesh()
120
  {
121
    Element::deletedDOFs.clear();
122
123
124
125
126

    for (std::deque<MacroElement*>::const_iterator it = macroElements.begin();
	 it != macroElements.end();
	 ++it) {
      (*it)->getElement()->deleteElementDOFs();
Thomas Witkowski's avatar
Thomas Witkowski committed
127
      delete *it;
128
129
    }    

130
    Element::deletedDOFs.clear();
131

132
    if (macroFileInfo != NULL) {
Thomas Witkowski's avatar
Thomas Witkowski committed
133
      delete macroFileInfo;
134
    }
135
    if (elementPrototype) {
Thomas Witkowski's avatar
Thomas Witkowski committed
136
      delete elementPrototype;
137
138
    }
    if (elementDataPrototype) {
Thomas Witkowski's avatar
Thomas Witkowski committed
139
      delete elementDataPrototype;
140
141
142
    }
    
    for (int i = 0; i < static_cast<int>(admin.size()); i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
143
      delete admin[i];
144
    }
145
  }
146
147

  Mesh& Mesh::operator=(const Mesh& m)
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
    FUNCNAME("Mesh::operator=()");

    if (this == &m)
      return *this;

    TEST_EXIT(dim == m.dim)("operator= works only on meshes with equal dim!\n");

    name = m.name;
    nVertices = m.nVertices;
    nEdges = m.nEdges;
    nLeaves = m.nLeaves;
    nElements = m.nElements;
    nFaces = m.nFaces;
    maxEdgeNeigh = m.maxEdgeNeigh;
    diam = m.diam;
    parametric = NULL;

    preserveCoarseDOFs = m.preserveCoarseDOFs;
    nDOFEl = m.nDOFEl;
    nDOF = m.nDOF;
    nNodeEl = m.nNodeEl;
    node = m.node;
    newDOF = m.newDOF;
    elementIndex = m.elementIndex;
    initialized = m.initialized;
    final_lambda = m.final_lambda;
    
    /* ====================== Create new DOFAdmins ================== */
    admin.resize(m.admin.size());
    for (int i = 0; i < static_cast<int>(admin.size()); i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
179
      admin[i] = new DOFAdmin(this);
180
      *(admin[i]) = *(m.admin[i]);
181
182
      admin[i]->setMesh(this);
    }
183

184

185
    /* ====================== Copy macro elements =================== */
186
  
187
188
189
190
191
192
193
194
195
    // mapIndex[i] is the index of the MacroElement element in the vector
    // macroElements, for which holds: element->getIndex() = i    
    std::map<int, int> mapIndex;

    // We use this map for coping the DOFs of the Elements within the
    // MacroElements objects.
    Mesh::serializedDOFs.clear();

    int insertCounter = 0;
196

197
198
    macroElements.clear();

199
200
201
202
203
    // Go through all MacroElements of mesh m, and create for every a new
    // MacroElement in this mesh.
    for (std::deque<MacroElement*>::const_iterator it = m.macroElements.begin();
	 it != m.macroElements.end();
	 ++it, insertCounter++) {
204

205
      // Create new MacroElement.
Thomas Witkowski's avatar
Thomas Witkowski committed
206
      MacroElement *el = new MacroElement(dim);
207

208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
      // Use copy operator to copy all the data to the new MacroElement.
      *el = **it;

      // Make a copy of the Element data, together with all DOFs
      el->setElement((*it)->getElement()->cloneWithDOFs());

      // Insert the new MacroElement in the vector of all MacroElements.
      macroElements.push_back(el);

      // Update the index map.
      mapIndex.insert(std::pair<int, int>(el->getIndex(), insertCounter));
    }

    // Now we have to go through all the new MacroElements, and update the neighbour
    // connections.
    insertCounter = 0;
    for (std::deque<MacroElement*>::const_iterator it = m.macroElements.begin();
	 it != m.macroElements.end();
	 ++it, insertCounter++) {
      // Go through all neighbours.
      for (int i = 0; i < dim; i++) {
	// 1. Get index of the old MacroElement for its i-th neighbour.
	// 2. Because the index in the new MacroElement is the same, search
	//    for the vector index the corresponding element is stored in.
	// 3. Get this element from macroElements, and set it as the i-th
	//    neighbour for the current element.
	macroElements[insertCounter]->
	  setNeighbour(i, macroElements[mapIndex[(*it)->getNeighbour(i)->getIndex()]]);
      }
    }

    // Cleanup
    Mesh::serializedDOFs.clear();

    /* ================== Things will be done when required ============ */
      
    TEST_EXIT(elementDataPrototype == NULL)("TODO\n");
    TEST_EXIT(m.parametric == NULL)("TODO\n");
    TEST_EXIT(periodicAssociations.size() == 0)("TODO\n");

    return *this;
  }

251
252
253
254
255
256
257
258
259
260
261
262
  void Mesh::updateNumberOfLeaves()
  {
    nLeaves = 0;

    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(this, -1, Mesh::CALL_LEAF_EL);
    while (elInfo) {
      nLeaves++;
      elInfo = stack.traverseNext(elInfo);
    }
  }

263
264
265
266
267
268
  void Mesh::addMacroElement(MacroElement* me) 
  {
    macroElements.push_back(me); 
    me->setIndex(macroElements.size());
  }

269
  void Mesh::removeMacroElements(std::vector<MacroElement*>& macros) 
270
271
272
273
274
  {
    FUNCNAME("Mesh::removeMacroElement()");

    TEST_EXIT(dim == 2)("Not yet implemented!\n");

275
276
277
278
279
280
281
282
283
284
285
    // Map that stores for each dof pointer (which may have a list of dofs)
    // all macro element indices that own the dof.
    std::map<const DegreeOfFreedom*, std::set<MacroElement*> > dofsOwner;
    
    // Determine all dof owner macro elements.
    for (std::deque<MacroElement*>::iterator macroIt = macroElements.begin();
	 macroIt != macroElements.end();
	 ++macroIt) {
      Element *el = (*macroIt)->getElement();
      for (int i = 0; i < 3; i++)
	dofsOwner[el->getDOF(i)].insert(*macroIt);      
286
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
287
   
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
    // Remove all the given macro elements.
    for (std::vector<MacroElement*>::iterator macroIt = macros.begin();
	 macroIt != macros.end();
	 ++macroIt) {
      bool found = false;

      // Remove the macro element from mesh's list of all macro elements.
      for (std::deque<MacroElement*>::iterator it = macroElements.begin();
	   it != macroElements.end();
	   ++it) {
	if (*it == *macroIt) {
	  macroElements.erase(it, it + 1);
	  found = true;
	  break;
	}
      }
      
      TEST_EXIT(found)("Cannot find MacroElement that should be removed!\n");
      
      // Go through all neighbours of the macro element and remove this macro element
      // to be neighbour of some other macro element.
      for (int i = 0; i <= dim; i++) {
	if ((*macroIt)->getNeighbour(i)) {
	  for (int j = 0; j <= dim; j++) {
	    if ((*macroIt)->getNeighbour(i)->getNeighbour(j) == *macroIt) {
	      (*macroIt)->getNeighbour(i)->setNeighbour(j, NULL);
	    }
315
	  }
316
317
318
319
	} else {
	  // There is no neighbour at this edge, so we have to decrease the number
	  // of edges in the mesh.
	  nEdges--;
320
	}
321
322
323
324
325
326
327
328
329
330
      }

      nLeaves--;
      nElements--;

      // Remove this macro element from the dof owner list.
      for (std::map<const DegreeOfFreedom*, std::set<MacroElement*> >::iterator dofsIt = dofsOwner.begin();
	   dofsIt != dofsOwner.end();
	   ++dofsIt) {
	std::set<MacroElement*>::iterator mIt = dofsIt->second.find(*macroIt);
Thomas Witkowski's avatar
Thomas Witkowski committed
331
	if (mIt != dofsIt->second.end()) {
332
	  dofsIt->second.erase(mIt);
Thomas Witkowski's avatar
Thomas Witkowski committed
333
	}
334
335
336
337
338
339
340
341
342
343
344
345
346
      }

      // And remove the macro element from memory
      delete *macroIt;
    }

    int nRemainDofs = 0;
    // Check now all the dofs, that have no owner anymore and therefore have to
    // be removed.
    for (std::map<const DegreeOfFreedom*, std::set<MacroElement*> >::iterator dofsIt = dofsOwner.begin();
	 dofsIt != dofsOwner.end();
	 ++dofsIt) {    
      if (dofsIt->second.size() == 0) {
Thomas Witkowski's avatar
Thomas Witkowski committed
347
	freeDOF(const_cast<DegreeOfFreedom*>(dofsIt->first), VERTEX);
348
      } else {
349
	nRemainDofs++;
350
351
352
      }
    }

353
    nVertices = nRemainDofs;
354
  }
355

356
  int Mesh::traverse(int level, Flag flag, int (*el_fct)(ElInfo*))
357
358
  {
    FUNCNAME("Mesh::traverse()");
Thomas Witkowski's avatar
Thomas Witkowski committed
359

Thomas Witkowski's avatar
Thomas Witkowski committed
360
    std::deque<MacroElement*>::iterator mel;
Thomas Witkowski's avatar
Thomas Witkowski committed
361
362
    ElInfoStack elInfoStack(this);
    ElInfo* elinfo = elInfoStack.getNextElement();
363
364
365
366
367
368
369
370
371
372
373
374
375
    Traverse tinfo(this, flag, level, el_fct);
    int sum = 0;
  
    elinfo->setFillFlag(flag);
  
    if (flag.isSet(Mesh::CALL_LEAF_EL_LEVEL) || 
	flag.isSet(Mesh::CALL_EL_LEVEL)      || 
	flag.isSet(Mesh::CALL_MG_LEVEL)) {
      TEST(level >= 0)("invalid level: %d\n", level);
    }
  
    for (mel = macroElements.begin(); mel != macroElements.end(); mel++) {
      elinfo->fillMacroInfo(*mel);
Thomas Witkowski's avatar
Thomas Witkowski committed
376
      sum += tinfo.recursive(&elInfoStack);
377
378
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
379
    elInfoStack.getBackElement();
380
381
382
383
384
385
386
387
388
389
    
    return (flag.isSet(Mesh::FILL_ADD_ALL)) ? sum : 0;
  }

  void Mesh::addDOFAdmin(DOFAdmin *localAdmin)
  {    
    FUNCNAME("Mesh::addDOFAdmin()");

    localAdmin->setMesh(this);

390
391
    std::vector<DOFAdmin*>::iterator dai = std::find(admin.begin(),admin.end(),localAdmin);

Thomas Witkowski's avatar
Thomas Witkowski committed
392
393
394
    TEST_EXIT(dai == admin.end())
      ("admin %s is already associated to mesh %s\n",
       localAdmin->getName().c_str(), this->getName().c_str());
395
396
397

    // ===== adding dofs to already existing elements ============================ 
    
Thomas Witkowski's avatar
Thomas Witkowski committed
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
    TEST_EXIT(!initialized)("Adding DOFAdmins to initilized meshes does not work yet!\n");

    if (initialized) {
      static bool pnd_1d_0[2] = {true, true};
      static bool pnd_1d_1[1] = {false};
      static bool pnd_2d_0[3] = {true, true, true};
      static bool pnd_2d_1[3] = {true, true, false};
      static bool pnd_2d_2[1] = {false};
      static bool pnd_3d_0[4] = {true, true, true, true};
      static bool pnd_3d_1[6] = {false, true, true, true, true, true};
      static bool pnd_3d_2[4] = {true, true, false, false};
      static bool pnd_3d_3[1] = {false};
      static bool *pnd_1d[2] = {pnd_1d_0, pnd_1d_1};
      static bool *pnd_2d[3] = {pnd_2d_0, pnd_2d_1, pnd_2d_2};
      static bool *pnd_3d[4] = {pnd_3d_0, pnd_3d_1, pnd_3d_2, pnd_3d_3};
      static bool **parentNeedsDOF[4] = {NULL, pnd_1d, pnd_2d, pnd_3d};

     
      std::list<struct delmem> delList;
      std::map< std::set<DegreeOfFreedom>, DegreeOfFreedom*> dofPtrMap;
      const DOFAdmin *vertexAdmin = getVertexAdmin();
      int vertexAdminPreDOFs = vertexAdmin->getNumberOfPreDOFs(VERTEX);

      // finding necessary node number for new admin

      int newNNode = 0;
      GeoIndex geoIndex;

      for (int d = 0; d < dim + 1; d++) {
	geoIndex = INDEX_OF_DIM(d, dim);
      
	if (localAdmin->getNumberOfDOFs(geoIndex)>0||nDOF[geoIndex]>0)
	  newNNode+=getGeo(geoIndex);
      }

      bool extendNodes = (newNNode>nNodeEl);
      nNodeEl = newNNode;

      TraverseStack stack;
      ElInfo *elInfo = NULL;
    
      elInfo = stack.traverseFirst(this, -1, CALL_EVERY_EL_PREORDER);
      while(elInfo) {
	Element *element = elInfo->getElement();
	DegreeOfFreedom *newDOF, **oldDOF, **dof = 
	  const_cast<DegreeOfFreedom**>(element->getDOF());

	int index = 0;

	if (extendNodes) {
	  oldDOF=dof;
	  element->setDOFPtrs();
	  dof=const_cast<DegreeOfFreedom**>(element->getDOF());
	  int index=0,oldIndex=0;
	  for (int d = 0; d < dim+1; d++) {
	    geoIndex = INDEX_OF_DIM(d, dim);
	    if (nDOF[geoIndex]>0) {
	      for(int i=0;i<getGeo(geoIndex);++i) 
		dof[index++]=oldDOF[oldIndex++];
	    }
	    else {
	      if (localAdmin->getNumberOfDOFs(geoIndex)>0) 
		index+=getGeo(geoIndex);
	    }
	  }
	
	  FREE_MEMORY(oldDOF, DegreeOfFreedom*, oldNNodes);

	  TEST_EXIT_DBG(index == nNodeEl)("ERROR: Number of entered nodes %f != number of nodes %f\n",index,nNodeEl);

	}


	index=0;

	// allocate new memory at elements
	for(int d = 0; d < dim+1; d++) {
	  geoIndex = INDEX_OF_DIM(d, dim);
      
	  int numberOfDOFs = localAdmin->getNumberOfDOFs(geoIndex);
	  int numberOfPreDOFs = nDOF[geoIndex];

	  if (numberOfDOFs>0||numberOfPreDOFs>0) {

	    // for all vertices/edges/...
	    for(int i = 0; i < getGeo(geoIndex); i++, index++) {
	      std::set<DegreeOfFreedom> dofSet;
	      for(int j = 0; j < d+1; j++) {
		dofSet.insert(dof[element->getVertexOfPosition(geoIndex, i, j)][vertexAdminPreDOFs]);
	      }
	    
	      if(element->isLeaf() || parentNeedsDOF[dim][d][i]) {
		if(dofPtrMap[dofSet] == NULL) {
		  if(localAdmin->getNumberOfDOFs(geoIndex)) {
		    newDOF = GET_MEMORY(DegreeOfFreedom, numberOfPreDOFs + numberOfDOFs);
		    // copy old dofs to new memory and free old memory
		    if(dof[index]) {
		      for(int j = 0; j < numberOfPreDOFs; j++) {
			newDOF[j] = dof[index][j];
		      }
		      //	  FREE_MEMORY(dof[index], DegreeOfFreedom, numberOfPreDOFs);
		      // Do not free memory. The information has to be used to identify the part in other elements.
		      // The memory is only marked for freeing.
		      struct delmem fm;
		      fm.ptr=dof[index];
		      fm.len=numberOfPreDOFs;
		      delList.push_back(fm);
		    }
		    for(int j = 0; j < numberOfDOFs; j++) {
		      newDOF[numberOfPreDOFs + j] = localAdmin->getDOFIndex();
		    }
		    dof[index] = newDOF;
		  }
		  dofPtrMap[dofSet] = dof[index];
		} else {
		  dof[index] = dofPtrMap[dofSet];
		}
	      }
	    }
	  }
	}
	elInfo = stack.traverseNext(elInfo);
      }
  
      // now free the old dof memory:

      std::list<struct delmem>::iterator it=delList.begin();
    
      while(it!=delList.end()) {
	FREE_MEMORY((*it).ptr, DegreeOfFreedom, (*it).len);
	it++;
      }

      delList.clear();

    }
534
535
536
537
538
539
540


    admin.push_back(localAdmin);

    nDOFEl = 0;

    localAdmin->setNumberOfPreDOFs(VERTEX,nDOF[VERTEX]);
541
    nDOF[VERTEX] += localAdmin->getNumberOfDOFs(VERTEX);
542
543
    nDOFEl += getGeo(VERTEX) * nDOF[VERTEX];

544
    if (dim > 1) {
545
      localAdmin->setNumberOfPreDOFs(EDGE,nDOF[EDGE]);
546
      nDOF[EDGE] += localAdmin->getNumberOfDOFs(EDGE);
547
548
549
550
551
552
553
      nDOFEl += getGeo(EDGE) * nDOF[EDGE];
    }

    localAdmin->setNumberOfPreDOFs(CENTER,nDOF[CENTER]);
    nDOF[CENTER]  += localAdmin->getNumberOfDOFs(CENTER);
    nDOFEl += nDOF[CENTER];

554
    TEST_EXIT_DBG(nDOF[VERTEX] > 0)("no vertex dofs\n");
555

556
557
    node[VERTEX] = 0;
    nNodeEl = getGeo(VERTEX);
558

559
560
    if (dim > 1) {
      node[EDGE] = nNodeEl;
561
562
      if (nDOF[EDGE] > 0) 
	nNodeEl += getGeo(EDGE);
563
564
    }

565
    if (dim == 3) {
566
      localAdmin->setNumberOfPreDOFs(FACE,nDOF[FACE]);
567
568
569
570
571
      nDOF[FACE] += localAdmin->getNumberOfDOFs(FACE);
      nDOFEl += getGeo(FACE) * nDOF[FACE];
      node[FACE] = nNodeEl;
      if (nDOF[FACE] > 0) 
	nNodeEl += getGeo(FACE);
572
573
    }

574
    node[CENTER] = nNodeEl;
575
    if (nDOF[CENTER] > 0) {
576
      nNodeEl += 1;
577
    }
578
579
580
581
  }

  void Mesh::dofCompress()
  {
582
583
    FUNCNAME("Mesh::dofCompress()");
    Flag fill_flag;
584

Thomas Witkowski's avatar
Thomas Witkowski committed
585
    for (int iadmin = 0; iadmin < static_cast<int>(admin.size()); iadmin++) {
586
587
588
      compressAdmin = admin[iadmin];

      TEST_EXIT_DBG(compressAdmin)("no admin[%d] in mesh\n", iadmin);
589
      
Thomas Witkowski's avatar
Thomas Witkowski committed
590
591
      int size = compressAdmin->getSize();
      if (size < 1) 
592
	continue;
Thomas Witkowski's avatar
Thomas Witkowski committed
593

594
595
      if (compressAdmin->getUsedDOFs() < 1)    
	continue;
Thomas Witkowski's avatar
Thomas Witkowski committed
596

597
598
      if (compressAdmin->getHoleCount() < 1)    
	continue;
Thomas Witkowski's avatar
Thomas Witkowski committed
599
  
600
601
602
603
604
605
606
607
      newDOF.resize(size);
      
      compressAdmin->compress(newDOF);
      
      if (preserveCoarseDOFs) {
	fill_flag = Mesh::CALL_EVERY_EL_PREORDER | Mesh::FILL_NOTHING;
      } else {
	fill_flag = Mesh::CALL_LEAF_EL | Mesh::FILL_NOTHING;
608
      }
609
      
Thomas Witkowski's avatar
Thomas Witkowski committed
610
611
      traverse(-1, fill_flag, newDOFFct1);
      traverse(-1, fill_flag, newDOFFct2);
612
613
614
      
      newDOF.resize(0);
    }   
615
616
617
618
619
  }


  DegreeOfFreedom *Mesh::getDOF(GeoIndex position)
  {
620
    FUNCNAME("Mesh::getDOF()");
621

622
    TEST_EXIT_DBG(position >= CENTER && position <= FACE)
623
      ("unknown position %d\n", position);
624

625
626
    int ndof = getNumberOfDOFs(position);
    if (ndof <= 0) 
627
      return NULL;
628

629
    DegreeOfFreedom *dof = new DegreeOfFreedom[ndof];
630

631
632
    for (int i = 0; i < getNumberOfDOFAdmin(); i++) {
      const DOFAdmin *localAdmin = &getDOFAdmin(i);
633
      TEST_EXIT_DBG(localAdmin)("no admin[%d]\n", i);
634
635
636
637
      
      int n  = localAdmin->getNumberOfDOFs(position);
      int n0 = localAdmin->getNumberOfPreDOFs(position);
      
638
      TEST_EXIT_DBG(n + n0 <= ndof)("n=%d, n0=%d too large: ndof=%d\n", n, n0, ndof);
639
      
640
      for (int j = 0; j < n; j++)
641
642
	dof[n0 + j] = const_cast<DOFAdmin*>(localAdmin)->getDOFIndex();
    }
643
  
644
    return dof;
645
646
647
648
649
  }


  DegreeOfFreedom **Mesh::createDOFPtrs()
  {
650
    FUNCNAME("Mesh::createDOFPtrs()");
651
652

    if (nNodeEl <= 0)
653
      return NULL;
654

655
    DegreeOfFreedom **ptrs = new DegreeOfFreedom*[nNodeEl];
656
    for (int i = 0; i < nNodeEl; i++)
657
658
      ptrs[i] = NULL;

659
    return ptrs;
660
661
662
663
  }

  void Mesh::freeDOFPtrs(DegreeOfFreedom **ptrs)
  {
664
    FUNCNAME("Mesh::freeDOFPtrs()");
665

666
    TEST_EXIT_DBG(ptrs)("ptrs=NULL\n");
667
668
669
670

    if (nNodeEl <= 0)
      return;
  
671
    delete [] ptrs;
672
673
674
  }


Thomas Witkowski's avatar
Thomas Witkowski committed
675
  const DOFAdmin *Mesh::createDOFAdmin(const std::string& lname, DimVec<int> lnDOF)
676
  {
677
    FUNCNAME("Mesh::createDOFAdmin()");
678

Thomas Witkowski's avatar
Thomas Witkowski committed
679
    DOFAdmin *localAdmin = new DOFAdmin(this, lname);
680

Thomas Witkowski's avatar
Thomas Witkowski committed
681
682
    for (int i = 0; i < dim + 1; i++)
      localAdmin->setNumberOfDOFs(i, lnDOF[i]);
683
684
685

    addDOFAdmin(localAdmin);

686
    return localAdmin;
687
688
689
690
691
692
693
  }


  const DOFAdmin* Mesh::getVertexAdmin() const
  {
    const DOFAdmin *localAdmin = NULL;

694
695
696
697
698
699
    for (int i = 0; i < static_cast<int>(admin.size()); i++) {
      if (admin[i]->getNumberOfDOFs(VERTEX)) {
	if (!localAdmin)  
	  localAdmin = admin[i];
	else if (admin[i]->getSize() < localAdmin->getSize())
	  localAdmin = admin[i];
700
      }
701
702
    }

703
    return localAdmin;
704
705
706
707
  }

  void Mesh::freeDOF(DegreeOfFreedom* dof, GeoIndex position)
  {
708
    FUNCNAME("Mesh::freeDOF()");
709

710
    TEST_EXIT_DBG(position >= CENTER && position <= FACE)
711
      ("unknown position %d\n", position);
712

713
714
715
716
    int ndof = nDOF[position];
    if (ndof) {
      if (!dof) {
	MSG("dof = NULL, but ndof=%d\n", ndof);
717
718
	return;
      }
719
720
721
722
723
724
    } else  {
      if (dof) {
	MSG("dof != NULL, but ndof=0\n");
      }
      return;
    }
725

726
    TEST_EXIT_DBG(ndof <= MAX_DOF)
727
      ("ndof too big: ndof=%d, MAX_DOF=%d\n", ndof, MAX_DOF);
728

729
730
731
732
733
734
735
736
737
738
    for (int i = 0; i < static_cast<int>(admin.size()); i++) {
      DOFAdmin *localAdmin = admin[i];
      int n = localAdmin->getNumberOfDOFs(position);
      int n0 = localAdmin->getNumberOfPreDOFs(position);
      
      TEST_EXIT_DBG(n + n0 <= ndof)("n=%d, n0=%d too large: ndof=%d\n", n, n0, ndof);
      
      for (int j = 0; j < n; j++)
	localAdmin->freeDOFIndex(dof[n0 + j]);
    }
739

740
    delete [] dof;
741
742
743
744
745
  }

  void Mesh::freeElement(Element* el)
  {
    freeDOFPtrs(const_cast<DegreeOfFreedom**>(el->getDOF()));
Thomas Witkowski's avatar
Thomas Witkowski committed
746
    delete el;
747
748
749
750
751
752
  }


  Element* Mesh::createNewElement(Element *parent)
  {
    FUNCNAME("Mesh::createNewElement()");
753
754

    TEST_EXIT_DBG(elementPrototype)("no element prototype\n");
755
756
757

    Element *el = parent ? parent->clone() : elementPrototype->clone();
  
758
    if (!parent && elementDataPrototype) {
759
760
761
762
763
764
765
766
      el->setElementData(elementDataPrototype->clone()); 
    } else {
      el->setElementData(NULL); // must be done in ElementData::refineElementData()
    }

    return el;
  }

767

768
769
770
771
  ElInfo* Mesh::createNewElInfo()
  {
    switch(dim) {
    case 1:
Thomas Witkowski's avatar
Thomas Witkowski committed
772
      return new ElInfo1d(this);
773
774
      break;
    case 2:
Thomas Witkowski's avatar
Thomas Witkowski committed
775
      return new ElInfo2d(this);
776
777
      break;
    case 3:
Thomas Witkowski's avatar
Thomas Witkowski committed
778
      return new ElInfo3d(this);
779
780
781
782
      break;
    default:
      ERROR_EXIT("invalid dim\n");
      return NULL;
783
    }
784
785
786
787
  }

  bool Mesh::findElInfoAtPoint(const WorldVector<double>& xy,
			       ElInfo *el_info,
788
789
			       DimVec<double>& bary,
			       const MacroElement *start_mel,
790
			       const WorldVector<double> *xy0,
791
			       double *sp)
792
793
794
795
796
797
798
799
800
801
  {
    static const MacroElement *mel = NULL;
    DimVec<double> lambda(dim, NO_INIT);
    ElInfo *mel_info = NULL;

    mel_info = createNewElInfo();

    if (start_mel != NULL)
      mel = start_mel;
    else
802
      if ((mel == NULL) || (mel->getElement()->getMesh() != this))
803
804
805
	mel = *(macroElements.begin());

    mel_info->setFillFlag(Mesh::FILL_COORDS);
806
    g_xy = &xy;
807
    g_xy0 = xy0;
808
    g_sp = sp;
809
810
811

    mel_info->fillMacroInfo(mel);

812
    int k;
813
814
815
816
817
818
819
820
821
822
    while ((k = mel_info->worldToCoord(xy, &lambda)) >= 0) {
      if (mel->getNeighbour(k)) {
	mel = mel->getNeighbour(k);
	mel_info->fillMacroInfo(mel);
	continue;
      }
      break;
    }

    /* now, descend in tree to find leaf element at point */
823
    bool inside = findElementAtPointRecursive(mel_info, lambda, k, el_info);
824
825
    for (int i = 0; i <= dim; i++)
      bary[i] = final_lambda[i];   
826
  
Thomas Witkowski's avatar
Thomas Witkowski committed
827
    delete mel_info;
828

829
    return inside;
830
831
832
  }

  bool Mesh::findElementAtPoint(const WorldVector<double>&  xy,
833
834
				Element **elp, 
				DimVec<double>& bary,
835
				const MacroElement *start_mel,
836
837
				const WorldVector<double> *xy0,
				double *sp)
838
  {
839
840
    ElInfo *el_info = createNewElInfo();
    int val = findElInfoAtPoint(xy, el_info, bary, start_mel, xy0, sp);
841
842
843

    *elp = el_info->getElement();

Thomas Witkowski's avatar
Thomas Witkowski committed
844
    delete el_info;
845

846
    return val;
847
848
  }

849
  bool Mesh::findElementAtPointRecursive(ElInfo *el_info,
850
					 const DimVec<double>& lambda,
851
					 int outside,
852
853
					 ElInfo* final_el_info)
  {
854
    FUNCNAME("Mesh::findElementAtPointRecursive()");
855
856
    Element *el = el_info->getElement();
    DimVec<double> c_lambda(dim, NO_INIT);
857
858
    int inside;
    int ichild, c_outside;
859
860
861
862

    if (el->isLeaf()) {
      *final_el_info = *el_info;
      if (outside < 0) {
863
864
865
866
	for (int i = 0; i <= dim; i++) {
	  final_lambda[i] = lambda[i];
	}

867
	return true;
868
869
870
871
872
873
874
875
876
877
878
879
      }  else {  /* outside */
	if (g_xy0) { /* find boundary point of [xy0, xy] */
	  el_info->worldToCoord(*(g_xy0), &c_lambda);
	  double s = lambda[outside] / (lambda[outside] - c_lambda[outside]);
	  for (int i = 0; i <= dim; i++) {
	    final_lambda[i] = s * c_lambda[i] + (1.0-s) * lambda[i];
	  }
	  if (g_sp) {
	    *(g_sp) = s;
	  }
	  if (dim == 3) 
	    MSG("outside finest level on el %d: s=%.3e\n", el->getIndex(), s);
880

881
	  return false;  /* ??? */
882
	} else {
883
	  return false;
884
	}
885
      }
886
887
    }

888
    ElInfo *c_el_info = createNewElInfo();
889

890
    if (dim == 1) {
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
      if (lambda[0] >= lambda[1]) {
	c_el_info->fillElInfo(0, el_info);
	if (outside >= 0) {
	  outside = el_info->worldToCoord(*(g_xy), &c_lambda);
	  if (outside >= 0) ERROR("point outside domain\n");
	} else {
	  c_lambda[0] = lambda[0] - lambda[1];
	  c_lambda[1] = 2.0 * lambda[1];
	}
      } else {
	c_el_info->fillElInfo(1, el_info);
	if (outside >= 0)  {
	  outside = el_info->worldToCoord(*(g_xy), &c_lambda);
	  if (outside >= 0) ERROR("point outside domain\n");
	} else {
	  c_lambda[1] = lambda[1] - lambda[0];
	  c_lambda[0] = 2.0 * lambda[0];
	}
      }
    } /* DIM == 1 */

912
    if (dim == 2) {
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
      if (lambda[0] >= lambda[1]) {
	c_el_info->fillElInfo(0, el_info);
	if (el->isNewCoordSet()) {
	  outside = c_el_info->worldToCoord(*(g_xy), &c_lambda);
	  if (outside >= 0) {
	    ERROR("outside curved boundary child 0\n");
	  }
	} else {
	  c_lambda[0] = lambda[2];
	  c_lambda[1] = lambda[0] - lambda[1];
	  c_lambda[2] = 2.0 * lambda[1];
	}
      } else {
	c_el_info->fillElInfo(1, el_info);
	if (el->isNewCoordSet()) {
	  outside = c_el_info->worldToCoord(*(g_xy), &c_lambda);
	  if (outside >= 0) {
	    ERROR("outside curved boundary child 1\n");
	  }
	} else {
	  c_lambda[0] = lambda[1] - lambda[0];
	  c_lambda[1] = lambda[2];
	  c_lambda[2] = 2.0 * lambda[0];
	}
      }
    } /* DIM == 2 */

940
    if (dim == 3) {
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
      if (el->isNewCoordSet()) {
	if (lambda[0] >= lambda[1])
	  ichild = 0;
	else
	  ichild = 1;
	c_el_info->fillElInfo(ichild, el_info);
	c_outside = c_el_info->worldToCoord(*(g_xy), &c_lambda);

	if (c_outside>=0) {  /* test is other child is better... */
	  DimVec<double> c_lambda2(dim, NO_INIT);
	  int c_outside2;
	  ElInfo *c_el_info2 = createNewElInfo();

	  c_el_info2->fillElInfo(1-ichild, el_info);
	  c_outside2 = c_el_info2->worldToCoord(*(g_xy), &c_lambda2);

	  MSG("new_coord CHILD %d: outside=%d, lambda=(%.2f %.2f %.2f %.2f)\n",
	      ichild, c_outside, c_lambda[0],c_lambda[1],c_lambda[2],c_lambda[3]);
	  MSG("new_coord CHILD %d: outside=%d, lambda=(%.2f %.2f %.2f %.2f)\n",
	      1-ichild, c_outside2, c_lambda2[0],c_lambda2[1],c_lambda2[2],
	      c_lambda2[3]);

	  if ((c_outside2 < 0) || (c_lambda2[c_outside2] > c_lambda[c_outside])) {
964
965
966
	    for (int i = 0; i <= dim; i++) {
	      c_lambda[i] = c_lambda2[i];
	    }
967
968
969
970
	    c_outside = c_outside2;
	    *c_el_info = *c_el_info2;
	    ichild = 1 - ichild;
	  }
Thomas Witkowski's avatar
Thomas Witkowski committed
971
	  delete c_el_info2;
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
	}
	outside = c_outside;
      } else {  /* no new_coord */
	if (lambda[0] >= lambda[1]) {
	  c_el_info->fillElInfo(0, el_info);
	  c_lambda[0] = lambda[0] - lambda[1];
	  c_lambda[1] = lambda[Tetrahedron::childVertex[(dynamic_cast<ElInfo3d*>(el_info))->
							getType()][0][1]];
	  c_lambda[2] = lambda[Tetrahedron::childVertex[(dynamic_cast<ElInfo3d*>(el_info))->
							getType()][0][2]];
	  c_lambda[3] = 2.0 * lambda[1];
	} else {
	  c_el_info->fillElInfo(1, el_info);
	  c_lambda[0] = lambda[1] - lambda[0];
	  c_lambda[1] = lambda[Tetrahedron::childVertex[(dynamic_cast<ElInfo3d*>(el_info))->
							getType()][1][1]];
	  c_lambda[2] = lambda[Tetrahedron::childVertex[(dynamic_cast<ElInfo3d*>(el_info))->
							getType()][1][2]];
	  c_lambda[3] = 2.0 * lambda[0];
	}
      }
    }  /* DIM == 3 */

    inside = findElementAtPointRecursive(c_el_info, c_lambda, outside, 
					 final_el_info);
Thomas Witkowski's avatar
Thomas Witkowski committed
997
    delete c_el_info;
998

999
    return inside; 
1000
1001
1002
  }


Thomas Witkowski's avatar
Thomas Witkowski committed
1003
1004
1005
1006
  void Mesh::setDiameter(const WorldVector<double>& w) 
  { 
    diam = w; 
  }
1007

Thomas Witkowski's avatar
Thomas Witkowski committed
1008
1009
1010
1011
  void Mesh::setDiameter(int i, double w) 
  { 
    diam[i] = w; 
  }
1012
1013
1014
1015
1016
1017
1018
1019
1020
1021
1022
1023


  int Mesh::newDOFFct1(ElInfo* ei) {
    ei->getElement()->newDOFFct1(compressAdmin);
    return 0;
  }

  int Mesh::newDOFFct2(ElInfo* ei) {
    ei->getElement()->newDOFFct2(compressAdmin);
    return 0;
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
1024
  void Mesh::serialize(std::ostream &out)
1025
1026
1027
1028
  {
    serializedDOFs.clear();

    // write name
Thomas Witkowski's avatar
Thomas Witkowski committed
1029
    out << name << "\n";
1030
1031
1032
1033
1034
1035
1036
1037
1038
1039
1040
1041
1042
1043
1044
1045
1046
1047
1048
1049
1050
1051
1052
1053
1054
1055
1056
1057
1058
1059
1060
1061
1062
1063
1064
1065
1066
1067
1068
1069
1070

    // write dim
    out.write(reinterpret_cast<const char*>(&dim), sizeof(int));

    // write nVertices
    out.write(reinterpret_cast<const char*>(&nVertices), sizeof(int));

    // write nEdges
    out.write(reinterpret_cast<const char*>(&nEdges), sizeof(int));

    // write nLeaves
    out.write(reinterpret_cast<const char*>(&nLeaves), sizeof(int));

    // write nElements
    out.write(reinterpret_cast<const char*>(&nElements), sizeof(int));

    // write nFaces
    out.write(reinterpret_cast<const char*>(&nFaces), sizeof(int));

    // write maxEdgeNeigh
    out.write(reinterpret_cast<const char*>(&maxEdgeNeigh), sizeof(int));

    // write diam
    diam.serialize(out);

    // write preserveCoarseDOFs
    out.write(reinterpret_cast<const char*>(&preserveCoarseDOFs), sizeof(bool));

    // write nDOFEl
    out.write(reinterpret_cast<const char*>(&nDOFEl), sizeof(int));

    // write nDOF
    nDOF.serialize(out);

    // write nNodeEl
    out.write(reinterpret_cast<const char*>(&nNodeEl), sizeof(int));

    // write node
    node.serialize(out);

    // write admins
1071
    int size = static_cast<int>(admin.size());
1072
    out.write(reinterpret_cast<const char*>(&size), sizeof(int));
1073
    for (int i = 0; i < size; i++)
1074
1075
1076
1077
1078
      admin[i]->serialize(out);

    // write macroElements
    size = static_cast<int>(macroElements.size());
    out.write(reinterpret_cast<const char*>(&size), sizeof(int));
1079
    for (int i = 0; i < size; i++)
1080
1081
1082
1083
1084
1085
1086
1087
1088
1089
1090
      macroElements[i]->serialize(out);

    // write elementIndex
    out.write(reinterpret_cast<const char*>(&elementIndex), sizeof(int));

    // write initialized
    out.write(reinterpret_cast<const char*>(&initialized), sizeof(bool));

    serializedDOFs.clear();
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
1091
  void Mesh::deserialize(std::istream &in)
1092
1093
1094
1095
1096
1097
1098
1099
1100
1101
  {
    serializedDOFs.clear();

    // read name
    in >> name;
    in.get();

    // read dim
    int oldVal = dim;
    in.read(reinterpret_cast<char*>(&dim), sizeof(int));
1102
    TEST_EXIT_DBG((oldVal == 0) || (dim == oldVal))("invalid dimension\n");
1103
1104
1105
1106
1107
1108
1109
1110
1111
1112
1113
1114
1115
1116
1117
1118
1119
1120
1121
1122
1123
1124
1125
1126
1127
1128
1129
1130

    // read nVertices
    in.read(reinterpret_cast<char*>(&nVertices), sizeof(int));

    // read nEdges
    in.read(reinterpret_cast<char*>(&nEdges), sizeof(int));

    // read nLeaves
    in.read(reinterpret_cast<char*>(&nLeaves), sizeof(int));

    // read nElements
    in.read(reinterpret_cast<char*>(&nElements), sizeof(int));

    // read nFaces
    in.read(reinterpret_cast<char*>(&nFaces), sizeof(int));

    // read maxEdgeNeigh
    in.read(reinterpret_cast<char*>(&maxEdgeNeigh), sizeof(int));

    // diam
    diam.deserialize(in);

    // read preserveCoarseDOFs
    in.read(reinterpret_cast<char*>(&preserveCoarseDOFs), sizeof(bool));

    // read nDOFEl
    oldVal = nDOFEl;
    in.read(reinterpret_cast<char*>(&nDOFEl), sizeof(int));
1131
    TEST_EXIT_DBG((oldVal == 0) || (nDOFEl == oldVal))("invalid nDOFEl\n");
1132
1133
1134
1135
1136
1137
1138

    // read nDOF
    nDOF.deserialize(in);

    // read nNodeEl
    oldVal = nNodeEl;
    in.read(reinterpret_cast<char*>(&nNodeEl), sizeof(int));
1139
    TEST_EXIT_DBG((oldVal == 0) || (nNodeEl == oldVal))("invalid nNodeEl\n");
1140
1141
1142
1143
1144

    // read node
    node.deserialize(in);

    // read admins
1145
    int size;
1146
1147
    in.read(reinterpret_cast<char*>(&size), sizeof(int));
    admin.resize(size, NULL);
1148
    for (int i = 0; i < size; i++) {
1149
      if (!admin[i])
Thomas Witkowski's avatar
Thomas Witkowski committed
1150
	admin[i] = new DOFAdmin(this);
1151

1152
1153
1154
1155
1156
1157
      admin[i]->deserialize(in);
    }

    // read macroElements
    in.read(reinterpret_cast<char*>(&size), sizeof(int));

Thomas Witkowski's avatar
Thomas Witkowski committed
1158
    std::vector< std::vector<int> > neighbourIndices(size);
1159

1160
    for (int i = 0; i < static_cast<int>(macroElements.size()); i++) {
1161
      if (macroElements[i])
Thomas Witkowski's avatar
Thomas Witkowski committed
1162
	delete macroElements[i];
1163
    }
1164

1165
    macroElements.resize(size);
1166
    for (int i = 0; i < size; i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
1167
      macroElements[i] = new MacroElement(dim);
1168
1169
1170
1171
1172
1173
1174
1175
1176
1177
1178
      macroElements[i]->writeNeighboursTo(&(neighbourIndices[i]));
      macroElements[i]->deserialize(in);
    }

    // read elementIndex
    in.read(reinterpret_cast<char*>(&elementIndex), sizeof(int));

    // read initialized
    in.read(reinterpret_cast<char*>(&initialized), sizeof(bool));

    // set neighbour pointer in macro elements
1179
1180
1181
    int neighs = getGeo(NEIGH);
    for (int i = 0; i < static_cast<int>(macroElements.size()); i++) {
      for (int j = 0; j < neighs; j++) {
1182
1183
1184
1185
1186
1187
1188
1189
1190
1191
1192
1193
	int index = neighbourIndices[i][j];
	if(index != -1) {
	  macroElements[i]->setNeighbour(j, macroElements[index]);
	} else {
	  macroElements[i]->setNeighbour(j, NULL);
	}
      }
    }

    // set mesh pointer in elements
    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(this, -1, CALL_EVERY_EL_PREORDER);
1194
    while (elInfo) {
1195
1196
1197
1198
1199
1200
1201
1202
1203
      elInfo->getElement()->setMesh(this);
      elInfo = stack.traverseNext(elInfo);
    }

    serializedDOFs.clear();
  }

  void Mesh::initialize() 
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
1204
1205
1206
    std::string macroFilename("");
    std::string valueFilename("");
    std::string periodicFile("");
1207
1208
1209
1210
1211
1212
1213
1214
1215
    int check = 1;

    GET_PARAMETER(0, name + "->macro file name",  &macroFilename);
    GET_PARAMETER(0, name + "->value file name",  &valueFilename);
    GET_PARAMETER(0, name + "->periodic file", &periodicFile);
    GET_PARAMETER(0, name + "->check", "%d", &check);
    GET_PARAMETER(0, name + "->preserve coarse dofs", "%d", &preserveCoarseDOFs);

    if (macroFilename.length()) {
1216
1217
1218
      macroFileInfo = MacroReader::readMacro(macroFilename.c_str(), this,
					     periodicFile == "" ? NULL : periodicFile.c_str(),
					     check);
1219
1220
1221
1222
1223
1224
1225
1226
1227
1228
1229
1230

      // If there is no value file which should be written, we can delete
      // the information of the macro file.
      if (!valueFilename.length()) {
	clearMacroFileInfo();
      }
    }

    initialized = true;
  }

  bool Mesh::associated(DegreeOfFreedom dof1, DegreeOfFreedom dof2) {
Thomas Witkowski's avatar
Thomas Witkowski committed
1231
1232
    std::map<BoundaryType, VertexVector*>::iterator it;
    std::map<BoundaryType, VertexVector*>::iterator end = periodicAssociations.end();
1233
1234
1235
1236
1237
1238
1239
1240
    for (it = periodicAssociations.begin(); it != end; ++it) {
      if ((*(it->second))[dof1] == dof2)
	return true;
    }
    return false;
  }

  bool Mesh::indirectlyAssociated(DegreeOfFreedom dof1, DegreeOfFreedom dof2) {
Thomas Witkowski's avatar
Thomas Witkowski committed
1241
1242
1243
    std::vector<DegreeOfFreedom> associatedToDOF1;
    std::map<BoundaryType, VertexVector*>::iterator it;
    std::map<BoundaryType, VertexVector*>::iterator end = periodicAssociations.end();
1244
1245
1246
    DegreeOfFreedom dof, assDOF;

    associatedToDOF1.push_back(dof1);
Thomas Witkowski's avatar
Thomas Witkowski committed
1247
1248
1249
    for (it = periodicAssociations.begin(); it != end; ++it) {
      int size = static_cast<int>(associatedToDOF1.size());
      for (int i = 0; i < size; i++) {
1250
1251
	dof = associatedToDOF1[i];
	assDOF = (*(it->second))[dof];
Thomas Witkowski's avatar
Thomas Witkowski committed
1252
	if (assDOF == dof2) {
1253
1254
	  return true;
	} else {
Thomas Witkowski's avatar
Thomas Witkowski committed
1255
1256
	  if (assDOF != dof) 
	    associatedToDOF1.push_back(assDOF);
1257
1258
1259
1260
1261
1262
1263
1264
	}
      }
    }
    return false;
  }

  void Mesh::clearMacroFileInfo()
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
1265
    macroFileInfo->clear();
Thomas Witkowski's avatar
Thomas Witkowski committed
1266
    delete macroFileInfo;
1267
    macroFileInfo = NULL;
1268
  }
1269
1270
1271

  int Mesh::calcMemoryUsage()
  {
1272
    int result = sizeof(Mesh);
1273

1274
1275
1276
1277
1278
1279
    result += nDOFEl;
    for (int i = 0; i < static_cast<int>(admin.size()); i++) {
      result += admin[i]->calcMemoryUsage();
      result += admin[i]->getUsedSize() * sizeof(DegreeOfFreedom);
    }
    
1280
1281
    for (int i = 0; i < static_cast<int>(macroElements.size()); i++)
      result += macroElements[i]->calcMemoryUsage();    
1282
1283
1284
    
    return result;
  }
1285
}