Mesh.cc 39.9 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 "io/MacroReader.h"
#include "io/MacroInfo.h"
#include "io/MacroWriter.h"

11
12
13
14
#include "AdaptStationary.h"
#include "AdaptInstationary.h"
#include "FiniteElemSpace.h"
#include "ElementData.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
15
#include "ElementDofIterator.h"
16
17
18
19
20
21
22
23
24
25
26
#include "MacroElement.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 "PeriodicMap.h"
#include "Projection.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
27
#include "ElInfoStack.h"
28
#include "Serializer.h"
29
#include "Lagrange.h"
30
31
32
33
34
35
36
37
38

namespace AMDiS {

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

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

39
40
41
42
43
44
45
46
47
48
49
50
51
52
  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);
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;
64
65
  const Flag Mesh::CALL_MG_LEVEL           = 0X4000L;  
  const Flag Mesh::CALL_REVERSE_MODE       = 0X8000L;
66
67


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

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


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

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

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

    elementPrototype->setIndex(-1);

119
120
    elementIndex = 0;
  }
121

122

123
  Mesh::~Mesh()
124
  {
125
    deleteMeshStructure();
126

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

138

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

175

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

188
189
    macroElements.clear();

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

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

241

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

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

254

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

261

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
    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");

273
274
    // === Determine to all DOFs in mesh the macro elements where the DOF  ===
    // === is part of.                                                     ===
275

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

281
282
283
284
285
    ElementDofIterator elDofIter(feSpace);
    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(this, -1, Mesh::CALL_LEAF_EL);
    while (elInfo) {
      elDofIter.reset(elInfo->getElement());
Thomas Witkowski's avatar
Thomas Witkowski committed
286
      do {
287
	dofsOwner[elDofIter.getDofPtr()].insert(elInfo->getMacroElement());
Thomas Witkowski's avatar
Thomas Witkowski committed
288
289
	dofsPosIndex[elDofIter.getDofPtr()] = elDofIter.getPosIndex();
      } while (elDofIter.nextStrict());
290
291
292
293

      elInfo = stack.traverseNext(elInfo);
    }		   

Thomas Witkowski's avatar
Thomas Witkowski committed
294

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
314

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

322
323
      // Go through all neighbours of the macro element and remove this macro element
      // to be neighbour of some other macro element.
324
325
      for (int i = 0; i <= dim; i++)
	if ((*macroIt)->getNeighbour(i))
Thomas Witkowski's avatar
Thomas Witkowski committed
326
327
	  for (int j = 0; j <= dim; j++)
	    if ((*macroIt)->getNeighbour(i)->getNeighbour(j) == *macroIt)
328
329
	      (*macroIt)->getNeighbour(i)->setNeighbour(j, NULL);

330
331
332
333
      // Delete element hierarchie
      if (!(*macroIt)->getElement()->isLeaf()) {
	delete (*macroIt)->getElement()->getChild(0);
	delete (*macroIt)->getElement()->getChild(1);
334
335
336
337
338

	(*macroIt)->getElement()->child[0] = NULL;
	(*macroIt)->getElement()->child[1] = NULL;

	(*macroIt)->getElement()->setElementData(elementDataPrototype->clone()); 
339
      }
340
    }     
341
342


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

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

365

366
    // === Update number of elements, vertices, etc. ===
367

368
369
370
371
372
373
374
375
    nLeaves = 0;
    nElements = 0;

    std::set<const DegreeOfFreedom*> allVertices;

    elInfo = stack.traverseFirst(this, -1, Mesh::CALL_EVERY_EL_PREORDER);
    while (elInfo) {
      nElements++;
376

377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
      if (elInfo->getElement()->isLeaf()) {
	nLeaves++;

	for (int i = 0; i < getGeo(VERTEX); i++)
	  allVertices.insert(elInfo->getElement()->getDof(i));
      }

      elInfo = stack.traverseNext(elInfo);
    }

    nVertices = allVertices.size();

    // === Note: Although the macro elements are removed from the mesh,   ===
    // === they are not deleted from memory. The macro elements are still ===
    // === stored in macroInfo structure. They are needed, if the mesh is ===
    // === redistributed between the ranks.                               ===
393
  }
394

395

