Mesh.cc 33.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
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

  //**************************************************************************
  //  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;
Thomas Witkowski's avatar
Thomas Witkowski committed
67
  std::vector<DegreeOfFreedom> Mesh::dof_used;
68
  const int Mesh::MAX_DOF = 100;
69
  std::map<std::pair<DegreeOfFreedom, int>, DegreeOfFreedom*> Mesh::serializedDOFs;
70
71
72

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


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

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

131
    Element::deletedDOFs.clear();
132

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

  Mesh& Mesh::operator=(const Mesh& m)
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
175
    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
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
  void Mesh::updateNumberOfLeaves()
  {
    nLeaves = 0;

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

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

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

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

278
279
    // 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
280
281
    DofElMap dofsOwner;
    DofPosMap dofsPosIndex;
282
283
284
    
    // Determine all dof owner macro elements.
    for (std::deque<MacroElement*>::iterator macroIt = macroElements.begin();
285
	 macroIt != macroElements.end(); ++macroIt) {
Thomas Witkowski's avatar
Thomas Witkowski committed
286
287
288
289
290
      elDofIter.reset((*macroIt)->getElement());
      do {
	dofsOwner[elDofIter.getDofPtr()].insert(*macroIt);
	dofsPosIndex[elDofIter.getDofPtr()] = elDofIter.getPosIndex();
      } while (elDofIter.nextStrict());
291
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
292

293
294
295
296
    // 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
297
298
299
300
301
302
303

      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);

304
305
306
307
      // 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
308
309
	  for (int j = 0; j <= dim; j++)
	    if ((*macroIt)->getNeighbour(i)->getNeighbour(j) == *macroIt)
310
311
312
313
314
	      (*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--;
315
	}
316
317
318
319
320
321
      }

      nLeaves--;
      nElements--;

      // Remove this macro element from the dof owner list.
Thomas Witkowski's avatar
Thomas Witkowski committed
322
323
      for (DofElMap::iterator dofsIt = dofsOwner.begin();
	   dofsIt != dofsOwner.end(); ++dofsIt) {
324
	std::set<MacroElement*>::iterator mIt = dofsIt->second.find(*macroIt);
Thomas Witkowski's avatar
Thomas Witkowski committed
325
	if (mIt != dofsIt->second.end())
326
327
328
329
330
331
332
333
334
335
	  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
336
337
338
    for (DofElMap::iterator dofsIt = dofsOwner.begin(); 
	 dofsIt != dofsOwner.end(); ++dofsIt) {    
      if (dofsIt->second.size() == 0)
339
	freeDOF(const_cast<DegreeOfFreedom*>(dofsIt->first), 
Thomas Witkowski's avatar
Thomas Witkowski committed
340
341
		dofsPosIndex[dofsIt->first]);
      else
342
	nRemainDofs++;
343
344
    }

345
    nVertices = nRemainDofs;
346
  }
347
348
349
350
351
352
353

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

    localAdmin->setMesh(this);

354
    std::vector<DOFAdmin*>::iterator dai = 
355
      std::find(admin.begin(), admin.end(), localAdmin);
356

Thomas Witkowski's avatar
Thomas Witkowski committed
357
358
359
    TEST_EXIT(dai == admin.end())
      ("admin %s is already associated to mesh %s\n",
       localAdmin->getName().c_str(), this->getName().c_str());
360

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

364
365
366
367
368
    admin.push_back(localAdmin);

    nDOFEl = 0;

    localAdmin->setNumberOfPreDOFs(VERTEX,nDOF[VERTEX]);
369
    nDOF[VERTEX] += localAdmin->getNumberOfDOFs(VERTEX);
370
371
    nDOFEl += getGeo(VERTEX) * nDOF[VERTEX];

372
    if (dim > 1) {
373
      localAdmin->setNumberOfPreDOFs(EDGE,nDOF[EDGE]);
374
      nDOF[EDGE] += localAdmin->getNumberOfDOFs(EDGE);
375
376
377
378
379
380
381
      nDOFEl += getGeo(EDGE) * nDOF[EDGE];
    }

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

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

384
385
    node[VERTEX] = 0;
    nNodeEl = getGeo(VERTEX);
386

387
388
    if (dim > 1) {
      node[EDGE] = nNodeEl;
389
390
      if (nDOF[EDGE] > 0) 
	nNodeEl += getGeo(EDGE);
391
392
    }

