Mesh.cc 39.5 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

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

94
    FUNCNAME("Mesh::Mesh()");
95
96
97
98

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

    elementPrototype->setIndex(-1);

113
114
    elementIndex = 0;
  }
115

116

117
  Mesh::~Mesh()
118
  {
119
    deleteMeshStructure();
120

121
122
123
124
125
126
    if (macroFileInfo != NULL)
      delete macroFileInfo;    
    if (elementPrototype)
      delete elementPrototype;    
    if (elementDataPrototype)
      delete elementDataPrototype;    
127
    
128
    for (unsigned int i = 0; i < admin.size(); i++)
129
      delete admin[i];    
130
  }
131

132

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

169

170
    /* ====================== Copy macro elements =================== */
171
  
172
173
174
175
176
177
178
179
180
    // 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;
181

182
183
    macroElements.clear();

184
185
186
    // 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();
187
	 it != m.macroElements.end(); ++it, insertCounter++) {
188

189
      // Create new MacroElement.
Thomas Witkowski's avatar
Thomas Witkowski committed
190
      MacroElement *el = new MacroElement(dim);
191

192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
      // 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;
  }

235

236
237
238
239
240
241
242
243
244
245
246
247
  void Mesh::updateNumberOfLeaves()
  {
    nLeaves = 0;

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

248

249
250
251
252
253
254
  void Mesh::addMacroElement(MacroElement* me) 
  {
    macroElements.push_back(me); 
    me->setIndex(macroElements.size());
  }

255

256
  void Mesh::removeMacroElements(std::set<MacroElement*>& macros,
Thomas Witkowski's avatar
Thomas Witkowski committed
257
				 const FiniteElemSpace *feSpace) 
258
259
260
  {
    FUNCNAME("Mesh::removeMacroElement()");

Thomas Witkowski's avatar
Thomas Witkowski committed
261
262
263
264
265
266
    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");

267
268
    // === Determine to all DOFs in mesh the macro elements where the DOF  ===
    // === is part of.                                                     ===
269

270
    // Map that stores for each dof pointer (which may have a list of dofs)
271
    // all macro element indices that own this dof.
Thomas Witkowski's avatar
Thomas Witkowski committed
272
273
    DofElMap dofsOwner;
    DofPosMap dofsPosIndex;
274

275
276
277
278
279
    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
280
      do {
281
	dofsOwner[elDofIter.getDofPtr()].insert(elInfo->getMacroElement());
Thomas Witkowski's avatar
Thomas Witkowski committed
282
283
	dofsPosIndex[elDofIter.getDofPtr()] = elDofIter.getPosIndex();
      } while (elDofIter.nextStrict());
284
285
286
287

      elInfo = stack.traverseNext(elInfo);
    }		   

Thomas Witkowski's avatar
Thomas Witkowski committed
288

289
290
291
292
293
294
295
296
297
298
299
300
301
    // === 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())
302
	newMacroElements.push_back(*elIter);     
303
    }
304

305
306
307
308
    // And replace the macro element list with the new one.
    macroElements.clear();
    macroElements = newMacroElements;

Thomas Witkowski's avatar
Thomas Witkowski committed
309

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

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

325
326
327
328
      // Delete element hierarchie
      if (!(*macroIt)->getElement()->isLeaf()) {
	delete (*macroIt)->getElement()->getChild(0);
	delete (*macroIt)->getElement()->getChild(1);
329
330
331
332
333

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

	(*macroIt)->getElement()->setElementData(elementDataPrototype->clone()); 
334
      }
335
336
337
338
339

      // 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);
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), 
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
    TEST_EXIT(std::find(admin.begin(), admin.end(), localAdmin) == admin.end())
Thomas Witkowski's avatar
Thomas Witkowski committed
403
404
      ("admin %s is already associated to mesh %s\n",
       localAdmin->getName().c_str(), this->getName().c_str());
405
406
407

    admin.push_back(localAdmin);

408
    nDofEl = 0;
409

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

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

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

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

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

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
449

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

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

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

465
466
      std::vector<int> newDofIndex(size);     
      compressAdmin->compress(newDofIndex);
Thomas Witkowski's avatar
Thomas Witkowski committed
467
468
469
470

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

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


487
  DegreeOfFreedom *Mesh::getDof(GeoIndex position)
488
  {
489
    FUNCNAME("Mesh::getDof()");
490

491
    TEST_EXIT_DBG(position >= CENTER && position <= FACE)
492
      ("unknown position %d\n", position);
493

494
    int ndof = getNumberOfDofs(position);
495
    if (ndof <= 0) 
496
      return NULL;
497

498
    DegreeOfFreedom *dof = new DegreeOfFreedom[ndof];
499

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


Thomas Witkowski's avatar
Thomas Witkowski committed
517
  DegreeOfFreedom **Mesh::createDofPtrs()
518
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
519
    FUNCNAME("Mesh::createDofPtrs()");
520
521

    if (nNodeEl <= 0)
522
      return NULL;
523

524
    DegreeOfFreedom **ptrs = new DegreeOfFreedom*[nNodeEl];
525
    for (int i = 0; i < nNodeEl; i++)
526
527
      ptrs[i] = NULL;

528
    return ptrs;
529
530
  }