396
397
398
399
400
401
  void Mesh::addDOFAdmin(DOFAdmin *localAdmin)
  {    
    FUNCNAME("Mesh::addDOFAdmin()");

    localAdmin->setMesh(this);

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

Thomas Witkowski's avatar
Thomas Witkowski committed
405
406
407
    TEST_EXIT(dai == admin.end())
      ("admin %s is already associated to mesh %s\n",
       localAdmin->getName().c_str(), this->getName().c_str());
408

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

412
413
    admin.push_back(localAdmin);

414
    nDofEl = 0;
415

416
417
418
    localAdmin->setNumberOfPreDOFs(VERTEX,nDof[VERTEX]);
    nDof[VERTEX] += localAdmin->getNumberOfDofs(VERTEX);
    nDofEl += getGeo(VERTEX) * nDof[VERTEX];
419

420
    if (dim > 1) {
421
422
423
      localAdmin->setNumberOfPreDOFs(EDGE,nDof[EDGE]);
      nDof[EDGE] += localAdmin->getNumberOfDofs(EDGE);
      nDofEl += getGeo(EDGE) * nDof[EDGE];
424
425
    }

426
427
428
    localAdmin->setNumberOfPreDOFs(CENTER,nDof[CENTER]);
    nDof[CENTER]  += localAdmin->getNumberOfDofs(CENTER);
    nDofEl += nDof[CENTER];
429

430
    TEST_EXIT_DBG(nDof[VERTEX] > 0)("no vertex dofs\n");
431

432
433
    node[VERTEX] = 0;
    nNodeEl = getGeo(VERTEX);
434

435
436
    if (dim > 1) {
      node[EDGE] = nNodeEl;
437
      if (nDof[EDGE] > 0) 
438
	nNodeEl += getGeo(EDGE);
439
440
    }

441
    if (dim == 3) {
442
443
444
      localAdmin->setNumberOfPreDOFs(FACE,nDof[FACE]);
      nDof[FACE] += localAdmin->getNumberOfDofs(FACE);
      nDofEl += getGeo(FACE) * nDof[FACE];
445
      node[FACE] = nNodeEl;
446
      if (nDof[FACE] > 0) 
447
	nNodeEl += getGeo(FACE);
448
449
    }

450
    node[CENTER] = nNodeEl;
451
    if (nDof[CENTER] > 0)
452
      nNodeEl += 1;
453
454
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
455

456
457
  void Mesh::dofCompress()
  {
458
    FUNCNAME("Mesh::dofCompress()");
459

Thomas Witkowski's avatar
Thomas Witkowski committed
460
461
    for (unsigned int iadmin = 0; iadmin < admin.size(); iadmin++) {
      DOFAdmin* compressAdmin = admin[iadmin];
462
463

      TEST_EXIT_DBG(compressAdmin)("no admin[%d] in mesh\n", iadmin);
464
      
Thomas Witkowski's avatar
Thomas Witkowski committed
465
      int size = compressAdmin->getSize();
Thomas Witkowski's avatar
Thomas Witkowski committed
466
      if (size < 1 || 
467
	  compressAdmin->getUsedDofs() < 1 || 
Thomas Witkowski's avatar
Thomas Witkowski committed
468
	  compressAdmin->getHoleCount() < 1)    
469
	continue;
Thomas Witkowski's avatar
Thomas Witkowski committed
470

471
472
      std::vector<int> newDofIndex(size);     
      compressAdmin->compress(newDofIndex);
Thomas Witkowski's avatar
Thomas Witkowski committed
473
474
475
476

      Flag fill_flag = (preserveCoarseDOFs ?  
			Mesh::CALL_EVERY_EL_PREORDER | Mesh::FILL_NOTHING :
			Mesh::CALL_LEAF_EL | Mesh::FILL_NOTHING);          
477
478
479
      TraverseStack stack;
      ElInfo *elInfo = stack.traverseFirst(this, -1, fill_flag);
      while (elInfo) {
480
	elInfo->getElement()->newDofFct1(compressAdmin, newDofIndex);
481
482
483
484
485
	elInfo = stack.traverseNext(elInfo);
      }

      elInfo = stack.traverseFirst(this, -1, fill_flag);
      while (elInfo) {
486
	elInfo->getElement()->newDofFct2(compressAdmin);
487
488
	elInfo = stack.traverseNext(elInfo);
      }
489
    }       
490
491
492
  }


493
  DegreeOfFreedom *Mesh::getDof(GeoIndex position)
494
  {
495
    FUNCNAME("Mesh::getDof()");
496

497
    TEST_EXIT_DBG(position >= CENTER && position <= FACE)
498
      ("unknown position %d\n", position);
499

500
    int ndof = getNumberOfDofs(position);
501
    if (ndof <= 0) 
502
      return NULL;
503

504
    DegreeOfFreedom *dof = new DegreeOfFreedom[ndof];
505

506
    for (int i = 0; i < getNumberOfDOFAdmin(); i++) {
507
      const DOFAdmin *localAdmin = &getDofAdmin(i);
508
      TEST_EXIT_DBG(localAdmin)("no admin[%d]\n", i);
509
      
510
      int n  = localAdmin->getNumberOfDofs(position);
511
512
      int n0 = localAdmin->getNumberOfPreDOFs(position);
      
513
      TEST_EXIT_DBG(n + n0 <= ndof)("n=%d, n0=%d too large: ndof=%d\n", n, n0, ndof);
514
      
515
      for (int j = 0; j < n; j++)
516
517
	dof[n0 + j] = const_cast<DOFAdmin*>(localAdmin)->getDOFIndex();
    }
518
 
519
    return dof;
520
521
522
  }


Thomas Witkowski's avatar
Thomas Witkowski committed
523
  DegreeOfFreedom **Mesh::createDofPtrs()
524
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
525
    FUNCNAME("Mesh::createDofPtrs()");
526
527

    if (nNodeEl <= 0)
528
      return NULL;
529

530
    DegreeOfFreedom **ptrs = new DegreeOfFreedom*[nNodeEl];
531
    for (int i = 0; i < nNodeEl; i++)
532
533
      ptrs[i] = NULL;

534
    return ptrs;
535
536
  }

