Mesh.cc 32.6 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
36
37
38
39
40
41
42
43
44
45
46

namespace AMDiS {

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

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

  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;


47
48
49
  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
62
63
64
65
66
67

  //**************************************************************************
  //  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 


  // const Flag Mesh::USE_PARAMETRIC          = 0X8000L ; // used in mg methods 

  DOFAdmin* Mesh::compressAdmin = NULL;
  Mesh* Mesh::traversePtr = NULL;
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
87
88
89
90
91
92
93
94
    : 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),
95
      macroFileInfo(NULL),
96
97
98
      final_lambda(dimension, DEFAULT_VALUE, 0.0)
  {

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

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

    elementPrototype->setIndex(-1);

118
119
    elementIndex = 0;
  }
120
121

  Mesh::~Mesh()
122
  {
123
    Element::deletedDOFs.clear();
124
125
126
127
128

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

132
    Element::deletedDOFs.clear();
133

134
    if (macroFileInfo != NULL) {
Thomas Witkowski's avatar
Thomas Witkowski committed
135
      delete macroFileInfo;
136
    }
137
    if (elementPrototype) {
Thomas Witkowski's avatar
Thomas Witkowski committed
138
      delete elementPrototype;
139
140
    }
    if (elementDataPrototype) {
Thomas Witkowski's avatar
Thomas Witkowski committed
141
      delete elementDataPrototype;
142
143
144
    }
    
    for (int i = 0; i < static_cast<int>(admin.size()); i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
145
      delete admin[i];
146
    }
147
  }
148
149

  Mesh& Mesh::operator=(const Mesh& m)
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
175
176
177
178
179
180
    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
181
      admin[i] = new DOFAdmin(this);
182
      *(admin[i]) = *(m.admin[i]);
183
184
      admin[i]->setMesh(this);
    }
185

186

187
    /* ====================== Copy macro elements =================== */
188
  
189
190
191
192
193
194
195
196
197
    // 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;
198

199
200
    macroElements.clear();

201
202
203
204
205
    // 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();
	 it != m.macroElements.end();
	 ++it, insertCounter++) {
206

207
      // Create new MacroElement.
Thomas Witkowski's avatar
Thomas Witkowski committed
208
      MacroElement *el = new MacroElement(dim);
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
247
248
249
250
251
252
      // 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;
  }

