Mesh.cc 39.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 "io/MacroReader.h"
#include "io/MacroInfo.h"
#include "io/MacroWriter.h"

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

namespace AMDiS {

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

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

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

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


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

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


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

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

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

    elementPrototype->setIndex(-1);

119
120
    elementIndex = 0;
  }
121

122

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

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

138

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

175

176
    /* ====================== Copy macro elements =================== */
177
  
178
179
180
181
182
183
184
185
186
    // mapIndex[i] is the index of the MacroElement element in the vector
    // macroElements, for which holds: element->getIndex() = i    
    std::map<int, int> mapIndex;

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

    int insertCounter = 0;
187

188
189
    macroElements.clear();

190
191
192
    // Go through all MacroElements of mesh m, and create for every a new
    // MacroElement in this mesh.
    for (std::deque<MacroElement*>::const_iterator it = m.macroElements.begin();
193
	 it != m.macroElements.end(); ++it, insertCounter++) {
194

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

198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
      // Use copy operator to copy all the data to the new MacroElement.
      *el = **it;

      // Make a copy of the Element data, together with all DOFs
      el->setElement((*it)->getElement()->cloneWithDOFs());

      // Insert the new MacroElement in the vector of all MacroElements.
      macroElements.push_back(el);

      // Update the index map.
      mapIndex.insert(std::pair<int, int>(el->getIndex(), insertCounter));
    }

    // Now we have to go through all the new MacroElements, and update the neighbour
    // connections.
    insertCounter = 0;
    for (std::deque<MacroElement*>::const_iterator it = m.macroElements.begin();
	 it != m.macroElements.end();
	 ++it, insertCounter++) {
      // Go through all neighbours.
      for (int i = 0; i < dim; i++) {
	// 1. Get index of the old MacroElement for its i-th neighbour.
	// 2. Because the index in the new MacroElement is the same, search
	//    for the vector index the corresponding element is stored in.
	// 3. Get this element from macroElements, and set it as the i-th
	//    neighbour for the current element.
	macroElements[insertCounter]->
	  setNeighbour(i, macroElements[mapIndex[(*it)->getNeighbour(i)->getIndex()]]);
      }
    }

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

    /* ================== Things will be done when required ============ */
      
    TEST_EXIT(elementDataPrototype == NULL)("TODO\n");
    TEST_EXIT(m.parametric == NULL)("TODO\n");
    TEST_EXIT(periodicAssociations.size() == 0)("TODO\n");

    return *this;
  }

241

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

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

254

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

261

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

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

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

273
274
    // === Determine to all DOFs in mesh the macro elements where the DOF  ===
    // === is part of.                                                     ===
275

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

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

      elInfo = stack.traverseNext(elInfo);
    }		   

Thomas Witkowski's avatar
Thomas Witkowski committed
294

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

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

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

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


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

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
		dofsPosIndex[dofsIt->first]);
358
359
    }

360

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

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

    std::set<const DegreeOfFreedom*> allVertices;

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

372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
      if (elInfo->getElement()->isLeaf()) {
	nLeaves++;

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

      elInfo = stack.traverseNext(elInfo);
    }

    nVertices = allVertices.size();

    // === Note: Although the macro elements are removed from the mesh,   ===
    // === they are not deleted from memory. The macro elements are still ===
    // === stored in macroInfo structure. They are needed, if the mesh is ===
    // === redistributed between the ranks.                               ===
388
  }
389

390

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

    localAdmin->setMesh(this);

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

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

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

407
408
    admin.push_back(localAdmin);

409
    nDofEl = 0;
410

411
412
413
    localAdmin->setNumberOfPreDOFs(VERTEX,nDof[VERTEX]);
    nDof[VERTEX] += localAdmin->getNumberOfDofs(VERTEX);
    nDofEl += getGeo(VERTEX) * nDof[VERTEX];
414

415
    if (dim > 1) {
416
417
418
      localAdmin->setNumberOfPreDOFs(EDGE,nDof[EDGE]);
      nDof[EDGE] += localAdmin->getNumberOfDofs(EDGE);
      nDofEl += getGeo(EDGE) * nDof[EDGE];
419
420
    }

421
422
423
    localAdmin->setNumberOfPreDOFs(CENTER,nDof[CENTER]);
    nDof[CENTER]  += localAdmin->getNumberOfDofs(CENTER);
    nDofEl += nDof[CENTER];
424

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

427
428
    node[VERTEX] = 0;
    nNodeEl = getGeo(VERTEX);
429

430
431
    if (dim > 1) {
      node[EDGE] = nNodeEl;
432
      if (nDof[EDGE] > 0) 
433
	nNodeEl += getGeo(EDGE);
434
435
    }

