Mesh.cc 34.5 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
#include "AdaptStationary.h"
#include "AdaptInstationary.h"
#include "FiniteElemSpace.h"
#include "ElementData.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
11
#include "ElementDofIterator.h"
12
13
14
15
16
17
18
19
20
21
22
23
24
#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
25
#include "ElInfoStack.h"
26
#include "Serializer.h"
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46

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;


47
48
49
  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);
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66

  //**************************************************************************
  //  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;
Thomas Witkowski's avatar
Thomas Witkowski committed
67
  std::vector<DegreeOfFreedom> Mesh::dof_used;
68
  const int Mesh::MAX_DOF = 100;
69
  std::map<std::pair<DegreeOfFreedom, int>, DegreeOfFreedom*> Mesh::serializedDOFs;
70
71
72

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


Thomas Witkowski's avatar
Thomas Witkowski committed
77
  Mesh::Mesh(std::string aName, int dimension) 
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
    : 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),
94
      macroFileInfo(NULL),
95
      changeIndex(0),
96
97
98
      final_lambda(dimension, DEFAULT_VALUE, 0.0)
  {

99
    FUNCNAME("Mesh::Mesh()");
100
101
102
103

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

    elementPrototype->setIndex(-1);

118
119
    elementIndex = 0;
  }
120
121

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

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

131
    Element::deletedDOFs.clear();
132

133
134
135
136
137
138
    if (macroFileInfo != NULL)
      delete macroFileInfo;    
    if (elementPrototype)
      delete elementPrototype;    
    if (elementDataPrototype)
      delete elementDataPrototype;    
139
    
140
141
    for (int i = 0; i < static_cast<int>(admin.size()); i++)
      delete admin[i];    
142
  }
143
144

  Mesh& Mesh::operator=(const Mesh& m)
145
  {
146
147
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
    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
176
      admin[i] = new DOFAdmin(this);
177
      *(admin[i]) = *(m.admin[i]);
178
179
      admin[i]->setMesh(this);
    }
180

181

182
    /* ====================== Copy macro elements =================== */
183
  
184
185
186
187
188
189
190
191
192
    // 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;
193

194
195
    macroElements.clear();

196
197
198
    // 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();
199
	 it != m.macroElements.end(); ++it, insertCounter++) {
200

201
      // Create new MacroElement.
Thomas Witkowski's avatar
Thomas Witkowski committed
202
      MacroElement *el = new MacroElement(dim);
203

204
205
206
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
      // 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;
  }