253
254
255
256
257
258
259
260
261
262
263
264
  void Mesh::updateNumberOfLeaves()
  {
    nLeaves = 0;

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

265
266
267
268
269
270
  void Mesh::addMacroElement(MacroElement* me) 
  {
    macroElements.push_back(me); 
    me->setIndex(macroElements.size());
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
271
272
  void Mesh::removeMacroElements(std::vector<MacroElement*>& macros,
				 const FiniteElemSpace *feSpace) 
273
274
275
  {
    FUNCNAME("Mesh::removeMacroElement()");

Thomas Witkowski's avatar
Thomas Witkowski committed
276
277
278
279
280
281
282
283
    typedef std::map<const DegreeOfFreedom*, std::set<MacroElement*> > DofElMap;
    typedef std::map<const DegreeOfFreedom*, GeoIndex> DofPosMap;

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

    ElementDofIterator elDofIter(feSpace);
284

285
286
    // Map that stores for each dof pointer (which may have a list of dofs)
    // all macro element indices that own the dof.
Thomas Witkowski's avatar
Thomas Witkowski committed
287
288
    DofElMap dofsOwner;
    DofPosMap dofsPosIndex;
289
290
291
292
293
    
    // Determine all dof owner macro elements.
    for (std::deque<MacroElement*>::iterator macroIt = macroElements.begin();
	 macroIt != macroElements.end();
	 ++macroIt) {
Thomas Witkowski's avatar
Thomas Witkowski committed
294
295
296
297
298
      elDofIter.reset((*macroIt)->getElement());
      do {
	dofsOwner[elDofIter.getDofPtr()].insert(*macroIt);
	dofsPosIndex[elDofIter.getDofPtr()] = elDofIter.getPosIndex();
      } while (elDofIter.nextStrict());
299
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
300

301
302
303
304
    // Remove all the given macro elements.
    for (std::vector<MacroElement*>::iterator macroIt = macros.begin();
	 macroIt != macros.end();
	 ++macroIt) {
Thomas Witkowski's avatar
Thomas Witkowski committed
305
306
307
308
309
310
311

      std::deque<MacroElement*>::iterator mEl = 
	find(macroElements.begin(), macroElements.end(), *macroIt);
      TEST_EXIT(mEl != macroElements.end())
	("Cannot find MacroElement that should be removed!\n");
      macroElements.erase(mEl);

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

      nLeaves--;
      nElements--;

      // Remove this macro element from the dof owner list.
Thomas Witkowski's avatar
Thomas Witkowski committed
330
331
      for (DofElMap::iterator dofsIt = dofsOwner.begin();
	   dofsIt != dofsOwner.end(); ++dofsIt) {
332
	std::set<MacroElement*>::iterator mIt = dofsIt->second.find(*macroIt);
Thomas Witkowski's avatar
Thomas Witkowski committed
333
	if (mIt != dofsIt->second.end())
334
335
336
337
338
339
340
341
342
343
	  dofsIt->second.erase(mIt);
      }

      // And remove the macro element from memory
      delete *macroIt;
    }

    int nRemainDofs = 0;
    // Check now all the dofs, that have no owner anymore and therefore have to
    // be removed.
Thomas Witkowski's avatar
Thomas Witkowski committed
344
345
346
    for (DofElMap::iterator dofsIt = dofsOwner.begin(); 
	 dofsIt != dofsOwner.end(); ++dofsIt) {    
      if (dofsIt->second.size() == 0)
347
	freeDOF(const_cast<DegreeOfFreedom*>(dofsIt->first), 
Thomas Witkowski's avatar
Thomas Witkowski committed
348
349
		dofsPosIndex[dofsIt->first]);
      else
350
	nRemainDofs++;
351
352
    }

353
    nVertices = nRemainDofs;
354
  }
355

356
  int Mesh::traverse(int level, Flag flag, int (*el_fct)(ElInfo*))
357
358
  {
    FUNCNAME("Mesh::traverse()");
Thomas Witkowski's avatar
Thomas Witkowski committed
359
360
361

    ElInfoStack elInfoStack(this);
    ElInfo* elinfo = elInfoStack.getNextElement();
362
363
364
365
366
    Traverse tinfo(this, flag, level, el_fct);
  
    elinfo->setFillFlag(flag);
  
    if (flag.isSet(Mesh::CALL_LEAF_EL_LEVEL) || 
367
	flag.isSet(Mesh::CALL_EL_LEVEL) || 
368
369
370
	flag.isSet(Mesh::CALL_MG_LEVEL)) {
      TEST(level >= 0)("invalid level: %d\n", level);
    }
371
372
373
374

    int sum = 0;  
    for (std::deque<MacroElement*>::iterator mel = macroElements.begin(); 
	 mel != macroElements.end(); mel++) {
375
      elinfo->fillMacroInfo(*mel);
Thomas Witkowski's avatar
Thomas Witkowski committed
376
      sum += tinfo.recursive(&elInfoStack);
377
378
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
379
    elInfoStack.getBackElement();
380
381
382
383
384
385
386
387
388
389
    
    return (flag.isSet(Mesh::FILL_ADD_ALL)) ? sum : 0;
  }

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

    localAdmin->setMesh(this);

390
    std::vector<DOFAdmin*>::iterator dai = 
391
      std::find(admin.begin(), admin.end(), localAdmin);
392

Thomas Witkowski's avatar
Thomas Witkowski committed
393
394
395
    TEST_EXIT(dai == admin.end())
      ("admin %s is already associated to mesh %s\n",
       localAdmin->getName().c_str(), this->getName().c_str());
396

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

400
401
402
403
404
    admin.push_back(localAdmin);

    nDOFEl = 0;

    localAdmin->setNumberOfPreDOFs(VERTEX,nDOF[VERTEX]);
405
    nDOF[VERTEX] += localAdmin->getNumberOfDOFs(VERTEX);
406
407
    nDOFEl += getGeo(VERTEX) * nDOF[VERTEX];

408
    if (dim > 1) {
409
      localAdmin->setNumberOfPreDOFs(EDGE,nDOF[EDGE]);
410
      nDOF[EDGE] += localAdmin->getNumberOfDOFs(EDGE);
411
412
413
414
415
416
417
      nDOFEl += getGeo(EDGE) * nDOF[EDGE];
    }

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

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

420
421
    node[VERTEX] = 0;
    nNodeEl = getGeo(VERTEX);
422

423
424
    if (dim > 1) {
      node[EDGE] = nNodeEl;
425
426
      if (nDOF[EDGE] > 0) 
	nNodeEl += getGeo(EDGE);
427
428
    }

429
    if (dim == 3) {
430
      localAdmin->setNumberOfPreDOFs(FACE,nDOF[FACE]);
431
432
433
434
435
      nDOF[FACE] += localAdmin->getNumberOfDOFs(FACE);
      nDOFEl += getGeo(FACE) * nDOF[FACE];
      node[FACE] = nNodeEl;
      if (nDOF[FACE] > 0) 
	nNodeEl += getGeo(FACE);
436
437
    }

438
    node[CENTER] = nNodeEl;
439
    if (nDOF[CENTER] > 0)
440
      nNodeEl += 1;
441
442
443
444
  }

  void Mesh::dofCompress()
  {
445
446
    FUNCNAME("Mesh::dofCompress()");
    Flag fill_flag;
447

Thomas Witkowski's avatar
Thomas Witkowski committed
448
    for (int iadmin = 0; iadmin < static_cast<int>(admin.size()); iadmin++) {
449
450
451
      compressAdmin = admin[iadmin];

      TEST_EXIT_DBG(compressAdmin)("no admin[%d] in mesh\n", iadmin);
452
      
Thomas Witkowski's avatar
Thomas Witkowski committed
453
454
      int size = compressAdmin->getSize();
      if (size < 1) 
455
	continue;
Thomas Witkowski's avatar
Thomas Witkowski committed
456

457
458
      if (compressAdmin->getUsedDOFs() < 1)    
	continue;
Thomas Witkowski's avatar
Thomas Witkowski committed
459

460
461
      if (compressAdmin->getHoleCount() < 1)    
	continue;
Thomas Witkowski's avatar
Thomas Witkowski committed
462
  
463
464
465
466
467
468
469
470
      newDOF.resize(size);
      
      compressAdmin->compress(newDOF);
      
      if (preserveCoarseDOFs) {
	fill_flag = Mesh::CALL_EVERY_EL_PREORDER | Mesh::FILL_NOTHING;
      } else {
	fill_flag = Mesh::CALL_LEAF_EL | Mesh::FILL_NOTHING;
471
      }
472
      
Thomas Witkowski's avatar
Thomas Witkowski committed
473
474
      traverse(-1, fill_flag, newDOFFct1);
      traverse(-1, fill_flag, newDOFFct2);
475
476
477
      
      newDOF.resize(0);
    }   
478
479
480
481
482
  }


  DegreeOfFreedom *Mesh::getDOF(GeoIndex position)
  {
483
    FUNCNAME("Mesh::getDOF()");
484

485
    TEST_EXIT_DBG(position >= CENTER && position <= FACE)
486
      ("unknown position %d\n", position);
487

488
489
    int ndof = getNumberOfDOFs(position);
    if (ndof <= 0) 
490
      return NULL;
491

492
    DegreeOfFreedom *dof = new DegreeOfFreedom[ndof];
493

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


  DegreeOfFreedom **Mesh::createDOFPtrs()
  {
513
    FUNCNAME("Mesh::createDOFPtrs()");
514
515

    if (nNodeEl <= 0)
516
      return NULL;
517

518
    DegreeOfFreedom **ptrs = new DegreeOfFreedom*[nNodeEl];
519
    for (int i = 0; i < nNodeEl; i++)
520
521
      ptrs[i] = NULL;

522
    return ptrs;
523
524
525
526
  }

  void Mesh::freeDOFPtrs(DegreeOfFreedom **ptrs)
  {
527
    FUNCNAME("Mesh::freeDOFPtrs()");
528

529
    TEST_EXIT_DBG(ptrs)("ptrs=NULL\n");
530
531
532
533

    if (nNodeEl <= 0)
      return;
  
534
    delete [] ptrs;
535
536
537
  }


Thomas Witkowski's avatar
Thomas Witkowski committed
538
  const DOFAdmin *Mesh::createDOFAdmin(std::string lname, DimVec<int> lnDOF)
539
  {
540
    FUNCNAME("Mesh::createDOFAdmin()");
541

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

Thomas Witkowski's avatar
Thomas Witkowski committed
544
545
    for (int i = 0; i < dim + 1; i++)
      localAdmin->setNumberOfDOFs(i, lnDOF[i]);
546
547
548

    addDOFAdmin(localAdmin);

549
    return localAdmin;
550
551
552
553
554
555
556
  }


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

557
558
559
560
561
562
    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];
563
      }
564
565
    }

566
    return localAdmin;
567
568
569
570
  }

  void Mesh::freeDOF(DegreeOfFreedom* dof, GeoIndex position)
  {
571
    FUNCNAME("Mesh::freeDOF()");
572

573
    TEST_EXIT_DBG(position >= CENTER && position <= FACE)
574
      ("unknown position %d\n", position);
575

576
577
578
    int ndof = nDOF[position];
    if (ndof) {
      if (!dof) {
579
	MSG("dof = NULL, but ndof = %d\n", ndof);
580
581
	return;
      }
582
    } else  {
583
584
585
      if (dof)
	MSG("dof != NULL, but ndof = 0\n");

586
587
      return;
    }
588

589
    TEST_EXIT_DBG(ndof <= MAX_DOF)
590
      ("ndof too big: ndof = %d, MAX_DOF = %d\n", ndof, MAX_DOF);
591

592
593
594
595
596
597
598
599
600
601
    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]);
    }
