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
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45

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;


46
47
48
  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);
49
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;
  Mesh* Mesh::traversePtr = 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
96
97
      final_lambda(dimension, DEFAULT_VALUE, 0.0)
  {

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

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

    elementPrototype->setIndex(-1);

117
118
    elementIndex = 0;
  }
119
120

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

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

131
    Element::deletedDOFs.clear();
132

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

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

185

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

198
199
    macroElements.clear();

200
201
202
203
204
    // 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++) {
205

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

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

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
275
276
277
278
279
280
281
282
    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);
283

284
285
    // 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
286
287
    DofElMap dofsOwner;
    DofPosMap dofsPosIndex;
288
289
290
291
292
    
    // 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
293
294
295
296
297
      elDofIter.reset((*macroIt)->getElement());
      do {
	dofsOwner[elDofIter.getDofPtr()].insert(*macroIt);
	dofsPosIndex[elDofIter.getDofPtr()] = elDofIter.getPosIndex();
      } while (elDofIter.nextStrict());
298
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
299

300
301
302
303
    // 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
304
305
306
307
308
309
310

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

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

      nLeaves--;
      nElements--;

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
340

341
342
343
    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
347
348
349
    for (DofElMap::iterator dofsIt = dofsOwner.begin(); 
	 dofsIt != dofsOwner.end(); ++dofsIt) {    
      if (dofsIt->second.size() == 0)
	freeDOF(const_cast<DegreeOfFreedom*>(dofsIt->first),
		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

Thomas Witkowski's avatar
Thomas Witkowski committed
360
    std::deque<MacroElement*>::iterator mel;
Thomas Witkowski's avatar
Thomas Witkowski committed
361
362
    ElInfoStack elInfoStack(this);
    ElInfo* elinfo = elInfoStack.getNextElement();
363
364
365
366
367
368
369
370
371
372
373
374
375
    Traverse tinfo(this, flag, level, el_fct);
    int sum = 0;
  
    elinfo->setFillFlag(flag);
  
    if (flag.isSet(Mesh::CALL_LEAF_EL_LEVEL) || 
	flag.isSet(Mesh::CALL_EL_LEVEL)      || 
	flag.isSet(Mesh::CALL_MG_LEVEL)) {
      TEST(level >= 0)("invalid level: %d\n", level);
    }
  
    for (mel = macroElements.begin(); mel != macroElements.end(); mel++) {
      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
391
    std::vector<DOFAdmin*>::iterator dai = 
      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
579
    int ndof = nDOF[position];
    if (ndof) {
      if (!dof) {
	MSG("dof = NULL, but ndof=%d\n", ndof);
580
581
	return;
      }
582
583
584
585
586
587
    } else  {
      if (dof) {
	MSG("dof != NULL, but ndof=0\n");
      }
      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
955
  {
    serializedDOFs.clear();

    // write name
Thomas Witkowski's avatar
Thomas Witkowski committed
956
    out << name << "\n";
957
958
959
960
961
962
963
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
989
990
991
992
993
994
995
996
997

    // write dim
    out.write(reinterpret_cast<const char*>(&dim), sizeof(int));

    // write nVertices
    out.write(reinterpret_cast<const char*>(&nVertices), sizeof(int));

    // write nEdges
    out.write(reinterpret_cast<const char*>(&nEdges), sizeof(int));

    // write nLeaves
    out.write(reinterpret_cast<const char*>(&nLeaves), sizeof(int));

    // write nElements
    out.write(reinterpret_cast<const char*>(&nElements), sizeof(int));

    // write nFaces
    out.write(reinterpret_cast<const char*>(&nFaces), sizeof(int));

    // write maxEdgeNeigh
    out.write(reinterpret_cast<const char*>(&maxEdgeNeigh), sizeof(int));

    // write diam
    diam.serialize(out);

    // write preserveCoarseDOFs
    out.write(reinterpret_cast<const char*>(&preserveCoarseDOFs), sizeof(bool));

    // write nDOFEl
    out.write(reinterpret_cast<const char*>(&nDOFEl), sizeof(int));

    // write nDOF
    nDOF.serialize(out);

    // write nNodeEl
    out.write(reinterpret_cast<const char*>(&nNodeEl), sizeof(int));

    // write node
    node.serialize(out);

    // write admins
998
    int size = static_cast<int>(admin.size());
999
    out.write(reinterpret_cast<const char*>(&size), sizeof(int));
1000
    for (int i = 0; i < size; i++)
1001
1002
1003
1004
1005
      admin[i]->serialize(out);

    // write macroElements
    size = static_cast<int>(macroElements.size());
    out.write(reinterpret_cast<const char*>(&size), sizeof(int));
1006
    for (int i = 0; i < size; i++)
1007
1008
1009
1010
1011
1012
1013
1014
1015
1016
1017
      macroElements[i]->serialize(out);

    // write elementIndex
    out.write(reinterpret_cast<const char*>(&elementIndex), sizeof(int));

    // write initialized
    out.write(reinterpret_cast<const char*>(&initialized), sizeof(bool));

    serializedDOFs.clear();
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
1018
  void Mesh::deserialize(std::istream &in)
1019
1020
1021
1022
1023
1024
1025
1026
1027
1028
  {
    serializedDOFs.clear();

    // read name
    in >> name;
    in.get();

    // read dim
    int oldVal = dim;
    in.read(reinterpret_cast<char*>(&dim), sizeof(int));
1029
    TEST_EXIT_DBG((oldVal == 0) || (dim == oldVal))("invalid dimension\n");
1030
1031
1032
1033
1034
1035
1036
1037
1038
1039
1040
1041
1042
1043
1044
1045
1046
1047
1048
1049
1050
1051
1052
1053
1054
1055
1056
1057

    // read nVertices
    in.read(reinterpret_cast<char*>(&nVertices), sizeof(int));

    // read nEdges
    in.read(reinterpret_cast<char*>(&nEdges), sizeof(int));

    // read nLeaves
    in.read(reinterpret_cast<char*>(&nLeaves), sizeof(int));

    // read nElements
    in.read(reinterpret_cast<char*>(&nElements), sizeof(int));

    // read nFaces
    in.read(reinterpret_cast<char*>(&nFaces), sizeof(int));

    // read maxEdgeNeigh
    in.read(reinterpret_cast<char*>(&maxEdgeNeigh), sizeof(int));

    // diam
    diam.deserialize(in);

    // read preserveCoarseDOFs
    in.read(reinterpret_cast<char*>(&preserveCoarseDOFs), sizeof(bool));

    // read nDOFEl
    oldVal = nDOFEl;
    in.read(reinterpret_cast<char*>(&nDOFEl), sizeof(int));
1058
    TEST_EXIT_DBG((oldVal == 0) || (nDOFEl == oldVal))("invalid nDOFEl\n");
1059
1060
1061
1062
1063
1064
1065

    // read nDOF
    nDOF.deserialize(in);

    // read nNodeEl
    oldVal = nNodeEl;
    in.read(reinterpret_cast<char*>(&nNodeEl), sizeof(int));
1066
    TEST_EXIT_DBG((oldVal == 0) || (nNodeEl == oldVal))("invalid nNodeEl\n");
1067
1068
1069
1070
1071

    // read node
    node.deserialize(in);

    // read admins
1072
    int size;
1073
1074
    in.read(reinterpret_cast<char*>(&size), sizeof(int));
    admin.resize(size, NULL);
1075
    for (int i = 0; i < size; i++) {
1076
      if (!admin[i])
Thomas Witkowski's avatar
Thomas Witkowski committed
1077
	admin[i] = new DOFAdmin(this);
1078

1079
1080
1081
1082
1083
1084
      admin[i]->deserialize(in);
    }

    // read macroElements
    in.read(reinterpret_cast<char*>(&size), sizeof(int));

Thomas Witkowski's avatar
Thomas Witkowski committed
1085
    std::vector< std::vector<int> > neighbourIndices(size);
1086

1087
    for (int i = 0; i < static_cast<int>(macroElements.size()); i++) {
1088
      if (macroElements[i])
Thomas Witkowski's avatar
Thomas Witkowski committed
1089
	delete macroElements[i];
1090
    }
1091

1092
    macroElements.resize(size);
1093
    for (int i = 0; i < size; i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
1094
      macroElements[i] = new MacroElement(dim);
1095
1096
1097
1098
1099
1100
1101
1102
1103
1104
1105
      macroElements[i]->writeNeighboursTo(&(neighbourIndices[i]));
      macroElements[i]->deserialize(in);
    }

    // read elementIndex
    in.read(reinterpret_cast<char*>(&elementIndex), sizeof(int));

    // read initialized
    in.read(reinterpret_cast<char*>(&initialized), sizeof(bool));

    // set neighbour pointer in macro elements
1106
1107
1108
    int neighs = getGeo(NEIGH);
    for (int i = 0; i < static_cast<int>(macroElements.size()); i++) {
      for (int j = 0; j < neighs; j++) {
1109
1110
1111
1112
1113
1114
1115
1116
1117
1118
1119
1120
	int index = neighbourIndices[i][j];
	if(index != -1) {
	  macroElements[i]->setNeighbour(j, macroElements[index]);
	} else {
	  macroElements[i]->setNeighbour(j, NULL);
	}
      }
    }

    // set mesh pointer in elements
    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(this, -1, CALL_EVERY_EL_PREORDER);
1121
    while (elInfo) {
1122
1123
1124
1125
1126
1127
1128
1129
1130
      elInfo->getElement()->setMesh(this);
      elInfo = stack.traverseNext(elInfo);
    }

    serializedDOFs.clear();
  }

  void Mesh::initialize() 
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
1131
1132
1133
    std::string macroFilename("");
    std::string valueFilename("");
    std::string periodicFile("");
1134
1135
1136
1137
1138
1139
1140
1141
1142
    int check = 1;

    GET_PARAMETER(0, name + "->macro file name",  &macroFilename);
    GET_PARAMETER(0, name + "->value file name",  &valueFilename);
    GET_PARAMETER(0, name + "->periodic file", &periodicFile);
    GET_PARAMETER(0, name + "->check", "%d", &check);
    GET_PARAMETER(0, name + "->preserve coarse dofs", "%d", &preserveCoarseDOFs);

    if (macroFilename.length()) {
Thomas Witkowski's avatar