247
248
249
250
251
252
253
254
255
256
257
258
  void Mesh::updateNumberOfLeaves()
  {
    nLeaves = 0;

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

259
260
261
262
263
264
  void Mesh::addMacroElement(MacroElement* me) 
  {
    macroElements.push_back(me); 
    me->setIndex(macroElements.size());
  }

265
  void Mesh::removeMacroElements(std::set<MacroElement*>& macros,
Thomas Witkowski's avatar
Thomas Witkowski committed
266
				 const FiniteElemSpace *feSpace) 
267
268
269
  {
    FUNCNAME("Mesh::removeMacroElement()");

Thomas Witkowski's avatar
Thomas Witkowski committed
270
271
272
273
274
275
276
    typedef std::map<const DegreeOfFreedom*, std::set<MacroElement*> > DofElMap;
    typedef std::map<const DegreeOfFreedom*, GeoIndex> DofPosMap;

    TEST_EXIT(admin.size() == 1)("Not yet implemented for multiple admins!\n");
    TEST_EXIT(admin[0])("There is something wrong!\n");

    ElementDofIterator elDofIter(feSpace);
277

278
    // Map that stores for each dof pointer (which may have a list of dofs)
279
    // all macro element indices that own this dof.
Thomas Witkowski's avatar
Thomas Witkowski committed
280
281
    DofElMap dofsOwner;
    DofPosMap dofsPosIndex;
282
283
284
285


    // === Determine to all dofs the macro elements where the dof is part of. ===

286
    for (std::deque<MacroElement*>::iterator macroIt = macroElements.begin();
287
	 macroIt != macroElements.end(); ++macroIt) {
Thomas Witkowski's avatar
Thomas Witkowski committed
288
289
290
291
292
      elDofIter.reset((*macroIt)->getElement());
      do {
	dofsOwner[elDofIter.getDofPtr()].insert(*macroIt);
	dofsPosIndex[elDofIter.getDofPtr()] = elDofIter.getPosIndex();
      } while (elDofIter.nextStrict());
293
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
294

295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
    // === Remove macro elements from mesh macro element list. ===

    // Removing arbitrary elements from an std::deque is very slow. Therefore, we
    // create a new deque with all macro elements that should not be deleted. The
    // macro element deque is than replaced by the new created one.

    std::deque<MacroElement*> newMacroElements;
    for (std::deque<MacroElement*>::iterator elIter = macroElements.begin();
	 elIter != macroElements.end(); ++elIter) {
      // If the current mesh macro element should not be deleted, i.e., it is not a
      // member of the list of macro elements to be deleted, is is inserted to the new
      // macro element list.
      if (macros.find(*elIter) == macros.end())
	newMacroElements.push_back(*elIter);      
    }
    // And replace the macro element list with the new one.
    macroElements.clear();
    macroElements = newMacroElements;

Thomas Witkowski's avatar
Thomas Witkowski committed
314

315
316
317
318
319
    // === For all macro elements to be deleted, delete them also to be neighbours ===
    // === of some other macro elements.                                           ===
    
    for (std::set<MacroElement*>::iterator macroIt = macros.begin();
	 macroIt != macros.end(); ++macroIt) {
Thomas Witkowski's avatar
Thomas Witkowski committed
320

321
322
323
324
      // 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)) {
Thomas Witkowski's avatar
Thomas Witkowski committed
325
326
	  for (int j = 0; j <= dim; j++)
	    if ((*macroIt)->getNeighbour(i)->getNeighbour(j) == *macroIt)
327
328
329
330
331
	      (*macroIt)->getNeighbour(i)->setNeighbour(j, NULL);
	} else {
	  // There is no neighbour at this edge, so we have to decrease the number
	  // of edges in the mesh.
	  nEdges--;
332
	}
333
334
      }

335
      // Decrease also the number of elements of the mesh.
336
337
      nLeaves--;
      nElements--;
338
    }     
339
340


341
342
    // === Check now all the dofs, that have no owner anymore and therefore have ===
    // === to be removed.                                                        ===
343
344

    int nRemainDofs = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
345
    for (DofElMap::iterator dofsIt = dofsOwner.begin(); 
346
347
348
349
350
351
352
353
354
355
356
357
358
359
	 dofsIt != dofsOwner.end(); ++dofsIt) {
      
      bool deleteDof = true;

      for (std::set<MacroElement*>::iterator elIter = dofsIt->second.begin();
	   elIter != dofsIt->second.end(); ++elIter) {
	std::set<MacroElement*>::iterator mIt = macros.find(*elIter);
	if (mIt == macros.end()) {
	  deleteDof = false;
	  break;
	}
      }

      if (deleteDof)
360
	freeDOF(const_cast<DegreeOfFreedom*>(dofsIt->first), 
Thomas Witkowski's avatar
Thomas Witkowski committed
361
362
		dofsPosIndex[dofsIt->first]);
      else
363
	nRemainDofs++;
364
365
    }

366
367
368
369
370
371
372

    // === Finally, remove the macro elements from memory. ===

    for (std::set<MacroElement*>::iterator macroIt = macros.begin();
	 macroIt != macros.end(); ++macroIt)
      delete *macroIt;

373
    nVertices = nRemainDofs;
374
  }
375
376
377
378
379
380
381

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

    localAdmin->setMesh(this);

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

Thomas Witkowski's avatar
Thomas Witkowski committed
385
386
387
    TEST_EXIT(dai == admin.end())
      ("admin %s is already associated to mesh %s\n",
       localAdmin->getName().c_str(), this->getName().c_str());
388

389
    // if this will be required, see the untested code in revision < 224
390
    //    TEST_EXIT(!initialized)("Adding DOFAdmins to initilized meshes does not work yet!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
391