602

603
    delete [] dof;
604
605
606
607
608
  }

  void Mesh::freeElement(Element* el)
  {
    freeDOFPtrs(const_cast<DegreeOfFreedom**>(el->getDOF()));
Thomas Witkowski's avatar
Thomas Witkowski committed
609
    delete el;
610
611
612
613
614
615
  }


  Element* Mesh::createNewElement(Element *parent)
  {
    FUNCNAME("Mesh::createNewElement()");
616
617

    TEST_EXIT_DBG(elementPrototype)("no element prototype\n");
618
619
620

    Element *el = parent ? parent->clone() : elementPrototype->clone();
  
621
    if (!parent && elementDataPrototype) {
622
623
624
625
626
627
628
629
      el->setElementData(elementDataPrototype->clone()); 
    } else {
      el->setElementData(NULL); // must be done in ElementData::refineElementData()
    }

    return el;
  }

630

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

  bool Mesh::findElInfoAtPoint(const WorldVector<double>& xy,
			       ElInfo *el_info,
651
652
			       DimVec<double>& bary,
			       const MacroElement *start_mel,
653
			       const WorldVector<double> *xy0,
654
			       double *sp)
655
656
657
658
659
660
661
662
663
664
  {
    static const MacroElement *mel = NULL;
    DimVec<double> lambda(dim, NO_INIT);
    ElInfo *mel_info = NULL;

    mel_info = createNewElInfo();

    if (start_mel != NULL)
      mel = start_mel;
    else
665
      if ((mel == NULL) || (mel->getElement()->getMesh() != this))
666
667
668
	mel = *(macroElements.begin());

    mel_info->setFillFlag(Mesh::FILL_COORDS);
669
    g_xy = &xy;
670
    g_xy0 = xy0;
671
    g_sp = sp;
672
673
674

    mel_info->fillMacroInfo(mel);

675
    int k;
676
677
678
679
680
681
682
683
684
685
    while ((k = mel_info->worldToCoord(xy, &lambda)) >= 0) {
      if (mel->getNeighbour(k)) {
	mel = mel->getNeighbour(k);
	mel_info->fillMacroInfo(mel);
	continue;
      }
      break;
    }

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

692
    return inside;
693
694
695
  }

  bool Mesh::findElementAtPoint(const WorldVector<double>&  xy,
696
697
				Element **elp, 
				DimVec<double>& bary,
698
				const MacroElement *start_mel,
699
700
				const WorldVector<double> *xy0,
				double *sp)
