Mesh.cc 27.8 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
11
12
13
14
15
16
17
18
19
20
21
22
23
#include "AdaptStationary.h"
#include "AdaptInstationary.h"
#include "FiniteElemSpace.h"
#include "ElementData.h"
#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
24
#include "ElInfoStack.h"
25
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
67

  //**************************************************************************
  //  flags for Mesh traversal                                                
  //**************************************************************************

  const Flag Mesh::CALL_EVERY_EL_PREORDER  = 0X0100L;
  const Flag Mesh::CALL_EVERY_EL_INORDER   = 0X0200L;
  const Flag Mesh::CALL_EVERY_EL_POSTORDER = 0X0400L;
  const Flag Mesh::CALL_LEAF_EL            = 0X0800L;
  const Flag Mesh::CALL_LEAF_EL_LEVEL      = 0X1000L;
  const Flag Mesh::CALL_EL_LEVEL           = 0X2000L;
  const Flag Mesh::CALL_MG_LEVEL           = 0X4000L ; // used in mg methods 


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

  DOFAdmin* Mesh::compressAdmin = NULL;
  Mesh* Mesh::traversePtr = NULL;
  int Mesh::iadmin = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
68
  std::vector<DegreeOfFreedom> Mesh::dof_used;
69
  const int Mesh::MAX_DOF = 100;
Thomas Witkowski's avatar
Thomas Witkowski committed
70
  std::map<DegreeOfFreedom, DegreeOfFreedom*> Mesh::serializedDOFs;
71
72
73

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


Thomas Witkowski's avatar
Thomas Witkowski committed
78
  Mesh::Mesh(const std::string& aName, int dimension) 
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
    : name(aName), 
      dim(dimension), 
      nVertices(0),
      nEdges(0),
      nLeaves(0), 
      nElements(0),
      parametric(NULL), 
      preserveCoarseDOFs(false),
      nDOFEl(0),
      nDOF(dimension, DEFAULT_VALUE, 0),
      nNodeEl(0),
      node(dimension, DEFAULT_VALUE, 0),
      elementPrototype(NULL),
      elementDataPrototype(NULL),
      elementIndex(-1),
      initialized(false),
95
      macroFileInfo(NULL),
96
97
98
      final_lambda(dimension, DEFAULT_VALUE, 0.0)
  {

99
    FUNCNAME("Mesh::Mesh()");
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117

    // set default element prototype
    switch(dim) {
    case 1:
      elementPrototype = NEW Line(this);
      break;
    case 2:
      elementPrototype = NEW Triangle(this);
      break;
    case 3:
      elementPrototype = NEW Tetrahedron(this);
      break;
    default:
      ERROR_EXIT("invalid dimension\n");
    }

    elementPrototype->setIndex(-1);

118
119
    elementIndex = 0;
  }
120
121

  Mesh::~Mesh()
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
  {
    serializedDOFs.empty();

    for (std::deque<MacroElement*>::const_iterator it = macroElements.begin();
	 it != macroElements.end();
	 ++it) {
      (*it)->getElement()->deleteElementDOFs();
      DELETE *it;
    }    

    serializedDOFs.empty();

    if (macroFileInfo) {
      DELETE macroFileInfo;
    }
  }
138
139

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

176
177
178
179
180
181
182
183
184
185
186
187
    /* ====================== Copy macro elements =================== */
    macroElements.clear();
    
    // 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;
188

189
190
    macroElements.clear();

191
192
193
194
195
    // 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++) {
196

197
198
      // Create new MacroElement.
      MacroElement *el = NEW MacroElement(dim);
199

200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
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;
  }

  void Mesh::addMacroElement(MacroElement* m) {
    macroElements.push_back(m); 
    m->setIndex(macroElements.size());
  }
247
248
249
250
251

  int Mesh::traverse(int level, Flag flag, 
		     int (*el_fct)(ElInfo*))
  {
    FUNCNAME("Mesh::traverse()");
Thomas Witkowski's avatar
Thomas Witkowski committed
252

Thomas Witkowski's avatar
Thomas Witkowski committed
253
    std::deque<MacroElement*>::iterator mel;
Thomas Witkowski's avatar
Thomas Witkowski committed
254
255
    ElInfoStack elInfoStack(this);
    ElInfo* elinfo = elInfoStack.getNextElement();
256
257
258
259
260
261
262
263
264
265
266
267
268
    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
269
      sum += tinfo.recursive(&elInfoStack);
270
271
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
272
    elInfoStack.getBackElement();
273
274
275
276
277
278
279
280
281
282
283
284
    
    return (flag.isSet(Mesh::FILL_ADD_ALL)) ? sum : 0;
  }



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

    localAdmin->setMesh(this);

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