392
393
394
395
396
    admin.push_back(localAdmin);

    nDOFEl = 0;

    localAdmin->setNumberOfPreDOFs(VERTEX,nDOF[VERTEX]);
397
    nDOF[VERTEX] += localAdmin->getNumberOfDOFs(VERTEX);
398
399
    nDOFEl += getGeo(VERTEX) * nDOF[VERTEX];

400
    if (dim > 1) {
401
      localAdmin->setNumberOfPreDOFs(EDGE,nDOF[EDGE]);
402
      nDOF[EDGE] += localAdmin->getNumberOfDOFs(EDGE);
403
404
405
406
407
408
409
      nDOFEl += getGeo(EDGE) * nDOF[EDGE];
    }

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

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

412
413
    node[VERTEX] = 0;
    nNodeEl = getGeo(VERTEX);
414

415
416
    if (dim > 1) {
      node[EDGE] = nNodeEl;
417
418
      if (nDOF[EDGE] > 0) 
	nNodeEl += getGeo(EDGE);
419
420
    }

421
    if (dim == 3) {
422
      localAdmin->setNumberOfPreDOFs(FACE,nDOF[FACE]);
423
424
425
426
427
      nDOF[FACE] += localAdmin->getNumberOfDOFs(FACE);
      nDOFEl += getGeo(FACE) * nDOF[FACE];
      node[FACE] = nNodeEl;
      if (nDOF[FACE] > 0) 
	nNodeEl += getGeo(FACE);
428
429
    }

430
    node[CENTER] = nNodeEl;
431
    if (nDOF[CENTER] > 0)
432
      nNodeEl += 1;
433
434
435
436
  }

  void Mesh::dofCompress()
  {
437
438
    FUNCNAME("Mesh::dofCompress()");
    Flag fill_flag;
439

Thomas Witkowski's avatar
Thomas Witkowski committed
440
    for (int iadmin = 0; iadmin < static_cast<int>(admin.size()); iadmin++) {
441
442
443
      compressAdmin = admin[iadmin];

      TEST_EXIT_DBG(compressAdmin)("no admin[%d] in mesh\n", iadmin);
444
      
Thomas Witkowski's avatar
Thomas Witkowski committed
445
446
      int size = compressAdmin->getSize();
      if (size < 1) 
447
	continue;
Thomas Witkowski's avatar
Thomas Witkowski committed
448

449
450
      if (compressAdmin->getUsedDOFs() < 1)    
	continue;
Thomas Witkowski's avatar
Thomas Witkowski committed
451

452
453
      if (compressAdmin->getHoleCount() < 1)    
	continue;
Thomas Witkowski's avatar
Thomas Witkowski committed
454
  
455
456
457
458
      newDOF.resize(size);
      
      compressAdmin->compress(newDOF);
      
459
      if (preserveCoarseDOFs)
460
	fill_flag = Mesh::CALL_EVERY_EL_PREORDER | Mesh::FILL_NOTHING;
461
      else
462
463
464
	fill_flag = Mesh::CALL_LEAF_EL | Mesh::FILL_NOTHING;
      
      
465
466
467
468
469
470
471
472
473
474
475
476
477
      TraverseStack stack;
      ElInfo *elInfo = stack.traverseFirst(this, -1, fill_flag);
      while (elInfo) {
	newDOFFct1(elInfo);
	elInfo = stack.traverseNext(elInfo);
      }

      elInfo = stack.traverseFirst(this, -1, fill_flag);
      while (elInfo) {
	newDOFFct2(elInfo);
	elInfo = stack.traverseNext(elInfo);
      }

478
      newDOF.resize(0);
479
    }       
480
481
482
483
484
  }


  DegreeOfFreedom *Mesh::getDOF(GeoIndex position)
  {
485
    FUNCNAME("Mesh::getDOF()");
486

487
    TEST_EXIT_DBG(position >= CENTER && position <= FACE)
488
      ("unknown position %d\n", position);
489

490
491
    int ndof = getNumberOfDOFs(position);
    if (ndof <= 0) 
492
      return NULL;
493

494
    DegreeOfFreedom *dof = new DegreeOfFreedom[ndof];
495

496
497
    for (int i = 0; i < getNumberOfDOFAdmin(); i++) {
      const DOFAdmin *localAdmin = &getDOFAdmin(i);
498
      TEST_EXIT_DBG(localAdmin)("no admin[%d]\n", i);
499
500
501
502
      
      int n  = localAdmin->getNumberOfDOFs(position);
      int n0 = localAdmin->getNumberOfPreDOFs(position);
      
503
      TEST_EXIT_DBG(n + n0 <= ndof)("n=%d, n0=%d too large: ndof=%d\n", n, n0, ndof);
504
      
505
      for (int j = 0; j < n; j++)
506
507
	dof[n0 + j] = const_cast<DOFAdmin*>(localAdmin)->getDOFIndex();
    }
508
 
509
    return dof;
510
511
512
513
514
  }


  DegreeOfFreedom **Mesh::createDOFPtrs()
  {
515
    FUNCNAME("Mesh::createDOFPtrs()");
516
517

    if (nNodeEl <= 0)
518
      return NULL;
519

520
    DegreeOfFreedom **ptrs = new DegreeOfFreedom*[nNodeEl];
521
    for (int i = 0; i < nNodeEl; i++)
522
523
      ptrs[i] = NULL;

524
    return ptrs;
525
526
527
528
  }

  void Mesh::freeDOFPtrs(DegreeOfFreedom **ptrs)
  {
529
    FUNCNAME("Mesh::freeDOFPtrs()");
530

531
    TEST_EXIT_DBG(ptrs)("ptrs=NULL\n");
532
533
534
535

    if (nNodeEl <= 0)
      return;
  
536
    delete [] ptrs;
537
538
539
  }


