Mesh.cc 34.6 KB
Newer Older
Thomas Witkowski's avatar
Thomas Witkowski committed
1
2
3
4
5
6
#include <algorithm>
#include <set>
#include <map>

#include "time.h"

7
8
9
10
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

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;


45
46
47
  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);
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65

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

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


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

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

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

    elementPrototype->setIndex(-1);

116
117
    elementIndex = 0;
  }
118
119

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

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

130
    Element::deletedDOFs.clear();
131

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

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

184

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

197
198
    macroElements.clear();

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

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

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

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

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

269
  void Mesh::removeMacroElements(std::vector<MacroElement*>& macros) 
270
271
272
273
274
  {
    FUNCNAME("Mesh::removeMacroElement()");

    TEST_EXIT(dim == 2)("Not yet implemented!\n");

275
276
277
278
279
280
281
282
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.
    std::map<const DegreeOfFreedom*, std::set<MacroElement*> > dofsOwner;
    
    // Determine all dof owner macro elements.
    for (std::deque<MacroElement*>::iterator macroIt = macroElements.begin();
	 macroIt != macroElements.end();
	 ++macroIt) {
      Element *el = (*macroIt)->getElement();
      for (int i = 0; i < 3; i++)
	dofsOwner[el->getDOF(i)].insert(*macroIt);      
286
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
287
   
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
    // Remove all the given macro elements.
    for (std::vector<MacroElement*>::iterator macroIt = macros.begin();
	 macroIt != macros.end();
	 ++macroIt) {
      bool found = false;

      // Remove the macro element from mesh's list of all macro elements.
      for (std::deque<MacroElement*>::iterator it = macroElements.begin();
	   it != macroElements.end();
	   ++it) {
	if (*it == *macroIt) {
	  macroElements.erase(it, it + 1);
	  found = true;
	  break;
	}
      }
      
      TEST_EXIT(found)("Cannot find MacroElement that should be removed!\n");
      
      // 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)) {
	  for (int j = 0; j <= dim; j++) {
	    if ((*macroIt)->getNeighbour(i)->getNeighbour(j) == *macroIt) {
	      (*macroIt)->getNeighbour(i)->setNeighbour(j, NULL);
	    }
315
	  }
316
317
318
319
	} else {
	  // There is no neighbour at this edge, so we have to decrease the number
	  // of edges in the mesh.
	  nEdges--;
320
	}
321
322
323
324
325
326
327
328
329
330
      }

      nLeaves--;
      nElements--;

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

      // 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.
    for (std::map<const DegreeOfFreedom*, std::set<MacroElement*> >::iterator dofsIt = dofsOwner.begin();
	 dofsIt != dofsOwner.end();
	 ++dofsIt) {    
      if (dofsIt->second.size() == 0) {
Thomas Witkowski's avatar
Thomas Witkowski committed
347
	freeDOF(const_cast<DegreeOfFreedom*>(dofsIt->first), VERTEX);
348
      } else {
349
	nRemainDofs++;
350
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);

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

    // ===== adding dofs to already existing elements ============================ 
    
Thomas Witkowski's avatar
Thomas Witkowski committed
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
    TEST_EXIT(!initialized)("Adding DOFAdmins to initilized meshes does not work yet!\n");

    if (initialized) {
      static bool pnd_1d_0[2] = {true, true};
      static bool pnd_1d_1[1] = {false};
      static bool pnd_2d_0[3] = {true, true, true};
      static bool pnd_2d_1[3] = {true, true, false};
      static bool pnd_2d_2[1] = {false};
      static bool pnd_3d_0[4] = {true, true, true, true};
      static bool pnd_3d_1[6] = {false, true, true, true, true, true};
      static bool pnd_3d_2[4] = {true, true, false, false};
      static bool pnd_3d_3[1] = {false};
      static bool *pnd_1d[2] = {pnd_1d_0, pnd_1d_1};
      static bool *pnd_2d[3] = {pnd_2d_0, pnd_2d_1, pnd_2d_2};
      static bool *pnd_3d[4] = {pnd_3d_0, pnd_3d_1, pnd_3d_2, pnd_3d_3};
      static bool **parentNeedsDOF[4] = {NULL, pnd_1d, pnd_2d, pnd_3d};

     
      std::list<struct delmem> delList;
      std::map< std::set<DegreeOfFreedom>, DegreeOfFreedom*> dofPtrMap;
      const DOFAdmin *vertexAdmin = getVertexAdmin();
      int vertexAdminPreDOFs = vertexAdmin->getNumberOfPreDOFs(VERTEX);

      // finding necessary node number for new admin

      int newNNode = 0;
      GeoIndex geoIndex;

      for (int d = 0; d < dim + 1; d++) {
	geoIndex = INDEX_OF_DIM(d, dim);
      
	if (localAdmin->getNumberOfDOFs(geoIndex)>0||nDOF[geoIndex]>0)
	  newNNode+=getGeo(geoIndex);
      }

      bool extendNodes = (newNNode>nNodeEl);
      nNodeEl = newNNode;

      TraverseStack stack;
      ElInfo *elInfo = NULL;
    
      elInfo = stack.traverseFirst(this, -1, CALL_EVERY_EL_PREORDER);
      while(elInfo) {
	Element *element = elInfo->getElement();
	DegreeOfFreedom *newDOF, **oldDOF, **dof = 
	  const_cast<DegreeOfFreedom**>(element->getDOF());

	int index = 0;

	if (extendNodes) {
	  oldDOF=dof;
	  element->setDOFPtrs();
	  dof=const_cast<DegreeOfFreedom**>(element->getDOF());
	  int index=0,oldIndex=0;
	  for (int d = 0; d < dim+1; d++) {
	    geoIndex = INDEX_OF_DIM(d, dim);
	    if (nDOF[geoIndex]>0) {
	      for(int i=0;i<getGeo(geoIndex);++i) 
		dof[index++]=oldDOF[oldIndex++];
	    }
	    else {
	      if (localAdmin->getNumberOfDOFs(geoIndex)>0) 
		index+=getGeo(geoIndex);
	    }
	  }
	
	  FREE_MEMORY(oldDOF, DegreeOfFreedom*, oldNNodes);

	  TEST_EXIT_DBG(index == nNodeEl)("ERROR: Number of entered nodes %f != number of nodes %f\n",index,nNodeEl);

	}


	index=0;

	// allocate new memory at elements
	for(int d = 0; d < dim+1; d++) {
	  geoIndex = INDEX_OF_DIM(d, dim);
      
	  int numberOfDOFs = localAdmin->getNumberOfDOFs(geoIndex);
	  int numberOfPreDOFs = nDOF[geoIndex];

	  if (numberOfDOFs>0||numberOfPreDOFs>0) {

	    // for all vertices/edges/...
	    for(int i = 0; i < getGeo(geoIndex); i++, index++) {
	      std::set<DegreeOfFreedom> dofSet;
	      for(int j = 0; j < d+1; j++) {
		dofSet.insert(dof[element->getVertexOfPosition(geoIndex, i, j)][vertexAdminPreDOFs]);
	      }
	    
	      if(element->isLeaf() || parentNeedsDOF[dim][d][i]) {
		if(dofPtrMap[dofSet] == NULL) {
		  if(localAdmin->getNumberOfDOFs(geoIndex)) {
		    newDOF = GET_MEMORY(DegreeOfFreedom, numberOfPreDOFs + numberOfDOFs);
		    // copy old dofs to new memory and free old memory
		    if(dof[index]) {
		      for(int j = 0; j < numberOfPreDOFs; j++) {
			newDOF[j] = dof[index][j];
		      }
		      //	  FREE_MEMORY(dof[index], DegreeOfFreedom, numberOfPreDOFs);
		      // Do not free memory. The information has to be used to identify the part in other elements.
		      // The memory is only marked for freeing.
		      struct delmem fm;
		      fm.ptr=dof[index];
		      fm.len=numberOfPreDOFs;
		      delList.push_back(fm);
		    }
		    for(int j = 0; j < numberOfDOFs; j++) {
		      newDOF[numberOfPreDOFs + j] = localAdmin->getDOFIndex();
		    }
		    dof[index] = newDOF;
		  }
		  dofPtrMap[dofSet] = dof[index];
		} else {
		  dof[index] = dofPtrMap[dofSet];
		}
	      }
	    }
	  }
	}
	elInfo = stack.traverseNext(elInfo);
      }
  
      // now free the old dof memory:

      std::list<struct delmem>::iterator it=delList.begin();
    
      while(it!=delList.end()) {
	FREE_MEMORY((*it).ptr, DegreeOfFreedom, (*it).len);
	it++;
      }

      delList.clear();

    }
534
535
536
537
538
539
540


    admin.push_back(localAdmin);

    nDOFEl = 0;

    localAdmin->setNumberOfPreDOFs(VERTEX,nDOF[VERTEX]);
541
    nDOF[VERTEX] += localAdmin->getNumberOfDOFs(VERTEX);
542
543
    nDOFEl += getGeo(VERTEX) * nDOF[VERTEX];

544
    if (dim > 1) {
545
      localAdmin->setNumberOfPreDOFs(EDGE,nDOF[EDGE]);
546
      nDOF[EDGE] += localAdmin->getNumberOfDOFs(EDGE);
547
548
549
550
551
552
553
      nDOFEl += getGeo(EDGE) * nDOF[EDGE];
    }

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

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

556
557
    node[VERTEX] = 0;
    nNodeEl = getGeo(VERTEX);
558

559
560
    if (dim > 1) {
      node[EDGE] = nNodeEl;
561
562
      if (nDOF[EDGE] > 0) 
	nNodeEl += getGeo(EDGE);
563
564
    }

565
    if (dim == 3) {
566
      localAdmin->setNumberOfPreDOFs(FACE,nDOF[FACE]);
567
568
569
570
571
      nDOF[FACE] += localAdmin->getNumberOfDOFs(FACE);
      nDOFEl += getGeo(FACE) * nDOF[FACE];
      node[FACE] = nNodeEl;
      if (nDOF[FACE] > 0) 
	nNodeEl += getGeo(FACE);
572
573
    }

574
    node[CENTER] = nNodeEl;
575
    if (nDOF[CENTER] > 0) {
576
      nNodeEl += 1;
577
    }
578
579
580
581
  }

  void Mesh::dofCompress()
  {
582
583
    FUNCNAME("Mesh::dofCompress()");
    Flag fill_flag;
584

Thomas Witkowski's avatar
Thomas Witkowski committed
585
    for (int iadmin = 0; iadmin < static_cast<int>(admin.size()); iadmin++) {
586
587
588
      compressAdmin = admin[iadmin];

      TEST_EXIT_DBG(compressAdmin)("no admin[%d] in mesh\n", iadmin);
589
      
Thomas Witkowski's avatar
Thomas Witkowski committed
590
591
      int size = compressAdmin->getSize();
      if (size < 1) 
592
	continue;
Thomas Witkowski's avatar
Thomas Witkowski committed
593

594
595
      if (compressAdmin->getUsedDOFs() < 1)    
	continue;
Thomas Witkowski's avatar
Thomas Witkowski committed
596

597
598
      if (compressAdmin->getHoleCount() < 1)    
	continue;
Thomas Witkowski's avatar
Thomas Witkowski committed
599
  
600
601
602
603
604
605
606
607
      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;
608
      }
609
      
Thomas Witkowski's avatar
Thomas Witkowski committed
610
611
      traverse(-1, fill_flag, newDOFFct1);
      traverse(-1, fill_flag, newDOFFct2);
612
613
614
      
      newDOF.resize(0);
    }   
615
616
617
618
619
  }


  DegreeOfFreedom *Mesh::getDOF(GeoIndex position)
  {
620
    FUNCNAME("Mesh::getDOF()");
621

622
    TEST_EXIT_DBG(position >= CENTER && position <= FACE)
623
      ("unknown position %d\n", position);
624

625
626
    int ndof = getNumberOfDOFs(position);
    if (ndof <= 0) 
627
      return NULL;
628

629
    DegreeOfFreedom *dof = new DegreeOfFreedom[ndof];
630

631
632
    for (int i = 0; i < getNumberOfDOFAdmin(); i++) {
      const DOFAdmin *localAdmin = &getDOFAdmin(i);
633
      TEST_EXIT_DBG(localAdmin)("no admin[%d]\n", i);
634
635
636
637
      
      int n  = localAdmin->getNumberOfDOFs(position);
      int n0 = localAdmin->getNumberOfPreDOFs(position);
      
638
      TEST_EXIT_DBG(n + n0 <= ndof)("n=%d, n0=%d too large: ndof=%d\n", n, n0, ndof);
639
      
640
      for (int j = 0; j < n; j++)
641
642
	dof[n0 + j] = const_cast<DOFAdmin*>(localAdmin)->getDOFIndex();
    }
643
  
644
    return dof;
645
646
647
648
649
  }


  DegreeOfFreedom **Mesh::createDOFPtrs()
  {
650
    FUNCNAME("Mesh::createDOFPtrs()");
651
652

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

655
    DegreeOfFreedom **ptrs = new DegreeOfFreedom*[nNodeEl];
656
    for (int i = 0; i < nNodeEl; i++)
657
658
      ptrs[i] = NULL;

659
    return ptrs;
660
661
662
663
  }

  void Mesh::freeDOFPtrs(DegreeOfFreedom **ptrs)
  {
664
    FUNCNAME("Mesh::freeDOFPtrs()");
665

666
    TEST_EXIT_DBG(ptrs)("ptrs=NULL\n");
667
668
669
670

    if (nNodeEl <= 0)
      return;
  
671
    delete [] ptrs;
672
673
674
  }


