Mesh.cc 39.8 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
    // === 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())
308
	newMacroElements.push_back(*elIter);     
309
    }
310

311
312
313
314
    // And replace the macro element list with the new one.
    macroElements.clear();
    macroElements = newMacroElements;

Thomas Witkowski's avatar
Thomas Witkowski committed
315

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

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

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

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

	(*macroIt)->getElement()->setElementData(elementDataPrototype->clone()); 
340
      }
341
342
343
344
345

      // We will delete at least some element DOFs in the next but will keep the
      // element object still in memory. Therefore the DOF pointer must be set to be
      // invalide.
      (*macroIt)->getElement()->setDofValid(false);
346
    }     
347
348


349
350
    // === Check now all the dofs, that have no owner anymore and therefore have ===
    // === to be removed.                                                        ===
351

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

371

372
    // === Update number of elements, vertices, etc. ===
373

374
375
376
377
378
379
380
381
    nLeaves = 0;
    nElements = 0;

    std::set<const DegreeOfFreedom*> allVertices;

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

383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
      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.                               ===
399
  }
400

401

402
403
404
405
406
407
  void Mesh::addDOFAdmin(DOFAdmin *localAdmin)
  {    
    FUNCNAME("Mesh::addDOFAdmin()");

    localAdmin->setMesh(this);

408
    std::vector<DOFAdmin*>::iterator dai = 
409
      std::find(admin.begin(), admin.end(), localAdmin);
410

Thomas Witkowski's avatar
Thomas Witkowski committed
411
412
413
    TEST_EXIT(dai == admin.end())
      ("admin %s is already associated to mesh %s\n",
       localAdmin->getName().c_str(), this->getName().c_str());
414

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

418
419
    admin.push_back(localAdmin);

420
    nDofEl = 0;
421

422
423
424
    localAdmin->setNumberOfPreDOFs(VERTEX,nDof[VERTEX]);
    nDof[VERTEX] += localAdmin->getNumberOfDofs(VERTEX);
    nDofEl += getGeo(VERTEX) * nDof[VERTEX];
425

426
    if (dim > 1) {
427
428
429
      localAdmin->setNumberOfPreDOFs(EDGE,nDof[EDGE]);
      nDof[EDGE] += localAdmin->getNumberOfDofs(EDGE);
      nDofEl += getGeo(EDGE) * nDof[EDGE];
430
431
    }

432
433
434
    localAdmin->setNumberOfPreDOFs(CENTER,nDof[CENTER]);
    nDof[CENTER]  += localAdmin->getNumberOfDofs(CENTER);
    nDofEl += nDof[CENTER];
435

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

438
439
    node[VERTEX] = 0;
    nNodeEl = getGeo(VERTEX);
440

441
442
    if (dim > 1) {
      node[EDGE] = nNodeEl;
443
      if (nDof[EDGE] > 0) 
444
	nNodeEl += getGeo(EDGE);
445
446
    }

447
    if (dim == 3) {
448
449
450
      localAdmin->setNumberOfPreDOFs(FACE,nDof[FACE]);
      nDof[FACE] += localAdmin->getNumberOfDofs(FACE);
      nDofEl += getGeo(FACE) * nDof[FACE];
451
      node[FACE] = nNodeEl;
452
      if (nDof[FACE] > 0) 
453
	nNodeEl += getGeo(FACE);
454
455
    }

456
    node[CENTER] = nNodeEl;
457
    if (nDof[CENTER] > 0)
458
      nNodeEl += 1;
459
460
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
461

462
463
  void Mesh::dofCompress()
  {
464
    FUNCNAME("Mesh::dofCompress()");
465

Thomas Witkowski's avatar
Thomas Witkowski committed
466
467
    for (unsigned int iadmin = 0; iadmin < admin.size(); iadmin++) {
      DOFAdmin* compressAdmin = admin[iadmin];
468
469

      TEST_EXIT_DBG(compressAdmin)("no admin[%d] in mesh\n", iadmin);
470
      
Thomas Witkowski's avatar
Thomas Witkowski committed
471
      int size = compressAdmin->getSize();
Thomas Witkowski's avatar
Thomas Witkowski committed
472
      if (size < 1 || 
473
	  compressAdmin->getUsedDofs() < 1 || 
Thomas Witkowski's avatar
Thomas Witkowski committed
474
	  compressAdmin->getHoleCount() < 1)    
475
	continue;
Thomas Witkowski's avatar
Thomas Witkowski committed
476

477
478
      std::vector<int> newDofIndex(size);     
      compressAdmin->compress(newDofIndex);
Thomas Witkowski's avatar
Thomas Witkowski committed
479
480
481
482

      Flag fill_flag = (preserveCoarseDOFs ?  
			Mesh::CALL_EVERY_EL_PREORDER | Mesh::FILL_NOTHING :
			Mesh::CALL_LEAF_EL | Mesh::FILL_NOTHING);          
483
484
485
      TraverseStack stack;
      ElInfo *elInfo = stack.traverseFirst(this, -1, fill_flag);
      while (elInfo) {
486
	elInfo->getElement()->newDofFct1(compressAdmin, newDofIndex);
487
488
489
490
491
	elInfo = stack.traverseNext(elInfo);
      }

      elInfo = stack.traverseFirst(this, -1, fill_flag);
      while (elInfo) {
492
	elInfo->getElement()->newDofFct2(compressAdmin);
493
494
	elInfo = stack.traverseNext(elInfo);
      }
495
    }       
496
497
498
  }


499
  DegreeOfFreedom *Mesh::getDof(GeoIndex position)
500
  {
501
    FUNCNAME("Mesh::getDof()");
502

503
    TEST_EXIT_DBG(position >= CENTER && position <= FACE)
504
      ("unknown position %d\n", position);
505

506
    int ndof = getNumberOfDofs(position);
507
    if (ndof <= 0) 
508
      return NULL;
509

510
    DegreeOfFreedom *dof = new DegreeOfFreedom[ndof];
511

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


Thomas Witkowski's avatar
Thomas Witkowski committed
529
  DegreeOfFreedom **Mesh::createDofPtrs()
530
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
531
    FUNCNAME("Mesh::createDofPtrs()");
532
533

    if (nNodeEl <= 0)
534
      return NULL;
535

536
    DegreeOfFreedom **ptrs = new DegreeOfFreedom*[nNodeEl];
537
    for (int i = 0; i < nNodeEl; i++)
538
539
      ptrs[i] = NULL;

540
    return ptrs;
541
542
  }

543

544
  void Mesh::freeDofPtrs(DegreeOfFreedom **ptrs)
545
  {
546
    FUNCNAME("Mesh::freeDofPtrs()");
547

548
    TEST_EXIT_DBG(ptrs)("ptrs is NULL!\n");
549
550
551
552

    if (nNodeEl <= 0)
      return;
  
553
    delete [] ptrs;
554
555
556
  }


557
  const DOFAdmin *Mesh::createDOFAdmin(std::string lname, DimVec<int> lnDof)
558
  {
559
    FUNCNAME("Mesh::createDOFAdmin()");
560

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

Thomas Witkowski's avatar
Thomas Witkowski committed
563
    for (int i = 0; i < dim + 1; i++)
564
      localAdmin->setNumberOfDOFs(i, lnDof[i]);
565
566
567

    addDOFAdmin(localAdmin);

568
    return localAdmin;
569
570
571
572
573
574
575
  }


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

576
    for (int i = 0; i < static_cast<int>(admin.size()); i++) {
577
      if (admin[i]->getNumberOfDofs(VERTEX)) {
578
579
580
581
	if (!localAdmin)  
	  localAdmin = admin[i];
	else if (admin[i]->getSize() < localAdmin->getSize())
	  localAdmin = admin[i];
582
      }
583
584
    }

585
    return localAdmin;
586
587
  }

