Mesh.cc 34.2 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
28
29
30
31
32
33
34
35

namespace AMDiS {

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

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

36
37
38
39
40
41
42
43
44
45
46
47
48
49
  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);
50
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;
  const Flag Mesh::CALL_MG_LEVEL           = 0X4000L ; // used in mg methods 
62
63
64


  // const Flag Mesh::USE_PARAMETRIC          = 0X8000L ; // used in mg methods 
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

  Mesh::~Mesh()
121
  {
122
    Element::deletedDOFs.clear();
123
124

    for (std::deque<MacroElement*>::const_iterator it = macroElements.begin();
125
	 it != macroElements.end(); ++it) {
126
      (*it)->getElement()->deleteElementDOFs();
Thomas Witkowski's avatar
Thomas Witkowski committed
127
      delete *it;
128
129
    }    

130
    Element::deletedDOFs.clear();
131

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

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

180

181
    /* ====================== Copy macro elements =================== */
182
  
183
184
185
186
187
188
189
190
191
    // 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;
192

193
194
    macroElements.clear();

195
196
197
    // 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();
198
	 it != m.macroElements.end(); ++it, insertCounter++) {
199

200
      // Create new MacroElement.
Thomas Witkowski's avatar
Thomas Witkowski committed
201
      MacroElement *el = new MacroElement(dim);
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
241
242
243
244
245
      // 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;
  }

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

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

258
259
260
261
262
263
  void Mesh::addMacroElement(MacroElement* me) 
  {
    macroElements.push_back(me); 
    me->setIndex(macroElements.size());
  }

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

Thomas Witkowski's avatar
Thomas Witkowski committed
269
270
271
272
273
274
275
    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);
276

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


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

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

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

314
315
316
317
318
    // === 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
319

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

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


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

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

365
366
367
368
369
370
371

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

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

372
    nVertices = nRemainDofs;
373
  }
374
375
376
377
378
379
380

  void Mesh::addDOFAdmin(DOFAdmin *localAdmin)
  {    
    FUNCNAME("Mesh::addDOFAdmin()");

    localAdmin->setMesh(this);

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

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

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

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

    nDOFEl = 0;

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

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

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

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

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

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
434

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

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

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

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

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

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

470
      newDOF.resize(0);
471
    }       
472
473
474
475
476
  }


  DegreeOfFreedom *Mesh::getDOF(GeoIndex position)
  {
477
    FUNCNAME("Mesh::getDOF()");
478

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

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

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

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


  DegreeOfFreedom **Mesh::createDOFPtrs()
  {
507
    FUNCNAME("Mesh::createDOFPtrs()");
508
509

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

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

516
    return ptrs;
517
518
519
520
  }

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

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

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


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

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

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

    addDOFAdmin(localAdmin);

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


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

551
552
553
554
555
556
    for (int i = 0; i < static_cast<int>(admin.size()); i++) {
      if (admin[i]->getNumberOfDOFs(VERTEX)) {
	if (!localAdmin)  
	  localAdmin = admin[i];
	else if (admin[i]->getSize() < localAdmin->getSize())
	  localAdmin = admin[i];
557
      }
558
559
    }

560
    return localAdmin;
561
562
563
564
  }

  void Mesh::freeDOF(DegreeOfFreedom* dof, GeoIndex position)
  {
565
    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
591
592
593
594
595
    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);
      
      TEST_EXIT_DBG(n + n0 <= ndof)("n=%d, n0=%d too large: ndof=%d\n", n, n0, ndof);
      
      for (int j = 0; j < n; j++)
	localAdmin->freeDOFIndex(dof[n0 + j]);
    }
596

597
    delete [] dof;
598
599
600
601
602
  }

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


  Element* Mesh::createNewElement(Element *parent)
  {
    FUNCNAME("Mesh::createNewElement()");
610
611

    TEST_EXIT_DBG(elementPrototype)("no element prototype\n");
612
613
614

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

623

624
625
  ElInfo* Mesh::createNewElInfo()
  {
626
627
    FUNCNAME("Mesh::createNewElInfo()");

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

644

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

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

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

    mel_info->fillMacroInfo(mel);

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

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

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

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

697
    return inside;
698
699
700
  }

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

    *elp = el_info->getElement();

Thomas Witkowski's avatar
Thomas Witkowski committed
712
    delete el_info;
713

714
    return val;
715
716
  }

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

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

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

749
	  return false;  /* ??? */
750
	} else {
751
	  return false;
752
	}
753
      }
754
755
    }

756
    ElInfo *c_el_info = createNewElInfo();
757

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

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

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

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

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

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

864
    return inside; 
865
866
  }

867
868

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

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

    return found;
  }
900

901
902
903

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

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

  }

924

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
935
  void Mesh::serialize(std::ostream &out)
936
937
938
  {
    serializedDOFs.clear();

Thomas Witkowski's avatar
Thomas Witkowski committed
939
    out << name << "\n";
940

941
942
943
944
945
946
947
    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);
948
949
950

    diam.serialize(out);

951
952
    SerUtil::serialize(out, preserveCoarseDOFs);
    SerUtil::serialize(out, nDOFEl);
953
954
955

    nDOF.serialize(out);

956
    SerUtil::serialize(out, nNodeEl);
957
958
959

    node.serialize(out);

960
961
962

    // === Write admins. ===

963
    int size = static_cast<int>(admin.size());
964
    SerUtil::serialize(out, size);
965
    for (int i = 0; i < size; i++)
966
967
      admin[i]->serialize(out);

968
969
970

    // === Write macroElements. ===

971
    size = static_cast<int>(macroElements.size());
972
    SerUtil::serialize(out, size);
973
    for (int i = 0; i < size; i++)
974
975
      macroElements[i]->serialize(out);

976
977
    SerUtil::serialize(out, elementIndex);
    SerUtil::serialize(out, initialized);
978

979
980
981
982
983
984
985
986
987