Mesh.cc 44.5 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
/******************************************************************************
 *
 * AMDiS - Adaptive multidimensional simulations
 *
 * Copyright (C) 2013 Dresden University of Technology. All Rights Reserved.
 * Web: https://fusionforge.zih.tu-dresden.de/projects/amdis
 *
 * Authors: 
 * Simon Vey, Thomas Witkowski, Andreas Naumann, Simon Praetorius, et al.
 *
 * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
 * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
 *
 *
 * This file is part of AMDiS
 *
 * See also license.opensource.txt in the distribution.
 * 
 ******************************************************************************/
20
21


Thomas Witkowski's avatar
Thomas Witkowski committed
22
23
24
25
26
27
#include <algorithm>
#include <set>
#include <map>

#include "time.h"

28
#include "io/Reader.h"
29
30
31
32
#include "io/MacroReader.h"
#include "io/MacroInfo.h"
#include "io/MacroWriter.h"

33
34
35
36
#include "AdaptStationary.h"
#include "AdaptInstationary.h"
#include "FiniteElemSpace.h"
#include "ElementData.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
37
#include "ElementDofIterator.h"
38
39
40
#include "MacroElement.h"
#include "Mesh.h"
#include "Traverse.h"
41
#include "Initfile.h"
42
43
44
45
46
47
#include "FixVec.h"
#include "DOFVector.h"
#include "CoarseningManager.h"
#include "DOFIterator.h"
#include "VertexVector.h"
#include "Projection.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
48
#include "ElInfoStack.h"
49
#include "Serializer.h"
50
#include "Lagrange.h"
51

52
using namespace std;
53

54
namespace AMDiS {
55
56
57
58
59

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

60
61
62
63
64
65
66
67
68
69
70
71
72
73
  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);
74
75
76
77
78
79
80
81
82
83
84

  //**************************************************************************
  //  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;
85
86
  const Flag Mesh::CALL_MG_LEVEL           = 0X4000L;  
  const Flag Mesh::CALL_REVERSE_MODE       = 0X8000L;
87
88


89
  vector<DegreeOfFreedom> Mesh::dof_used;
90
  const int Mesh::MAX_DOF = 100;
91
  map<pair<DegreeOfFreedom, int>, DegreeOfFreedom*> Mesh::serializedDOFs;
92
  std::set<string> Mesh::refinedMeshNames;
93

94
  Mesh::Mesh(string aName, int dimension) 
95
96
97
98
99
100
    : name(aName), 
      dim(dimension), 
      nVertices(0),
      nEdges(0),
      nLeaves(0), 
      nElements(0),
101
      parametric(NULL), 
102
      preserveCoarseDOFs(false),
103
104
      nDofEl(0),
      nDof(dimension, DEFAULT_VALUE, 0),
105
106
      nNodeEl(0),
      node(dimension, DEFAULT_VALUE, 0),
107
108
      elementPrototype(NULL),
      elementDataPrototype(NULL),
109
110
      elementIndex(-1),
      initialized(false),
111
      macroFileInfo(NULL),
112
      changeIndex(0),
113
114
      final_lambda(dimension, DEFAULT_VALUE, 0.0)
  {
115
    FUNCNAME("Mesh::Mesh()");
116

117
118
119
120
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    nParallelPreRefinements = 0;
#endif

121
122
123
    // set default element prototype
    switch(dim) {
    case 1:
Thomas Witkowski's avatar
Thomas Witkowski committed
124
      elementPrototype = new Line(this);
125
126
      break;
    case 2:
Thomas Witkowski's avatar
Thomas Witkowski committed
127
      elementPrototype = new Triangle(this);
128
129
      break;
    case 3:
Thomas Witkowski's avatar
Thomas Witkowski committed
130
      elementPrototype = new Tetrahedron(this);
131
132
133
134
135
136
137
      break;
    default:
      ERROR_EXIT("invalid dimension\n");
    }

    elementPrototype->setIndex(-1);

138
139
    elementIndex = 0;
  }
140

141

142
  Mesh::~Mesh()
143
  {
144
    deleteMeshStructure();
145

146
    if (macroFileInfo != NULL)
147
148
149
150
      delete macroFileInfo;    
    if (elementPrototype)
      delete elementPrototype;    
    if (elementDataPrototype)
151
      delete elementDataPrototype;        
152
  }
