Mesh.cc 38.7 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
14
15
16
17
18
19
20
21
22
23
24
#include "MacroElement.h"
#include "MacroReader.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 "MacroWriter.h"
#include "PeriodicMap.h"
#include "Projection.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
25
#include "ElInfoStack.h"
26
#include "Serializer.h"
27
#include "Lagrange.h"
28
29
30
31
32
33
34
35
36

namespace AMDiS {

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

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

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

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


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

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


Thomas Witkowski's avatar
Thomas Witkowski committed
76
  Mesh::Mesh(std::string aName, int dimension) 
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
    : 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),
93
      macroFileInfo(NULL),
94
      changeIndex(0),
95
96
97
      final_lambda(dimension, DEFAULT_VALUE, 0.0)
  {

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

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

    elementPrototype->setIndex(-1);

117
118
    elementIndex = 0;
  }
119

120

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

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

136

137
  Mesh& Mesh::operator=(const Mesh& m)
138
  {
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
    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
169
      admin[i] = new DOFAdmin(this);
170
      *(admin[i]) = *(m.admin[i]);
171
172
      admin[i]->setMesh(this);
    }
173

174

175
    /* ====================== Copy macro elements =================== */
176
  
177
178
179
180
181
182
183
184
185
    // 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;
186

187
188
    macroElements.clear();

189
190
191
    // 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();
192
	 it != m.macroElements.end(); ++it, insertCounter++) {
193

194
      // Create new MacroElement.
Thomas Witkowski's avatar
Thomas Witkowski committed
195
      MacroElement *el = new MacroElement(dim);
196

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

240

241
242
243
244
245
246
247
248
249
250
251
252
  void Mesh::updateNumberOfLeaves()
  {
    nLeaves = 0;

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

253

254
255
256
257
258
259
  void Mesh::addMacroElement(MacroElement* me) 
  {
    macroElements.push_back(me); 
    me->setIndex(macroElements.size());
  }

260

261
  void Mesh::removeMacroElements(std::set<MacroElement*>& macros,
Thomas Witkowski's avatar
Thomas Witkowski committed
262
				 const FiniteElemSpace *feSpace) 
263
264
265
  {
    FUNCNAME("Mesh::removeMacroElement()");

Thomas Witkowski's avatar
Thomas Witkowski committed
266
267
268
269
270
271
272
    typedef std::map<const DegreeOfFreedom*, std::set<MacroElement*> > DofElMap;
    typedef std::map<const DegreeOfFreedom*, GeoIndex> DofPosMap;

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

    ElementDofIterator elDofIter(feSpace);
273

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


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

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

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

311
312
313
314
315
    // === 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
316

317
318
319
320
      // 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
321
322
	  for (int j = 0; j <= dim; j++)
	    if ((*macroIt)->getNeighbour(i)->getNeighbour(j) == *macroIt)
323
324
325
326
327
	      (*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--;
328
	}
329
330
      }

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


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

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

362
363
364
365
366
367
368

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

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

369
    nVertices = nRemainDofs;
370
  }
371

372

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

    localAdmin->setMesh(this);

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

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

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

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

    nDOFEl = 0;

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

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

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

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

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

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
432

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

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

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

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

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

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

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


  DegreeOfFreedom *Mesh::getDOF(GeoIndex position)
  {
475
    FUNCNAME("Mesh::getDOF()");
476

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

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

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

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


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

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

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

514
    return ptrs;
515
516
  }

517

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

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

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


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

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

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

    addDOFAdmin(localAdmin);

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


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

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

559
    return localAdmin;
560
561
  }

562

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

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

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

580
581
      return;
    }
582

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

586
587
588
589
590
    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);
      
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
	localAdmin->freeDofIndex(dof[n0 + j]);      
596
    }
597

598
    delete [] dof;
599
600
  }

601

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


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

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

    Element *el = parent ? parent->clone() : elementPrototype->clone();
  
Thomas Witkowski's avatar
Thomas Witkowski committed
617
    if (!parent && elementDataPrototype)
618
      el->setElementData(elementDataPrototype->clone()); 
Thomas Witkowski's avatar
Thomas Witkowski committed
619
    else
620
      el->setElementData(NULL); // must be done in ElementData::refineElementData()
Thomas Witkowski's avatar
Thomas Witkowski committed
621
    
622
623
624
    return el;
  }

625

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

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

646

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

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

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

    mel_info->fillMacroInfo(mel);

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

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

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

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

699
    return inside;
700
701
702
  }

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

    *elp = el_info->getElement();

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

716
    return val;
717
718
  }

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

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

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

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

758
    ElInfo *c_el_info = createNewElInfo();
759

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

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

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

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

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

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

866
    return inside; 
867
868
  }

869
870

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

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

    return found;
  }
902

903
904
905

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

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

  }

926

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

932

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

938

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

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

945
946
947
948
949
950
951
    SerUtil::serialize(out, dim);
    SerUtil::serialize(out, nVertices);
    SerUtil::serialize(out, nEdges);
    SerUtil::serialize(out, nLeaves);
    SerUtil::serialize(out, nElements);
    SerUtil::serialize(out, nFaces);
    SerUtil::serialize(out, maxEdgeNeigh);
952
953
954

    diam.serialize(out);

955
956
    SerUtil::serialize(out, preserveCoarseDOFs);
    SerUtil::serialize(out, nDOFEl);
957
958
959

    nDOF.serialize(out);

960
    SerUtil::serialize(out, nNodeEl);