Mesh.cc 40.3 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
//
// Software License for AMDiS
//
// Copyright (c) 2010 Dresden University of Technology 
// All rights reserved.
// Authors: Simon Vey, Thomas Witkowski et al.
//
// This file is part of AMDiS
//
// See also license.opensource.txt in the distribution.


Thomas Witkowski's avatar
Thomas Witkowski committed
13
14
15
16
17
18
#include <algorithm>
#include <set>
#include <map>

#include "time.h"

19
20
21
22
#include "io/MacroReader.h"
#include "io/MacroInfo.h"
#include "io/MacroWriter.h"

23
24
25
26
#include "AdaptStationary.h"
#include "AdaptInstationary.h"
#include "FiniteElemSpace.h"
#include "ElementData.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
27
#include "ElementDofIterator.h"
28
29
30
31
32
33
34
35
36
37
38
#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
39
#include "ElInfoStack.h"
40
#include "Serializer.h"
41
#include "Lagrange.h"
42
43
44
45
46
47
48
49
50

namespace AMDiS {

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

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

51
52
53
54
55
56
57
58
59
60
61
62
63
64
  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);
65
66
67
68
69
70
71
72
73
74
75

  //**************************************************************************
  //  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;
76
77
  const Flag Mesh::CALL_MG_LEVEL           = 0X4000L;  
  const Flag Mesh::CALL_REVERSE_MODE       = 0X8000L;
78
79


Thomas Witkowski's avatar
Thomas Witkowski committed
80
  std::vector<DegreeOfFreedom> Mesh::dof_used;
81
  const int Mesh::MAX_DOF = 100;
82
  std::map<std::pair<DegreeOfFreedom, int>, DegreeOfFreedom*> Mesh::serializedDOFs;
83

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

106
    FUNCNAME("Mesh::Mesh()");
107
108
109
110

    // set default element prototype
    switch(dim) {
    case 1:
Thomas Witkowski's avatar
Thomas Witkowski committed
111
      elementPrototype = new Line(this);
112
113
      break;
    case 2:
Thomas Witkowski's avatar
Thomas Witkowski committed
114
      elementPrototype = new Triangle(this);
115
116
      break;
    case 3:
Thomas Witkowski's avatar
Thomas Witkowski committed
117
      elementPrototype = new Tetrahedron(this);
118
119
120
121
122
123
124
      break;
    default:
      ERROR_EXIT("invalid dimension\n");
    }

    elementPrototype->setIndex(-1);

125
126
    elementIndex = 0;
  }
127

128

129
  Mesh::~Mesh()
130
  {
131
    deleteMeshStructure();
132

133
134
135
136
137
138
    if (macroFileInfo != NULL)
      delete macroFileInfo;    
    if (elementPrototype)
      delete elementPrototype;    
    if (elementDataPrototype)
      delete elementDataPrototype;    
139
    
140
    for (unsigned int i = 0; i < admin.size(); i++)
141
      delete admin[i];    
142
  }
143

144

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

181

182
    /* ====================== Copy macro elements =================== */
183
  
184
185
186
187
188
189
190
191
192
    // 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;
193

194
195
    macroElements.clear();

196
197
198
    // 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();
199
	 it != m.macroElements.end(); ++it, insertCounter++) {
200

201
      // Create new MacroElement.
Thomas Witkowski's avatar
Thomas Witkowski committed
202
      MacroElement *el = new MacroElement(dim);
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
      // 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.
230
	if((*it)->getNeighbour(i)!=NULL) {
231
232
233
	macroElements[insertCounter]->
	  setNeighbour(i, macroElements[mapIndex[(*it)->getNeighbour(i)->getIndex()]]);
      }
234
      }
235
236
237
238
239
240
241
242
243
244
245
246
247
248
    }

    // 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;
  }

249