153

154

155
  Mesh& Mesh::operator=(const Mesh& m)
156
  {
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
    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;
172
    parametric = NULL;
173
174

    preserveCoarseDOFs = m.preserveCoarseDOFs;
175
176
    nDofEl = m.nDofEl;
    nDof = m.nDof;
177
178
179
180
181
182
183
184
185
    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
186
      admin[i] = new DOFAdmin(this);
187
      *(admin[i]) = *(m.admin[i]);
188
189
      admin[i]->setMesh(this);
    }
190

191

192
    /* ====================== Copy macro elements =================== */
193
  
194
195
    // mapIndex[i] is the index of the MacroElement element in the vector
    // macroElements, for which holds: element->getIndex() = i    
196
    map<int, int> mapIndex;
197
198
199
200
201
202

    // We use this map for coping the DOFs of the Elements within the
    // MacroElements objects.
    Mesh::serializedDOFs.clear();

    int insertCounter = 0;
203

204
205
    macroElements.clear();

206
207
    // Go through all MacroElements of mesh m, and create for every a new
    // MacroElement in this mesh.
208
    for (deque<MacroElement*>::const_iterator it = m.macroElements.begin();
209
	 it != m.macroElements.end(); ++it, insertCounter++) {
210

211
      // Create new MacroElement.
Thomas Witkowski's avatar
Thomas Witkowski committed
212
      MacroElement *el = new MacroElement(dim);
213

214
215
216
217
218
219
220
221
222
223
      // 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.
224
      mapIndex.insert(pair<int, int>(el->getIndex(), insertCounter));
225
226
227
228
229
    }

    // Now we have to go through all the new MacroElements, and update the neighbour
    // connections.
    insertCounter = 0;
230
    for (deque<MacroElement*>::const_iterator it = m.macroElements.begin();
231
232
233
234
235
236
237
238
239
	 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.
240
	if ((*it)->getNeighbour(i)!=NULL) {
241
242
243
	macroElements[insertCounter]->
	  setNeighbour(i, macroElements[mapIndex[(*it)->getNeighbour(i)->getIndex()]]);
      }
244
      }
245
246
247
248
249
250
251
    }

    // Cleanup
    Mesh::serializedDOFs.clear();

    /* ================== Things will be done when required ============ */
      
252
253
    TEST_EXIT(elementDataPrototype == NULL)("TODO\n");
    TEST_EXIT(m.parametric == NULL)("TODO\n");
254
255
256
257
258
    TEST_EXIT(periodicAssociations.size() == 0)("TODO\n");

    return *this;
  }

259

