Mesh.cc 34.1 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

namespace AMDiS {

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

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

36
37
38
39
40
41
42
43
44
45
46
47
48
49
  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;


  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

  //**************************************************************************
  //  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;
61
62
  const Flag Mesh::CALL_MG_LEVEL           = 0X4000L;  
  const Flag Mesh::CALL_REVERSE_MODE       = 0X8000L;
63
64


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

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


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

120
  Mesh::~Mesh()
121
  {
122
    deleteMeshStructure();
123

124
125
126
127
128
129
    if (macroFileInfo != NULL)
      delete macroFileInfo;    
    if (elementPrototype)
      delete elementPrototype;    
    if (elementDataPrototype)
      delete elementDataPrototype;    
130
    
131
132
    for (int i = 0; i < static_cast<int>(admin.size()); i++)
      delete admin[i];    
133
  }
134

135

136
  Mesh& Mesh::operator=(const Mesh& m)
137
  {
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
    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
168
      admin[i] = new DOFAdmin(this);
169
      *(admin[i]) = *(m.admin[i]);
170
171
      admin[i]->setMesh(this);
    }
172

173

174
    /* ====================== Copy macro elements =================== */
175
  
176
177
178
179
180
181
182
183
184
    // 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;
185

186
187
    macroElements.clear();

188
189
190
    // 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();
191
	 it != m.macroElements.end(); ++it, insertCounter++) {
192

193
      // Create new MacroElement.
Thomas Witkowski's avatar
Thomas Witkowski committed
194
      MacroElement *el = new MacroElement(dim);
195

196
197
198
199
200
201
202
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
      // 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;
  }

239

240
241
242
243
244
245
246
247
248
249
250
251
  void Mesh::updateNumberOfLeaves()
  {
    nLeaves = 0;

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

252

253
254
255
256
257
258
  void Mesh::addMacroElement(MacroElement* me) 
  {
    macroElements.push_back(me); 
    me->setIndex(macroElements.size());
  }

259

260
  void Mesh::removeMacroElements(std::set<MacroElement*>& macros,
Thomas Witkowski's avatar
Thomas Witkowski committed
261
				 const FiniteElemSpace *feSpace) 
262
263
264
  {
    FUNCNAME("Mesh::removeMacroElement()");

Thomas Witkowski's avatar
Thomas Witkowski committed
265
266
267
268
269
270
271
    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);
272

273
    // Map that stores for each dof pointer (which may have a list of dofs)
274
    // all macro element indices that own this dof.
Thomas Witkowski's avatar
Thomas Witkowski committed
275
276
    DofElMap dofsOwner;
    DofPosMap dofsPosIndex;
277
278
279
280


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

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

290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
    // === 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
309

310
311
312
313
314
    // === 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
315

316
317
318
319
      // 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
320
321
	  for (int j = 0; j <= dim; j++)
	    if ((*macroIt)->getNeighbour(i)->getNeighbour(j) == *macroIt)
322
323
324
325
326
	      (*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--;
327
	}
328
329
      }

330
      // Decrease also the number of elements of the mesh.
331
332
      nLeaves--;
      nElements--;
333
    }     
334
335


336
337
    // === Check now all the dofs, that have no owner anymore and therefore have ===
    // === to be removed.                                                        ===
338
339

    int nRemainDofs = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
340
    for (DofElMap::iterator dofsIt = dofsOwner.begin(); 
341
342
343
344
345
346
347
348
349
350
351
352
353
354
	 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)
355
	freeDOF(const_cast<DegreeOfFreedom*>(dofsIt->first), 
Thomas Witkowski's avatar
Thomas Witkowski committed
356
357
		dofsPosIndex[dofsIt->first]);
      else
358
	nRemainDofs++;
359
360
    }

361
362
363
364
365
366
367

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

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

368
    nVertices = nRemainDofs;
369
  }
370

371