436
    if (dim == 3) {
437
438
439
      localAdmin->setNumberOfPreDOFs(FACE,nDof[FACE]);
      nDof[FACE] += localAdmin->getNumberOfDofs(FACE);
      nDofEl += getGeo(FACE) * nDof[FACE];
440
      node[FACE] = nNodeEl;
441
      if (nDof[FACE] > 0) 
442
	nNodeEl += getGeo(FACE);
443
444
    }

445
    node[CENTER] = nNodeEl;
446
    if (nDof[CENTER] > 0)
447
      nNodeEl += 1;
448
449
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
450

451
452
  void Mesh::dofCompress()
  {
453
    FUNCNAME("Mesh::dofCompress()");
454

Thomas Witkowski's avatar
Thomas Witkowski committed
455
456
    for (unsigned int iadmin = 0; iadmin < admin.size(); iadmin++) {
      DOFAdmin* compressAdmin = admin[iadmin];
457
458

      TEST_EXIT_DBG(compressAdmin)("no admin[%d] in mesh\n", iadmin);
459
      
Thomas Witkowski's avatar
Thomas Witkowski committed
460
      int size = compressAdmin->getSize();
Thomas Witkowski's avatar
Thomas Witkowski committed
461
      if (size < 1 || 
462
	  compressAdmin->getUsedDofs() < 1 || 
Thomas Witkowski's avatar
Thomas Witkowski committed
463
	  compressAdmin->getHoleCount() < 1)    
464
	continue;
Thomas Witkowski's avatar
Thomas Witkowski committed
465

466
467
      std::vector<int> newDofIndex(size);     
      compressAdmin->compress(newDofIndex);
Thomas Witkowski's avatar
Thomas Witkowski committed
468
469
470
471

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

      elInfo = stack.traverseFirst(this, -1, fill_flag);
      while (elInfo) {
481
	elInfo->getElement()->newDofFct2(compressAdmin);
482
483
	elInfo = stack.traverseNext(elInfo);
      }
484
    }       
485
486
487
  }


488
  DegreeOfFreedom *Mesh::getDof(GeoIndex position)
489
  {
490
    FUNCNAME("Mesh::getDof()");
491

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

495
    int ndof = getNumberOfDofs(position);
496
    if (ndof <= 0) 
497
      return NULL;
498

499
    DegreeOfFreedom *dof = new DegreeOfFreedom[ndof];
500

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


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

    if (nNodeEl <= 0)
523
      return NULL;
524

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

529
    return ptrs;
530
531
  }

532

533
  void Mesh::freeDofPtrs(DegreeOfFreedom **ptrs)
534
  {
535
    FUNCNAME("Mesh::freeDofPtrs()");
536

537
    TEST_EXIT_DBG(ptrs)("ptrs is NULL!\n");
538
539
540
541

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


546
  const DOFAdmin *Mesh::createDOFAdmin(std::string lname, DimVec<int> lnDof)
547
  {
548
    FUNCNAME("Mesh::createDOFAdmin()");
549

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

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

    addDOFAdmin(localAdmin);

557
    return localAdmin;
558
559
560
561
562
563
564
  }


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

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

574
    return localAdmin;
575
576
  }

577

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

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

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

595
596
      return;
    }
597

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

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

613
    delete [] dof;
614
    dof = NULL;
615
616
  }

617

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


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

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

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

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

638
639
640
    return el;
  }

641

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

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

662

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

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

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

    mel_info->fillMacroInfo(mel);

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

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

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

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

715
    return inside;
716
717
718
  }

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

    *elp = el_info->getElement();

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

732
    return val;
733
734
  }

735
  bool Mesh::findElementAtPointRecursive(ElInfo *el_info,
736
					 const DimVec<double>& lambda,
737
					 int outside,
738
739
					 ElInfo* final_el_info)
  {
740
    FUNCNAME("Mesh::findElementAtPointRecursive()");
Thomas Witkowski's avatar
Thomas Witkowski committed
741

742
743
    Element *el = el_info->getElement();
    DimVec<double> c_lambda(dim, NO_INIT);
744
745
    int inside;
    int ichild, c_outside;
746
747
748
749

    if (el->isLeaf()) {
      *final_el_info = *el_info;
      if (outside < 0) {
Thomas Witkowski's avatar
Thomas Witkowski committed
750
	for (int i = 0; i <= dim; i++)
751
	  final_lambda[i] = lambda[i];
Thomas Witkowski's avatar
Thomas Witkowski committed
752
	
753
	return true;
754
755
756
757
      }  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
758
759
760
761
	  for (int i = 0; i <= dim; i++) 
	    final_lambda[i] = s * c_lambda[i] + (1.0 - s) * lambda[i];
	  
	  if (g_sp)
762
	    *(g_sp) = s;
Thomas Witkowski's avatar
Thomas Witkowski committed
763
	  
764
	  if (dim == 3) 
765
	    MSG("Outside finest level on el %d: s = %.3e\n", el->getIndex(), s);
766

767
	  return false;  /* ??? */
768
	} else {
769
	  return false;
770
	}
771
      }