393
    if (dim == 3) {
394
      localAdmin->setNumberOfPreDOFs(FACE,nDOF[FACE]);
395
396
397
398
399
      nDOF[FACE] += localAdmin->getNumberOfDOFs(FACE);
      nDOFEl += getGeo(FACE) * nDOF[FACE];
      node[FACE] = nNodeEl;
      if (nDOF[FACE] > 0) 
	nNodeEl += getGeo(FACE);
400
401
    }

402
    node[CENTER] = nNodeEl;
403
    if (nDOF[CENTER] > 0)
404
      nNodeEl += 1;
405
406
407
408
  }

  void Mesh::dofCompress()
  {
409
410
    FUNCNAME("Mesh::dofCompress()");
    Flag fill_flag;
411

Thomas Witkowski's avatar
Thomas Witkowski committed
412
    for (int iadmin = 0; iadmin < static_cast<int>(admin.size()); iadmin++) {
413
414
415
      compressAdmin = admin[iadmin];

      TEST_EXIT_DBG(compressAdmin)("no admin[%d] in mesh\n", iadmin);
416
      
Thomas Witkowski's avatar
Thomas Witkowski committed
417
418
      int size = compressAdmin->getSize();
      if (size < 1) 
419
	continue;
Thomas Witkowski's avatar
Thomas Witkowski committed
420

421
422
      if (compressAdmin->getUsedDOFs() < 1)    
	continue;
Thomas Witkowski's avatar
Thomas Witkowski committed
423

424
425
      if (compressAdmin->getHoleCount() < 1)    
	continue;
Thomas Witkowski's avatar
Thomas Witkowski committed
426
  
427
428
429
430
      newDOF.resize(size);
      
      compressAdmin->compress(newDOF);
      
431
      if (preserveCoarseDOFs)
432
	fill_flag = Mesh::CALL_EVERY_EL_PREORDER | Mesh::FILL_NOTHING;
433
      else
434
435
436
	fill_flag = Mesh::CALL_LEAF_EL | Mesh::FILL_NOTHING;
      
      
437
438
439
440
441
442
443
444
445
446
447
448
449
      TraverseStack stack;
      ElInfo *elInfo = stack.traverseFirst(this, -1, fill_flag);
      while (elInfo) {
	newDOFFct1(elInfo);
	elInfo = stack.traverseNext(elInfo);
      }

      elInfo = stack.traverseFirst(this, -1, fill_flag);
      while (elInfo) {
	newDOFFct2(elInfo);
	elInfo = stack.traverseNext(elInfo);
      }

450
      newDOF.resize(0);
451
    }       
452
453
454
455
456
  }


  DegreeOfFreedom *Mesh::getDOF(GeoIndex position)
  {
457
    FUNCNAME("Mesh::getDOF()");
458

459
    TEST_EXIT_DBG(position >= CENTER && position <= FACE)
460
      ("unknown position %d\n", position);
461

462
463
    int ndof = getNumberOfDOFs(position);
    if (ndof <= 0) 
464
      return NULL;
465

466
    DegreeOfFreedom *dof = new DegreeOfFreedom[ndof];
467

468
469
    for (int i = 0; i < getNumberOfDOFAdmin(); i++) {
      const DOFAdmin *localAdmin = &getDOFAdmin(i);
470
      TEST_EXIT_DBG(localAdmin)("no admin[%d]\n", i);
471
472
473
474
      
      int n  = localAdmin->getNumberOfDOFs(position);
      int n0 = localAdmin->getNumberOfPreDOFs(position);
      
475
      TEST_EXIT_DBG(n + n0 <= ndof)("n=%d, n0=%d too large: ndof=%d\n", n, n0, ndof);
476
      
477
      for (int j = 0; j < n; j++)
478
479
	dof[n0 + j] = const_cast<DOFAdmin*>(localAdmin)->getDOFIndex();
    }
480
 
481
    return dof;
482
483
484
485
486
  }


  DegreeOfFreedom **Mesh::createDOFPtrs()
  {
487
    FUNCNAME("Mesh::createDOFPtrs()");
488
489

    if (nNodeEl <= 0)
490
      return NULL;
491

492
    DegreeOfFreedom **ptrs = new DegreeOfFreedom*[nNodeEl];
493
    for (int i = 0; i < nNodeEl; i++)
494
495
      ptrs[i] = NULL;

496
    return ptrs;
497
498
499
500
  }

  void Mesh::freeDOFPtrs(DegreeOfFreedom **ptrs)
  {
501
    FUNCNAME("Mesh::freeDOFPtrs()");
502

503
    TEST_EXIT_DBG(ptrs)("ptrs=NULL\n");
504
505
506
507

    if (nNodeEl <= 0)
      return;
  
508
    delete [] ptrs;
509
510
511
  }


