Mesh.cc 33.1 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();
  
Thomas Witkowski's avatar
Thomas Witkowski committed
595
    if (!parent && elementDataPrototype)
596
      el->setElementData(elementDataPrototype->clone()); 
Thomas Witkowski's avatar
Thomas Witkowski committed
597
    else
598
      el->setElementData(NULL); // must be done in ElementData::refineElementData()
Thomas Witkowski's avatar
Thomas Witkowski committed
599
    
600
601
602
    return el;
  }

603

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

  bool Mesh::findElInfoAtPoint(const WorldVector<double>& xy,
			       ElInfo *el_info,
624
625
			       DimVec<double>& bary,
			       const MacroElement *start_mel,
626
			       const WorldVector<double> *xy0,
627
			       double *sp)
628
629
630
  {
    static const MacroElement *mel = NULL;
    DimVec<double> lambda(dim, NO_INIT);
Thomas Witkowski's avatar
Thomas Witkowski committed
631
    ElInfo *mel_info = createNewElInfo();
632
633
634
635

    if (start_mel != NULL)
      mel = start_mel;
    else
Thomas Witkowski's avatar
Thomas Witkowski committed
636
      if (mel == NULL || mel->getElement()->getMesh() != this)
637
638
639
	mel = *(macroElements.begin());

    mel_info->setFillFlag(Mesh::FILL_COORDS);
640
    g_xy = &xy;
641
    g_xy0 = xy0;
642
    g_sp = sp;
643
644
645

    mel_info->fillMacroInfo(mel);

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

663
    return inside;
664
665
666
  }

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

    *elp = el_info->getElement();

Thomas Witkowski's avatar
Thomas Witkowski committed
678
    delete el_info;
679

680
    return val;
681
682
  }

683
  bool Mesh::findElementAtPointRecursive(ElInfo *el_info,
684
					 const DimVec<double>& lambda,
685
					 int outside,
686
687
					 ElInfo* final_el_info)
  {
688
    FUNCNAME("Mesh::findElementAtPointRecursive()");
Thomas Witkowski's avatar
Thomas Witkowski committed
689

690
691
    Element *el = el_info->getElement();
    DimVec<double> c_lambda(dim, NO_INIT);
692
693
    int inside;
    int ichild, c_outside;
694
695
696
697

    if (el->isLeaf()) {
      *final_el_info = *el_info;
      if (outside < 0) {
Thomas Witkowski's avatar
Thomas Witkowski committed
698
	for (int i = 0; i <= dim; i++)
699
	  final_lambda[i] = lambda[i];
Thomas Witkowski's avatar
Thomas Witkowski committed
700
	
701
	return true;
702
703
704
705
      }  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
706
707
708
709
	  for (int i = 0; i <= dim; i++) 
	    final_lambda[i] = s * c_lambda[i] + (1.0 - s) * lambda[i];
	  
	  if (g_sp)
710
	    *(g_sp) = s;
Thomas Witkowski's avatar
Thomas Witkowski committed
711
	  
712
713
	  if (dim == 3) 
	    MSG("outside finest level on el %d: s=%.3e\n", el->getIndex(), s);
714

715
	  return false;  /* ??? */
716
	} else {
717
	  return false;
718
	}
719
      }
720
721
    }

722
    ElInfo *c_el_info = createNewElInfo();
723

724
    if (dim == 1) {
725
726
727
728
      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
729
	  TEST_EXIT(outside == 0)("point outside domain\n");
730
731
732
733
734
735
736
737
	} 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
738
	  TEST_EXIT(outside == 0)("point outside domain\n");
739
740
741
742
743
744
745
	} else {
	  c_lambda[1] = lambda[1] - lambda[0];
	  c_lambda[0] = 2.0 * lambda[0];
	}
      }
    } /* DIM == 1 */

746
    if (dim == 2) {
747
748
749
750
      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
751
	  TEST_EXIT(outside == 0)("outside curved boundary child 0\n");	  
752
753
754
755
756
757
758
759
760
	} 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
761
	  TEST_EXIT(outside == 0)("outside curved boundary child 1\n");	  
762
763
764
765
766
767
768
769
	} else {
	  c_lambda[0] = lambda[1] - lambda[0];
	  c_lambda[1] = lambda[2];
	  c_lambda[2] = 2.0 * lambda[0];
	}
      }
    } /* DIM == 2 */

770
    if (dim == 3) {
771
772
773
774
775
776
777
778
779
780
781
782
783
      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
784
	  int c_outside2 = c_el_info2->worldToCoord(*(g_xy), &c_lambda2);
785
786
787
788
789
790
791
792

	  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])) {
793
794
	    for (int i = 0; i <= dim; i++) 
	      c_lambda[i] = c_lambda2[i];	    
795
796
797
798
	    c_outside = c_outside2;
	    *c_el_info = *c_el_info2;
	    ichild = 1 - ichild;
	  }
