Mesh.cc 39.6 KB
Newer Older
Thomas Witkowski's avatar
Thomas Witkowski committed
1
2
3
4
5
6
#include <algorithm>
#include <set>
#include <map>

#include "time.h"

7
8
9
10
#include "AdaptStationary.h"
#include "AdaptInstationary.h"
#include "FiniteElemSpace.h"
#include "ElementData.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
11
#include "ElementDofIterator.h"
12
13
#include "MacroElement.h"
#include "MacroReader.h"
14
#include "MacroInfo.h"
15
16
17
18
19
20
21
22
23
24
25
#include "Mesh.h"
#include "Traverse.h"
#include "Parameters.h"
#include "FixVec.h"
#include "DOFVector.h"
#include "CoarseningManager.h"
#include "DOFIterator.h"
#include "VertexVector.h"
#include "MacroWriter.h"
#include "PeriodicMap.h"
#include "Projection.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
26
#include "ElInfoStack.h"
27
#include "Serializer.h"
28
#include "Lagrange.h"
29
30
31
32
33
34
35
36
37

namespace AMDiS {

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

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

38
39
40
41
42
43
44
45
46
47
48
49
50
51
  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);
52
53
54
55
56
57
58
59
60
61
62

  //**************************************************************************
  //  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;
63
64
  const Flag Mesh::CALL_MG_LEVEL           = 0X4000L;  
  const Flag Mesh::CALL_REVERSE_MODE       = 0X8000L;
65
66


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

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


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

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

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

    elementPrototype->setIndex(-1);

118
119
    elementIndex = 0;
  }
120

121

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

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

137

138
  Mesh& Mesh::operator=(const Mesh& m)
139
  {
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
    FUNCNAME("Mesh::operator=()");

    if (this == &m)
      return *this;

    TEST_EXIT(dim == m.dim)("operator= works only on meshes with equal dim!\n");

    name = m.name;
    nVertices = m.nVertices;
    nEdges = m.nEdges;
    nLeaves = m.nLeaves;
    nElements = m.nElements;
    nFaces = m.nFaces;
    maxEdgeNeigh = m.maxEdgeNeigh;
    diam = m.diam;
    parametric = NULL;

    preserveCoarseDOFs = m.preserveCoarseDOFs;
    nDOFEl = m.nDOFEl;
    nDOF = m.nDOF;
    nNodeEl = m.nNodeEl;
    node = m.node;
    newDOF = m.newDOF;
    elementIndex = m.elementIndex;
    initialized = m.initialized;
    final_lambda = m.final_lambda;
    
    /* ====================== Create new DOFAdmins ================== */
    admin.resize(m.admin.size());
    for (int i = 0; i < static_cast<int>(admin.size()); i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
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
    // === Determine to all dofs the macro elements where the dof is part of. ===
274

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

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

      elInfo = stack.traverseNext(elInfo);
    }		   

Thomas Witkowski's avatar
Thomas Witkowski committed
293

294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
    // === 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
313

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

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

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


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

Thomas Witkowski's avatar
Thomas Witkowski committed
340
    for (DofElMap::iterator dofsIt = dofsOwner.begin(); 
341
342
343
344
345
346
347
348
349
350
351
352
353
354
	 dofsIt != dofsOwner.end(); ++dofsIt) {
      
      bool deleteDof = true;

      for (std::set<MacroElement*>::iterator elIter = dofsIt->second.begin();
	   elIter != dofsIt->second.end(); ++elIter) {
	std::set<MacroElement*>::iterator mIt = macros.find(*elIter);
	if (mIt == macros.end()) {
	  deleteDof = false;
	  break;
	}
      }

      if (deleteDof)
355
	freeDof(const_cast<DegreeOfFreedom*>(dofsIt->first), 
Thomas Witkowski's avatar
Thomas Witkowski committed
356
		dofsPosIndex[dofsIt->first]);
357
358
    }

359

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

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

    std::set<const DegreeOfFreedom*> allVertices;

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

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

389

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

    localAdmin->setMesh(this);

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

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

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

406
407
408
409
410
    admin.push_back(localAdmin);

    nDOFEl = 0;

    localAdmin->setNumberOfPreDOFs(VERTEX,nDOF[VERTEX]);
411
    nDOF[VERTEX] += localAdmin->getNumberOfDOFs(VERTEX);
412
413
    nDOFEl += getGeo(VERTEX) * nDOF[VERTEX];

414
    if (dim > 1) {
415
      localAdmin->setNumberOfPreDOFs(EDGE,nDOF[EDGE]);
416
      nDOF[EDGE] += localAdmin->getNumberOfDOFs(EDGE);
417
418
419
420
421
422
423
      nDOFEl += getGeo(EDGE) * nDOF[EDGE];
    }

    localAdmin->setNumberOfPreDOFs(CENTER,nDOF[CENTER]);
    nDOF[CENTER]  += localAdmin->getNumberOfDOFs(CENTER);
    nDOFEl += nDOF[CENTER];

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
432
      if (nDOF[EDGE] > 0) 
	nNodeEl += getGeo(EDGE);
433
434
    }