Thomas Witkowski's avatar
Thomas Witkowski committed
675
  const DOFAdmin *Mesh::createDOFAdmin(const std::string& lname, DimVec<int> lnDOF)
676
  {
677
    FUNCNAME("Mesh::createDOFAdmin()");
678

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

Thomas Witkowski's avatar
Thomas Witkowski committed
681
682
    for (int i = 0; i < dim + 1; i++)
      localAdmin->setNumberOfDOFs(i, lnDOF[i]);
683
684
685

    addDOFAdmin(localAdmin);

686
    return localAdmin;
687
688
689
690
691
692
693
  }


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

694
695
696
697
698
699
    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];
700
      }
701
702
    }

703
    return localAdmin;
704
705
706
707
  }

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

710
    TEST_EXIT_DBG(position >= CENTER && position <= FACE)
711
      ("unknown position %d\n", position);
712

713
714
715
716
    int ndof = nDOF[position];
    if (ndof) {
      if (!dof) {
	MSG("dof = NULL, but ndof=%d\n", ndof);
717
718
	return;
      }
719
720
721
722
723
724
    } else  {
      if (dof) {
	MSG("dof != NULL, but ndof=0\n");
      }
      return;
    }
725

726
    TEST_EXIT_DBG(ndof <= MAX_DOF)