537

538
  void Mesh::freeDofPtrs(DegreeOfFreedom **ptrs)
539
  {
540
    FUNCNAME("Mesh::freeDofPtrs()");
541

542
    TEST_EXIT_DBG(ptrs)("ptrs is NULL!\n");
543
544
545
546

    if (nNodeEl <= 0)
      return;
  
547
    delete [] ptrs;
548
549
550
  }


551
  const DOFAdmin *Mesh::createDOFAdmin(std::string lname, DimVec<int> lnDof)
552
  {
553
    FUNCNAME("Mesh::createDOFAdmin()");
554

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

Thomas Witkowski's avatar
Thomas Witkowski committed
557
    for (int i = 0; i < dim + 1; i++)
558
      localAdmin->setNumberOfDOFs(i, lnDof[i]);
559
560
561

    addDOFAdmin(localAdmin);

562
    return localAdmin;
563
564
565
566
567
568
569
  }


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

570
    for (int i = 0; i < static_cast<int>(admin.size()); i++) {
571
      if (admin[i]->getNumberOfDofs(VERTEX)) {
572
573
574
575
	if (!localAdmin)  
	  localAdmin = admin[i];
	else if (admin[i]->getSize() < localAdmin->getSize())
	  localAdmin = admin[i];
576
      }
577
578
    }

579
    return localAdmin;
580
581
  }

582

583
  void Mesh::freeDof(DegreeOfFreedom* dof, GeoIndex position)
584
  {
585
    FUNCNAME("Mesh::freeDof()");
586

587
    TEST_EXIT_DBG(position >= CENTER && position <= FACE)
588
      ("unknown position %d\n", position);
589

590
    int ndof = nDof[position];
591
592
    if (ndof) {
      if (!dof) {
593
	MSG("dof = NULL, but ndof = %d\n", ndof);
594
595
	return;
      }
596
    } else  {
597
598
599
      if (dof)
	MSG("dof != NULL, but ndof = 0\n");

600
601
      return;
    }
602

603
    TEST_EXIT_DBG(ndof <= MAX_DOF)
604
      ("ndof too big: ndof = %d, MAX_DOF = %d\n", ndof, MAX_DOF);
605

606
    for (unsigned int i = 0; i < admin.size(); i++) {
607
      DOFAdmin *localAdmin = admin[i];
608
      int n = localAdmin->getNumberOfDofs(position);
609
610
      int n0 = localAdmin->getNumberOfPreDOFs(position);
      
611
612
      TEST_EXIT_DBG(n + n0 <= ndof)
	("n = %d, n0 = %d too large: ndof = %d\n", n, n0, ndof);
613
614
      
      for (int j = 0; j < n; j++)
615
	localAdmin->freeDofIndex(dof[n0 + j]);      
616
    }
617

618
    delete [] dof;
619
    dof = NULL;
620
621
  }

622

