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


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

#include "time.h"

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

23
24
25
26
#include "AdaptStationary.h"
#include "AdaptInstationary.h"
#include "FiniteElemSpace.h"
#include "ElementData.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
27
#include "ElementDofIterator.h"
28
29
30
31
32
33
34
35
36
37
38
#include "MacroElement.h"
#include "Mesh.h"
#include "Traverse.h"
#include "Parameters.h"
#include "FixVec.h"
#include "DOFVector.h"
#include "CoarseningManager.h"
#include "DOFIterator.h"
#include "VertexVector.h"
#include "PeriodicMap.h"
#include "Projection.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
39
#include "ElInfoStack.h"
40
#include "Serializer.h"
41
#include "Lagrange.h"
42
43
44
45
46
47
48
49
50

namespace AMDiS {

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

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

51
52
53
54
55
56
57
58
59
60
61
62
63
64
  const Flag Mesh::FILL_NOTHING    = 0X00L;
  const Flag Mesh::FILL_COORDS     = 0X01L;
  const Flag Mesh::FILL_BOUND      = 0X02L;
  const Flag Mesh::FILL_NEIGH      = 0X04L;
  const Flag Mesh::FILL_OPP_COORDS = 0X08L;
  const Flag Mesh::FILL_ORIENTATION= 0X10L;
  const Flag Mesh::FILL_DET        = 0X20L;
  const Flag Mesh::FILL_GRD_LAMBDA = 0X40L;
  const Flag Mesh::FILL_ADD_ALL    = 0X80L;


  const Flag Mesh::FILL_ANY_1D = (0X01L|0X02L|0X04L|0X08L|0x20L|0X40L|0X80L);
  const Flag Mesh::FILL_ANY_2D = (0X01L|0X02L|0X04L|0X08L|0x20L|0X40L|0X80L);
  const Flag Mesh::FILL_ANY_3D = (0X01L|0X02L|0X04L|0X08L|0X10L|0x20L|0X40L|0X80L);
65
66
67
68
69
70
71
72
73
74
75

  //**************************************************************************
  //  flags for Mesh traversal                                                
  //**************************************************************************

  const Flag Mesh::CALL_EVERY_EL_PREORDER  = 0X0100L;
  const Flag Mesh::CALL_EVERY_EL_INORDER   = 0X0200L;
  const Flag Mesh::CALL_EVERY_EL_POSTORDER = 0X0400L;
  const Flag Mesh::CALL_LEAF_EL            = 0X0800L;
  const Flag Mesh::CALL_LEAF_EL_LEVEL      = 0X1000L;
  const Flag Mesh::CALL_EL_LEVEL           = 0X2000L;
76
77
  const Flag Mesh::CALL_MG_LEVEL           = 0X4000L;  
  const Flag Mesh::CALL_REVERSE_MODE       = 0X8000L;
78
79


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

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

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

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

    elementPrototype->setIndex(-1);

125
126
    elementIndex = 0;
  }
127

128

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

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

144

145
  Mesh& Mesh::operator=(const Mesh& m)
146
  {
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
    FUNCNAME("Mesh::operator=()");

    if (this == &m)
      return *this;

    TEST_EXIT(dim == m.dim)("operator= works only on meshes with equal dim!\n");

    name = m.name;
    nVertices = m.nVertices;
    nEdges = m.nEdges;
    nLeaves = m.nLeaves;
    nElements = m.nElements;
    nFaces = m.nFaces;
    maxEdgeNeigh = m.maxEdgeNeigh;
    diam = m.diam;
    parametric = NULL;

    preserveCoarseDOFs = m.preserveCoarseDOFs;
165
166
    nDofEl = m.nDofEl;
    nDof = m.nDof;
167
168
169
170
171
172
173
174
175
    nNodeEl = m.nNodeEl;
    node = m.node;
    elementIndex = m.elementIndex;
    initialized = m.initialized;
    final_lambda = m.final_lambda;
    
    /* ====================== Create new DOFAdmins ================== */
    admin.resize(m.admin.size());
    for (int i = 0; i < static_cast<int>(admin.size()); i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
176
      admin[i] = new DOFAdmin(this);
177
      *(admin[i]) = *(m.admin[i]);
178
179
      admin[i]->setMesh(this);
    }
180

181

182
    /* ====================== Copy macro elements =================== */
183
  
184
185
186
187
188
189
190
191
192
    // mapIndex[i] is the index of the MacroElement element in the vector
    // macroElements, for which holds: element->getIndex() = i    
    std::map<int, int> mapIndex;

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

    int insertCounter = 0;
193

194
195
    macroElements.clear();

196
197
198
    // Go through all MacroElements of mesh m, and create for every a new
    // MacroElement in this mesh.
    for (std::deque<MacroElement*>::const_iterator it = m.macroElements.begin();
199
	 it != m.macroElements.end(); ++it, insertCounter++) {
200

201
      // Create new MacroElement.
Thomas Witkowski's avatar
Thomas Witkowski committed
202
      MacroElement *el = new MacroElement(dim);
203

204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
      // 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;
  }

247

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

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

260

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

267

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

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

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

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

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

      elInfo = stack.traverseNext(elInfo);
    }		   