701
  {
702
703
    ElInfo *el_info = createNewElInfo();
    int val = findElInfoAtPoint(xy, el_info, bary, start_mel, xy0, sp);
704
705
706

    *elp = el_info->getElement();

Thomas Witkowski's avatar
Thomas Witkowski committed
707
    delete el_info;
708

709
    return val;
710
711
  }

712
  bool Mesh::findElementAtPointRecursive(ElInfo *el_info,
713
					 const DimVec<double>& lambda,
714
					 int outside,
715
716
					 ElInfo* final_el_info)
  {
717
    FUNCNAME("Mesh::findElementAtPointRecursive()");
718
719
    Element *el = el_info->getElement();
    DimVec<double> c_lambda(dim, NO_INIT);
720
721
    int inside;
    int ichild, c_outside;
722
723
724
725

    if (el->isLeaf()) {
      *final_el_info = *el_info;
      if (outside < 0) {
726
727
728
729
	for (int i = 0; i <= dim; i++) {
	  final_lambda[i] = lambda[i];
	}

730
	return true;
731
732
733
734
735
736
737
738
739
740
741
742
      }  else {  /* outside */
	if (g_xy0) { /* find boundary point of [xy0, xy] */
	  el_info->worldToCoord(*(g_xy0), &c_lambda);
	  double s = lambda[outside] / (lambda[outside] - c_lambda[outside]);
	  for (int i = 0; i <= dim; i++) {
	    final_lambda[i] = s * c_lambda[i] + (1.0-s) * lambda[i];
	  }
	  if (g_sp) {
	    *(g_sp) = s;
	  }
	  if (dim == 3) 
	    MSG("outside finest level on el %d: s=%.3e\n", el->getIndex(), s);
743

744
	  return false;  /* ??? */
745
	} else {
746
	  return false;
747
	}
748
      }
749
750
    }

751
    ElInfo *c_el_info = createNewElInfo();
752

753
    if (dim == 1) {
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
      if (lambda[0] >= lambda[1]) {
	c_el_info->fillElInfo(0, el_info);
	if (outside >= 0) {
	  outside = el_info->worldToCoord(*(g_xy), &c_lambda);
	  if (outside >= 0) ERROR("point outside domain\n");
	} 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);
	  if (outside >= 0) ERROR("point outside domain\n");
	} else {
	  c_lambda[1] = lambda[1] - lambda[0];
	  c_lambda[0] = 2.0 * lambda[0];
	}
      }
    } /* DIM == 1 */

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