260
261
262
263
264
265
266
267
268
269
270
271
  void Mesh::updateNumberOfLeaves()
  {
    nLeaves = 0;

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

272

273
274
275
276
277
278
  void Mesh::addMacroElement(MacroElement* me) 
  {
    macroElements.push_back(me); 
    me->setIndex(macroElements.size());
  }

279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
  void Mesh::removeAllMacroElements()
  {
    // Delete all the dofs
    Element::deletedDOFs.clear();
    for (deque<MacroElement*>::const_iterator it = macroElements.begin();
	 it != macroElements.end(); ++it) {
      (*it)->getElement()->deleteElementDOFs();
    }
    Element::deletedDOFs.clear();
    
    // Set all neighbors null
    for (deque<MacroElement*>::const_iterator macroIt = macroElements.begin();
	 macroIt != macroElements.end(); ++macroIt) {

      for (int i = 0; i < getGeo(NEIGH); i++)
	(*macroIt)->setNeighbour(i, NULL);

      Element *mel = (*macroIt)->getElement();
      // Delete element hierarchie
      if (!(mel->isLeaf())) {
	delete mel->getChild(0);
	delete mel->getChild(1);

	mel->child[0] = NULL;
	mel->child[1] = NULL;

	mel->setElementData(elementDataPrototype->clone()); 
      }
    }
    
    macroElements.clear();
    nLeaves = 0;
    nElements = 0;
    nVertices = 0;
    
    for (size_t i = 0; i < admin.size(); i++)
      {
	TEST_EXIT_DBG(admin[i]->getUsedSize() == admin[i]->getHoleCount())
	  ("All macro elements has been removed. But not all dofs are cleaned. (UsedSize = %d, HoleCount = %d)\n", 
	    admin[i]->getUsedSize(), admin[i]->getHoleCount());
	  
	admin[i]->reset();
      }
  }
323
  
324
  void Mesh::removeMacroElements(std::set<MacroElement*>& delMacros,
325
				 vector<const FiniteElemSpace*>& feSpaces) 
326
327
328
  {
    FUNCNAME("Mesh::removeMacroElement()");

329
330
    typedef map<const DegreeOfFreedom*, std::set<MacroElement*> > DofElMap;
    typedef map<const DegreeOfFreedom*, GeoIndex> DofPosMap;
Thomas Witkowski's avatar
Thomas Witkowski committed
331

332
    TEST_EXIT(feSpaces.size() > 0)("Should not happen!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
333

334
335
336
337
    // Search for the FE space with the highest degree of polynomials. Using this
    // FE space ensures that deleting DOFs defined on it, also DOFs of lower
    // order FE spaces will be deleted correctly.
    const FiniteElemSpace *feSpace = FiniteElemSpace::getHighest(feSpaces);
338
339

    
340
341
    // === Determine to all DOFs in mesh the macro elements where the DOF  ===
    // === is part of.                                                     ===
342

343
344
    // Map that stores for each DOF pointer (which may have a list of DOFs)
    // all macro element indices that own this DOF.
Thomas Witkowski's avatar
Thomas Witkowski committed
345
346
    DofElMap dofsOwner;
    DofPosMap dofsPosIndex;
347

348
349
350
351
352
    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
353
      do {
354
355
	dofsOwner[elDofIter.getBaseDof()].insert(elInfo->getMacroElement());
	dofsPosIndex[elDofIter.getBaseDof()] = elDofIter.getPosIndex();
Thomas Witkowski's avatar
Thomas Witkowski committed
356
      } while (elDofIter.nextStrict());
357
358
359
360

      elInfo = stack.traverseNext(elInfo);
    }		   

Thomas Witkowski's avatar
Thomas Witkowski committed
361

362
363
    // === Remove macro elements from mesh macro element list. ===

364
    // Removing arbitrary elements from an deque is very slow. Therefore, we
365
366
367
    // 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.

368
369
    deque<MacroElement*> newMacroElements;
    for (deque<MacroElement*>::iterator elIter = macroElements.begin();
370
	 elIter != macroElements.end(); ++elIter) {
371
372
373
      // 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.
374
      if (delMacros.find(*elIter) == delMacros.end())
375
	newMacroElements.push_back(*elIter);     
376
    }
377

378
379
380
381
    // And replace the macro element list with the new one.
    macroElements.clear();
    macroElements = newMacroElements;

Thomas Witkowski's avatar
Thomas Witkowski committed
382

383
384
385
    // === For all macro elements to be deleted, delete them also to be       ===
    // === neighbours of some other macro elements. Furtheremore, delete the  ===
    // === whole element hierarchie structure of the macro element.           ===
386
    
387
388
    for (std::set<MacroElement*>::iterator macroIt = delMacros.begin();
	 macroIt != delMacros.end(); ++macroIt) {
Thomas Witkowski's avatar
Thomas Witkowski committed
389

390
391
      // Go through all neighbours of the macro element and remove this macro
      // element to be neighbour of some other macro element.
392
      for (int i = 0; i < getGeo(NEIGH); i++)
393
	if ((*macroIt)->getNeighbour(i))
394
	  for (int j = 0; j < getGeo(NEIGH); j++)
Thomas Witkowski's avatar
Thomas Witkowski committed
395
	    if ((*macroIt)->getNeighbour(i)->getNeighbour(j) == *macroIt)
396
	      (*macroIt)->getNeighbour(i)->setNeighbour(j, NULL);
397

398

399
      Element *mel = (*macroIt)->getElement();
400
      // Delete element hierarchie
401
402
403
      if (!(mel->isLeaf())) {
	delete mel->getChild(0);
	delete mel->getChild(1);
404

405
406
	mel->child[0] = NULL;
	mel->child[1] = NULL;
407

408
	mel->setElementData(elementDataPrototype->clone()); 
409
      }
410

411
      mel->delDofPtr();
412
    }     
413
414


415
416
    // === Check now all the DOFs that have no owner anymore and therefore  ===
    // === have to be removed.                                              ===
417

Thomas Witkowski's avatar
Thomas Witkowski committed
418
    for (DofElMap::iterator dofsIt = dofsOwner.begin(); 
419
420
421
422
423
424
	 dofsIt != dofsOwner.end(); ++dofsIt) {
      
      bool deleteDof = true;

      for (std::set<MacroElement*>::iterator elIter = dofsIt->second.begin();
	   elIter != dofsIt->second.end(); ++elIter) {
425
426
	std::set<MacroElement*>::iterator mIt = delMacros.find(*elIter);
	if (mIt == delMacros.end()) {
427
428
429
430
431
432
	  deleteDof = false;
	  break;
	}
      }

      if (deleteDof)
433
	freeDof(const_cast<DegreeOfFreedom*>(dofsIt->first), 
434
		dofsPosIndex[dofsIt->first]);      
435
436
    }

437

438
    // === Update number of elements, vertices, etc. ===
439

440
441
    nLeaves = 0;
    nElements = 0;
442
    nVertices = 0;
443

444
445
    if (!macroElements.empty()) {
      std::set<const DegreeOfFreedom*> allVertices;
446

447
448
449
      elInfo = stack.traverseFirst(this, -1, Mesh::CALL_EVERY_EL_PREORDER);
      while (elInfo) {
	nElements++;
450

451
452
	if (elInfo->getElement()->isLeaf()) {
	  nLeaves++;
453

454
455
456
457
458
	  for (int i = 0; i < getGeo(VERTEX); i++)
	    allVertices.insert(elInfo->getElement()->getDof(i));
	}

	elInfo = stack.traverseNext(elInfo);
459
460
      }

461
462
463
464
465
466
467
468
469
470
471
      nVertices = allVertices.size();
    } else {
      
      for (size_t i = 0; i < admin.size(); i++)
      {
	TEST_EXIT_DBG(admin[i]->getUsedSize() == admin[i]->getHoleCount())
	  ("All macro elements has been removed. But not all dofs are cleaned. (UsedSize = %d, HoleCount = %d)\n", 
	    admin[i]->getUsedSize(), admin[i]->getHoleCount());
	  
	admin[i]->reset();
      }
472
473
474
475
476
    }
    // === 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.                               ===
477
  }
478

479

480
481
482
483
484
485
  void Mesh::addDOFAdmin(DOFAdmin *localAdmin)
  {    
    FUNCNAME("Mesh::addDOFAdmin()");

    localAdmin->setMesh(this);

486
    TEST_EXIT(find(admin.begin(), admin.end(), localAdmin) == admin.end())
Thomas Witkowski's avatar
Thomas Witkowski committed
487
488
      ("admin %s is already associated to mesh %s\n",
       localAdmin->getName().c_str(), this->getName().c_str());
489
490
491

    admin.push_back(localAdmin);

492
    nDofEl = 0;
493

494
    localAdmin->setNumberOfPreDofs(VERTEX, nDof[VERTEX]);
495
496
    nDof[VERTEX] += localAdmin->getNumberOfDofs(VERTEX);
    nDofEl += getGeo(VERTEX) * nDof[VERTEX];
497

498
    if (dim > 1) {
499
      localAdmin->setNumberOfPreDofs(EDGE, nDof[EDGE]);
500
501
      nDof[EDGE] += localAdmin->getNumberOfDofs(EDGE);
      nDofEl += getGeo(EDGE) * nDof[EDGE];
502
503
    }

504
505
    localAdmin->setNumberOfPreDofs(CENTER, nDof[CENTER]);
    nDof[CENTER] += localAdmin->getNumberOfDofs(CENTER);
506
    nDofEl += nDof[CENTER];
507

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

510
511
    node[VERTEX] = 0;
    nNodeEl = getGeo(VERTEX);
512

513
514
    if (dim > 1) {
      node[EDGE] = nNodeEl;
515
      if (nDof[EDGE] > 0) 
516
	nNodeEl += getGeo(EDGE);
517
518
    }

519
    if (dim == 3) {
520
      localAdmin->setNumberOfPreDofs(FACE, nDof[FACE]);
521
522
      nDof[FACE] += localAdmin->getNumberOfDofs(FACE);
      nDofEl += getGeo(FACE) * nDof[FACE];
523
      node[FACE] = nNodeEl;
524
      if (nDof[FACE] > 0) 
525
	nNodeEl += getGeo(FACE);
526
527
    }

528
    node[CENTER] = nNodeEl;
529
    if (nDof[CENTER] > 0)
530
      nNodeEl += 1;
531
532
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
533

534
535
  void Mesh::dofCompress()
  {
536
    FUNCNAME_DBG("Mesh::dofCompress()");
537

Thomas Witkowski's avatar
Thomas Witkowski committed
538
539
    for (unsigned int iadmin = 0; iadmin < admin.size(); iadmin++) {
      DOFAdmin* compressAdmin = admin[iadmin];
540
541

      TEST_EXIT_DBG(compressAdmin)("no admin[%d] in mesh\n", iadmin);
542
      
Thomas Witkowski's avatar
Thomas Witkowski committed
543
      int size = compressAdmin->getSize();
Thomas Witkowski's avatar
Thomas Witkowski committed
544
      if (size < 1 || 
545
	  compressAdmin->getUsedDofs() < 1 || 
Thomas Witkowski's avatar
Thomas Witkowski committed
546
	  compressAdmin->getHoleCount() < 1)    
547
	continue;
Thomas Witkowski's avatar
Thomas Witkowski committed
548

549
      vector<DegreeOfFreedom> newDofIndex(size);     
550
      compressAdmin->compress(newDofIndex);
Thomas Witkowski's avatar
Thomas Witkowski committed
551
552
553
554

      Flag fill_flag = (preserveCoarseDOFs ?  
			Mesh::CALL_EVERY_EL_PREORDER | Mesh::FILL_NOTHING :
			Mesh::CALL_LEAF_EL | Mesh::FILL_NOTHING);          
555
556
557
      TraverseStack stack;
      ElInfo *elInfo = stack.traverseFirst(this, -1, fill_flag);
      while (elInfo) {
558
	elInfo->getElement()->newDofFct1(compressAdmin, newDofIndex);
559
560
561
562
563
	elInfo = stack.traverseNext(elInfo);
      }

      elInfo = stack.traverseFirst(this, -1, fill_flag);
      while (elInfo) {
564
	elInfo->getElement()->newDofFct2(compressAdmin);
565
566
	elInfo = stack.traverseNext(elInfo);
      }
567
    }       
568
569
570
  }


571
  DegreeOfFreedom *Mesh::getDof(GeoIndex position)
572
  {
573
    FUNCNAME_DBG("Mesh::getDof()");
574

575
    TEST_EXIT_DBG(position >= CENTER && position <= FACE)
576
      ("unknown position %d\n", position);
577

578
    int ndof = getNumberOfDofs(position);
579
    if (ndof <= 0) 
580
      return NULL;
581

582
    DegreeOfFreedom *dof = new DegreeOfFreedom[ndof];
583

584
    for (int i = 0; i < getNumberOfDOFAdmin(); i++) {
585
      const DOFAdmin *localAdmin = &getDofAdmin(i);
586
      TEST_EXIT_DBG(localAdmin)("no admin[%d]\n", i);
587
      
588
      int n  = localAdmin->getNumberOfDofs(position);
589
      int n0 = localAdmin->getNumberOfPreDofs(position);
590
      
591
592
      TEST_EXIT_DBG(n + n0 <= ndof)
	("n = %d, n0 = %d too large: ndof = %d\n", n, n0, ndof);
593
      
594
      for (int j = 0; j < n; j++)
595
596
	dof[n0 + j] = const_cast<DOFAdmin*>(localAdmin)->getDOFIndex();
    }
597
    
598
    return dof;
599
600
601
  }


Thomas Witkowski's avatar
Thomas Witkowski committed
602
  DegreeOfFreedom **Mesh::createDofPtrs()
603
604
  {
    if (nNodeEl <= 0)
605
      return NULL;
606

607
    DegreeOfFreedom **ptrs = new DegreeOfFreedom*[nNodeEl];
608
    for (int i = 0; i < nNodeEl; i++)
609
      ptrs[i] = NULL;
610

611
    return ptrs;
612
613
  }

614

615
  void Mesh::freeDofPtrs(DegreeOfFreedom **ptrs)
616
  {
617
    FUNCNAME_DBG("Mesh::freeDofPtrs()");
618

619
    TEST_EXIT_DBG(ptrs)("ptrs is NULL!\n");
620
621
622
623

    if (nNodeEl <= 0)
      return;
  
624
    delete [] ptrs;
625
    ptrs = NULL;
626
627
628
  }


629
  const DOFAdmin *Mesh::createDOFAdmin(string lname, DimVec<int> lnDof)
630
  {    
Thomas Witkowski's avatar
Thomas Witkowski committed
631
    DOFAdmin *localAdmin = new DOFAdmin(this, lname);
632

Thomas Witkowski's avatar
Thomas Witkowski committed
633
    for (int i = 0; i < dim + 1; i++)
634
      localAdmin->setNumberOfDofs(i, lnDof[i]);
635
636
637

    addDOFAdmin(localAdmin);

638
    return localAdmin;
639
640
641
642
643
  }


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

Thomas Witkowski's avatar
Thomas Witkowski committed
646
    for (unsigned int i = 0; i < admin.size(); i++) {
647
      if (admin[i]->getNumberOfDofs(VERTEX)) {
648
649
650
651
	if (!localAdmin)  
	  localAdmin = admin[i];
	else if (admin[i]->getSize() < localAdmin->getSize())
	  localAdmin = admin[i];
652
      }
653
654
    }

655
    return localAdmin;
656
657
  }

658

659
  void Mesh::freeDof(DegreeOfFreedom* dof, GeoIndex position)
660
  {
661
    FUNCNAME_DBG("Mesh::freeDof()");
662

663
    TEST_EXIT_DBG(position >= CENTER && position <= FACE)
664
      ("unknown position %d\n", position);
665

666
    if (nDof[position]) {
667
      TEST_EXIT_DBG(dof != NULL)("dof = NULL, but ndof = %d\n", nDof[position]);
668
    } else  {
669
      TEST_EXIT_DBG(dof == NULL)("dof != NULL, but ndof = 0\n");
670
    }
671

672
673
    TEST_EXIT_DBG(nDof[position] <= MAX_DOF)
      ("ndof too big: ndof = %d, MAX_DOF = %d\n", nDof[position], MAX_DOF);
674

675
    for (unsigned int i = 0; i < admin.size(); i++) {
676
677
      int n = admin[i]->getNumberOfDofs(position);
      int n0 = admin[i]->getNumberOfPreDofs(position);
678
      
679
680
      TEST_EXIT_DBG(n + n0 <= nDof[position])
	("n = %d, n0 = %d too large: ndof = %d\n", n, n0, nDof[position]);
681
682
      
      for (int j = 0; j < n; j++)
683
	admin[i]->freeDofIndex(dof[n0 + j]);      
684
    }
685
    delete [] dof;
686
687
  }

688

689
690
  void Mesh::freeElement(Element* el)
  {
691
    freeDofPtrs(const_cast<DegreeOfFreedom**>(el->getDof()));
Thomas Witkowski's avatar
Thomas Witkowski committed
692
    delete el;
693
694
695
696
697
  }


  Element* Mesh::createNewElement(Element *parent)
  {
698
    FUNCNAME_DBG("Mesh::createNewElement()");
699
700

    TEST_EXIT_DBG(elementPrototype)("no element prototype\n");
701
702

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

Thomas Witkowski's avatar
Thomas Witkowski committed
704
    if (!parent && elementDataPrototype)
705
      el->setElementData(elementDataPrototype->clone()); 
Thomas Witkowski's avatar
Thomas Witkowski committed
706
    else
707
      el->setElementData(NULL); // must be done in ElementData::refineElementData()
708

709
710
711
    return el;
  }

712

713
714
  ElInfo* Mesh::createNewElInfo()
  {
715
716
    FUNCNAME("Mesh::createNewElInfo()");

717
    switch (dim) {
718
    case 1:
Thomas Witkowski's avatar
Thomas Witkowski committed
719
      return new ElInfo1d(this);
720
721
      break;
    case 2:
Thomas Witkowski's avatar
Thomas Witkowski committed
722
      return new ElInfo2d(this);
723
724
      break;
    case 3:
Thomas Witkowski's avatar
Thomas Witkowski committed
725
      return new ElInfo3d(this);
726
727
      break;
    default:
728
      ERROR_EXIT("invalid dim [%d]\n",dim);
729
      return NULL;
730
    }
731
732
  }

733

734
735
  bool Mesh::findElInfoAtPoint(const WorldVector<double>& xy,
			       ElInfo *el_info,
736
737
			       DimVec<double>& bary,
			       const MacroElement *start_mel,
738
			       const WorldVector<double> *xy0,
739
			       double *sp)
740
  {
741
    static const MacroElement *mel = NULL;
742
    DimVec<double> lambda(dim, NO_INIT);
Thomas Witkowski's avatar
Thomas Witkowski committed
743
    ElInfo *mel_info = createNewElInfo();
744

745
    if (start_mel != NULL)
746
747
      mel = start_mel;
    else
748
      if (mel == NULL || mel->getElement()->getMesh() != this)
749
750
751
	mel = *(macroElements.begin());

    mel_info->setFillFlag(Mesh::FILL_COORDS);
752
    g_xy = &xy;
753
    g_xy0 = xy0;
754
    g_sp = sp;
755
756
757

    mel_info->fillMacroInfo(mel);

758
759
760
761
762
    // 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;
763
764
765
    std::stack<MacroElement*> active;
    
//     macrosVisited.insert(mel->getIndex());
766

767
    int k;
768
    while ((k = mel_info->worldToCoord(xy, &lambda)) >= 0) {
769
      macrosVisited.insert(mel->getIndex());
770
      
771
      if (mel->getNeighbour(k) && !macrosVisited.count(mel->getNeighbour(k)->getIndex())) {
772
	// look for next macro-element in the direction of the coordinates xy
773
	mel = mel->getNeighbour(k);
774
775
	mel_info->fillMacroInfo(mel);
	continue;
776
      } else {
777
	// consider all neighbors of the current macro-element to visit next
778
	for (int i = 0; i < dim + 1; ++i) {
779
	  if (i !=  k && mel->getNeighbour(i) && !macrosVisited.count(mel->getNeighbour(i)->getIndex()))
780
781
	    active.push(mel->getNeighbour(i));    
	}
782
783
784
785
	// if all neighbors are visited already
	if (active.empty()) { 	  
	  if (macrosVisited.size() == static_cast<size_t>(getNumberOfMacros())) {
	    // if all macro-elements are visited -> no element found!
786
787
788
	    delete mel_info;
	    return false;
	  } else {
789
	    // go to an arbitrary macro-element to continue the search
790
791
792
	    deque<MacroElement*>::iterator it;
	    bool found = false;
	    for (it = firstMacroElement(); it != endOfMacroElements(); it++) {
793
	      if (!macrosVisited.count((*it)->getIndex())) {
794
795
796
797
		active.push(*it);
		found = true;
	      }
	    }
798
799
800
801
	    if (!found) {
	      delete mel_info;
	      return false;
	    }
802
	  }
803
	}
804
805
	mel = active.top();
	active.pop();
806
807

	mel_info->fillMacroInfo(mel);
808
809
810
811
      }
    }

    /* now, descend in tree to find leaf element at point */
812
    bool inside = findElementAtPointRecursive(mel_info, lambda, k, el_info);
813
814
    for (int i = 0; i <= dim; i++)
      bary[i] = final_lambda[i];   
815
  
Thomas Witkowski's avatar
Thomas Witkowski committed
816
    delete mel_info;
817

818
    return inside;
819
820
  }

821

822
  bool Mesh::findElementAtPoint(const WorldVector<double>&  xy,
823
824
				Element **elp, 
				DimVec<double>& bary,
825
				const MacroElement *start_mel,
826
827
				const WorldVector<double> *xy0,
				double *sp)
828
  {
829
830
    ElInfo *el_info = createNewElInfo();
    int val = findElInfoAtPoint(xy, el_info, bary, start_mel, xy0, sp);
831
832
833

    *elp = el_info->getElement();

Thomas Witkowski's avatar
Thomas Witkowski committed
834
    delete el_info;
835

836
    return val;
837
838
  }

839

840
  bool Mesh::findElementAtPointRecursive(ElInfo *el_info,
841
					 const DimVec<double>& lambda,