Thomas Witkowski's avatar
Thomas Witkowski committed
512
  const DOFAdmin *Mesh::createDOFAdmin(std::string lname, DimVec<int> lnDOF)
513
  {
514
    FUNCNAME("Mesh::createDOFAdmin()");
515

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

Thomas Witkowski's avatar
Thomas Witkowski committed
518
519
    for (int i = 0; i < dim + 1; i++)
      localAdmin->setNumberOfDOFs(i, lnDOF[i]);
520
521
522

    addDOFAdmin(localAdmin);

523
    return localAdmin;
524
525
526
527
528
529
530
  }


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

531
532
533
534
535
536
    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];
537
      }
538
539
    }

540
    return localAdmin;
541
542
543
544
  }

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

547
    TEST_EXIT_DBG(position >= CENTER && position <= FACE)
548
      ("unknown position %d\n", position);
549

550
551
552
    int ndof = nDOF[position];
    if (ndof) {
      if (!dof) {
553
	MSG("dof = NULL, but ndof = %d\n", ndof);
554
555
	return;
      }
556
    } else  {
557
558
559
      if (dof)
	MSG("dof != NULL, but ndof = 0\n");

560
561
      return;
    }
562

563
    TEST_EXIT_DBG(ndof <= MAX_DOF)
564
      ("ndof too big: ndof = %d, MAX_DOF = %d\n", ndof, MAX_DOF);
565

566
567
568
569
570
571
572
573
574
575
    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]);
    }
576

577
    delete [] dof;
578
579
580
581
582
  }

  void Mesh::freeElement(Element* el)
  {
    freeDOFPtrs(const_cast<DegreeOfFreedom**>(el->getDOF()));
Thomas Witkowski's avatar
Thomas Witkowski committed
583
    delete el;
584
585
586
587
588
589
  }


  Element* Mesh::createNewElement(Element *parent)
  {
    FUNCNAME("Mesh::createNewElement()");
590
591

    TEST_EXIT_DBG(elementPrototype)("no element prototype\n");
592
593
594

    Element *el = parent ? parent->clone() : elementPrototype->clone();
  
595
    if (!parent && elementDataPrototype) {
596
597
598
599
600
601
602
603
      el->setElementData(elementDataPrototype->clone()); 
    } else {
      el->setElementData(NULL); // must be done in ElementData::refineElementData()
    }

    return el;
  }

604

605
606
607
608
  ElInfo* Mesh::createNewElInfo()
  {
    switch(dim) {
    case 1:
Thomas Witkowski's avatar
Thomas Witkowski committed
609
      return new ElInfo1d(this);
610
611
      break;
    case 2:
Thomas Witkowski's avatar
Thomas Witkowski committed
612
      return new ElInfo2d(this);
613
614
      break;
    case 3:
Thomas Witkowski's avatar
Thomas Witkowski committed
615
      return new ElInfo3d(this);
616
617
618
619
      break;
    default:
      ERROR_EXIT("invalid dim\n");
      return NULL;
620
    }
621
622
623
624
  }

  bool Mesh::findElInfoAtPoint(const WorldVector<double>& xy,
			       ElInfo *el_info,
625
626
			       DimVec<double>& bary,
			       const MacroElement *start_mel,
627
			       const WorldVector<double> *xy0,
628
			       double *sp)
629
630
631
632
633
634
635
636
637
638
  {
    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
639
      if ((mel == NULL) || (mel->getElement()->getMesh() != this))
640
641
642
	mel = *(macroElements.begin());

    mel_info->setFillFlag(Mesh::FILL_COORDS);
643
    g_xy = &xy;
644
    g_xy0 = xy0;
645
    g_sp = sp;
646
647
648

    mel_info->fillMacroInfo(mel);

649
    int k;
650
651
652
653
654
655
656
657
658
659
    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 */
660
    bool inside = findElementAtPointRecursive(mel_info, lambda, k, el_info);
661
662
    for (int i = 0; i <= dim; i++)
      bary[i] = final_lambda[i];   
663
  
Thomas Witkowski's avatar
Thomas Witkowski committed
664
    delete mel_info;
665

666
    return inside;
667
668
669
  }

  bool Mesh::findElementAtPoint(const WorldVector<double>&  xy,
670
671
				Element **elp, 
				DimVec<double>& bary,
672
				const MacroElement *start_mel,
673
674
				const WorldVector<double> *xy0,
				double *sp)