727
      ("ndof too big: ndof=%d, MAX_DOF=%d\n", ndof, MAX_DOF);
728

729
730
731
732
733
734
735
736
737
738
    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]);
    }
739

740
    delete [] dof;
741
742
743
744
745
  }

  void Mesh::freeElement(Element* el)
  {
    freeDOFPtrs(const_cast<DegreeOfFreedom**>(el->getDOF()));
Thomas Witkowski's avatar
Thomas Witkowski committed
746
    delete el;
747
748
749
750
751
752
  }


  Element* Mesh::createNewElement(Element *parent)
  {
    FUNCNAME("Mesh::createNewElement()");
753
754

    TEST_EXIT_DBG(elementPrototype)("no element prototype\n");
755
756
757

    Element *el = parent ? parent->clone() : elementPrototype->clone();
  
758
    if (!parent && elementDataPrototype) {
759
760
761
762
763
764
765
766
      el->setElementData(elementDataPrototype->clone()); 
    } else {
      el->setElementData(NULL); // must be done in ElementData::refineElementData()
    }

    return el;
  }

767

768
769
770
771
  ElInfo* Mesh::createNewElInfo()
  {
    switch(dim) {
    case 1:
Thomas Witkowski's avatar
Thomas Witkowski committed
772
      return new ElInfo1d(this);
773
774
      break;
    case 2:
Thomas Witkowski's avatar
Thomas Witkowski committed
775
      return new ElInfo2d(this);
776
777
      break;
    case 3:
Thomas Witkowski's avatar
Thomas Witkowski committed
778
      return new ElInfo3d(this);
779
780
781
782
      break;
    default:
      ERROR_EXIT("invalid dim\n");
      return NULL;
783
    }
784
785
786
787
  }

  bool Mesh::findElInfoAtPoint(const WorldVector<double>& xy,
			       ElInfo *el_info,
788
789
			       DimVec<double>& bary,
			       const MacroElement *start_mel,
790
			       const WorldVector<double> *xy0,
791
			       double *sp)