Thomas Witkowski's avatar
Thomas Witkowski committed
540
  const DOFAdmin *Mesh::createDOFAdmin(std::string lname, DimVec<int> lnDOF)
541
  {
542
    FUNCNAME("Mesh::createDOFAdmin()");
543

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

Thomas Witkowski's avatar
Thomas Witkowski committed
546
547
    for (int i = 0; i < dim + 1; i++)
      localAdmin->setNumberOfDOFs(i, lnDOF[i]);
548
549
550

    addDOFAdmin(localAdmin);

551
    return localAdmin;
552
553
554
555
556
557
558
  }


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

559
560
561
562
563
564
    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];
565
      }
566
567
    }

568
    return localAdmin;
569
570
571
572
  }

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

575
    TEST_EXIT_DBG(position >= CENTER && position <= FACE)
576
      ("unknown position %d\n", position);
577

578
579
580
    int ndof = nDOF[position];
    if (ndof) {
      if (!dof) {
581
	MSG("dof = NULL, but ndof = %d\n", ndof);
582
583
	return;
      }
584
    } else  {
585
586
587
      if (dof)
	MSG("dof != NULL, but ndof = 0\n");

588
589
      return;
    }
590

591
    TEST_EXIT_DBG(ndof <= MAX_DOF)
592
      ("ndof too big: ndof = %d, MAX_DOF = %d\n", ndof, MAX_DOF);
593

594
595
596
597
598
599
600
601
602
603
    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]);
    }
604

605
    delete [] dof;
606
607
608
609
610
  }

  void Mesh::freeElement(Element* el)
  {
    freeDOFPtrs(const_cast<DegreeOfFreedom**>(el->getDOF()));
Thomas Witkowski's avatar
Thomas Witkowski committed
611
    delete el;
612
613
614
615
616
617
  }


  Element* Mesh::createNewElement(Element *parent)
  {
    FUNCNAME("Mesh::createNewElement()");
618
619

    TEST_EXIT_DBG(elementPrototype)("no element prototype\n");
620
621
622

    Element *el = parent ? parent->clone() : elementPrototype->clone();
  
Thomas Witkowski's avatar
Thomas Witkowski committed
623
    if (!parent && elementDataPrototype)
624
      el->setElementData(elementDataPrototype->clone()); 
Thomas Witkowski's avatar
Thomas Witkowski committed
625
    else
626
      el->setElementData(NULL); // must be done in ElementData::refineElementData()
Thomas Witkowski's avatar
Thomas Witkowski committed
627
    
628
629
630
    return el;
  }