287
288
289
290
291
292
293
    if (dai!= admin.end()) {
      ERROR("admin %s is already associated to mesh %s\n",
	    localAdmin->getName().c_str(), this->getName().c_str());
    }

    // ===== adding dofs to already existing elements ============================ 
    
294
295
296
    // If adding DOFAdmins to already initilized meshes is required, see older
    // AMDiS version (revision < 244) at the same code position.
    TEST_EXIT(!initialized)("Adding DOFAdmins to initilized meshes does not work!\n");
297
298
299
300
301
302
303


    admin.push_back(localAdmin);

    nDOFEl = 0;

    localAdmin->setNumberOfPreDOFs(VERTEX,nDOF[VERTEX]);
304
    nDOF[VERTEX] += localAdmin->getNumberOfDOFs(VERTEX);
305
306
    nDOFEl += getGeo(VERTEX) * nDOF[VERTEX];

307
    if (dim > 1) {
308
      localAdmin->setNumberOfPreDOFs(EDGE,nDOF[EDGE]);
309
      nDOF[EDGE] += localAdmin->getNumberOfDOFs(EDGE);
310
311
312
313
314
315
316
      nDOFEl += getGeo(EDGE) * nDOF[EDGE];
    }

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

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

319
320
    node[VERTEX] = 0;
    nNodeEl = getGeo(VERTEX);
321

322
323
    if (dim > 1) {
      node[EDGE] = nNodeEl;
324
325
      if (nDOF[EDGE] > 0) 
	nNodeEl += getGeo(EDGE);
326
327
    }

328
    if (dim == 3) {
329
      localAdmin->setNumberOfPreDOFs(FACE,nDOF[FACE]);
330
331
332
333
334
      nDOF[FACE] += localAdmin->getNumberOfDOFs(FACE);
      nDOFEl += getGeo(FACE) * nDOF[FACE];
      node[FACE] = nNodeEl;
      if (nDOF[FACE] > 0) 
	nNodeEl += getGeo(FACE);
335
336
    }

337
    node[CENTER] = nNodeEl;
338
    if (nDOF[CENTER] > 0) {
339
      nNodeEl += 1;
340
    }
341
342
343
344
345
346
347
348
349
  }


  /****************************************************************************/
  /*  dofCompress: remove holes in dof vectors                                */
  /****************************************************************************/

  void Mesh::dofCompress()
  {
350
351
352
    FUNCNAME("Mesh::dofCompress()");
    int size;
    Flag fill_flag;
353

354
    for (iadmin = 0; iadmin < static_cast<int>(admin.size()); iadmin++) {
355
356
357
      compressAdmin = admin[iadmin];

      TEST_EXIT_DBG(compressAdmin)("no admin[%d] in mesh\n", iadmin);
358
359
360
      
      if ((size = compressAdmin->getSize()) < 1) 
	continue;
Thomas Witkowski's avatar
Thomas Witkowski committed
361

362
363
      if (compressAdmin->getUsedDOFs() < 1)    
	continue;
Thomas Witkowski's avatar
Thomas Witkowski committed
364

365
366
      if (compressAdmin->getHoleCount() < 1)    
	continue;
Thomas Witkowski's avatar
Thomas Witkowski committed
367
  
368
369
370
371
372
373
374
375
      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;
376
      }
377
378
379
380
381
382
      
      traverse( -1, fill_flag, newDOFFct1);
      traverse( -1, fill_flag, newDOFFct2);
      
      newDOF.resize(0);
    }   