Thomas Witkowski's avatar
Thomas Witkowski committed
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())
314
	newMacroElements.push_back(*elIter);     
315
    }
316

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

Thomas Witkowski's avatar
Thomas Witkowski committed
321

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

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

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

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

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

      // We will delete at least some element DOFs in the next but will keep the
      // element object still in memory. Therefore the DOF pointer must be set to be
      // invalide.
      (*macroIt)->getElement()->setDofValid(false);
352
    }     
353
354


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

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

377

378
    // === Update number of elements, vertices, etc. ===
379

380
381
382
383
384
385
386
387
    nLeaves = 0;
    nElements = 0;

    std::set<const DegreeOfFreedom*> allVertices;

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

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

407

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

    localAdmin->setMesh(this);

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

    admin.push_back(localAdmin);

420
    nDofEl = 0;
421

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

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

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

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

438
439
    node[VERTEX] = 0;
    nNodeEl = getGeo(VERTEX);
440

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
461

462
463
  void Mesh::dofCompress()
  {
464
    FUNCNAME("Mesh::dofCompress()");
465

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

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

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

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

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


499
  DegreeOfFreedom *Mesh::getDof(GeoIndex position)
500
  {
501
    FUNCNAME("Mesh::getDof()");
502

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

506
    int ndof = getNumberOfDofs(position);
507
    if (ndof <= 0) 
508
      return NULL;
509

510
    DegreeOfFreedom *dof = new DegreeOfFreedom[ndof];
511

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


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

    if (nNodeEl <= 0)
534
      return NULL;
535

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

540
    return ptrs;
541
542
  }

543

544
  void Mesh::freeDofPtrs(DegreeOfFreedom **ptrs)
545
  {
546
    FUNCNAME("Mesh::freeDofPtrs()");
547

548
    TEST_EXIT_DBG(ptrs)("ptrs is NULL!\n");
549
550
551
552

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


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

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

    addDOFAdmin(localAdmin);

568
    return localAdmin;
569
570
571
572
573
574
575
  }


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

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

585
    return localAdmin;
586
587
  }

588

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

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

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

606
607
      return;
    }
608

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

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

624
    delete [] dof;
625
    dof = NULL;
626
627
  }

628

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


  Element* Mesh::createNewElement(Element *parent)
  {
    FUNCNAME("Mesh::createNewElement()");
639
640

    TEST_EXIT_DBG(elementPrototype)("no element prototype\n");
641
642

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

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

649
650
651
    return el;
  }

652

653
654
  ElInfo* Mesh::createNewElInfo()
  {
655
656
    FUNCNAME("Mesh::createNewElInfo()");

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

673

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

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

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

    mel_info->fillMacroInfo(mel);

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

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

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

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

726
    return inside;
727
728
  }

729

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

    *elp = el_info->getElement();

Thomas Witkowski's avatar
Thomas Witkowski committed
742
    delete el_info;
743

744
    return val;
745
746
  }

747

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

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

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

780
	  return false;  /* ??? */
781
	} else {
782
	  return false;
783
	}
784
      }
785
786
    }

787
    ElInfo *c_el_info = createNewElInfo();
788

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

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

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

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

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

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

895
    return inside; 
896
897
  }

898
899

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

908
    ElInfo *elInfo = stack.traverseFirst(this, -1, 
Thomas Witkowski's avatar
Thomas Witkowski committed
909
					 Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS);    
910
    while (elInfo) {