631

632
633
  ElInfo* Mesh::createNewElInfo()
  {
634
635
    FUNCNAME("Mesh::createNewElInfo()");

636
637
    switch(dim) {
    case 1:
Thomas Witkowski's avatar
Thomas Witkowski committed
638
      return new ElInfo1d(this);
639
640
      break;
    case 2:
Thomas Witkowski's avatar
Thomas Witkowski committed
641
      return new ElInfo2d(this);
642
643
      break;
    case 3:
Thomas Witkowski's avatar
Thomas Witkowski committed
644
      return new ElInfo3d(this);
645
646
647
648
      break;
    default:
      ERROR_EXIT("invalid dim\n");
      return NULL;
649
    }
650
651
  }

652

653
654
  bool Mesh::findElInfoAtPoint(const WorldVector<double>& xy,
			       ElInfo *el_info,
655
656
			       DimVec<double>& bary,
			       const MacroElement *start_mel,
657
			       const WorldVector<double> *xy0,
658
			       double *sp)
659
660
661
  {
    static const MacroElement *mel = NULL;
    DimVec<double> lambda(dim, NO_INIT);
Thomas Witkowski's avatar
Thomas Witkowski committed
662
    ElInfo *mel_info = createNewElInfo();
663
664
665
666

    if (start_mel != NULL)
      mel = start_mel;
    else
Thomas Witkowski's avatar
Thomas Witkowski committed
667
      if (mel == NULL || mel->getElement()->getMesh() != this)
668
669
670
	mel = *(macroElements.begin());

    mel_info->setFillFlag(Mesh::FILL_COORDS);
671
    g_xy = &xy;
672
    g_xy0 = xy0;
673
    g_sp = sp;
674
675
676

    mel_info->fillMacroInfo(mel);

677
678
679
680
681
682
683
    // We have the care about not to visite a macro element twice. In this case the
    // function would end up in an infinite loop. If a macro element is visited a 
    // second time, what can happen with periodic boundary conditions, the point is
    // not within the mesh!
    std::set<int> macrosVisited;
    macrosVisited.insert(mel->getIndex());

684
    int k;
685
686
687
    while ((k = mel_info->worldToCoord(xy, &lambda)) >= 0) {
      if (mel->getNeighbour(k)) {
	mel = mel->getNeighbour(k);
688
689
690
691
692
	mel = mel->getNeighbour(k);
	if (macrosVisited.count(mel->getIndex()))
	  return false;

	macrosVisited.insert(mel->getIndex());
693
694
695
696
697
698
699
	mel_info->fillMacroInfo(mel);
	continue;
      }
      break;
    }

    /* now, descend in tree to find leaf element at point */
700
    bool inside = findElementAtPointRecursive(mel_info, lambda, k, el_info);
701
702
    for (int i = 0; i <= dim; i++)
      bary[i] = final_lambda[i];   
703
  
Thomas Witkowski's avatar
Thomas Witkowski committed
704
    delete mel_info;
705

706
    return inside;
707
708
709
  }

  bool Mesh::findElementAtPoint(const WorldVector<double>&  xy,
710
711
				Element **elp, 
				DimVec<double>& bary,
712
				const MacroElement *start_mel,
713
714
				const WorldVector<double> *xy0,
				double *sp)
715
  {
716
717
    ElInfo *el_info = createNewElInfo();
    int val = findElInfoAtPoint(xy, el_info, bary, start_mel, xy0, sp);
718
719
720

    *elp = el_info->getElement();

Thomas Witkowski's avatar
Thomas Witkowski committed
721
    delete el_info;
722

723
    return val;
724
725
  }

726
  bool Mesh::findElementAtPointRecursive(ElInfo *el_info,
727
					 const DimVec<double>& lambda,
728
					 int outside,
729
730
					 ElInfo* final_el_info)
  {
731
    FUNCNAME("Mesh::findElementAtPointRecursive()");
Thomas Witkowski's avatar
Thomas Witkowski committed
732

733
734
    Element *el = el_info->getElement();
    DimVec<double> c_lambda(dim, NO_INIT);
735
736
    int inside;
    int ichild, c_outside;
737
738
739
740

    if (el->isLeaf()) {
      *final_el_info = *el_info;
      if (outside < 0) {
Thomas Witkowski's avatar
Thomas Witkowski committed
741
	for (int i = 0; i <= dim; i++)
742
	  final_lambda[i] = lambda[i];
Thomas Witkowski's avatar
Thomas Witkowski committed
743
	
744
	return true;
745
746
747
748
      }  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]);