772
773
    }

774
    ElInfo *c_el_info = createNewElInfo();
775

776
    if (dim == 1) {
777
778
779
780
      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
781
	  TEST_EXIT(outside == 0)("point outside domain\n");
782
783
784
785
786
787
788
789
	} 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
790
	  TEST_EXIT(outside == 0)("point outside domain\n");
791
792
793
794
795
796
797
	} else {
	  c_lambda[1] = lambda[1] - lambda[0];
	  c_lambda[0] = 2.0 * lambda[0];
	}
      }
    } /* DIM == 1 */

798
    if (dim == 2) {
799
800
801
802
      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
803
	  TEST_EXIT(outside == 0)("outside curved boundary child 0\n");	  
804
805
806
807
808
809
810
811
812
	} 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
813
	  TEST_EXIT(outside == 0)("outside curved boundary child 1\n");	  
814
815
816
817
818
819
820
821
	} else {
	  c_lambda[0] = lambda[1] - lambda[0];
	  c_lambda[1] = lambda[2];
	  c_lambda[2] = 2.0 * lambda[0];
	}
      }
    } /* DIM == 2 */

822
    if (dim == 3) {
823
824
825
826
827
828
829
830
831
832
833
834
835
      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
836
	  int c_outside2 = c_el_info2->worldToCoord(*(g_xy), &c_lambda2);
837
838

	  MSG("new_coord CHILD %d: outside=%d, lambda=(%.2f %.2f %.2f %.2f)\n",
839
	      ichild, c_outside, c_lambda[0], c_lambda[1], c_lambda[2], c_lambda[3]);
840
	  MSG("new_coord CHILD %d: outside=%d, lambda=(%.2f %.2f %.2f %.2f)\n",
841
	      1 - ichild, c_outside2, c_lambda2[0], c_lambda2[1], c_lambda2[2],
842
843
844
	      c_lambda2[3]);

	  if ((c_outside2 < 0) || (c_lambda2[c_outside2] > c_lambda[c_outside])) {
845
846
	    for (int i = 0; i <= dim; i++) 
	      c_lambda[i] = c_lambda2[i];	    
847
848
849
850
	    c_outside = c_outside2;
	    *c_el_info = *c_el_info2;
	    ichild = 1 - ichild;
	  }
Thomas Witkowski's avatar
Thomas Witkowski committed
851
	  delete c_el_info2;
852
853
854
855
856
857
	}
	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];
858
859
860
861
862
863
	  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]];
864
865
866
867
	  c_lambda[3] = 2.0 * lambda[1];
	} else {
	  c_el_info->fillElInfo(1, el_info);
	  c_lambda[0] = lambda[1] - lambda[0];
868
869
870
871
872
873
	  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]];
874
875
876
877
878
	  c_lambda[3] = 2.0 * lambda[0];
	}
      }
    }  /* DIM == 3 */

Thomas Witkowski's avatar
Thomas Witkowski committed
879
    inside = findElementAtPointRecursive(c_el_info, c_lambda, outside, final_el_info);
Thomas Witkowski's avatar
Thomas Witkowski committed
880
    delete c_el_info;
881

882
    return inside; 
883
884
  }

885
886

  bool Mesh::getDofIndexCoords(DegreeOfFreedom dof, 
887
888
889
890
891
892
			       const FiniteElemSpace* feSpace,
			       WorldVector<double>& coords)
  {
    DimVec<double>* baryCoords;
    bool found = false;
    TraverseStack stack;
893
    std::vector<DegreeOfFreedom> dofVec(feSpace->getBasisFcts()->getNumber());
Thomas Witkowski's avatar
Thomas Witkowski committed
894

895
    ElInfo *elInfo = stack.traverseFirst(this, -1, 
Thomas Witkowski's avatar
Thomas Witkowski committed
896
					 Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS);    
897
    while (elInfo) {
898
899
900
      feSpace->getBasisFcts()->getLocalIndices(elInfo->getElement(),
					       feSpace->getAdmin(),
					       dofVec);
Thomas Witkowski's avatar
Thomas Witkowski committed
901
      for (int i = 0; i < feSpace->getBasisFcts()->getNumber(); i++) {
902
	if (dofVec[i] == dof) {
903
904
905
	  baryCoords = feSpace->getBasisFcts()->getCoords(i);
	  elInfo->coordToWorld(*baryCoords, coords);
	  found = true;
Thomas Witkowski's avatar
Thomas Witkowski committed
906
	  break;	  
907
	}
Thomas Witkowski's avatar
Thomas Witkowski committed
908
909
      }
      
910
911
      if (found)
	break;
Thomas Witkowski's avatar
Thomas Witkowski committed
912
      
913
914
915
916
917
      elInfo = stack.traverseNext(elInfo);
    }

    return found;
  }
918