435
    if (dim == 3) {
436
      localAdmin->setNumberOfPreDOFs(FACE,nDOF[FACE]);
437
438
439
440
441
      nDOF[FACE] += localAdmin->getNumberOfDOFs(FACE);
      nDOFEl += getGeo(FACE) * nDOF[FACE];
      node[FACE] = nNodeEl;
      if (nDOF[FACE] > 0) 
	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
461
462
      if (size < 1 || 
	  compressAdmin->getUsedDOFs() < 1 || 
	  compressAdmin->getHoleCount() < 1)    
463
	continue;
Thomas Witkowski's avatar
Thomas Witkowski committed
464

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

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

      elInfo = stack.traverseFirst(this, -1, fill_flag);
      while (elInfo) {
Thomas Witkowski's avatar
Thomas Witkowski committed
481
	elInfo->getElement()->newDOFFct2(compressAdmin);
482
483
484
	elInfo = stack.traverseNext(elInfo);
      }

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


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

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

497
498
    int ndof = getNumberOfDOFs(position);
    if (ndof <= 0) 
499
      return NULL;
500

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

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


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

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

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

531
    return ptrs;
532
533
  }

534

535
536
  void Mesh::freeDOFPtrs(DegreeOfFreedom **ptrs)
  {
537
    FUNCNAME("Mesh::freeDOFPtrs()");
538

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

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


Thomas Witkowski's avatar
Thomas Witkowski committed
548
  const DOFAdmin *Mesh::createDOFAdmin(std::string lname, DimVec<int> lnDOF)
549
  {
550
    FUNCNAME("Mesh::createDOFAdmin()");
551

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

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

    addDOFAdmin(localAdmin);

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


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

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

576
    return localAdmin;
577
578
  }

579

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

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

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

597
598
      return;
    }
599

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

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

615
    delete [] dof;
616
617
  }

618

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


  Element* Mesh::createNewElement(Element *parent)
  {
    FUNCNAME("Mesh::createNewElement()");
629
630

    TEST_EXIT_DBG(elementPrototype)("no element prototype\n");
631
632

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

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

639
640
641
    return el;
  }

642

643
644
  ElInfo* Mesh::createNewElInfo()
  {
645
646
    FUNCNAME("Mesh::createNewElInfo()");

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

663

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

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

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

    mel_info->fillMacroInfo(mel);

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

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

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

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

716
    return inside;
717
718
719
  }

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

    *elp = el_info->getElement();

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

733
    return val;
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->getBasisFcts()->getCoords(i);
	  elInfo->coordToWorld(*baryCoords, coords);
	  found = true;
Thomas Witkowski's avatar
Thomas Witkowski committed
907
	  break;	  
908
	}
Thomas Witkowski's avatar
Thomas Witkowski committed
909
910
      }
      
911
912
      if (found)
	break;
Thomas Witkowski's avatar
Thomas Witkowski committed
913
      
914
915
916
917
918
      elInfo = stack.traverseNext(elInfo);
    }

    return found;
  }
919

920
921
922

  void Mesh::getDofIndexCoords(const FiniteElemSpace* feSpace,
			       DOFVector<WorldVector<double> >& coords)
923
  {
924
925
926
    const BasisFunction* basFcts = feSpace->getBasisFcts();
    int nBasFcts = basFcts->getNumber();
    std::vector<DegreeOfFreedom> dofVec(nBasFcts);
927

928
929
930
    TraverseStack stack;
    ElInfo *elInfo = 
      stack.traverseFirst(this, -1, Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS);
931
    while (elInfo) {
932
933
934
935
      basFcts->getLocalIndices(elInfo->getElement(), feSpace->getAdmin(), dofVec);
      for (int i = 0; i < nBasFcts; i++) {
	DimVec<double> *baryCoords = basFcts->getCoords(i);
	elInfo->coordToWorld(*baryCoords, coords[dofVec[i]]);
936
      }
937
         
938
939
940
941
942
      elInfo = stack.traverseNext(elInfo);
    }

  }

943

Thomas Witkowski's avatar
Thomas Witkowski committed
944
945
946
947
  void Mesh::setDiameter(const WorldVector<double>& w) 
  { 
    diam = w; 
  }
948