Thomas Witkowski's avatar
Thomas Witkowski committed
749
750
751
752
	  for (int i = 0; i <= dim; i++) 
	    final_lambda[i] = s * c_lambda[i] + (1.0 - s) * lambda[i];
	  
	  if (g_sp)
753
	    *(g_sp) = s;
Thomas Witkowski's avatar
Thomas Witkowski committed
754
	  
755
756
	  if (dim == 3) 
	    MSG("outside finest level on el %d: s=%.3e\n", el->getIndex(), s);
757

758
	  return false;  /* ??? */
759
	} else {
760
	  return false;
761
	}
762
      }
763
764
    }

765
    ElInfo *c_el_info = createNewElInfo();
766

767
    if (dim == 1) {
768
769
770
771
      if (lambda[0] >= lambda[1]) {
	c_el_info->fillElInfo(0, el_info);
	if (outside >= 0) {
	  outside = el_info->worldToCoord(*(g_xy), &c_lambda);
Thomas Witkowski's avatar
Thomas Witkowski committed
772
	  TEST_EXIT(outside == 0)("point outside domain\n");
773
774
775
776
777
778
779
780
	} 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);
Thomas Witkowski's avatar
Thomas Witkowski committed
781
	  TEST_EXIT(outside == 0)("point outside domain\n");
782
783
784
785
786
787
788
	} else {
	  c_lambda[1] = lambda[1] - lambda[0];
	  c_lambda[0] = 2.0 * lambda[0];
	}
      }
    } /* DIM == 1 */

789
    if (dim == 2) {
790
791
792
793
      if (lambda[0] >= lambda[1]) {
	c_el_info->fillElInfo(0, el_info);
	if (el->isNewCoordSet()) {
	  outside = c_el_info->worldToCoord(*(g_xy), &c_lambda);
Thomas Witkowski's avatar
Thomas Witkowski committed
794
	  TEST_EXIT(outside == 0)("outside curved boundary child 0\n");	  
795
796
797
798
799
800
801
802
803
	} 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);
Thomas Witkowski's avatar
Thomas Witkowski committed
804
	  TEST_EXIT(outside == 0)("outside curved boundary child 1\n");	  
805
806
807
808
809
810
811
812
	} else {
	  c_lambda[0] = lambda[1] - lambda[0];
	  c_lambda[1] = lambda[2];
	  c_lambda[2] = 2.0 * lambda[0];
	}
      }
    } /* DIM == 2 */

813
    if (dim == 3) {
814
815
816
817
818
819
820
821
822
823
824
825
826
      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);
	  ElInfo *c_el_info2 = createNewElInfo();

	  c_el_info2->fillElInfo(1-ichild, el_info);
Thomas Witkowski's avatar
Thomas Witkowski committed
827
	  int c_outside2 = c_el_info2->worldToCoord(*(g_xy), &c_lambda2);
828
829

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

	  if ((c_outside2 < 0) || (c_lambda2[c_outside2] > c_lambda[c_outside])) {
836
837
	    for (int i = 0; i <= dim; i++) 
	      c_lambda[i] = c_lambda2[i];	    
838
839
840
841
	    c_outside = c_outside2;
	    *c_el_info = *c_el_info2;
	    ichild = 1 - ichild;
	  }
Thomas Witkowski's avatar
Thomas Witkowski committed
842
	  delete c_el_info2;
843
844
845
846
847
848
	}
	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];
849
850
851
852
853
854
	  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]];
855
856
857
858
	  c_lambda[3] = 2.0 * lambda[1];
	} else {
	  c_el_info->fillElInfo(1, el_info);
	  c_lambda[0] = lambda[1] - lambda[0];
859
860
861
862
863
864
	  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]];