372
373
374
375
376
377
  void Mesh::addDOFAdmin(DOFAdmin *localAdmin)
  {    
    FUNCNAME("Mesh::addDOFAdmin()");

    localAdmin->setMesh(this);

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

Thomas Witkowski's avatar
Thomas Witkowski committed
381
382
383
    TEST_EXIT(dai == admin.end())
      ("admin %s is already associated to mesh %s\n",
       localAdmin->getName().c_str(), this->getName().c_str());
384

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

388
389
390
391
392
    admin.push_back(localAdmin);

    nDOFEl = 0;

    localAdmin->setNumberOfPreDOFs(VERTEX,nDOF[VERTEX]);
393
    nDOF[VERTEX] += localAdmin->getNumberOfDOFs(VERTEX);
394
395
    nDOFEl += getGeo(VERTEX) * nDOF[VERTEX];

396
    if (dim > 1) {
397
      localAdmin->setNumberOfPreDOFs(EDGE,nDOF[EDGE]);
398
      nDOF[EDGE] += localAdmin->getNumberOfDOFs(EDGE);
399
400
401
402
403
404
405
      nDOFEl += getGeo(EDGE) * nDOF[EDGE];
    }

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

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

408
409
    node[VERTEX] = 0;
    nNodeEl = getGeo(VERTEX);
410

411
412
    if (dim > 1) {
      node[EDGE] = nNodeEl;
413
414
      if (nDOF[EDGE] > 0) 
	nNodeEl += getGeo(EDGE);
415
416
    }

417
    if (dim == 3) {
418
      localAdmin->setNumberOfPreDOFs(FACE,nDOF[FACE]);
419
420
421
422
423
      nDOF[FACE] += localAdmin->getNumberOfDOFs(FACE);
      nDOFEl += getGeo(FACE) * nDOF[FACE];
      node[FACE] = nNodeEl;
      if (nDOF[FACE] > 0) 
	nNodeEl += getGeo(FACE);
424
425
    }

426
    node[CENTER] = nNodeEl;
427
    if (nDOF[CENTER] > 0)
428
      nNodeEl += 1;
429
430
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
431

432
433
  void Mesh::dofCompress()
  {
434
    FUNCNAME("Mesh::dofCompress()");
435

Thomas Witkowski's avatar
Thomas Witkowski committed
436
437
    for (unsigned int iadmin = 0; iadmin < admin.size(); iadmin++) {
      DOFAdmin* compressAdmin = admin[iadmin];
438
439

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

447
448
449
      newDOF.resize(size);
      
      compressAdmin->compress(newDOF);
Thomas Witkowski's avatar
Thomas Witkowski committed
450
451
452
453

      Flag fill_flag = (preserveCoarseDOFs ?  
			Mesh::CALL_EVERY_EL_PREORDER | Mesh::FILL_NOTHING :
			Mesh::CALL_LEAF_EL | Mesh::FILL_NOTHING);          
454
455
456
      TraverseStack stack;
      ElInfo *elInfo = stack.traverseFirst(this, -1, fill_flag);
      while (elInfo) {
Thomas Witkowski's avatar
Thomas Witkowski committed
457
	elInfo->getElement()->newDOFFct1(compressAdmin);
458
459
460
461
462
	elInfo = stack.traverseNext(elInfo);
      }

      elInfo = stack.traverseFirst(this, -1, fill_flag);
      while (elInfo) {
Thomas Witkowski's avatar
Thomas Witkowski committed
463
	elInfo->getElement()->newDOFFct2(compressAdmin);
464
465
466
	elInfo = stack.traverseNext(elInfo);
      }

467
      newDOF.resize(0);
468
    }       
469
470
471
472
473
  }


  DegreeOfFreedom *Mesh::getDOF(GeoIndex position)
  {
474
    FUNCNAME("Mesh::getDOF()");
475

476
    TEST_EXIT_DBG(position >= CENTER && position <= FACE)
477
      ("unknown position %d\n", position);
478

479
480
    int ndof = getNumberOfDOFs(position);
    if (ndof <= 0) 
481
      return NULL;
482

483
    DegreeOfFreedom *dof = new DegreeOfFreedom[ndof];
484

485
486
    for (int i = 0; i < getNumberOfDOFAdmin(); i++) {
      const DOFAdmin *localAdmin = &getDOFAdmin(i);
487
      TEST_EXIT_DBG(localAdmin)("no admin[%d]\n", i);
488
489
490
491
      
      int n  = localAdmin->getNumberOfDOFs(position);
      int n0 = localAdmin->getNumberOfPreDOFs(position);
      
492
      TEST_EXIT_DBG(n + n0 <= ndof)("n=%d, n0=%d too large: ndof=%d\n", n, n0, ndof);
493
      
494
      for (int j = 0; j < n; j++)
495
496
	dof[n0 + j] = const_cast<DOFAdmin*>(localAdmin)->getDOFIndex();
    }
497
 
498
    return dof;
499
500
501
502
503
  }


  DegreeOfFreedom **Mesh::createDOFPtrs()
  {
504
    FUNCNAME("Mesh::createDOFPtrs()");
505
506

    if (nNodeEl <= 0)
507
      return NULL;
508

509
    DegreeOfFreedom **ptrs = new DegreeOfFreedom*[nNodeEl];
510
    for (int i = 0; i < nNodeEl; i++)
511
512
      ptrs[i] = NULL;

513
    return ptrs;
514
515
  }

516

517
518
  void Mesh::freeDOFPtrs(DegreeOfFreedom **ptrs)
  {
519
    FUNCNAME("Mesh::freeDOFPtrs()");
520

521
    TEST_EXIT_DBG(ptrs)("ptrs=NULL\n");
522
523
524
525

    if (nNodeEl <= 0)
      return;
  
526
    delete [] ptrs;
527
528
529
  }


Thomas Witkowski's avatar
Thomas Witkowski committed
530
  const DOFAdmin *Mesh::createDOFAdmin(std::string lname, DimVec<int> lnDOF)
531
  {
532
    FUNCNAME("Mesh::createDOFAdmin()");
533

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

Thomas Witkowski's avatar
Thomas Witkowski committed
536
537
    for (int i = 0; i < dim + 1; i++)
      localAdmin->setNumberOfDOFs(i, lnDOF[i]);
538
539
540

    addDOFAdmin(localAdmin);

541
    return localAdmin;
542
543
544
545
546
547
548
  }


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

549
550
551
552
553
554
    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];
