Mesh.cc 38.8 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
273
    typedef std::map<const DegreeOfFreedom*, std::set<MacroElement*> > DofElMap;
    typedef std::map<const DegreeOfFreedom*, GeoIndex> DofPosMap;

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

    ElementDofIterator elDofIter(feSpace);
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


    // === Determine to all dofs the macro elements where the dof is part of. ===

283
    for (std::deque<MacroElement*>::iterator macroIt = macroElements.begin();
284
	 macroIt != macroElements.end(); ++macroIt) {
Thomas Witkowski's avatar
Thomas Witkowski committed
285
286
287
288
289
      elDofIter.reset((*macroIt)->getElement());
      do {
	dofsOwner[elDofIter.getDofPtr()].insert(*macroIt);
	dofsPosIndex[elDofIter.getDofPtr()] = elDofIter.getPosIndex();
      } while (elDofIter.nextStrict());
290
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
291

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

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

318
319
320
321
      // Go through all neighbours of the macro element and remove this macro element
      // to be neighbour of some other macro element.
      for (int i = 0; i <= dim; i++) {
	if ((*macroIt)->getNeighbour(i)) {
Thomas Witkowski's avatar
Thomas Witkowski committed
322
323
	  for (int j = 0; j <= dim; j++)
	    if ((*macroIt)->getNeighbour(i)->getNeighbour(j) == *macroIt)
324
325
326
327
328
	      (*macroIt)->getNeighbour(i)->setNeighbour(j, NULL);
	} else {
	  // There is no neighbour at this edge, so we have to decrease the number
	  // of edges in the mesh.
	  nEdges--;
329
	}
330
331
      }

332
      // Decrease also the number of elements of the mesh.
333
334
      nLeaves--;
      nElements--;
335
    }     
336
337


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

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

363
364
365
366
367
368
369

    // === Finally, remove the macro elements from memory. ===

    for (std::set<MacroElement*>::iterator macroIt = macros.begin();
	 macroIt != macros.end(); ++macroIt)
      delete *macroIt;

370
    nVertices = nRemainDofs;
371
  }
372

373