250
251
252
253
254
255
256
257
258
259
260
261
  void Mesh::updateNumberOfLeaves()
  {
    nLeaves = 0;

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

262

263
264
265
266
267
268
  void Mesh::addMacroElement(MacroElement* me) 
  {
    macroElements.push_back(me); 
    me->setIndex(macroElements.size());
  }

269

270
  void Mesh::removeMacroElements(std::set<MacroElement*>& macros,
Thomas Witkowski's avatar
Thomas Witkowski committed
271
				 const FiniteElemSpace *feSpace) 
272
273
274
  {
    FUNCNAME("Mesh::removeMacroElement()");

Thomas Witkowski's avatar
Thomas Witkowski committed
275
276
277
278
279
280
    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");

281
282
    // === Determine to all DOFs in mesh the macro elements where the DOF  ===
    // === is part of.                                                     ===
283

284
    // Map that stores for each dof pointer (which may have a list of dofs)
285
    // all macro element indices that own this dof.
Thomas Witkowski's avatar
Thomas Witkowski committed
286
287
    DofElMap dofsOwner;
    DofPosMap dofsPosIndex;
288

289
290
291
292
293
    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
294
      do {
295
	dofsOwner[elDofIter.getDofPtr()].insert(elInfo->getMacroElement());
Thomas Witkowski's avatar
Thomas Witkowski committed
296
297
	dofsPosIndex[elDofIter.getDofPtr()] = elDofIter.getPosIndex();
      } while (elDofIter.nextStrict());
298
299
300
301

      elInfo = stack.traverseNext(elInfo);
    }		   

Thomas Witkowski's avatar
Thomas Witkowski committed
302

303
304
305
306
307
308
309
310
311
312
313
314
315
    // === 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())
316
	newMacroElements.push_back(*elIter);     
317
    }
318

319
320
321
322
    // And replace the macro element list with the new one.
    macroElements.clear();
    macroElements = newMacroElements;

Thomas Witkowski's avatar
Thomas Witkowski committed
323

324
    // === For all macro elements to be deleted, delete them also to be neighbours ===
325
326
    // === of some other macro elements. Furtheremore, delete the whole element    ===
    // === hierarchie structure of the macro element.                              ===
327
328
329
    
    for (std::set<MacroElement*>::iterator macroIt = macros.begin();
	 macroIt != macros.end(); ++macroIt) {
Thomas Witkowski's avatar
Thomas Witkowski committed
330

331
332
      // Go through all neighbours of the macro element and remove this macro element
      // to be neighbour of some other macro element.
333
334
      for (int i = 0; i <= dim; i++)
	if ((*macroIt)->getNeighbour(i))
Thomas Witkowski's avatar
Thomas Witkowski committed
335
336
	  for (int j = 0; j <= dim; j++)
	    if ((*macroIt)->getNeighbour(i)->getNeighbour(j) == *macroIt)
337
338
	      (*macroIt)->getNeighbour(i)->setNeighbour(j, NULL);

339
340
341
342
      // Delete element hierarchie
      if (!(*macroIt)->getElement()->isLeaf()) {
	delete (*macroIt)->getElement()->getChild(0);
	delete (*macroIt)->getElement()->getChild(1);
343
344
345
346
347

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

	(*macroIt)->getElement()->setElementData(elementDataPrototype->clone()); 
348
      }
349
350
351
352
353

      // 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);
354
    }     
355
356


357
358
    // === Check now all the dofs, that have no owner anymore and therefore have ===
    // === to be removed.                                                        ===
359

Thomas Witkowski's avatar
Thomas Witkowski committed
360
    for (DofElMap::iterator dofsIt = dofsOwner.begin(); 
361
362
363
364
365
366
367
368
369
370
371
372
373
374
	 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)
375
	freeDof(const_cast<DegreeOfFreedom*>(dofsIt->first), 
376
		dofsPosIndex[dofsIt->first]);      
377
378
    }

379

380
    // === Update number of elements, vertices, etc. ===
381

382
383
384
385
386
387
388
389
    nLeaves = 0;
    nElements = 0;

    std::set<const DegreeOfFreedom*> allVertices;

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