675
  {
676
677
    ElInfo *el_info = createNewElInfo();
    int val = findElInfoAtPoint(xy, el_info, bary, start_mel, xy0, sp);
678
679
680

    *elp = el_info->getElement();

Thomas Witkowski's avatar
Thomas Witkowski committed
681
    delete el_info;
682

683
    return val;
684
685
  }

686
  bool Mesh::findElementAtPointRecursive(ElInfo *el_info,
687
					 const DimVec<double>& lambda,
688
					 int outside,
689
690
					 ElInfo* final_el_info)
  {
691
    FUNCNAME("Mesh::findElementAtPointRecursive()");
692
693
    Element *el = el_info->getElement();
    DimVec<double> c_lambda(dim, NO_INIT);
694
695
    int inside;
    int ichild, c_outside;
696
697
698
699

    if (el->isLeaf()) {
      *final_el_info = *el_info;
      if (outside < 0) {
700
701
702
703
	for (int i = 0; i <= dim; i++) {
	  final_lambda[i] = lambda[i];
	}

704
	return true;
705
706
707
708
709
710
711
712
713
714
715
716
      }  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);
717

718
	  return false;  /* ??? */
719
	} else {
720
	  return false;
721
	}
722
      }
723
724
    }

725
    ElInfo *c_el_info = createNewElInfo();
726

727
    if (dim == 1) {
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
      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 */

749
    if (dim == 2) {
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
      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 */

777
    if (dim == 3) {
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
      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])) {
801
802
	    for (int i = 0; i <= dim; i++) 
	      c_lambda[i] = c_lambda2[i];	    
803
804
805
806
	    c_outside = c_outside2;
	    *c_el_info = *c_el_info2;
	    ichild = 1 - ichild;
	  }
Thomas Witkowski's avatar
Thomas Witkowski committed
807
	  delete c_el_info2;
808
809
810
811
812
813
	}
	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];
814
815
816
817
818
819
	  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]];
820
821
822
823
	  c_lambda[3] = 2.0 * lambda[1];
	} else {
	  c_el_info->fillElInfo(1, el_info);
	  c_lambda[0] = lambda[1] - lambda[0];
824
825
826
827
828
829
	  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]];
830
831
832
833
834
835
836
	  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
837
    delete c_el_info;
838

839
    return inside; 
840
841
  }

842
843
844
845
846
847
848
  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
849
850
    Vector<DegreeOfFreedom> dofVec(feSpace->getBasisFcts()->getNumber());

851
    ElInfo *elInfo = stack.traverseFirst(this, -1, 
Thomas Witkowski's avatar
Thomas Witkowski committed
852
					 Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS);    
853
    while (elInfo) {
Thomas Witkowski's avatar
Thomas Witkowski committed
854
855
856
857
858
      feSpace->getBasisFcts()->getLocalIndicesVec(elInfo->getElement(),
						  feSpace->getAdmin(),
						  &dofVec);
      for (int i = 0; i < feSpace->getBasisFcts()->getNumber(); i++) {
	if (dofVec[i] == *dof) {
859
860
861
	  baryCoords = feSpace->getBasisFcts()->getCoords(i);
	  elInfo->coordToWorld(*baryCoords, coords);
	  found = true;
Thomas Witkowski's avatar
Thomas Witkowski committed
862
	  break;	  
863
	}
Thomas Witkowski's avatar
Thomas Witkowski committed
864
865
      }
      
866
867
      if (found)
	break;
Thomas Witkowski's avatar
Thomas Witkowski committed
868
      
869
870
871
872
873
      elInfo = stack.traverseNext(elInfo);
    }

    return found;
  }
874

875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
  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
908
909
910
911
  void Mesh::setDiameter(const WorldVector<double>& w) 
  { 
    diam = w; 
  }
912

Thomas Witkowski's avatar
Thomas Witkowski committed
913
914
915
916
  void Mesh::setDiameter(int i, double w) 
  { 
    diam[i] = w; 
  }
917

918
919
  int Mesh::newDOFFct1(ElInfo* ei) 
  {
920
921
922
923
    ei->getElement()->newDOFFct1(compressAdmin);
    return 0;
  }