374
375
376
377
378
379
  void Mesh::addDOFAdmin(DOFAdmin *localAdmin)
  {    
    FUNCNAME("Mesh::addDOFAdmin()");

    localAdmin->setMesh(this);

380
    std::vector<DOFAdmin*>::iterator dai = 
381
      std::find(admin.begin(), admin.end(), localAdmin);
382

Thomas Witkowski's avatar
Thomas Witkowski committed
383
384
385
    TEST_EXIT(dai == admin.end())
      ("admin %s is already associated to mesh %s\n",
       localAdmin->getName().c_str(), this->getName().c_str());
386

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

390
391
392
393
394
    admin.push_back(localAdmin);

    nDOFEl = 0;

    localAdmin->setNumberOfPreDOFs(VERTEX,nDOF[VERTEX]);
395
    nDOF[VERTEX] += localAdmin->getNumberOfDOFs(VERTEX);
396
397
    nDOFEl += getGeo(VERTEX) * nDOF[VERTEX];

398
    if (dim > 1) {
399
      localAdmin->setNumberOfPreDOFs(EDGE,nDOF[EDGE]);
400
      nDOF[EDGE] += localAdmin->getNumberOfDOFs(EDGE);
401
402
403
404
405
406
407
      nDOFEl += getGeo(EDGE) * nDOF[EDGE];
    }

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

408
    TEST_EXIT_DBG(nDOF[VERTEX] > 0)("no vertex dofs\n");
409

410
411
    node[VERTEX] = 0;
    nNodeEl = getGeo(VERTEX);
412

413
414
    if (dim > 1) {
      node[EDGE] = nNodeEl;
415
416
      if (nDOF[EDGE] > 0) 
	nNodeEl += getGeo(EDGE);
417
418
    }

419
    if (dim == 3) {
420
      localAdmin->setNumberOfPreDOFs(FACE,nDOF[FACE]);
421
422
423
424
425
      nDOF[FACE] += localAdmin->getNumberOfDOFs(FACE);
      nDOFEl += getGeo(FACE) * nDOF[FACE];
      node[FACE] = nNodeEl;
      if (nDOF[FACE] > 0) 
	nNodeEl += getGeo(FACE);
426
427
    }

428
    node[CENTER] = nNodeEl;
429
    if (nDOF[CENTER] > 0)
430
      nNodeEl += 1;
431
432
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
433

434
435
  void Mesh::dofCompress()
  {
436
    FUNCNAME("Mesh::dofCompress()");
437

Thomas Witkowski's avatar
Thomas Witkowski committed
438
439
    for (unsigned int iadmin = 0; iadmin < admin.size(); iadmin++) {
      DOFAdmin* compressAdmin = admin[iadmin];
440
441

      TEST_EXIT_DBG(compressAdmin)("no admin[%d] in mesh\n", iadmin);
442
      
Thomas Witkowski's avatar
Thomas Witkowski committed
443
      int size = compressAdmin->getSize();
Thomas Witkowski's avatar
Thomas Witkowski committed
444
445
446
      if (size < 1 || 
	  compressAdmin->getUsedDOFs() < 1 || 
	  compressAdmin->getHoleCount() < 1)    
447
	continue;
Thomas Witkowski's avatar
Thomas Witkowski committed
448

449
450
451
      newDOF.resize(size);
      
      compressAdmin->compress(newDOF);
Thomas Witkowski's avatar
Thomas Witkowski committed
452
453
454
455

      Flag fill_flag = (preserveCoarseDOFs ?  
			Mesh::CALL_EVERY_EL_PREORDER | Mesh::FILL_NOTHING :
			Mesh::CALL_LEAF_EL | Mesh::FILL_NOTHING);          
456
457
458
      TraverseStack stack;
      ElInfo *elInfo = stack.traverseFirst(this, -1, fill_flag);
      while (elInfo) {
Thomas Witkowski's avatar
Thomas Witkowski committed
459
	elInfo->getElement()->newDOFFct1(compressAdmin);
460
461
462
463
464
	elInfo = stack.traverseNext(elInfo);
      }

      elInfo = stack.traverseFirst(this, -1, fill_flag);
      while (elInfo) {
Thomas Witkowski's avatar
Thomas Witkowski committed
465
	elInfo->getElement()->newDOFFct2(compressAdmin);
466
467
468
	elInfo = stack.traverseNext(elInfo);
      }

469
      newDOF.resize(0);
470
    }       
471
472
473
  }


474
  DegreeOfFreedom *Mesh::getDof(GeoIndex position)
475
  {
476
    FUNCNAME("Mesh::getDof()");
477

478
    TEST_EXIT_DBG(position >= CENTER && position <= FACE)
479
      ("unknown position %d\n", position);
480

481
482
    int ndof = getNumberOfDOFs(position);
    if (ndof <= 0) 
483
      return NULL;
484

485
    DegreeOfFreedom *dof = new DegreeOfFreedom[ndof];
486

487
488
    for (int i = 0; i < getNumberOfDOFAdmin(); i++) {
      const DOFAdmin *localAdmin = &getDOFAdmin(i);
489
      TEST_EXIT_DBG(localAdmin)("no admin[%d]\n", i);
490
491
492
493
      
      int n  = localAdmin->getNumberOfDOFs(position);
      int n0 = localAdmin->getNumberOfPreDOFs(position);
      
494
      TEST_EXIT_DBG(n + n0 <= ndof)("n=%d, n0=%d too large: ndof=%d\n", n, n0, ndof);
495
      
496
      for (int j = 0; j < n; j++)
497
498
	dof[n0 + j] = const_cast<DOFAdmin*>(localAdmin)->getDOFIndex();
    }
499
 
500
    return dof;
501
502
503
  }


Thomas Witkowski's avatar
Thomas Witkowski committed
504
  DegreeOfFreedom **Mesh::createDofPtrs()
505
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
506
    FUNCNAME("Mesh::createDofPtrs()");
507
508

    if (nNodeEl <= 0)
509
      return NULL;
510

511
    DegreeOfFreedom **ptrs = new DegreeOfFreedom*[nNodeEl];
512
    for (int i = 0; i < nNodeEl; i++)
513
514
      ptrs[i] = NULL;

515
    return ptrs;
516
517
  }

518

519
520
  void Mesh::freeDOFPtrs(DegreeOfFreedom **ptrs)
  {
521
    FUNCNAME("Mesh::freeDOFPtrs()");
522

523
    TEST_EXIT_DBG(ptrs)("ptrs is NULL!\n");
524
525
526
527

    if (nNodeEl <= 0)
      return;
  
528
    delete [] ptrs;
529
530
531
  }


Thomas Witkowski's avatar
Thomas Witkowski committed
532
  const DOFAdmin *Mesh::createDOFAdmin(std::string lname, DimVec<int> lnDOF)