865
866
867
868
869
	  c_lambda[3] = 2.0 * lambda[0];
	}
      }
    }  /* DIM == 3 */

Thomas Witkowski's avatar
Thomas Witkowski committed
870
    inside = findElementAtPointRecursive(c_el_info, c_lambda, outside, final_el_info);
Thomas Witkowski's avatar
Thomas Witkowski committed
871
    delete c_el_info;
872

873
    return inside; 
874
875
  }

876
877

  bool Mesh::getDofIndexCoords(DegreeOfFreedom dof, 
878
879
880
881
882
883
			       const FiniteElemSpace* feSpace,
			       WorldVector<double>& coords)
  {
    DimVec<double>* baryCoords;
    bool found = false;
    TraverseStack stack;
884
    std::vector<DegreeOfFreedom> dofVec(feSpace->getBasisFcts()->getNumber());
Thomas Witkowski's avatar
Thomas Witkowski committed
885

886
    ElInfo *elInfo = stack.traverseFirst(this, -1, 
Thomas Witkowski's avatar
Thomas Witkowski committed
887
					 Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS);    
888
    while (elInfo) {
889
890
891
      feSpace->getBasisFcts()->getLocalIndices(elInfo->getElement(),
					       feSpace->getAdmin(),
					       dofVec);
Thomas Witkowski's avatar
Thomas Witkowski committed
892
      for (int i = 0; i < feSpace->getBasisFcts()->getNumber(); i++) {
893
	if (dofVec[i] == dof) {
894
895
896
	  baryCoords = feSpace->getBasisFcts()->getCoords(i);
	  elInfo->coordToWorld(*baryCoords, coords);
	  found = true;
Thomas Witkowski's avatar
Thomas Witkowski committed
897
	  break;	  
898
	}
Thomas Witkowski's avatar
Thomas Witkowski committed
899
900
      }
      
901
902
      if (found)
	break;
Thomas Witkowski's avatar
Thomas Witkowski committed
903
      
904
905
906
907
908
      elInfo = stack.traverseNext(elInfo);
    }

    return found;
  }
909

910
911
912

  void Mesh::getDofIndexCoords(const FiniteElemSpace* feSpace,
			       DOFVector<WorldVector<double> >& coords)
913
  {
914
915
916
    const BasisFunction* basFcts = feSpace->getBasisFcts();
    int nBasFcts = basFcts->getNumber();
    std::vector<DegreeOfFreedom> dofVec(nBasFcts);
917

918
919
920
    TraverseStack stack;
    ElInfo *elInfo = 
      stack.traverseFirst(this, -1, Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS);
921
    while (elInfo) {
922
923
924
925
      basFcts->getLocalIndices(elInfo->getElement(), feSpace->getAdmin(), dofVec);
      for (int i = 0; i < nBasFcts; i++) {
	DimVec<double> *baryCoords = basFcts->getCoords(i);
	elInfo->coordToWorld(*baryCoords, coords[dofVec[i]]);
926
      }
927
         
928
929
930
931
932
      elInfo = stack.traverseNext(elInfo);
    }

  }

933

Thomas Witkowski's avatar
Thomas Witkowski committed
934
935
936
937
  void Mesh::setDiameter(const WorldVector<double>& w) 
  { 
    diam = w; 
  }
938

Thomas Witkowski's avatar
Thomas Witkowski committed
939
940
941
942
  void Mesh::setDiameter(int i, double w) 
  { 
    diam[i] = w; 
  }
943

944
945
  int Mesh::newDOFFct1(ElInfo* ei) 
  {
946
947
948
949
    ei->getElement()->newDOFFct1(compressAdmin);
    return 0;
  }