623
624
  void Mesh::freeElement(Element* el)
  {
625
    freeDofPtrs(const_cast<DegreeOfFreedom**>(el->getDof()));
Thomas Witkowski's avatar
Thomas Witkowski committed
626
    delete el;
627
628
629
630
631
632
  }


  Element* Mesh::createNewElement(Element *parent)
  {
    FUNCNAME("Mesh::createNewElement()");
633
634

    TEST_EXIT_DBG(elementPrototype)("no element prototype\n");
635
636

    Element *el = parent ? parent->clone() : elementPrototype->clone();
637

Thomas Witkowski's avatar
Thomas Witkowski committed
638
    if (!parent && elementDataPrototype)
639
      el->setElementData(elementDataPrototype->clone()); 
Thomas Witkowski's avatar
Thomas Witkowski committed
640
    else
641
      el->setElementData(NULL); // must be done in ElementData::refineElementData()
642

643
644
645
    return el;
  }

646

647
648
  ElInfo* Mesh::createNewElInfo()
  {
649
650
    FUNCNAME("Mesh::createNewElInfo()");

651
652
    switch(dim) {
    case 1:
Thomas Witkowski's avatar
Thomas Witkowski committed
653
      return new ElInfo1d(this);
654
655
      break;
    case 2:
Thomas Witkowski's avatar
Thomas Witkowski committed
656
      return new ElInfo2d(this);
657
658
      break;
    case 3:
Thomas Witkowski's avatar
Thomas Witkowski committed
659
      return new ElInfo3d(this);
660
661
662
663
      break;
    default:
      ERROR_EXIT("invalid dim\n");
      return NULL;
664
    }
665
666
  }

667

668
669
  bool Mesh::findElInfoAtPoint(const WorldVector<double>& xy,
			       ElInfo *el_info,
670
671
			       DimVec<double>& bary,
			       const MacroElement *start_mel,
672
			       const WorldVector<double> *xy0,
673
			       double *sp)
674
675
676
  {
    static const MacroElement *mel = NULL;
    DimVec<double> lambda(dim, NO_INIT);
Thomas Witkowski's avatar
Thomas Witkowski committed
677
    ElInfo *mel_info = createNewElInfo();
678
679
680
681

    if (start_mel != NULL)
      mel = start_mel;
    else
Thomas Witkowski's avatar
Thomas Witkowski committed
682
      if (mel == NULL || mel->getElement()->getMesh() != this)
683
684
685
	mel = *(macroElements.begin());

    mel_info->setFillFlag(Mesh::FILL_COORDS);
686
    g_xy = &xy;
687
    g_xy0 = xy0;
688
    g_sp = sp;
689
690
691

    mel_info->fillMacroInfo(mel);

692
693
694
695
696
697
698
    // 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());

699
    int k;
700
701
    while ((k = mel_info->worldToCoord(xy, &lambda)) >= 0) {
      if (mel->getNeighbour(k)) {
702
703
704
705
706
	mel = mel->getNeighbour(k);
	if (macrosVisited.count(mel->getIndex()))
	  return false;

	macrosVisited.insert(mel->getIndex());
707
708
709
710
711
712
713
	mel_info->fillMacroInfo(mel);
	continue;
      }
      break;
    }

    /* now, descend in tree to find leaf element at point */
714
    bool inside = findElementAtPointRecursive(mel_info, lambda, k, el_info);
715
716
    for (int i = 0; i <= dim; i++)
      bary[i] = final_lambda[i];   
717
  
Thomas Witkowski's avatar
Thomas Witkowski committed
718
    delete mel_info;
719

720
    return inside;
721
722
723
  }

  bool Mesh::findElementAtPoint(const WorldVector<double>&  xy,
724
725
				Element **elp, 
				DimVec<double>& bary,
726
				const MacroElement *start_mel,
727
728
				const WorldVector<double> *xy0,
				double *sp)
729
  {
730
731
    ElInfo *el_info = createNewElInfo();
    int val = findElInfoAtPoint(xy, el_info, bary, start_mel, xy0, sp);
732
733
734

    *elp = el_info->getElement();

Thomas Witkowski's avatar
Thomas Witkowski committed
735
    delete el_info;
736

737
    return val;
738
739
  }