555
      }
556
557
    }

558
    return localAdmin;
559
560
  }

561

562
563
  void Mesh::freeDOF(DegreeOfFreedom* dof, GeoIndex position)
  {
564
    FUNCNAME("Mesh::freeDOF()");
565

566
    TEST_EXIT_DBG(position >= CENTER && position <= FACE)
567
      ("unknown position %d\n", position);
568

569
570
571
    int ndof = nDOF[position];
    if (ndof) {
      if (!dof) {
572
	MSG("dof = NULL, but ndof = %d\n", ndof);
573
574
	return;
      }
575
    } else  {
576
577
578
      if (dof)
	MSG("dof != NULL, but ndof = 0\n");

579
580
      return;
    }
581

582
    TEST_EXIT_DBG(ndof <= MAX_DOF)
583
      ("ndof too big: ndof = %d, MAX_DOF = %d\n", ndof, MAX_DOF);
584

585
586
587
588
589
590
591
592
593
594
    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]);
    }
595

596
    delete [] dof;
597
598
  }

599

600
601
602
  void Mesh::freeElement(Element* el)
  {
    freeDOFPtrs(const_cast<DegreeOfFreedom**>(el->getDOF()));
Thomas Witkowski's avatar
Thomas Witkowski committed
603
    delete el;
604
605
606
607
608
609
  }


  Element* Mesh::createNewElement(Element *parent)
  {
    FUNCNAME("Mesh::createNewElement()");
610
611

    TEST_EXIT_DBG(elementPrototype)("no element prototype\n");
612
613
614

    Element *el = parent ? parent->clone() : elementPrototype->clone();
  
Thomas Witkowski's avatar
Thomas Witkowski committed
615
    if (!parent && elementDataPrototype)
616
      el->setElementData(elementDataPrototype->clone()); 
Thomas Witkowski's avatar
Thomas Witkowski committed
617
    else
618
      el->setElementData(NULL); // must be done in ElementData::refineElementData()
Thomas Witkowski's avatar
Thomas Witkowski committed
619
    
620
621
622
    return el;
  }