533
  {
534
    FUNCNAME("Mesh::createDOFAdmin()");
535

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

Thomas Witkowski's avatar
Thomas Witkowski committed
538
539
    for (int i = 0; i < dim + 1; i++)
      localAdmin->setNumberOfDOFs(i, lnDOF[i]);
540
541
542

    addDOFAdmin(localAdmin);

543
    return localAdmin;
544
545
546
547
548
549
550
  }


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

551
552
553
554
555
556
    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];
557
      }
558
559
    }

560
    return localAdmin;
561
562
  }

563

564
  void Mesh::freeDof(DegreeOfFreedom* dof, GeoIndex position)
565
  {
566
    FUNCNAME("Mesh::freeDof()");
567

568
    TEST_EXIT_DBG(position >= CENTER && position <= FACE)
569
      ("unknown position %d\n", position);
570

571
572
573
    int ndof = nDOF[position];
    if (ndof) {
      if (!dof) {
574
	MSG("dof = NULL, but ndof = %d\n", ndof);
575
576
	return;
      }
577
    } else  {
578
579
580
      if (dof)
	MSG("dof != NULL, but ndof = 0\n");

581
582
      return;
    }
583

584
    TEST_EXIT_DBG(ndof <= MAX_DOF)
585
      ("ndof too big: ndof = %d, MAX_DOF = %d\n", ndof, MAX_DOF);
586

587
588
589
590
591
    for (int i = 0; i < static_cast<int>(admin.size()); i++) {
      DOFAdmin *localAdmin = admin[i];
      int n = localAdmin->getNumberOfDOFs(position);
      int n0 = localAdmin->getNumberOfPreDOFs(position);
      
592
593
      TEST_EXIT_DBG(n + n0 <= ndof)
	("n = %d, n0 = %d too large: ndof = %d\n", n, n0, ndof);
594
595
      
      for (int j = 0; j < n; j++)
596
	localAdmin->freeDofIndex(dof[n0 + j]);      
597
    }
598

599
    delete [] dof;
600
601
  }

602

603
604
605
  void Mesh::freeElement(Element* el)
  {
    freeDOFPtrs(const_cast<DegreeOfFreedom**>(el->getDOF()));
Thomas Witkowski's avatar
Thomas Witkowski committed
606
    delete el;
607
608
609
610
611
612
  }


  Element* Mesh::createNewElement(Element *parent)
  {
    FUNCNAME("Mesh::createNewElement()");
613
614

    TEST_EXIT_DBG(elementPrototype)("no element prototype\n");
615
616

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

Thomas Witkowski's avatar
Thomas Witkowski committed
618
    if (!parent && elementDataPrototype)
619
      el->setElementData(elementDataPrototype->clone()); 
Thomas Witkowski's avatar
Thomas Witkowski committed
620
    else
621
      el->setElementData(NULL); // must be done in ElementData::refineElementData()
622

623
624
625
    return el;
  }

626

627
628
  ElInfo* Mesh::createNewElInfo()
  {
629
630
    FUNCNAME("Mesh::createNewElInfo()");

631
632
    switch(dim) {
    case 1:
Thomas Witkowski's avatar
Thomas Witkowski committed
633
      return new ElInfo1d(this);
634
635
      break;
    case 2:
Thomas Witkowski's avatar
Thomas Witkowski committed
636
      return new ElInfo2d(this);
637
638
      break;
    case 3:
Thomas Witkowski's avatar
Thomas Witkowski committed
639
      return new ElInfo3d(this);
640
641
642
643
      break;
    default:
      ERROR_EXIT("invalid dim\n");
      return NULL;
644
    }
645
646
  }

647

648
649
  bool Mesh::findElInfoAtPoint(const WorldVector<double>& xy,
			       ElInfo *el_info,
650
651
			       DimVec<double>& bary,
			       const MacroElement *start_mel,
652
			       const WorldVector<double> *xy0,
653
			       double *sp)