792
793
794
795
796
797
798
799
800
801
  {
    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
802
      if ((mel == NULL) || (mel->getElement()->getMesh() != this))
803
804
805
	mel = *(macroElements.begin());

    mel_info->setFillFlag(Mesh::FILL_COORDS);
806
    g_xy = &xy;
807
    g_xy0 = xy0;
808
    g_sp = sp;
809
810
811

    mel_info->fillMacroInfo(mel);

812
    int k;
813
814
815
816
817
818
819
820
821
822
    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 */
823
    bool inside = findElementAtPointRecursive(mel_info, lambda, k, el_info);
824
825
    for (int i = 0; i <= dim; i++)
      bary[i] = final_lambda[i];   
826
  
Thomas Witkowski's avatar
Thomas Witkowski committed
827
    delete mel_info;
828

829
    return inside;
830
831
832
  }

  bool Mesh::findElementAtPoint(const WorldVector<double>&  xy,
833
834
				Element **elp, 
				DimVec<double>& bary,
835
				const MacroElement *start_mel,
836
837
				const WorldVector<double> *xy0,
				double *sp)
838
  {
839
840
    ElInfo *el_info = createNewElInfo();
    int val = findElInfoAtPoint(xy, el_info, bary, start_mel, xy0, sp);
841
842
843

    *elp = el_info->getElement();

Thomas Witkowski's avatar
Thomas Witkowski committed
844
    delete el_info;
845

846
    return val;
847
848
  }