588

589
  void Mesh::freeDof(DegreeOfFreedom* dof, GeoIndex position)
590
  {
591
    FUNCNAME("Mesh::freeDof()");
592

593
    TEST_EXIT_DBG(position >= CENTER && position <= FACE)
594
      ("unknown position %d\n", position);
595

596
    int ndof = nDof[position];
597
598
    if (ndof) {
      if (!dof) {
599
	MSG("dof = NULL, but ndof = %d\n", ndof);
600
601
	return;
      }
602
    } else  {
603
604
605
      if (dof)
	MSG("dof != NULL, but ndof = 0\n");

606
607
      return;
    }
608

609
    TEST_EXIT_DBG(ndof <= MAX_DOF)
610
      ("ndof too big: ndof = %d, MAX_DOF = %d\n", ndof, MAX_DOF);
611

612
    for (unsigned int i = 0; i < admin.size(); i++) {
613
      DOFAdmin *localAdmin = admin[i];
614
      int n = localAdmin->getNumberOfDofs(position);
615
616
      int n0 = localAdmin->getNumberOfPreDOFs(position);
      
617
618
      TEST_EXIT_DBG(n + n0 <= ndof)
	("n = %d, n0 = %d too large: ndof = %d\n", n, n0, ndof);
619
620
      
      for (int j = 0; j < n; j++)
621
	localAdmin->freeDofIndex(dof[n0 + j]);      
622
    }
623

624
    delete [] dof;
625
    dof = NULL;
626
627
  }