623

624
625
  ElInfo* Mesh::createNewElInfo()
  {
626
627
    FUNCNAME("Mesh::createNewElInfo()");

628
629
    switch(dim) {
    case 1:
Thomas Witkowski's avatar
Thomas Witkowski committed
630
      return new ElInfo1d(this);
631
632
      break;
    case 2:
Thomas Witkowski's avatar
Thomas Witkowski committed
633
      return new ElInfo2d(this);
634
635
      break;
    case 3:
Thomas Witkowski's avatar
Thomas Witkowski committed
636
      return new ElInfo3d(this);
637
638
639
640
      break;
    default:
      ERROR_EXIT("invalid dim\n");
      return NULL;
641
    }
642
643
  }

644

645
646
  bool Mesh::findElInfoAtPoint(const WorldVector<double>& xy,
			       ElInfo *el_info,
647
648
			       DimVec<double>& bary,
			       const MacroElement *start_mel,
649
			       const WorldVector<double> *xy0,
650
			       double *sp)
651
652
653
  {
    static const MacroElement *mel = NULL;
    DimVec<double> lambda(dim, NO_INIT);
Thomas Witkowski's avatar
Thomas Witkowski committed
654
    ElInfo *mel_info = createNewElInfo();
655
656
657
658

    if (start_mel != NULL)
      mel = start_mel;
    else
Thomas Witkowski's avatar
Thomas Witkowski committed
659
      if (mel == NULL || mel->getElement()->getMesh() != this)
660
661
662
	mel = *(macroElements.begin());

    mel_info->setFillFlag(Mesh::FILL_COORDS);
663
    g_xy = &xy;
664
    g_xy0 = xy0;
665
    g_sp = sp;
666
667
668

    mel_info->fillMacroInfo(mel);

669
670
671
672
673
674
675
    // 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());

676
    int k;
677
678
    while ((k = mel_info->worldToCoord(xy, &lambda)) >= 0) {
      if (mel->getNeighbour(k)) {
679
680
681
682
683
	mel = mel->getNeighbour(k);
	if (macrosVisited.count(mel->getIndex()))
	  return false;

	macrosVisited.insert(mel->getIndex());
684
685
686
687
688
689
690
	mel_info->fillMacroInfo(mel);
	continue;
      }
      break;
    }

    /* now, descend in tree to find leaf element at point */
691
    bool inside = findElementAtPointRecursive(mel_info, lambda, k, el_info);
692
693
    for (int i = 0; i <= dim; i++)
      bary[i] = final_lambda[i];   
694
  
Thomas Witkowski's avatar
Thomas Witkowski committed
695
    delete mel_info;
696

697
    return inside;
698
699
700
  }

  bool Mesh::findElementAtPoint(const WorldVector<double>&  xy,
701
702
				Element **elp, 
				DimVec<double>& bary,
703
				const MacroElement *start_mel,
704
705
				const WorldVector<double> *xy0,
				double *sp)
706
  {
707
708
    ElInfo *el_info = createNewElInfo();
    int val = findElInfoAtPoint(xy, el_info, bary, start_mel, xy0, sp);
709
710
711

    *elp = el_info->getElement();

Thomas Witkowski's avatar
Thomas Witkowski committed
712
    delete el_info;
713

714
    return val;
715
716
  }