654
655
656
  {
    static const MacroElement *mel = NULL;
    DimVec<double> lambda(dim, NO_INIT);
Thomas Witkowski's avatar
Thomas Witkowski committed
657
    ElInfo *mel_info = createNewElInfo();
658
659
660
661

    if (start_mel != NULL)
      mel = start_mel;
    else
Thomas Witkowski's avatar
Thomas Witkowski committed
662
      if (mel == NULL || mel->getElement()->getMesh() != this)
663
664
665
	mel = *(macroElements.begin());

    mel_info->setFillFlag(Mesh::FILL_COORDS);
666
    g_xy = &xy;
667
    g_xy0 = xy0;
668
    g_sp = sp;
669
670
671

    mel_info->fillMacroInfo(mel);

672
673
674
675
676
677
678
    // 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());

679
    int k;
680
681
    while ((k = mel_info->worldToCoord(xy, &lambda)) >= 0) {
      if (mel->getNeighbour(k)) {
682
683
684
685
686
	mel = mel->getNeighbour(k);
	if (macrosVisited.count(mel->getIndex()))
	  return false;

	macrosVisited.insert(mel->getIndex());
687
688
689
690
691
692
693
	mel_info->fillMacroInfo(mel);
	continue;
      }
      break;
    }

    /* now, descend in tree to find leaf element at point */
694
    bool inside = findElementAtPointRecursive(mel_info, lambda, k, el_info);
695
696
    for (int i = 0; i <= dim; i++)
      bary[i] = final_lambda[i];   
697
  
Thomas Witkowski's avatar
Thomas Witkowski committed
698
    delete mel_info;
699

700
    return inside;
701
702
703
  }

  bool Mesh::findElementAtPoint(const WorldVector<double>&  xy,
704
705
				Element **elp, 
				DimVec<double>& bary,
706
				const MacroElement *start_mel,
707
708
				const WorldVector<double> *xy0,
				double *sp)
709
  {
710
711
    ElInfo *el_info = createNewElInfo();
    int val = findElInfoAtPoint(xy, el_info, bary, start_mel, xy0, sp);
712
713
714

    *elp = el_info->getElement();

Thomas Witkowski's avatar
Thomas Witkowski committed
715
    delete el_info;
716

717
    return val;
718
719
  }

720
  bool Mesh::findElementAtPointRecursive(ElInfo *el_info,
721
					 const DimVec<double>& lambda,
722
					 int outside,
723
724
					 ElInfo* final_el_info)
  {
725
    FUNCNAME("Mesh::findElementAtPointRecursive()");
Thomas Witkowski's avatar
Thomas Witkowski committed
726

727
728
    Element *el = el_info->getElement();
    DimVec<double> c_lambda(dim, NO_INIT);
729
730
    int inside;
    int ichild, c_outside;
731
732
733
734

    if (el->isLeaf()) {
      *final_el_info = *el_info;
      if (outside < 0) {
Thomas Witkowski's avatar
Thomas Witkowski committed
735
	for (int i = 0; i <= dim; i++)
736
	  final_lambda[i] = lambda[i];
Thomas Witkowski's avatar
Thomas Witkowski committed
737
	
738
	return true;
739
740
741
742
      }  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
743
744
745
746
	  for (int i = 0; i <= dim; i++) 
	    final_lambda[i] = s * c_lambda[i] + (1.0 - s) * lambda[i];
	  
	  if (g_sp)
747
	    *(g_sp) = s;
Thomas Witkowski's avatar
Thomas Witkowski committed
748
	  
749
	  if (dim == 3) 
750
	    MSG("Outside finest level on el %d: s = %.3e\n", el->getIndex(), s);
751

752
	  return false;  /* ??? */
753
	} else {
754
	  return false;
755
	}
756
      }
757
758
    }

759
    ElInfo *c_el_info = createNewElInfo();
760

761
    if (dim == 1) {
762
763
764
765
      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
766
	  TEST_EXIT(outside == 0)("point outside domain\n");
767
768
769
770
771
772
773
774
	} 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
775
	  TEST_EXIT(outside == 0)("point outside domain\n");
776
777
778
779
780
781
782
	} else {
	  c_lambda[1] = lambda[1] - lambda[0];
	  c_lambda[0] = 2.0 * lambda[0];
	}
      }
    } /* DIM == 1 */

783
    if (dim == 2) {
784
785
786
787
      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
788
	  TEST_EXIT(outside == 0)("outside curved boundary child 0\n");	  
789
790
791
792
793
794
795
796
797
	} 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
798
	  TEST_EXIT(outside == 0)("outside curved boundary child 1\n");	  
799
800
801
802
803
804
805
806
	} else {
	  c_lambda[0] = lambda[1] - lambda[0];
	  c_lambda[1] = lambda[2];
	  c_lambda[2] = 2.0 * lambda[0];
	}
      }
    } /* DIM == 2 */