Thomas Witkowski's avatar
Thomas Witkowski committed
799
	  delete c_el_info2;
800
801
802
803
804
805
	}
	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];
806
807
808
809
810
811
	  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]];
812
813
814
815
	  c_lambda[3] = 2.0 * lambda[1];
	} else {
	  c_el_info->fillElInfo(1, el_info);
	  c_lambda[0] = lambda[1] - lambda[0];
816
817
818
819
820
821
	  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]];
822
823
824
825
826
	  c_lambda[3] = 2.0 * lambda[0];
	}
      }
    }  /* DIM == 3 */

Thomas Witkowski's avatar
Thomas Witkowski committed
827
    inside = findElementAtPointRecursive(c_el_info, c_lambda, outside, final_el_info);
Thomas Witkowski's avatar
Thomas Witkowski committed
828
    delete c_el_info;
829

830
    return inside; 
831
832
  }

833
834
835
836
837
838
839
  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
840
841
    Vector<DegreeOfFreedom> dofVec(feSpace->getBasisFcts()->getNumber());

842
    ElInfo *elInfo = stack.traverseFirst(this, -1, 
Thomas Witkowski's avatar
Thomas Witkowski committed
843
					 Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS);    
844
    while (elInfo) {
Thomas Witkowski's avatar
Thomas Witkowski committed
845
846
847
848
849
      feSpace->getBasisFcts()->getLocalIndicesVec(elInfo->getElement(),
						  feSpace->getAdmin(),
						  &dofVec);
      for (int i = 0; i < feSpace->getBasisFcts()->getNumber(); i++) {
	if (dofVec[i] == *dof) {
850
851
852
	  baryCoords = feSpace->getBasisFcts()->getCoords(i);
	  elInfo->coordToWorld(*baryCoords, coords);
	  found = true;
Thomas Witkowski's avatar
Thomas Witkowski committed
853
	  break;	  
854
	}
Thomas Witkowski's avatar
Thomas Witkowski committed
855
856
      }
      
857
858
      if (found)
	break;
Thomas Witkowski's avatar
Thomas Witkowski committed
859
      
860
861
862
863
864
      elInfo = stack.traverseNext(elInfo);
    }

    return found;
  }
865

866
867
868
869
870
871
872
873
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
  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
899
900
901
902
  void Mesh::setDiameter(const WorldVector<double>& w) 
  { 
    diam = w; 
  }
903

Thomas Witkowski's avatar
Thomas Witkowski committed
904
905
906
907
  void Mesh::setDiameter(int i, double w) 
  { 
    diam[i] = w; 
  }
908

909
910
  int Mesh::newDOFFct1(ElInfo* ei) 
  {
911
912
913
914
    ei->getElement()->newDOFFct1(compressAdmin);
    return 0;
  }

915
916
  int Mesh::newDOFFct2(ElInfo* ei) 
  {
917
918
919
920
    ei->getElement()->newDOFFct2(compressAdmin);
    return 0;
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
921
  void Mesh::serialize(std::ostream &out)
922
923
924
  {
    serializedDOFs.clear();

Thomas Witkowski's avatar
Thomas Witkowski committed
925
    out << name << "\n";
926

927
928
929
930
931
932
933
    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);
934
935
936

    diam.serialize(out);

937
938
    SerUtil::serialize(out, preserveCoarseDOFs);
    SerUtil::serialize(out, nDOFEl);
939
940
941

    nDOF.serialize(out);

942
    SerUtil::serialize(out, nNodeEl);
943
944
945

    node.serialize(out);

946
947
948

    // === Write admins. ===

949
    int size = static_cast<int>(admin.size());
950
    SerUtil::serialize(out, size);
951
    for (int i = 0; i < size; i++)
952
953
      admin[i]->serialize(out);

954
955
956

    // === Write macroElements. ===

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

962
963
    SerUtil::serialize(out, elementIndex);
    SerUtil::serialize(out, initialized);
964

965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988

    // === 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);
    }

989
990
991
    serializedDOFs.clear();
  }

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

996
997
998
999
1000
1001
    serializedDOFs.clear();

    in >> name;
    in.get();

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

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

    diam.deserialize(in);

1014
    SerUtil::deserialize(in, preserveCoarseDOFs);
1015
1016

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

    nDOF.deserialize(in);

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

    node.deserialize(in);

1028
1029
1030

    // === Read admins. ===

1031
    int size;
1032
    SerUtil::deserialize(in, size);