803
    if (dim == 3) {
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
      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);
	  int c_outside2;
	  ElInfo *c_el_info2 = createNewElInfo();

	  c_el_info2->fillElInfo(1-ichild, el_info);
	  c_outside2 = c_el_info2->worldToCoord(*(g_xy), &c_lambda2);

	  MSG("new_coord CHILD %d: outside=%d, lambda=(%.2f %.2f %.2f %.2f)\n",
	      ichild, c_outside, c_lambda[0],c_lambda[1],c_lambda[2],c_lambda[3]);
	  MSG("new_coord CHILD %d: outside=%d, lambda=(%.2f %.2f %.2f %.2f)\n",
	      1-ichild, c_outside2, c_lambda2[0],c_lambda2[1],c_lambda2[2],
	      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
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
	}
	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];
	  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]];
	  c_lambda[3] = 2.0 * lambda[1];
	} else {
	  c_el_info->fillElInfo(1, el_info);
	  c_lambda[0] = lambda[1] - lambda[0];
	  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]];
	  c_lambda[3] = 2.0 * lambda[0];
	}
      }
    }  /* DIM == 3 */

    inside = findElementAtPointRecursive(c_el_info, c_lambda, outside, 
					 final_el_info);
Thomas Witkowski's avatar
Thomas Witkowski committed
859
    delete c_el_info;
