Mesh.cc 39.7 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
170
    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
171
      admin[i] = new DOFAdmin(this);
172
      *(admin[i]) = *(m.admin[i]);
173
174
      admin[i]->setMesh(this);
    }
175

176

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

189
190
    macroElements.clear();

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

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

242

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

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

255

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

262

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

Thomas Witkowski's avatar
Thomas Witkowski committed
268
269
270
271
272
273
    typedef std::map<const DegreeOfFreedom*, std::set<MacroElement*> > DofElMap;
    typedef std::map<const DegreeOfFreedom*, GeoIndex> DofPosMap;

    TEST_EXIT(admin.size() == 1)("Not yet implemented for multiple admins!\n");
    TEST_EXIT(admin[0])("There is something wrong!\n");

274
    // === Determine to all dofs 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
334
      // Delete element hierarchie
      if (!(*macroIt)->getElement()->isLeaf()) {
	delete (*macroIt)->getElement()->getChild(0);
	delete (*macroIt)->getElement()->getChild(1);
      }
335
    }     
336
337


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

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

360

361
    // === Update number of elements, vertices, etc. ===
362

363
364
365
366
367
368
369
370
    nLeaves = 0;
    nElements = 0;

    std::set<const DegreeOfFreedom*> allVertices;

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

372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
      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.                               ===
388
  }
389

390

391
392
393
394
395
396
  void Mesh::addDOFAdmin(DOFAdmin *localAdmin)
  {    
    FUNCNAME("Mesh::addDOFAdmin()");

    localAdmin->setMesh(this);

397
    std::vector<DOFAdmin*>::iterator dai = 
398
      std::find(admin.begin(), admin.end(), localAdmin);
399

Thomas Witkowski's avatar
Thomas Witkowski committed
400
401
402
    TEST_EXIT(dai == admin.end())
      ("admin %s is already associated to mesh %s\n",
       localAdmin->getName().c_str(), this->getName().c_str());
403

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

407
408
    admin.push_back(localAdmin);

409
    nDofEl = 0;
410

411
412
413
    localAdmin->setNumberOfPreDOFs(VERTEX,nDof[VERTEX]);
    nDof[VERTEX] += localAdmin->getNumberOfDofs(VERTEX);
    nDofEl += getGeo(VERTEX) * nDof[VERTEX];
414

415
    if (dim > 1) {
416
417
418
      localAdmin->setNumberOfPreDOFs(EDGE,nDof[EDGE]);
      nDof[EDGE] += localAdmin->getNumberOfDofs(EDGE);
      nDofEl += getGeo(EDGE) * nDof[EDGE];
419
420
    }

421
422
423
    localAdmin->setNumberOfPreDOFs(CENTER,nDof[CENTER]);
    nDof[CENTER]  += localAdmin->getNumberOfDofs(CENTER);
    nDofEl += nDof[CENTER];
424

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

427
428
    node[VERTEX] = 0;
    nNodeEl = getGeo(VERTEX);
429

430
431
    if (dim > 1) {
      node[EDGE] = nNodeEl;
432
      if (nDof[EDGE] > 0) 
433
	nNodeEl += getGeo(EDGE);
434
435
    }

436
    if (dim == 3) {
437
438
439
      localAdmin->setNumberOfPreDOFs(FACE,nDof[FACE]);
      nDof[FACE] += localAdmin->getNumberOfDofs(FACE);
      nDofEl += getGeo(FACE) * nDof[FACE];
440
      node[FACE] = nNodeEl;
441
      if (nDof[FACE] > 0) 
442
	nNodeEl += getGeo(FACE);
443
444
    }

445
    node[CENTER] = nNodeEl;
446
    if (nDof[CENTER] > 0)
447
      nNodeEl += 1;
448
449
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
450

451
452
  void Mesh::dofCompress()
  {
453
    FUNCNAME("Mesh::dofCompress()");
454

Thomas Witkowski's avatar
Thomas Witkowski committed
455
456
    for (unsigned int iadmin = 0; iadmin < admin.size(); iadmin++) {
      DOFAdmin* compressAdmin = admin[iadmin];
457
458

      TEST_EXIT_DBG(compressAdmin)("no admin[%d] in mesh\n", iadmin);
459
      
Thomas Witkowski's avatar
Thomas Witkowski committed
460
      int size = compressAdmin->getSize();
Thomas Witkowski's avatar
Thomas Witkowski committed
461
      if (size < 1 || 
462
	  compressAdmin->getUsedDofs() < 1 || 
Thomas Witkowski's avatar
Thomas Witkowski committed
463
	  compressAdmin->getHoleCount() < 1)    
464
	continue;
Thomas Witkowski's avatar
Thomas Witkowski committed
465

466
467
468
      newDOF.resize(size);
      
      compressAdmin->compress(newDOF);
Thomas Witkowski's avatar
Thomas Witkowski committed
469
470
471
472

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

      elInfo = stack.traverseFirst(this, -1, fill_flag);
      while (elInfo) {
482
	elInfo->getElement()->newDofFct2(compressAdmin);
483
484
485
	elInfo = stack.traverseNext(elInfo);
      }

486
      newDOF.resize(0);
487
    }       
488
489
490
  }


491
  DegreeOfFreedom *Mesh::getDof(GeoIndex position)
492
  {
493
    FUNCNAME("Mesh::getDof()");
494

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

498
    int ndof = getNumberOfDofs(position);
499
    if (ndof <= 0) 
500
      return NULL;
501

502
    DegreeOfFreedom *dof = new DegreeOfFreedom[ndof];
503

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


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

    if (nNodeEl <= 0)
526
      return NULL;
527

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

532
    return ptrs;
533
534
  }

535

536
  void Mesh::freeDofPtrs(DegreeOfFreedom **ptrs)
537
  {
538
    FUNCNAME("Mesh::freeDofPtrs()");
539

540
    TEST_EXIT_DBG(ptrs)("ptrs is NULL!\n");
541
542
543
544

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


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

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

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

    addDOFAdmin(localAdmin);

560
    return localAdmin;
561
562
563
564
565
566
567
  }


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

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

577
    return localAdmin;
578
579
  }

580

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

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

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

598
599
      return;
    }
600

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

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

616
    delete [] dof;
617
    dof = NULL;
618
619
  }

620

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


  Element* Mesh::createNewElement(Element *parent)
  {
    FUNCNAME("Mesh::createNewElement()");
631
632

    TEST_EXIT_DBG(elementPrototype)("no element prototype\n");
633
634

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

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

641
642
643
    return el;
  }

644

645
646
  ElInfo* Mesh::createNewElInfo()
  {
647
648
    FUNCNAME("Mesh::createNewElInfo()");

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

665

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

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

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

    mel_info->fillMacroInfo(mel);

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

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

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

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

718
    return inside;
719
720
721
  }

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

    *elp = el_info->getElement();

Thomas Witkowski's avatar
Thomas Witkowski committed
733
    delete el_info;
734

735
    return val;
736
737
  }

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

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

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

770
	  return false;  /* ??? */
771
	} else {
772
	  return false;
773
	}
774
      }
775
776
    }

777
    ElInfo *c_el_info = createNewElInfo();
778

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

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

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

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

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

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

885
    return inside; 
886
887
  }

888
889

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

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

    return found;
  }