391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
      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.                               ===
407
  }
408

409

410
411
412
413
414
415
  void Mesh::addDOFAdmin(DOFAdmin *localAdmin)
  {    
    FUNCNAME("Mesh::addDOFAdmin()");

    localAdmin->setMesh(this);

416
    TEST_EXIT(std::find(admin.begin(), admin.end(), localAdmin) == admin.end())
Thomas Witkowski's avatar
Thomas Witkowski committed
417
418
      ("admin %s is already associated to mesh %s\n",
       localAdmin->getName().c_str(), this->getName().c_str());
419
420
421

    admin.push_back(localAdmin);

422
    nDofEl = 0;
423

424
    localAdmin->setNumberOfPreDofs(VERTEX, nDof[VERTEX]);
425
426
    nDof[VERTEX] += localAdmin->getNumberOfDofs(VERTEX);
    nDofEl += getGeo(VERTEX) * nDof[VERTEX];
427

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

434
435
    localAdmin->setNumberOfPreDofs(CENTER, nDof[CENTER]);
    nDof[CENTER] += localAdmin->getNumberOfDofs(CENTER);
436
    nDofEl += nDof[CENTER];
437

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

440
441
    node[VERTEX] = 0;
    nNodeEl = getGeo(VERTEX);
442

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
463

464
465
  void Mesh::dofCompress()
  {
466
    FUNCNAME("Mesh::dofCompress()");
467

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

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

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

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

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


501
  DegreeOfFreedom *Mesh::getDof(GeoIndex position)
502
  {
503
    FUNCNAME("Mesh::getDof()");
504

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

508
    int ndof = getNumberOfDofs(position);
509
    if (ndof <= 0) 
510
      return NULL;
511

512
    DegreeOfFreedom *dof = new DegreeOfFreedom[ndof];
513

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


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

    if (nNodeEl <= 0)
536
      return NULL;
537

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

542
    return ptrs;
543
544
  }

545

546
  void Mesh::freeDofPtrs(DegreeOfFreedom **ptrs)
547
  {
548
    FUNCNAME("Mesh::freeDofPtrs()");
549

550
    TEST_EXIT_DBG(ptrs)("ptrs is NULL!\n");
551
552
553
554

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


559
  const DOFAdmin *Mesh::createDOFAdmin(std::string lname, DimVec<int> lnDof)
560
  {
561
    FUNCNAME("Mesh::createDOFAdmin()");
562
    
Thomas Witkowski's avatar
Thomas Witkowski committed
563
    DOFAdmin *localAdmin = new DOFAdmin(this, lname);
564

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

    addDOFAdmin(localAdmin);

570
    return localAdmin;
571
572
573
574
575
576
577
  }


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

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

587
    return localAdmin;
588
589
  }

590

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

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

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

608
609
      return;
    }
610

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

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

626
    delete [] dof;
627
    dof = NULL;
628
629
  }

630

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


  Element* Mesh::createNewElement(Element *parent)
  {
    FUNCNAME("Mesh::createNewElement()");
641
642

    TEST_EXIT_DBG(elementPrototype)("no element prototype\n");
643
644

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

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

651
652
653
    return el;
  }

654

655
656
  ElInfo* Mesh::createNewElInfo()
  {
657
658
    FUNCNAME("Mesh::createNewElInfo()");

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

675

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

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

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

    mel_info->fillMacroInfo(mel);

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

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

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

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

728
    return inside;
729
730
  }

731

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

    *elp = el_info->getElement();

Thomas Witkowski's avatar
Thomas Witkowski committed
744
    delete el_info;
745

746
    return val;
747
748
  }

749

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

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

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

782
	  return false;  /* ??? */
783
	} else {
784
	  return false;
785
	}
786
      }
787
788
    }

789
    ElInfo *c_el_info = createNewElInfo();
790

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

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

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

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

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

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

897
    return inside; 
898
899
  }

900
901

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