383
384
385
386
387
  }


  DegreeOfFreedom *Mesh::getDOF(GeoIndex position)
  {
388
    FUNCNAME("Mesh::getDOF()");
389

390
    TEST_EXIT_DBG(position >= CENTER && position <= FACE)
391
      ("unknown position %d\n", position);
392

393
394
395
    int ndof = getNumberOfDOFs(position);
    if (ndof <= 0) 
      return(NULL);
396

397
    DegreeOfFreedom *dof = GET_MEMORY(DegreeOfFreedom, ndof);
398

399
400
    for (int i = 0; i < getNumberOfDOFAdmin(); i++) {
      const DOFAdmin *localAdmin = &getDOFAdmin(i);
401
      TEST_EXIT_DBG(localAdmin)("no admin[%d]\n", i);
402
403
404
405
      
      int n  = localAdmin->getNumberOfDOFs(position);
      int n0 = localAdmin->getNumberOfPreDOFs(position);
      
406
      TEST_EXIT_DBG(n + n0 <= ndof)("n=%d, n0=%d too large: ndof=%d\n", n, n0, ndof);
407
408
409
      
      for (int j = 0; j < n; j++) {
	dof[n0 + j] = const_cast<DOFAdmin*>(localAdmin)->getDOFIndex();
410
      }
411
    }
412
413
414
415
416
417
418
  
    return(dof);
  }


  DegreeOfFreedom **Mesh::createDOFPtrs()
  {
419
    FUNCNAME("Mesh::createDOFPtrs()");
420
421
422
423

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

424
425
    DegreeOfFreedom **ptrs = GET_MEMORY(DegreeOfFreedom*, nNodeEl);
    for (int i = 0; i < nNodeEl; i++)
426
427
428
429
430
431
432
      ptrs[i] = NULL;

    return(ptrs);
  }

  void Mesh::freeDOFPtrs(DegreeOfFreedom **ptrs)
  {
433
    FUNCNAME("Mesh::freeDOFPtrs()");
434

435
    TEST_EXIT_DBG(ptrs)("ptrs=NULL\n");
436
437
438
439
440
441
442
443

    if (nNodeEl <= 0)
      return;
  
    FREE_MEMORY(ptrs, DegreeOfFreedom*, nNodeEl);
  }


Thomas Witkowski's avatar
Thomas Witkowski committed
444
  const DOFAdmin *Mesh::createDOFAdmin(const std::string& lname,DimVec<int> lnDOF)
445
  {
446
    FUNCNAME("Mesh::createDOFAdmin()");
447

448
    DOFAdmin *localAdmin = NEW DOFAdmin(this, lname);
449

450
    for (int i = 0; i < dim+1; i++)
451
452
453
454
455
456
457
458
459
460
461
462
      localAdmin->setNumberOfDOFs(i,lnDOF[i]);

    addDOFAdmin(localAdmin);

    return(localAdmin);
  }


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

463
464
465
466
467
468
    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];
469
      }
470
471
    }

472
473
474
475
476
    return(localAdmin);
  }

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

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

482
483
484
485
    int ndof = nDOF[position];
    if (ndof) {
      if (!dof) {
	MSG("dof = NULL, but ndof=%d\n", ndof);
486
487
	return;
      }
488
489
490
491
492
493
    } else  {
      if (dof) {
	MSG("dof != NULL, but ndof=0\n");
      }
      return;
    }
494

495
    TEST_EXIT_DBG(ndof <= MAX_DOF)
496
497
      ("ndof too big: ndof=%d, MAX_DOF=%d\n",ndof,MAX_DOF);

498
499
500
501
502
503
504
505
506
507
    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]);
    }
508
509
510
511
512
513
514
515
516
517
518
519
520
521

    FREE_MEMORY(dof, DegreeOfFreedom, ndof);
  }

  void Mesh::freeElement(Element* el)
  {
    freeDOFPtrs(const_cast<DegreeOfFreedom**>(el->getDOF()));
    DELETE el;
  }


  Element* Mesh::createNewElement(Element *parent)
  {
    FUNCNAME("Mesh::createNewElement()");
522
523

    TEST_EXIT_DBG(elementPrototype)("no element prototype\n");
524
525
526

    Element *el = parent ? parent->clone() : elementPrototype->clone();
  
527
    if (!parent && elementDataPrototype) {
528
529
530
531
532
533
534
535
      el->setElementData(elementDataPrototype->clone()); 
    } else {
      el->setElementData(NULL); // must be done in ElementData::refineElementData()
    }

    return el;
  }

536

537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
  ElInfo* Mesh::createNewElInfo()
  {
    switch(dim) {
    case 1:
      return NEW ElInfo1d(this);
      break;
    case 2:
      return NEW ElInfo2d(this);
      break;
    case 3:
      return NEW ElInfo3d(this);
      break;
    default:
      ERROR_EXIT("invalid dim\n");
      return NULL;
    };
  }



  bool Mesh::findElInfoAtPoint(const WorldVector<double>& xy,
			       ElInfo *el_info,
559
560
			       DimVec<double>& bary,
			       const MacroElement *start_mel,
561
			       const WorldVector<double> *xy0,
562
			       double *sp)