628

629
630
  void Mesh::freeElement(Element* el)
  {
631
    freeDofPtrs(const_cast<DegreeOfFreedom**>(el->getDof()));
Thomas Witkowski's avatar
Thomas Witkowski committed
632
    delete el;
633
634
635
636
637
638
  }


  Element* Mesh::createNewElement(Element *parent)
  {
    FUNCNAME("Mesh::createNewElement()");
639
640

    TEST_EXIT_DBG(elementPrototype)("no element prototype\n");
641
642

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

Thomas Witkowski's avatar
Thomas Witkowski committed
644
    if (!parent && elementDataPrototype)
645
      el->setElementData(elementDataPrototype->clone()); 
Thomas Witkowski's avatar
Thomas Witkowski committed
646
    else
647
      el->setElementData(NULL); // must be done in ElementData::refineElementData()
648

649
650
651
    return el;
  }

652

653
654
  ElInfo* Mesh::createNewElInfo()
  {
655
656
    FUNCNAME("Mesh::createNewElInfo()");

657
658
    switch(dim) {
    case 1:
Thomas Witkowski's avatar
Thomas Witkowski committed
659
      return new ElInfo1d(this);
660
661
      break;
    case 2:
Thomas Witkowski's avatar
Thomas Witkowski committed
662
      return new ElInfo2d(this);
663
664
      break;
    case 3:
Thomas Witkowski's avatar
Thomas Witkowski committed
665
      return new ElInfo3d(this);
666
667
668
669
      break;
    default:
      ERROR_EXIT("invalid dim\n");
      return NULL;
670
    }
671
672
  }

673

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

    if (start_mel != NULL)
      mel = start_mel;
    else
Thomas Witkowski's avatar
Thomas Witkowski committed
688
      if (mel == NULL || mel->getElement()->getMesh() != this)
689
690
691
	mel = *(macroElements.begin());

    mel_info->setFillFlag(Mesh::FILL_COORDS);
692
    g_xy = &xy;
693
    g_xy0 = xy0;
694
    g_sp = sp;
695
696
697

    mel_info->fillMacroInfo(mel);

698
699
700
701
702
703
704
    // 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());

705
    int k;
706
707
    while ((k = mel_info->worldToCoord(xy, &lambda)) >= 0) {
      if (mel->getNeighbour(k)) {
708
709
710
711
712
	mel = mel->getNeighbour(k);
	if (macrosVisited.count(mel->getIndex()))
	  return false;

	macrosVisited.insert(mel->getIndex());
713
714
715
716
717
718
719
	mel_info->fillMacroInfo(mel);
	continue;
      }
      break;
    }

    /* now, descend in tree to find leaf element at point */
720
    bool inside = findElementAtPointRecursive(mel_info, lambda, k, el_info);
721
722
    for (int i = 0; i <= dim; i++)
      bary[i] = final_lambda[i];   
723
  
Thomas Witkowski's avatar
Thomas Witkowski committed
724
    delete mel_info;
725

726
    return inside;
727
728
  }

729

730
  bool Mesh::findElementAtPoint(const WorldVector<double>&  xy,
731
732
				Element **elp, 
				DimVec<double>& bary,
733
				const MacroElement *start_mel,
734
735
				const WorldVector<double> *xy0,
				double *sp)