924
925
  int Mesh::newDOFFct2(ElInfo* ei) 
  {
926
927
928
929
    ei->getElement()->newDOFFct2(compressAdmin);
    return 0;
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
930
  void Mesh::serialize(std::ostream &out)
931
932
933
  {
    serializedDOFs.clear();

Thomas Witkowski's avatar
Thomas Witkowski committed
934
    out << name << "\n";
935

936
937
938
939
940
941
942
    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);
943
944
945

    diam.serialize(out);

946
947
    SerUtil::serialize(out, preserveCoarseDOFs);
    SerUtil::serialize(out, nDOFEl);
948
949
950

    nDOF.serialize(out);

951
    SerUtil::serialize(out, nNodeEl);
952
953
954

    node.serialize(out);

955
956
957

    // === Write admins. ===

958
    int size = static_cast<int>(admin.size());
959
    SerUtil::serialize(out, size);
960
    for (int i = 0; i < size; i++)
961
962
      admin[i]->serialize(out);

963
964
965

    // === Write macroElements. ===

966
    size = static_cast<int>(macroElements.size());
967
    SerUtil::serialize(out, size);
968
    for (int i = 0; i < size; i++)
969
970
      macroElements[i]->serialize(out);

971
972
    SerUtil::serialize(out, elementIndex);
    SerUtil::serialize(out, initialized);
973

974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997

    // === Write periodic associations. ===
    int mapSize = periodicAssociations.size();
    SerUtil::serialize(out, mapSize);
    for (std::map<BoundaryType, VertexVector*>::iterator it = periodicAssociations.begin();
	 it != periodicAssociations.end(); ++it) {
      BoundaryType b = it->first;

      // Check which DOFAdmin is used for the current VertexVector we want to serialize.
      int ithAdmin = -1;
      for (int i = 0; i < static_cast<int>(admin.size()); i++) {
	if (admin[i] == it->second->getAdmin()) {
	  ithAdmin = i;
	  break;
	}
      }
      TEST_EXIT(ithAdmin >= 0)
	("No DOFAdmin found for serialization of periodic associations!\n");

      SerUtil::serialize(out, b);
      SerUtil::serialize(out, ithAdmin);
      it->second->serialize(out);
    }

998
999
1000
    serializedDOFs.clear();
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
1001
  void Mesh::deserialize(std::istream &in)
1002
  {
1003
1004
    FUNCNAME("Mesh::deserialize()");

1005
1006
1007
1008
1009
1010
    serializedDOFs.clear();

    in >> name;
    in.get();

    int oldVal = dim;
1011
1012
    SerUtil::deserialize(in, dim);
    TEST_EXIT_DBG(oldVal == 0 || dim == oldVal)("Invalid dimension!\n");
1013

1014
1015
1016
1017
1018
1019
    SerUtil::deserialize(in, nVertices);
    SerUtil::deserialize(in, nEdges);
    SerUtil::deserialize(in, nLeaves);
    SerUtil::deserialize(in, nElements);
    SerUtil::deserialize(in, nFaces);
    SerUtil::deserialize(in, maxEdgeNeigh);
1020
1021
1022

    diam.deserialize(in);

1023
    SerUtil::deserialize(in, preserveCoarseDOFs);
1024
1025

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

    nDOF.deserialize(in);

    oldVal = nNodeEl;
1032
1033
    SerUtil::deserialize(in, nNodeEl);
    TEST_EXIT_DBG(oldVal == 0 || nNodeEl == oldVal)("Invalid nNodeEl!\n");
1034
1035
1036

    node.deserialize(in);

1037
1038
1039

    // === Read admins. ===

1040
    int size;
1041
    SerUtil::deserialize(in, size);
1042
    admin.resize(size, NULL);
1043
    for (int i = 0; i < size; i++) {
1044
      if (!admin[i])
Thomas Witkowski's avatar
Thomas Witkowski committed
1045
	admin[i] = new DOFAdmin(this);
1046

1047
1048
1049
      admin[i]->deserialize(in);
    }

1050
    SerUtil::deserialize(in, size);
Thomas Witkowski's avatar
Thomas Witkowski committed
1051
    std::vector< std::vector<int> > neighbourIndices(size);
1052

1053
    for (int i = 0; i < static_cast<int>(macroElements.size()); i++)
1054
      if (macroElements[i])
Thomas Witkowski's avatar
Thomas Witkowski committed
1055
	delete macroElements[i];
1056
1057
1058
1059
1060
1061
1062
1063

    // 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;
1064

1065
    macroElements.resize(size);