807
    if (dim == 3) {
808
809
810
811
812
813
814
815
816
817
818
819
820
      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
821
	  int c_outside2 = c_el_info2->worldToCoord(*(g_xy), &c_lambda2);
822
823

	  MSG("new_coord CHILD %d: outside=%d, lambda=(%.2f %.2f %.2f %.2f)\n",
824
	      ichild, c_outside, c_lambda[0], c_lambda[1], c_lambda[2], c_lambda[3]);
825
	  MSG("new_coord CHILD %d: outside=%d, lambda=(%.2f %.2f %.2f %.2f)\n",
826
	      1 - ichild, c_outside2, c_lambda2[0], c_lambda2[1], c_lambda2[2],
827
828
829
	      c_lambda2[3]);

	  if ((c_outside2 < 0) || (c_lambda2[c_outside2] > c_lambda[c_outside])) {
830
831
	    for (int i = 0; i <= dim; i++) 
	      c_lambda[i] = c_lambda2[i];	    
832
833
834
835
	    c_outside = c_outside2;
	    *c_el_info = *c_el_info2;
	    ichild = 1 - ichild;
	  }
Thomas Witkowski's avatar
Thomas Witkowski committed
836
	  delete c_el_info2;
837
838
839
840
841
842
	}
	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];
843
844
845
846
847
848
	  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]];
849
850
851
852
	  c_lambda[3] = 2.0 * lambda[1];
	} else {
	  c_el_info->fillElInfo(1, el_info);
	  c_lambda[0] = lambda[1] - lambda[0];
853
854
855
856
857
858
	  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]];
859
860
861
862
863
	  c_lambda[3] = 2.0 * lambda[0];
	}
      }
    }  /* DIM == 3 */

Thomas Witkowski's avatar
Thomas Witkowski committed
864
    inside = findElementAtPointRecursive(c_el_info, c_lambda, outside, final_el_info);
Thomas Witkowski's avatar
Thomas Witkowski committed
865
    delete c_el_info;
866

867
    return inside; 
868
869
  }

870
871

  bool Mesh::getDofIndexCoords(DegreeOfFreedom dof, 
872
873
874
875
876
877
			       const FiniteElemSpace* feSpace,
			       WorldVector<double>& coords)
  {
    DimVec<double>* baryCoords;
    bool found = false;
    TraverseStack stack;
878
    std::vector<DegreeOfFreedom> dofVec(feSpace->getBasisFcts()->getNumber());
Thomas Witkowski's avatar
Thomas Witkowski committed
879

880
    ElInfo *elInfo = stack.traverseFirst(this, -1, 
Thomas Witkowski's avatar
Thomas Witkowski committed
881
					 Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS);    
882
    while (elInfo) {
883
884
885
      feSpace->getBasisFcts()->getLocalIndices(elInfo->getElement(),
					       feSpace->getAdmin(),
					       dofVec);
Thomas Witkowski's avatar
Thomas Witkowski committed
886
      for (int i = 0; i < feSpace->getBasisFcts()->getNumber(); i++) {
887
	if (dofVec[i] == dof) {
888
889
890
	  baryCoords = feSpace->getBasisFcts()->getCoords(i);
	  elInfo->coordToWorld(*baryCoords, coords);
	  found = true;
Thomas Witkowski's avatar
Thomas Witkowski committed
891
	  break;	  
892
	}
Thomas Witkowski's avatar
Thomas Witkowski committed
893
894
      }
      
895
896
      if (found)
	break;
Thomas Witkowski's avatar
Thomas Witkowski committed
897
      
898
899
900
901
902
      elInfo = stack.traverseNext(elInfo);
    }

    return found;
  }
903

904
905
906

  void Mesh::getDofIndexCoords(const FiniteElemSpace* feSpace,
			       DOFVector<WorldVector<double> >& coords)
907
  {
908
909
910
    const BasisFunction* basFcts = feSpace->getBasisFcts();
    int nBasFcts = basFcts->getNumber();
    std::vector<DegreeOfFreedom> dofVec(nBasFcts);
911

912
913
914
    TraverseStack stack;
    ElInfo *elInfo = 
      stack.traverseFirst(this, -1, Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS);
915
    while (elInfo) {
916
917
918
919
      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]]);
920
      }
921
         
922
923
924
925
926
      elInfo = stack.traverseNext(elInfo);
    }

  }

927

Thomas Witkowski's avatar
Thomas Witkowski committed
928
929
930
931
  void Mesh::setDiameter(const WorldVector<double>& w) 
  { 
    diam = w; 
  }
932

933

Thomas Witkowski's avatar
Thomas Witkowski committed
934
935
936
937
  void Mesh::setDiameter(int i, double w) 
  { 
    diam[i] = w; 
  }
938

939

Thomas Witkowski's avatar
Thomas Witkowski committed
940
  void Mesh::serialize(std::ostream &out)
941
942
943
  {
    serializedDOFs.clear();

Thomas Witkowski's avatar
Thomas Witkowski committed
944
    out << name << "\n";
945

946
947