563
564
565
566
567
568
569
570
571
572
  {
    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
573
      if ((mel == NULL) || (mel->getElement()->getMesh() != this))
574
575
576
	mel = *(macroElements.begin());

    mel_info->setFillFlag(Mesh::FILL_COORDS);
577
    g_xy = &xy;
578
    g_xy0 = xy0;
579
    g_sp = sp;
580
581
582

    mel_info->fillMacroInfo(mel);

583
    int k;
584
585
586
587
588
589
590
591
592
593
    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 */
594
595
596
597
    bool inside = findElementAtPointRecursive(mel_info, lambda, k, el_info);
    for (int i = 0; i <= dim; i++) {
      bary[i] = final_lambda[i];
    }
598
599
600
601
602
603
604
  
    DELETE mel_info;

    return(inside);
  }

  bool Mesh::findElementAtPoint(const WorldVector<double>&  xy,
605
606
				Element **elp, 
				DimVec<double>& bary,
607
				const MacroElement *start_mel,
608
609
				const WorldVector<double> *xy0,
				double *sp)
610
  {
611
612
    ElInfo *el_info = createNewElInfo();
    int val = findElInfoAtPoint(xy, el_info, bary, start_mel, xy0, sp);
613
614
615
616
617
618
619
620
621
622

    *elp = el_info->getElement();

    DELETE el_info;

    return(val);
  }



623
  bool Mesh::findElementAtPointRecursive(ElInfo *el_info,
624
					 const DimVec<double>& lambda,
625
					 int outside,
626
627
					 ElInfo* final_el_info)
  {
628
    FUNCNAME("Mesh::findElementAtPointRecursive()");
629
630
    Element *el = el_info->getElement();
    DimVec<double> c_lambda(dim, NO_INIT);
631
632
    int inside;
    int ichild, c_outside;
633
634
635
636

    if (el->isLeaf()) {
      *final_el_info = *el_info;
      if (outside < 0) {
637
638
639
640
	for (int i = 0; i <= dim; i++) {
	  final_lambda[i] = lambda[i];
	}

641
	return(true);
642
643
644
645
646
647
648
649
650
651
652
653
      }  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);
654

655
	  return(false);  /* ??? */
656
657
	} else {
	  return(false);
658
	}
659
      }
660
661
    }

662
    ElInfo *c_el_info = createNewElInfo();
663

664
    if (dim == 1) {
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
      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 */

686
    if (dim == 2) {
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
      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 */

714
    if (dim == 3) {
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
      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])) {
738
739
740
	    for (int i = 0; i <= dim; i++) {
	      c_lambda[i] = c_lambda2[i];
	    }
741
742
743
744
745
746
747
748
749
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
	    c_outside = c_outside2;
	    *c_el_info = *c_el_info2;
	    ichild = 1 - ichild;
	  }
	  DELETE c_el_info2;
	}
	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);
    DELETE c_el_info;

    return(inside); 
  }


Thomas Witkowski's avatar
Thomas Witkowski committed
777
778
779
780
  void Mesh::setDiameter(const WorldVector<double>& w) 
  { 
    diam = w; 
  }
781

Thomas Witkowski's avatar
Thomas Witkowski committed
782
783
784
785
  void Mesh::setDiameter(int i, double w) 
  { 
    diam[i] = w; 
  }
786
787
788
789
790
791
792
793
794
795
796
797


  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
798
  void Mesh::serialize(std::ostream &out)