717
  bool Mesh::findElementAtPointRecursive(ElInfo *el_info,
718
					 const DimVec<double>& lambda,
719
					 int outside,
720
721
					 ElInfo* final_el_info)
  {
722
    FUNCNAME("Mesh::findElementAtPointRecursive()");
Thomas Witkowski's avatar
Thomas Witkowski committed
723

724
725
    Element *el = el_info->getElement();
    DimVec<double> c_lambda(dim, NO_INIT);
726
727
    int inside;
    int ichild, c_outside;
728
729
730
731

    if (el->isLeaf()) {
      *final_el_info = *el_info;
      if (outside < 0) {
Thomas Witkowski's avatar
Thomas Witkowski committed
732
	for (int i = 0; i <= dim; i++)
733
	  final_lambda[i] = lambda[i];
Thomas Witkowski's avatar
Thomas Witkowski committed
734
	
735
	return true;
736
737
738
739
      }  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
740
741
742
743
	  for (int i = 0; i <= dim; i++) 
	    final_lambda[i] = s * c_lambda[i] + (1.0 - s) * lambda[i];
	  
	  if (g_sp)
744
	    *(g_sp) = s;
Thomas Witkowski's avatar
Thomas Witkowski committed
745
	  
746
747
	  if (dim == 3) 
	    MSG("outside finest level on el %d: s=%.3e\n", el->getIndex(), s);
748

749
	  return false;  /* ??? */
750
	} else {
751
	  return false;
752
	}
753
      }
754
755
    }

756
    ElInfo *c_el_info = createNewElInfo();
757

758
    if (dim == 1) {
759
760
761
762
      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
763
	  TEST_EXIT(outside == 0)("point outside domain\n");
764
765
766
767
768
769
770
771
	} 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
772
	  TEST_EXIT(outside == 0)("point outside domain\n");
773
774
775
776
777
778
779
	} else {
	  c_lambda[1] = lambda[1] - lambda[0];
	  c_lambda[0] = 2.0 * lambda[0];
	}
      }
    } /* DIM == 1 */

780
    if (dim == 2) {
781
782
783
784
      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
785
	  TEST_EXIT(outside == 0)("outside curved boundary child 0\n");	  
786
787
788
789
790
791
792
793
794
	} 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
795
	  TEST_EXIT(outside == 0)("outside curved boundary child 1\n");	  
796
797
798
799
800
801
802
803
	} else {
	  c_lambda[0] = lambda[1] - lambda[0];
	  c_lambda[1] = lambda[2];
	  c_lambda[2] = 2.0 * lambda[0];
	}
      }
    } /* DIM == 2 */

804
    if (dim == 3) {
805
806
807
808
809
810
811
812
813
814
815
816
817
      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
818
	  int c_outside2 = c_el_info2->worldToCoord(*(g_xy), &c_lambda2);
819
820

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
861
    inside = findElementAtPointRecursive(c_el_info, c_lambda, outside, final_el_info);
Thomas Witkowski's avatar
Thomas Witkowski committed
862
    delete c_el_info;
863

864
    return inside; 
865
866
  }

867
868

  bool Mesh::getDofIndexCoords(DegreeOfFreedom dof, 
869
870
871
872
873
874
			       const FiniteElemSpace* feSpace,
			       WorldVector<double>& coords)
  {
    DimVec<double>* baryCoords;
    bool found = false;
    TraverseStack stack;
875
    std::vector<DegreeOfFreedom> dofVec(feSpace->getBasisFcts()->getNumber());
Thomas Witkowski's avatar
Thomas Witkowski committed
876

877
    ElInfo *elInfo = stack.traverseFirst(this, -1, 
Thomas Witkowski's avatar
Thomas Witkowski committed
878
					 Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS);    
879
    while (elInfo) {
880
881
882
      feSpace->getBasisFcts()->getLocalIndices(elInfo->getElement(),
					       feSpace->getAdmin(),
					       dofVec);
Thomas Witkowski's avatar
Thomas Witkowski committed
883
      for (int i = 0; i < feSpace->getBasisFcts()->getNumber(); i++) {
884
	if (dofVec[i] == dof) {
885
886
887
	  baryCoords = feSpace->getBasisFcts()->getCoords(i);
	  elInfo->coordToWorld(*baryCoords, coords);
	  found = true;
Thomas Witkowski's avatar
Thomas Witkowski committed
888
	  break;	  
889
	}
Thomas Witkowski's avatar
Thomas Witkowski committed
890
891
      }
      
892
893
      if (found)
	break;
Thomas Witkowski's avatar
Thomas Witkowski committed
894
      
895
896
897
898
899
      elInfo = stack.traverseNext(elInfo);
    }

    return found;
  }
900

901
902
903

  void Mesh::getDofIndexCoords(const FiniteElemSpace* feSpace,
			       DOFVector<WorldVector<double> >& coords)