849
  bool Mesh::findElementAtPointRecursive(ElInfo *el_info,
850
					 const DimVec<double>& lambda,
851
					 int outside,
852
853
					 ElInfo* final_el_info)
  {
854
    FUNCNAME("Mesh::findElementAtPointRecursive()");
855
856
    Element *el = el_info->getElement();
    DimVec<double> c_lambda(dim, NO_INIT);
857
858
    int inside;
    int ichild, c_outside;
859
860
861
862

    if (el->isLeaf()) {
      *final_el_info = *el_info;
      if (outside < 0) {
863
864
865
866
	for (int i = 0; i <= dim; i++) {
	  final_lambda[i] = lambda[i];
	}

867
	return true;
868
869
870
871
872
873
874
875
876
877
878
879
      }  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);
880

881
	  return false;  /* ??? */
882
	} else {
883
	  return false;
884
	}
885
      }
886
887
    }

888
    ElInfo *c_el_info = createNewElInfo();
889

890
    if (dim == 1) {
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
      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 */

912
    if (dim == 2) {
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
      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 */

940
    if (dim == 3) {
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
      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])) {
964
965
966
	    for (int i = 0; i <= dim; i++) {
	      c_lambda[i] = c_lambda2[i];
	    }
967
968
969
970
	    c_outside = c_outside2;
	    *c_el_info = *c_el_info2;
	    ichild = 1 - ichild;
	  }
Thomas Witkowski's avatar
Thomas Witkowski committed
971
	  delete c_el_info2;
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
	}
	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
997
    delete c_el_info;
998

999
    return inside; 
1000
  }
For faster browsing, not all history is shown. View entire blame