950
951
  int Mesh::newDOFFct2(ElInfo* ei) 
  {
952
953
954
955
    ei->getElement()->newDOFFct2(compressAdmin);
    return 0;
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
956
  void Mesh::serialize(std::ostream &out)
957
958
959
  {
    serializedDOFs.clear();

Thomas Witkowski's avatar
Thomas Witkowski committed
960
    out << name << "\n";
961

962
963
964
965
966
967
968
    SerUtil::serialize(out, dim);
    SerUtil::serialize(out, nVertices);
    SerUtil::serialize(out, nEdges);
    SerUtil::serialize(out, nLeaves);
    SerUtil::serialize(out, nElements);
    SerUtil::serialize(out, nFaces);
    SerUtil::serialize(out, maxEdgeNeigh);
969
970
971

    diam.serialize(out);

972
973
    SerUtil::serialize(out, preserveCoarseDOFs);
    SerUtil::serialize(out, nDOFEl);
974
975
976

    nDOF.serialize(out);

977
    SerUtil::serialize(out, nNodeEl);
978
979
980

    node.serialize(out);

981
982
983

    // === Write admins. ===

984
    int size = static_cast<int>(admin.size());
985
    SerUtil::serialize(out, size);
986
    for (int i = 0; i < size; i++)
987
988
      admin[i]->serialize(out);

989
990
991

    // === Write macroElements. ===

992
    size = static_cast<int>(macroElements.size());
993
    SerUtil::serialize(out, size);
994
    for (int i = 0; i < size; i++)
995
996
      macroElements[i]->serialize(out);

997
998
    SerUtil::serialize(out, elementIndex);
    SerUtil::serialize(out, initialized);
999

1000
1001
1002
1003
1004
1005
1006
1007
1008
1009
1010
1011
1012
1013
1014
1015
1016
1017
1018
1019
1020
1021
1022
1023

    // === Write periodic associations. ===
    int mapSize = periodicAssociations.size();
    SerUtil::serialize(out, mapSize);
    for (std::map<BoundaryType, VertexVector*>::iterator it = periodicAssociations.begin();
	 it != periodicAssociations.end(); ++it) {
      BoundaryType b = it->first;

      // Check which DOFAdmin is used for the current VertexVector we want to serialize.
      int ithAdmin = -1;
      for (int i = 0; i < static_cast<int>(admin.size()); i++) {
	if (admin[i] == it->second->getAdmin()) {
	  ithAdmin = i;
	  break;
	}
      }
      TEST_EXIT(ithAdmin >= 0)
	("No DOFAdmin found for serialization of periodic associations!\n");

      SerUtil::serialize(out, b);
      SerUtil::serialize(out, ithAdmin);
      it->second->serialize(out);
    }

1024
1025
1026
    serializedDOFs.clear();
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
1027
  void Mesh::deserialize(std::istream &in)
1028
  {
1029
1030
    FUNCNAME("Mesh::deserialize()");

1031
1032
1033
1034
1035
1036
    serializedDOFs.clear();

    in >> name;
    in.get();

    int oldVal = dim;
1037
1038
    SerUtil::deserialize(in, dim);
    TEST_EXIT_DBG(oldVal == 0 || dim == oldVal)("Invalid dimension!\n");
1039

1040
1041
1042
1043
1044
1045
    SerUtil::deserialize(in, nVertices);
    SerUtil::deserialize(in, nEdges);
    SerUtil::deserialize(in, nLeaves);
    SerUtil::deserialize(in, nElements);
    SerUtil::deserialize(in, nFaces);
    SerUtil::deserialize(in, maxEdgeNeigh);
1046
1047
1048

    diam.deserialize(in);

1049
    SerUtil::deserialize(in, preserveCoarseDOFs);
1050
1051

    oldVal = nDOFEl;
1052
1053
    SerUtil::deserialize(in, nDOFEl);
    TEST_EXIT_DBG(oldVal == 0 || nDOFEl == oldVal)("Invalid nDOFEl!\n");
1054
1055
1056
1057

    nDOF.deserialize(in);

    oldVal = nNodeEl;
1058
1059
    SerUtil::deserialize(in, nNodeEl);
    TEST_EXIT_DBG(oldVal == 0 || nNodeEl == oldVal)("Invalid nNodeEl!\n");
1060
1061
1062

    node.deserialize(in);

1063
1064
1065

    // === Read admins. ===

1066
    int size;
1067
    SerUtil::deserialize(in, size);
1068
    admin.resize(size, NULL);
1069
    for (int i = 0; i < size; i++) {
1070
      if (!admin[i])
Thomas Witkowski's avatar
Thomas Witkowski committed
1071
	admin[i] = new DOFAdmin(this);