740
  bool Mesh::findElementAtPointRecursive(ElInfo *el_info,
741
					 const DimVec<double>& lambda,
742
					 int outside,
743
744
					 ElInfo* final_el_info)
  {
745
    FUNCNAME("Mesh::findElementAtPointRecursive()");
Thomas Witkowski's avatar
Thomas Witkowski committed
746

747
748
    Element *el = el_info->getElement();
    DimVec<double> c_lambda(dim, NO_INIT);
749
750
    int inside;
    int ichild, c_outside;
751
752
753
754

    if (el->isLeaf()) {
      *final_el_info = *el_info;
      if (outside < 0) {
Thomas Witkowski's avatar
Thomas Witkowski committed
755
	for (int i = 0; i <= dim; i++)
756
	  final_lambda[i] = lambda[i];
Thomas Witkowski's avatar
Thomas Witkowski committed
757
	
758
	return true;
759
760
761
762
      }  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
763
764
765
766
	  for (int i = 0; i <= dim; i++) 
	    final_lambda[i] = s * c_lambda[i] + (1.0 - s) * lambda[i];
	  
	  if (g_sp)
767
	    *(g_sp) = s;
Thomas Witkowski's avatar
Thomas Witkowski committed
768
	  
769
	  if (dim == 3) 
770
	    MSG("Outside finest level on el %d: s = %.3e\n", el->getIndex(), s);
771

772
	  return false;  /* ??? */
773
	} else {
774
	  return false;
775
	}
776
      }
777
778
    }

779
    ElInfo *c_el_info = createNewElInfo();
780

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

803
    if (dim == 2) {
804
805
806
807
      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
808
	  TEST_EXIT(outside == 0)("outside curved boundary child 0\n");	  
809
810
811
812
813
814
815
816
817
	} 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
818
	  TEST_EXIT(outside == 0)("outside curved boundary child 1\n");	  
819
820
821
822
823
824
825
826
	} else {
	  c_lambda[0] = lambda[1] - lambda[0];
	  c_lambda[1] = lambda[2];
	  c_lambda[2] = 2.0 * lambda[0];
	}
      }
    } /* DIM == 2 */

827
    if (dim == 3) {
828
829
830
831
832
833
834
835
836
837
838
839
840
      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
841
	  int c_outside2 = c_el_info2->worldToCoord(*(g_xy), &c_lambda2);
842
843

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

	  if ((c_outside2 < 0) || (c_lambda2[c_outside2] > c_lambda[c_outside])) {
850
851
	    for (int i = 0; i <= dim; i++) 
	      c_lambda[i] = c_lambda2[i];	    
852
853
854
855
	    c_outside = c_outside2;
	    *c_el_info = *c_el_info2;
	    ichild = 1 - ichild;
	  }
Thomas Witkowski's avatar
Thomas Witkowski committed
856
	  delete c_el_info2;
857
858
859
860
861
862
	}
	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];
863
864
865
866
867
868
	  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]];
869
870
871
872
	  c_lambda[3] = 2.0 * lambda[1];
	} else {
	  c_el_info->fillElInfo(1, el_info);
	  c_lambda[0] = lambda[1] - lambda[0];
873
874
875
876
877
878
	  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]];
879
880
881
882
883
	  c_lambda[3] = 2.0 * lambda[0];
	}
      }
    }  /* DIM == 3 */

Thomas Witkowski's avatar
Thomas Witkowski committed
884
    inside = findElementAtPointRecursive(c_el_info, c_lambda, outside, final_el_info);
Thomas Witkowski's avatar
Thomas Witkowski committed
885
    delete c_el_info;
886

887
    return inside; 
888
889
  }

890
891

  bool Mesh::getDofIndexCoords(DegreeOfFreedom dof, 
892
893
894
895
896
897
			       const FiniteElemSpace* feSpace,
			       WorldVector<double>& coords)
  {
    DimVec<double>* baryCoords;
    bool found = false;
    TraverseStack stack;
898
    std::vector<DegreeOfFreedom> dofVec(feSpace->getBasisFcts()->getNumber());
Thomas Witkowski's avatar
Thomas Witkowski committed
899

900
    ElInfo *elInfo = stack.traverseFirst(this, -1, 
Thomas Witkowski's avatar
Thomas Witkowski committed
901
					 Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS);    
902
    while (elInfo) {
903
904
905
      feSpace->getBasisFcts()->getLocalIndices(elInfo->getElement(),
					       feSpace->getAdmin(),
					       dofVec);
Thomas Witkowski's avatar
Thomas Witkowski committed
906
      for (int i = 0; i < feSpace->getBasisFcts()->getNumber(); i++) {
907
	if (dofVec[i] == dof) {
908
909
910
	  baryCoords = feSpace->getBasisFcts()->getCoords(i);
	  elInfo->coordToWorld(*baryCoords, coords);
	  found = true;
Thomas Witkowski's avatar
Thomas Witkowski committed
911
	  break;	  
912
	}
Thomas Witkowski's avatar
Thomas Witkowski committed
913
914
      }
      
Thomas Witkowski's avatar