531

532
  void Mesh::freeDofPtrs(DegreeOfFreedom **ptrs)
533
  {
534
    FUNCNAME("Mesh::freeDofPtrs()");
535

536
    TEST_EXIT_DBG(ptrs)("ptrs is NULL!\n");
537
538
539
540

    if (nNodeEl <= 0)
      return;
  
541
    delete [] ptrs;
542
543
544
  }


545
  const DOFAdmin *Mesh::createDOFAdmin(std::string lname, DimVec<int> lnDof)
546
  {
547
    FUNCNAME("Mesh::createDOFAdmin()");
548
    
Thomas Witkowski's avatar
Thomas Witkowski committed
549
    DOFAdmin *localAdmin = new DOFAdmin(this, lname);
550

Thomas Witkowski's avatar
Thomas Witkowski committed
551
    for (int i = 0; i < dim + 1; i++)
552
      localAdmin->setNumberOfDofs(i, lnDof[i]);
553
554
555

    addDOFAdmin(localAdmin);

556
    return localAdmin;
557
558
559
560
561
562
563
  }


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

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

573
    return localAdmin;
574
575
  }

576

577
  void Mesh::freeDof(DegreeOfFreedom* dof, GeoIndex position)
578
  {
579
    FUNCNAME("Mesh::freeDof()");
580

581
    TEST_EXIT_DBG(position >= CENTER && position <= FACE)
582
      ("unknown position %d\n", position);
583

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

594
595
      return;
    }
596

597
    TEST_EXIT_DBG(ndof <= MAX_DOF)
598
      ("ndof too big: ndof = %d, MAX_DOF = %d\n", ndof, MAX_DOF);
599

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

612
    delete [] dof;
613
    dof = NULL;
614
615
  }

616

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


  Element* Mesh::createNewElement(Element *parent)
  {
    FUNCNAME("Mesh::createNewElement()");
627
628

    TEST_EXIT_DBG(elementPrototype)("no element prototype\n");
629
630

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

Thomas Witkowski's avatar
Thomas Witkowski committed
632
    if (!parent && elementDataPrototype)
633
      el->setElementData(elementDataPrototype->clone()); 
Thomas Witkowski's avatar
Thomas Witkowski committed
634
    else
635
      el->setElementData(NULL); // must be done in ElementData::refineElementData()
636

637
638
639
    return el;
  }

640

641
642
  ElInfo* Mesh::createNewElInfo()
  {
643
644
    FUNCNAME("Mesh::createNewElInfo()");

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

661

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

    if (start_mel != NULL)
      mel = start_mel;
    else
Thomas Witkowski's avatar
Thomas Witkowski committed
676
      if (mel == NULL || mel->getElement()->getMesh() != this)
677
678
679
	mel = *(macroElements.begin());

    mel_info->setFillFlag(Mesh::FILL_COORDS);
680
    g_xy = &xy;
681
    g_xy0 = xy0;
682
    g_sp = sp;
683
684
685

    mel_info->fillMacroInfo(mel);

686
687
688
689
690
691
692
    // 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());

693
    int k;
694
695
    while ((k = mel_info->worldToCoord(xy, &lambda)) >= 0) {
      if (mel->getNeighbour(k)) {
696
697
698
699
700
	mel = mel->getNeighbour(k);
	if (macrosVisited.count(mel->getIndex()))
	  return false;

	macrosVisited.insert(mel->getIndex());
701
702
703
704
705
706
707
	mel_info->fillMacroInfo(mel);
	continue;
      }
      break;
    }

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

714
    return inside;
715
716
  }

717

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

    *elp = el_info->getElement();

Thomas Witkowski's avatar
Thomas Witkowski committed
730
    delete el_info;
731

732
    return val;
733
734
  }

735

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

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

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

768
	  return false;  /* ??? */
769
	} else {
770
	  return false;
771
	}
772
      }
773
774
    }

775
    ElInfo *c_el_info = createNewElInfo();
776

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

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

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

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

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

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

883
    return inside; 
884
885
  }

886
887

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

896
    ElInfo *elInfo = stack.traverseFirst(this, -1, 
Thomas Witkowski's avatar
Thomas Witkowski committed
897
					 Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS);    
898
    while (elInfo) {
899
900
901
      feSpace->getBasisFcts()->getLocalIndices(elInfo->getElement(),
					       feSpace->getAdmin(),
					       dofVec);
Thomas Witkowski's avatar
Thomas Witkowski committed
902
      for (int i = 0; i < feSpace->getBasisFcts()->getNumber(); i++) {
903
	if (dofVec[i] == dof) {
904
905
906
	  baryCoords = feSpace->getBasisFct