860

861
    return inside; 
862
863
  }

864
865
866
867
868
869
870
  bool Mesh::getDofIndexCoords(const DegreeOfFreedom* dof, 
			       const FiniteElemSpace* feSpace,
			       WorldVector<double>& coords)
  {
    DimVec<double>* baryCoords;
    bool found = false;
    TraverseStack stack;
Thomas Witkowski's avatar
Thomas Witkowski committed
871
872
    Vector<DegreeOfFreedom> dofVec(feSpace->getBasisFcts()->getNumber());

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

    return found;
  }
896

897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
  bool Mesh::getDofIndexCoords(DegreeOfFreedom dof, 
			       const FiniteElemSpace* feSpace,
			       WorldVector<double>& coords)
  {
    DimVec<double>* baryCoords;
    bool found = false;
    TraverseStack stack;
    Vector<DegreeOfFreedom> dofVec(feSpace->getBasisFcts()->getNumber());

    ElInfo *elInfo = stack.traverseFirst(this, -1, 
					 Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS);    
    while (elInfo) {
      feSpace->getBasisFcts()->getLocalIndicesVec(elInfo->getElement(),
						  feSpace->getAdmin(),
						  &dofVec);
      for (int i = 0; i < feSpace->getBasisFcts()->getNumber(); i++) {
	if (dofVec[i] == dof) {
	  baryCoords = feSpace->getBasisFcts()->getCoords(i);
	  elInfo->coordToWorld(*baryCoords, coords);
	  found = true;
	  break;	  
	}
      }
      
      if (found)
	break;
      
      elInfo = stack.traverseNext(elInfo);
    }

    return found;
  }

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

Thomas Witkowski's avatar
Thomas Witkowski committed
935
936
937
938
  void Mesh::setDiameter(int i, double w) 
  { 
    diam[i] = w; 
  }