736
  {
737
738
    ElInfo *el_info = createNewElInfo();
    int val = findElInfoAtPoint(xy, el_info, bary, start_mel, xy0, sp);
739
740
741

    *elp = el_info->getElement();

Thomas Witkowski's avatar
Thomas Witkowski committed
742
    delete el_info;
743

744
    return val;
745
746
  }

747

748
  bool Mesh::findElementAtPointRecursive(ElInfo *el_info,
749
					 const DimVec<double>& lambda,
750
					 int outside,
751
752
					 ElInfo* final_el_info)
  {
753
    FUNCNAME("Mesh::findElementAtPointRecursive()");
Thomas Witkowski's avatar
Thomas Witkowski committed
754

755
756
    Element *el = el_info->getElement();
    DimVec<double> c_lambda(dim, NO_INIT);
757
758
    int inside;
    int ichild, c_outside;
759
760
761
762

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

780
	  return false;  /* ??? */
781
	} else {
782
	  return false;
783
	}
784
      }
785
786
    }

787
    ElInfo *c_el_info = createNewElInfo();
788

789
    if (dim == 1) {
790
791
792
793
      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
794
	  TEST_EXIT(outside == 0)("point outside domain\n");
795
796
797
798
799
800
801
802
	} 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
803
	  TEST_EXIT(outside == 0)("point outside domain\n");
804
805
806
807
808
809
810
	} else {
	  c_lambda[1] = lambda[1] - lambda[0];
	  c_lambda[0] = 2.0 * lambda[0];
	}
      }
    } /* DIM == 1 */

811
    if (dim == 2) {
812
813
814
815
      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
816
	  TEST_EXIT(outside == 0)("outside curved boundary child 0\n");	  
817
818
819
820
821
822
823
824
825
	} 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
826
	  TEST_EXIT(outside == 0)("outside curved boundary child 1\n");	  
827
828
829
830
831
832
833
834
	} else {
	  c_lambda[0] = lambda[1] - lambda[0];
	  c_lambda[1] = lambda[2];
	  c_lambda[2] = 2.0 * lambda[0];
	}
      }
    } /* DIM == 2 */

835
    if (dim == 3) {
836
837
838
839
840
841
842
843
844
845
846
847
848
      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
849
	  int c_outside2 = c_el_info2->worldToCoord(*(g_xy), &c_lambda2);
850
851

	  MSG("new_coord CHILD %d: outside=%d, lambda=(%.2f %.2f %.2f %.2f)\n",
852
	      ichild, c_outside, c_lambda[0], c_lambda[1], c_lambda[2], c_lambda[3]);
853
	  MSG("new_coord CHILD %d: outside=%d, lambda=(%.2f %.2f %.2f %.2f)\n",
854
	      1 - ichild, c_outside2, c_lambda2[0], c_lambda2[1], c_lambda2[2],
855
856
857
	      c_lambda2[3]);

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

Thomas Witkowski's avatar
Thomas Witkowski committed
892
    inside = findElementAtPointRecursive(c_el_info, c_lambda, outside, final_el_info);
Thomas Witkowski's avatar
Thomas Witkowski committed
893
    delete c_el_info;
894

895
    return inside; 
896
897
  }

898
899

  bool Mesh::getDofIndexCoords(DegreeOfFreedom dof, 
900
901
902
903
904
905
			       const FiniteElemSpace* feSpace,
			       WorldVector<double>& coords)
  {
    DimVec<double>* baryCoords;
    bool found = false;
    TraverseStack stack;
906
    std::vector<DegreeOfFreedom> dofVec(feSpace->getBasisFcts()->getNumber());
Thomas Witkowski's avatar
Thomas Witkowski committed
907

908
    ElInfo *elInfo = stack.traverseFirst(this, -1, 
Thomas Witkowski's avatar
Thomas Witkowski committed
909
					 Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS);    
910
    while (elInfo) {
911
912
913
      feSpace->getBasisFcts()->getLocalIndices(elInfo->getElement(),
					       feSpace->getAdmin(),
					       dofVec);
Thomas Witkowski's avatar
Thomas Witkowski committed