Mesh.cc 34.2 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
50
  const Flag Mesh::FILL_NOTHING       = 0X000L;
  const Flag Mesh::FILL_COORDS        = 0X001L;
  const Flag Mesh::FILL_BOUND         = 0X002L;
  const Flag Mesh::FILL_NEIGH         = 0X004L;
  const Flag Mesh::FILL_OPP_COORDS    = 0X008L;
  const Flag Mesh::FILL_ORIENTATION   = 0X010L;
  const Flag Mesh::FILL_DET           = 0X020L;
  const Flag Mesh::FILL_GRD_LAMBDA    = 0X040L;
  const Flag Mesh::FILL_LOCAL_INDICES = 0X080L;
  const Flag Mesh::FILL_ADD_ALL       = 0X100L;


  const Flag Mesh::FILL_ANY_1D = (0X01L|0X02L|0X04L|0X08L|0x20L|0X40L|0X100L);
  const Flag Mesh::FILL_ANY_2D = (0X01L|0X02L|0X04L|0X08L|0x20L|0X40L|0X100L);
  const Flag Mesh::FILL_ANY_3D = (0X01L|0X02L|0X04L|0X08L|0X10L|0x20L|0X40L|0X100L);
51
52
53
54
55
56
57
58
59
60
61
62
63

  //**************************************************************************
  //  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 

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

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


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

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

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

    elementPrototype->setIndex(-1);

115
116
    elementIndex = 0;
  }
117
118

  Mesh::~Mesh()
119
  {
120
    Element::deletedDOFs.clear();
121
122

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

128
    Element::deletedDOFs.clear();
129

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

  Mesh& Mesh::operator=(const Mesh& m)
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
168
169
170
171
172
    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
173
      admin[i] = new DOFAdmin(this);
174
      *(admin[i]) = *(m.admin[i]);
175
176
      admin[i]->setMesh(this);
    }
177

178

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

191
192
    macroElements.clear();

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

198
      // Create new MacroElement.
Thomas Witkowski's avatar
Thomas Witkowski committed
199
      MacroElement *el = new MacroElement(dim);
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
239
240
241
242
243
      // 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;
  }

244
245
246
247
248
249
250
251
252
253
254
255
  void Mesh::updateNumberOfLeaves()
  {
    nLeaves = 0;

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

256
257
258
259
260
261
  void Mesh::addMacroElement(MacroElement* me) 
  {
    macroElements.push_back(me); 
    me->setIndex(macroElements.size());
  }

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

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

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


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

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

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

312
313
314
315
316
    // === 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
317

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

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


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

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

363
364
365
366
367
368
369

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

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

370
    nVertices = nRemainDofs;
371
  }
372
373
374
375
376
377
378

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

    localAdmin->setMesh(this);

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

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

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

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

    nDOFEl = 0;

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

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

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

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

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

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
432

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

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

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

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

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

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

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


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

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

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

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

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


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

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

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

514
    return ptrs;
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
  }

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

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

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

578
579
      return;
    }
580

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

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

595
    delete [] dof;
596
597
598
599
600
  }

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


  Element* Mesh::createNewElement(Element *parent)
  {
    FUNCNAME("Mesh::createNewElement()");
608
609

    TEST_EXIT_DBG(elementPrototype)("no element prototype\n");
610
611
612

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

621

622
623
  ElInfo* Mesh::createNewElInfo()
  {
624
625
    FUNCNAME("Mesh::createNewElInfo()");

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

642

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

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

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

    mel_info->fillMacroInfo(mel);

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

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

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

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

695
    return inside;
696
697
698
  }

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

    *elp = el_info->getElement();

Thomas Witkowski's avatar
Thomas Witkowski committed
710
    delete el_info;
711

712
    return val;
713
714
  }

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

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

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

747
	  return false;  /* ??? */
748
	} else {
749
	  return false;
750
	}
751
      }
752
753
    }

754
    ElInfo *c_el_info = createNewElInfo();
755

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

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

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

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

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

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

862
    return inside; 
863
864
  }

865
866

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

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

    return found;
  }
898

899
900
901

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

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

  }

922

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
933
  void Mesh::serialize(std::ostream &out)
934
935
936
  {
    serializedDOFs.clear();

Thomas Witkowski's avatar
Thomas Witkowski committed
937
    out << name << "\n";
938

939
940
941
942
943
944
945
    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);
946
947
948

    diam.serialize(out);

949
950
    SerUtil::serialize(out, preserveCoarseDOFs);
    SerUtil::serialize(out, nDOFEl);
951
952
953

    nDOF.serialize(out);

954
    SerUtil::serialize(out, nNodeEl);
955
956
957

    node.serialize(out);

958
959
960

    // === Write admins. ===

961
    int size = static_cast<int>(admin.size());
962
    SerUtil::serialize(out, size);
963
    for (int i = 0; i < size; i++)
964
965
      admin[i]->serialize(out);

966
967
968

    // === Write macroElements. ===

969
    size = static_cast<int>(macroElements.size());
970
    SerUtil::serialize(out, size);
971
    for (int i = 0; i < size; i++)
972
973
      macroElements[i]->serialize(out);

974
975
    SerUtil::serialize(out, elementIndex);
    SerUtil::serialize(out, initialized);
976

977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000

    // === Write periodic associations. ===
    int mapSize = periodicAssociations.size();
    SerUtil::serialize(out, mapSize);
    for (std::map<BoundaryType, VertexVector*>::iterator