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
  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 is 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
  {
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
    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);
      
590
591
      TEST_EXIT_DBG(n + n0 <= ndof)
	("n = %d, n0 = %d too large: ndof = %d\n", n, n0, ndof);
592
593
      
      for (int j = 0; j < n; j++)
594
	localAdmin->freeDofIndex(dof[n0 + j]);      
595
    }
596

597
    delete [] dof;
598
599
  }

600

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


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

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

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

624

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

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

645

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

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

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

    mel_info->fillMacroInfo(mel);

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

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

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

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

698
    return inside;
699
700
701
  }

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

    *elp = el_info->getElement();

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

715
    return val;
716
717
  }

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

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

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

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

757
    ElInfo *c_el_info = createNewElInfo();
758

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

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

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

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

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

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

865
    return inside; 
866
867
  }

868
869

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

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

    return found;
  }
901

902
903
904

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

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

  }

925

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

931

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

937

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

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

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

    diam.serialize(out);

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

    nDOF.serialize(out);

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

    node.serialize(out);

963
964
965

    // === Write admins. ===

966
    int size =