939
940
941
942
943
944
945
946
947
948
949
950


  int Mesh::newDOFFct1(ElInfo* ei) {
    ei->getElement()->newDOFFct1(compressAdmin);
    return 0;
  }

  int Mesh::newDOFFct2(ElInfo* ei) {
    ei->getElement()->newDOFFct2(compressAdmin);
    return 0;
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
951
  void Mesh::serialize(std::ostream &out)
952
953
954
  {
    serializedDOFs.clear();

Thomas Witkowski's avatar
Thomas Witkowski committed
955
    out << name << "\n";
956

957
958
959
960
961
962
963
    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);
964
965
966

    diam.serialize(out);

967
968
    SerUtil::serialize(out, preserveCoarseDOFs);
    SerUtil::serialize(out, nDOFEl);
969
970
971

    nDOF.serialize(out);

972
    SerUtil::serialize(out, nNodeEl);
973
974
975
976

    node.serialize(out);

    // write admins
977
    int size = static_cast<int>(admin.size());
978
    SerUtil::serialize(out, size);
979
    for (int i = 0; i < size; i++)
980
981
982
983
      admin[i]->serialize(out);

    // write macroElements
    size = static_cast<int>(macroElements.size());
984
    SerUtil::serialize(out, size);
985
    for (int i = 0; i < size; i++)
986
987
      macroElements[i]->serialize(out);

988
989
    SerUtil::serialize(out, elementIndex);
    SerUtil::serialize(out, initialized);
990
991
992
993

    serializedDOFs.clear();
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
994
  void Mesh::deserialize(std::istream &in)
995
  {
996
997
    FUNCNAME("Mesh::deserialize()");

998
999
1000
1001
1002
1003
    serializedDOFs.clear();

    in >> name;
    in.get();

    int oldVal = dim;
1004
1005
    SerUtil::deserialize(in, dim);
    TEST_EXIT_DBG(oldVal == 0 || dim == oldVal)("Invalid dimension!\n");
1006

1007
1008
1009
1010
1011
1012
    SerUtil::deserialize(in, nVertices);
    SerUtil::deserialize(in, nEdges);
    SerUtil::deserialize(in, nLeaves);
    SerUtil::deserialize(in, nElements);
    SerUtil::deserialize(in, nFaces);
    SerUtil::deserialize(in, maxEdgeNeigh);
1013
1014
1015

    diam.deserialize(in);

1016
    SerUtil::deserialize(in, preserveCoarseDOFs);
1017
1018

    oldVal = nDOFEl;
1019
1020
    SerUtil::deserialize(in, nDOFEl);
    TEST_EXIT_DBG(oldVal == 0 || nDOFEl == oldVal)("Invalid nDOFEl!\n");
1021
1022
1023
1024

    nDOF.deserialize(in);

    oldVal = nNodeEl;
1025
1026
    SerUtil::deserialize(in, nNodeEl);
    TEST_EXIT_DBG(oldVal == 0 || nNodeEl == oldVal)("Invalid nNodeEl!\n");
1027
1028
1029
1030

    node.deserialize(in);

    // read admins
1031
    int size;
1032
    SerUtil::deserialize(in, size);
1033
    admin.resize(size, NULL);
1034
    for (int i = 0; i < size; i++) {
1035
      if (!admin[i])
Thomas Witkowski's avatar
Thomas Witkowski committed
1036
	admin[i] = new DOFAdmin(this);
1037

1038
1039
1040
      admin[i]->deserialize(in);
    }

1041
    SerUtil::deserialize(in, size);
Thomas Witkowski's avatar
Thomas Witkowski committed
1042
    std::vector< std::vector<int> > neighbourIndices(size);
1043

1044
    for (int i = 0; i < static_cast<int>(macroElements.size()); i++)
1045
      if (macroElements[i])
Thomas Witkowski's avatar
Thomas Witkowski committed
1046
	delete macroElements[i];
1047
1048
1049
1050
1051
1052
1053
1054

    // All macro elements are stored in the list \ref macroElements with continous 
    // index from 0 to n - 1. But the macro element index may be different and may
    // contain holes (for example when macro elements were removed because of domain
    // decomposition based parallelization. Therefore we create a temporary map
    // from macro element indices to the continous index of \ref macroElements. This
    // will be used later to find the correct neighbours of the macro elements.
    std::map<int, int> elIndexVecIndex;
1055

1056
    macroElements.resize(size);
1057
    for (int i = 0; i < size; i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
1058
      macroElements[i] = new MacroElement(dim);
1059
1060
      macroElements[i]->writeNeighboursTo(&(neighbourIndices[i]));
      macroElements[i]->deserialize(in);
1061
      elIndexVecIndex[macroElements[i]->getIndex()] = i;
1062
1063
    }