799
800
801
802
  {
    serializedDOFs.clear();

    // write name
Thomas Witkowski's avatar
Thomas Witkowski committed
803
    out << name << "\n";
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
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
859
860
861
862
863
864
865
866

    // 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
    int i, size = static_cast<int>(admin.size());
    out.write(reinterpret_cast<const char*>(&size), sizeof(int));
    for (i = 0; i < size; i++) {
      admin[i]->serialize(out);
    }

    // write macroElements
    size = static_cast<int>(macroElements.size());
    out.write(reinterpret_cast<const char*>(&size), sizeof(int));
    for (i = 0; i < size; i++) {
      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
867
  void Mesh::deserialize(std::istream &in)
868
869
870
871
872
873
874
875
876
877
  {
    serializedDOFs.clear();

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

    // read dim
    int oldVal = dim;
    in.read(reinterpret_cast<char*>(&dim), sizeof(int));
878
    TEST_EXIT_DBG((oldVal == 0) || (dim == oldVal))("invalid dimension\n");
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

    // 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));
907
    TEST_EXIT_DBG((oldVal == 0) || (nDOFEl == oldVal))("invalid nDOFEl\n");
908
909
910
911
912
913
914

    // read nDOF
    nDOF.deserialize(in);

    // read nNodeEl
    oldVal = nNodeEl;
    in.read(reinterpret_cast<char*>(&nNodeEl), sizeof(int));
915
    TEST_EXIT_DBG((oldVal == 0) || (nNodeEl == oldVal))("invalid nNodeEl\n");
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933

    // read node
    node.deserialize(in);

    // read admins
    int i, size;
    in.read(reinterpret_cast<char*>(&size), sizeof(int));
    admin.resize(size, NULL);
    for (i = 0; i < size; i++) {
      if (!admin[i]) {
	admin[i] = NEW DOFAdmin(this);
      }
      admin[i]->deserialize(in);
    }

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

Thomas Witkowski's avatar
Thomas Witkowski committed
934
    std::vector< std::vector<int> > neighbourIndices(size);
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979

    for (i = 0; i < static_cast<int>(macroElements.size()); i++) {
      if (macroElements[i]) {
	DELETE macroElements[i];
      }
    }
    macroElements.resize(size);
    for(i = 0; i < size; i++) {
      macroElements[i] = NEW MacroElement(dim);
      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
    int j, neighs = getGeo(NEIGH);
    for(i = 0; i < static_cast<int>(macroElements.size()); i++) {
      for(j = 0; j < neighs; j++) {
	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);
    while(elInfo) {
      elInfo->getElement()->setMesh(this);
      elInfo = stack.traverseNext(elInfo);
    }

    serializedDOFs.clear();
  }

  void Mesh::initialize() 
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
980
981
982
    std::string macroFilename("");
    std::string valueFilename("");
    std::string periodicFile("");
983
984
985
986
987
988
989
990
991
    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()) {
992
993
994
      macroFileInfo = MacroReader::readMacro(macroFilename.c_str(), this,
					     periodicFile == "" ? NULL : periodicFile.c_str(),
					     check);
995
996
997
998
999
1000
1001
1002
1003
1004
1005
1006

      // If there is no value file which should be written, we can delete
      // the information of the macro file.
      if (!valueFilename.length()) {
	clearMacroFileInfo();
      }
    }

    initialized = true;
  }

  bool Mesh::associated(DegreeOfFreedom dof1, DegreeOfFreedom dof2) {
Thomas Witkowski's avatar
Thomas Witkowski committed
1007
1008
    std::map<BoundaryType, VertexVector*>::iterator it;
    std::map<BoundaryType, VertexVector*>::iterator end = periodicAssociations.end();
1009
1010
1011
1012
1013
1014
1015
1016
    for (it = periodicAssociations.begin(); it != end; ++it) {
      if ((*(it->second))[dof1] == dof2)
	return true;
    }
    return false;
  }

  bool Mesh::indirectlyAssociated(DegreeOfFreedom dof1, DegreeOfFreedom dof2) {
Thomas Witkowski's avatar
Thomas Witkowski committed
1017
1018
1019
    std::vector<DegreeOfFreedom> associatedToDOF1;
    std::map<BoundaryType, VertexVector*>::iterator it;
    std::map<BoundaryType, VertexVector*>::iterator end = periodicAssociations.end();
1020
1021
1022
    DegreeOfFreedom dof, assDOF;

    associatedToDOF1.push_back(dof1);
Thomas Witkowski's avatar
Thomas Witkowski committed
1023
1024
1025
    for (it = periodicAssociations.begin(); it != end; ++it) {
      int size = static_cast<int>(associatedToDOF1.size());
      for (int i = 0; i < size; i++) {
1026
1027
	dof = associatedToDOF1[i];
	assDOF = (*(it->second))[dof];
Thomas Witkowski's avatar
Thomas Witkowski committed
1028
	if (assDOF == dof2) {
1029
1030
	  return true;
	} else {
Thomas Witkowski's avatar
Thomas Witkowski committed
1031
1032
	  if (assDOF != dof) 
	    associatedToDOF1.push_back(assDOF);
1033
1034
1035
1036
1037
1038
1039
1040
	}
      }
    }
    return false;
  }

  void Mesh::clearMacroFileInfo()
  {
1041
1042
1043
    macroFileInfo->clear(getNumberOfEdges(),
			 getNumberOfVertices());
    DELETE macroFileInfo;
1044
  }
1045
1046
1047
1048
1049
1050
1051
1052
1053
1054
1055
1056

  int Mesh::calcMemoryUsage()
  {
    int result = 0;

    result += sizeof(Mesh);
    for (int i = 0; i < static_cast<int>(macroElements.size()); i++) {
      result += macroElements[i]->calcMemoryUsage();
    }
    
    return result;
  }
1057
}