904
  {
905
906
907
    const BasisFunction* basFcts = feSpace->getBasisFcts();
    int nBasFcts = basFcts->getNumber();
    std::vector<DegreeOfFreedom> dofVec(nBasFcts);
908

909
910
911
    TraverseStack stack;
    ElInfo *elInfo = 
      stack.traverseFirst(this, -1, Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS);
912
    while (elInfo) {
913
914
915
916
      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]]);
917
      }
918
         
919
920
921
922
923
      elInfo = stack.traverseNext(elInfo);
    }

  }

924

Thomas Witkowski's avatar
Thomas Witkowski committed
925
926
927
928
  void Mesh::setDiameter(const WorldVector<double>& w) 
  { 
    diam = w; 
  }
929

930

Thomas Witkowski's avatar
Thomas Witkowski committed
931
932
933
934
  void Mesh::setDiameter(int i, double w) 
  { 
    diam[i] = w; 
  }
935

936

Thomas Witkowski's avatar
Thomas Witkowski committed
937
  void Mesh::serialize(std::ostream &out)
938
939
940
  {
    serializedDOFs.clear();

Thomas Witkowski's avatar
Thomas Witkowski committed
941
    out << name << "\n";
942

943
944
945
946
947
948
949
    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);
950
951
952

    diam.serialize(out);

953
954
    SerUtil::serialize(out, preserveCoarseDOFs);
    SerUtil::serialize(out, nDOFEl);
955
956
957

    nDOF.serialize(out);

958
    SerUtil::serialize(out, nNodeEl);
959
960
961

    node.serialize(out);

962
963
964

    // === Write admins. ===

965
    int size = static_cast<int>(admin.size());
966
    SerUtil::serialize(out, size);
967
    for (int i = 0; i < size; i++)
968
969
      admin[i]->serialize(out);

970
971
972

    // === Write macroElements. ===

973
    size = static_cast<int>(macroElements.size());
974
    SerUtil::serialize(out, size);
975
    for (int i = 0; i < size; i++)
976
977
      macroElements[i]->serialize(out);

978
979
    SerUtil::serialize(out, elementIndex);
    SerUtil::serialize(out, initialized);
980

981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
1001
1002
1003
1004

    // === 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);
    }

1005
1006
1007
    serializedDOFs.clear();
  }

1008

Thomas Witkowski's avatar
Thomas Witkowski committed
1009
  void Mesh::deserialize(std::istream &in)
1010
  {
1011
1012
    FUNCNAME("Mesh::deserialize()");

1013
1014
1015
1016
1017
1018
    serializedDOFs.clear();

    in >> name;
    in.get();

    int oldVal = dim;
1019
1020
    SerUtil::deserialize(in, dim);
    TEST_EXIT_DBG(oldVal == 0 || dim == oldVal)("Invalid dimension!\n");
1021

1022
1023
1024
1025
1026
1027
    SerUtil::deserialize(in, nVertices);
    SerUtil::deserialize(in, nEdges);
    SerUtil::deserialize(in, nLeaves);
    SerUtil::deserialize(in, nElements);
    SerUtil::deserialize(in, nFaces);
    SerUtil::deserialize(in, maxEdgeNeigh);
1028
1029
1030

    diam.deserialize(in);

1031
    SerUtil::deserialize(in, preserveCoarseDOFs);
1032
1033

    oldVal = nDOFEl;
1034
1035
    SerUtil::deserialize(in, nDOFEl);
    TEST_EXIT_DBG(oldVal == 0 || nDOFEl == oldVal)("Invalid nDOFEl!\n");
1036
1037
1038
1039

    nDOF.deserialize(in);

    oldVal = nNodeEl;
1040
1041
    SerUtil::deserialize(in, nNodeEl);
    TEST_EXIT_DBG(oldVal == 0 || nNodeEl == oldVal)("Invalid nNodeEl!\n");
1042
1043
1044

    node.deserialize(in);

1045
1046
1047

    // === Read admins. ===

1048
    int size;
1049
    SerUtil::deserialize(in, size);