Mesh.cc 44 KB
Newer Older
1
2
3
4
5
6
7
/******************************************************************************
 *
 * AMDiS - Adaptive multidimensional simulations
 *
 * Copyright (C) 2013 Dresden University of Technology. All Rights Reserved.
 * Web: https://fusionforge.zih.tu-dresden.de/projects/amdis
 *
8
 * Authors:
9
10
11
12
13
14
15
16
17
 * Simon Vey, Thomas Witkowski, Andreas Naumann, Simon Praetorius, et al.
 *
 * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
 * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
 *
 *
 * This file is part of AMDiS
 *
 * See also license.opensource.txt in the distribution.
18
 *
19
 ******************************************************************************/
20
21


Thomas Witkowski's avatar
Thomas Witkowski committed
22
23
24
25
26
27
#include <algorithm>
#include <set>
#include <map>

#include "time.h"

28
#include "io/Reader.h"
29
30
31
32
#include "io/MacroReader.h"
#include "io/MacroInfo.h"
#include "io/MacroWriter.h"

33
34
35
36
#include "AdaptStationary.h"
#include "AdaptInstationary.h"
#include "FiniteElemSpace.h"
#include "ElementData.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
37
#include "ElementDofIterator.h"
38
39
40
#include "MacroElement.h"
#include "Mesh.h"
#include "Traverse.h"
41
#include "Initfile.h"
42
43
44
45
46
47
#include "FixVec.h"
#include "DOFVector.h"
#include "CoarseningManager.h"
#include "DOFIterator.h"
#include "VertexVector.h"
#include "Projection.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
48
#include "ElInfoStack.h"
49
#include "Serializer.h"
50
#include "Lagrange.h"
51

52
using namespace std;
53

54
namespace AMDiS {
55
56

  //**************************************************************************
57
  //  flags, which information should be present in the elInfo structure
58
59
  //**************************************************************************

60
61
62
63
64
65
66
67
68
69
70
71
72
73
  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;


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

  //**************************************************************************
76
  //  flags for Mesh traversal
77
78
79
80
81
82
83
84
  //**************************************************************************

  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;
85
  const Flag Mesh::CALL_MG_LEVEL           = 0X4000L;
86
  const Flag Mesh::CALL_REVERSE_MODE       = 0X8000L;
87
88


89
  vector<DegreeOfFreedom> Mesh::dof_used;
90
  const int Mesh::MAX_DOF = 100;
91
  map<pair<DegreeOfFreedom, int>, DegreeOfFreedom*> Mesh::serializedDOFs;
92
  std::set<string> Mesh::refinedMeshNames;
93

94
95
96
  Mesh::Mesh(string aName, int dimension)
    : name(aName),
      dim(dimension),
97
98
      nVertices(0),
      nEdges(0),
99
      nLeaves(0),
100
      nElements(0),
101
      parametric(NULL),
102
      preserveCoarseDOFs(false),
103
104
      nDofEl(0),
      nDof(dimension, DEFAULT_VALUE, 0),
105
106
      nNodeEl(0),
      node(dimension, DEFAULT_VALUE, 0),
107
108
      elementPrototype(NULL),
      elementDataPrototype(NULL),
109
110
      elementIndex(-1),
      initialized(false),
111
      macroFileInfo(NULL),
112
      changeIndex(0),
113
114
      final_lambda(dimension, DEFAULT_VALUE, 0.0)
  {
115
    FUNCNAME("Mesh::Mesh()");
116

117
118
119
120
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    nParallelPreRefinements = 0;
#endif

121
122
123
    // set default element prototype
    switch(dim) {
    case 1:
Thomas Witkowski's avatar
Thomas Witkowski committed
124
      elementPrototype = new Line(this);
125
126
      break;
    case 2:
Thomas Witkowski's avatar
Thomas Witkowski committed
127
      elementPrototype = new Triangle(this);
128
129
      break;
    case 3:
Thomas Witkowski's avatar
Thomas Witkowski committed
130
      elementPrototype = new Tetrahedron(this);
131
132
133
134
135
136
137
      break;
    default:
      ERROR_EXIT("invalid dimension\n");
    }

    elementPrototype->setIndex(-1);

138
139
    elementIndex = 0;
  }
140

141

142
  Mesh::~Mesh()
143
  {
144
    deleteMeshStructure();
145

146
    if (macroFileInfo != NULL)
147
      delete macroFileInfo;
148
    if (elementPrototype)
149
      delete elementPrototype;
150
    if (elementDataPrototype)
151
      delete elementDataPrototype;
152
  }
153

154

155
  Mesh& Mesh::operator=(const Mesh& m)
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;
172
    parametric = NULL;
173
174

    preserveCoarseDOFs = m.preserveCoarseDOFs;
175
176
    nDofEl = m.nDofEl;
    nDof = m.nDof;
177
178
179
180
181
    nNodeEl = m.nNodeEl;
    node = m.node;
    elementIndex = m.elementIndex;
    initialized = m.initialized;
    final_lambda = m.final_lambda;
182

183
184
185
    /* ====================== 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
186
      admin[i] = new DOFAdmin(this);
187
      *(admin[i]) = *(m.admin[i]);
188
189
      admin[i]->setMesh(this);
    }
190

191

192
    /* ====================== Copy macro elements =================== */
193

194
    // mapIndex[i] is the index of the MacroElement element in the vector
195
    // macroElements, for which holds: element->getIndex() = i
196
    map<int, int> mapIndex;
197
198
199
200
201
202

    // We use this map for coping the DOFs of the Elements within the
    // MacroElements objects.
    Mesh::serializedDOFs.clear();

    int insertCounter = 0;
203

204
205
    macroElements.clear();

206
207
    // Go through all MacroElements of mesh m, and create for every a new
    // MacroElement in this mesh.
208
    for (deque<MacroElement*>::const_iterator it = m.macroElements.begin();
209
	 it != m.macroElements.end(); ++it, insertCounter++) {
210

211
      // Create new MacroElement.
Thomas Witkowski's avatar
Thomas Witkowski committed
212
      MacroElement *el = new MacroElement(dim);
213

214
215
216
217
218
219
220
221
222
223
      // 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.
224
      mapIndex.insert(pair<int, int>(el->getIndex(), insertCounter));
225
226
227
228
229
    }

    // Now we have to go through all the new MacroElements, and update the neighbour
    // connections.
    insertCounter = 0;
230
    for (deque<MacroElement*>::const_iterator it = m.macroElements.begin();
231
232
233
234
235
236
237
238
239
	 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.
240
	if ((*it)->getNeighbour(i)!=NULL) {
241
242
243
	macroElements[insertCounter]->
	  setNeighbour(i, macroElements[mapIndex[(*it)->getNeighbour(i)->getIndex()]]);
      }
244
      }
245
246
247
248
249
250
    }

    // Cleanup
    Mesh::serializedDOFs.clear();

    /* ================== Things will be done when required ============ */
251

252
253
    TEST_EXIT(elementDataPrototype == NULL)("TODO\n");
    TEST_EXIT(m.parametric == NULL)("TODO\n");
254
255
256
257
258
    TEST_EXIT(periodicAssociations.size() == 0)("TODO\n");

    return *this;
  }

259

260
261
262
263
264
265
266
267
268
269
270
271
  void Mesh::updateNumberOfLeaves()
  {
    nLeaves = 0;

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

272

273
  void Mesh::addMacroElement(MacroElement* me)
274
  {
275
    macroElements.push_back(me);
276
277
278
    me->setIndex(macroElements.size());
  }

279
280
281
282
283
284
285
286
287
  void Mesh::removeAllMacroElements()
  {
    // Delete all the dofs
    Element::deletedDOFs.clear();
    for (deque<MacroElement*>::const_iterator it = macroElements.begin();
	 it != macroElements.end(); ++it) {
      (*it)->getElement()->deleteElementDOFs();
    }
    Element::deletedDOFs.clear();
288

289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
    // Set all neighbors null
    for (deque<MacroElement*>::const_iterator macroIt = macroElements.begin();
	 macroIt != macroElements.end(); ++macroIt) {

      for (int i = 0; i < getGeo(NEIGH); i++)
	(*macroIt)->setNeighbour(i, NULL);

      Element *mel = (*macroIt)->getElement();
      // Delete element hierarchie
      if (!(mel->isLeaf())) {
	delete mel->getChild(0);
	delete mel->getChild(1);

	mel->child[0] = NULL;
	mel->child[1] = NULL;

305
	mel->setElementData(elementDataPrototype->clone());
306
307
      }
    }
308

309
310
311
312
    macroElements.clear();
    nLeaves = 0;
    nElements = 0;
    nVertices = 0;
313

314
315
316
    for (size_t i = 0; i < admin.size(); i++)
      {
	TEST_EXIT_DBG(admin[i]->getUsedSize() == admin[i]->getHoleCount())
317
	  ("All macro elements has been removed. But not all dofs are cleaned. (UsedSize = %d, HoleCount = %d)\n",
318
	    admin[i]->getUsedSize(), admin[i]->getHoleCount());
319

320
321
322
	admin[i]->reset();
      }
  }
323

324
  void Mesh::removeMacroElements(std::set<MacroElement*>& delMacros,
325
				 vector<const FiniteElemSpace*>& feSpaces)
326
327
328
  {
    FUNCNAME("Mesh::removeMacroElement()");

329
330
    typedef map<const DegreeOfFreedom*, std::set<MacroElement*> > DofElMap;
    typedef map<const DegreeOfFreedom*, GeoIndex> DofPosMap;
Thomas Witkowski's avatar
Thomas Witkowski committed
331

332
    TEST_EXIT(feSpaces.size() > 0)("Should not happen!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
333

334
335
336
337
    // Search for the FE space with the highest degree of polynomials. Using this
    // FE space ensures that deleting DOFs defined on it, also DOFs of lower
    // order FE spaces will be deleted correctly.
    const FiniteElemSpace *feSpace = FiniteElemSpace::getHighest(feSpaces);
338

339

340
341
    // === Determine to all DOFs in mesh the macro elements where the DOF  ===
    // === is part of.                                                     ===
342

343
344
    // Map that stores for each DOF pointer (which may have a list of DOFs)
    // all macro element indices that own this DOF.
Thomas Witkowski's avatar
Thomas Witkowski committed
345
346
    DofElMap dofsOwner;
    DofPosMap dofsPosIndex;
347

348
349
350
351
352
    ElementDofIterator elDofIter(feSpace);
    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(this, -1, Mesh::CALL_LEAF_EL);
    while (elInfo) {
      elDofIter.reset(elInfo->getElement());
Thomas Witkowski's avatar
Thomas Witkowski committed
353
      do {
354
355
	dofsOwner[elDofIter.getBaseDof()].insert(elInfo->getMacroElement());
	dofsPosIndex[elDofIter.getBaseDof()] = elDofIter.getPosIndex();
Thomas Witkowski's avatar
Thomas Witkowski committed
356
      } while (elDofIter.nextStrict());
357
358

      elInfo = stack.traverseNext(elInfo);
359
    }
360

Thomas Witkowski's avatar
Thomas Witkowski committed
361

362
363
    // === Remove macro elements from mesh macro element list. ===

364
    // Removing arbitrary elements from an deque is very slow. Therefore, we
365
366
367
    // create a new deque with all macro elements that should not be deleted. The
    // macro element deque is than replaced by the new created one.

368
369
    deque<MacroElement*> newMacroElements;
    for (deque<MacroElement*>::iterator elIter = macroElements.begin();
370
	 elIter != macroElements.end(); ++elIter) {
371
372
373
      // If the current mesh macro element should not be deleted, i.e., it is not
      // a member of the list of macro elements to be deleted, is is inserted to
      // the new macro element list.
374
      if (delMacros.find(*elIter) == delMacros.end())
375
	newMacroElements.push_back(*elIter);
376
    }
377

378
379
380
381
    // And replace the macro element list with the new one.
    macroElements.clear();
    macroElements = newMacroElements;

Thomas Witkowski's avatar
Thomas Witkowski committed
382

383
384
385
    // === For all macro elements to be deleted, delete them also to be       ===
    // === neighbours of some other macro elements. Furtheremore, delete the  ===
    // === whole element hierarchie structure of the macro element.           ===
386

387
388
    for (std::set<MacroElement*>::iterator macroIt = delMacros.begin();
	 macroIt != delMacros.end(); ++macroIt) {
Thomas Witkowski's avatar
Thomas Witkowski committed
389

390
391
      // Go through all neighbours of the macro element and remove this macro
      // element to be neighbour of some other macro element.
392
      for (int i = 0; i < getGeo(NEIGH); i++)
393
	if ((*macroIt)->getNeighbour(i))
394
	  for (int j = 0; j < getGeo(NEIGH); j++)
Thomas Witkowski's avatar
Thomas Witkowski committed
395
	    if ((*macroIt)->getNeighbour(i)->getNeighbour(j) == *macroIt)
396
	      (*macroIt)->getNeighbour(i)->setNeighbour(j, NULL);
397

398

399
      Element *mel = (*macroIt)->getElement();
400
      // Delete element hierarchie
401
402
403
      if (!(mel->isLeaf())) {
	delete mel->getChild(0);
	delete mel->getChild(1);
404

405
406
	mel->child[0] = NULL;
	mel->child[1] = NULL;
407

408
	mel->setElementData(elementDataPrototype->clone());
409
      }
410

411
      mel->delDofPtr();
412
    }
413
414


415
416
    // === Check now all the DOFs that have no owner anymore and therefore  ===
    // === have to be removed.                                              ===
417

418
    for (DofElMap::iterator dofsIt = dofsOwner.begin();
419
	 dofsIt != dofsOwner.end(); ++dofsIt) {
420

421
422
423
424
      bool deleteDof = true;

      for (std::set<MacroElement*>::iterator elIter = dofsIt->second.begin();
	   elIter != dofsIt->second.end(); ++elIter) {
425
426
	std::set<MacroElement*>::iterator mIt = delMacros.find(*elIter);
	if (mIt == delMacros.end()) {
427
428
429
430
431
432
	  deleteDof = false;
	  break;
	}
      }

      if (deleteDof)
433
434
	freeDof(const_cast<DegreeOfFreedom*>(dofsIt->first),
		dofsPosIndex[dofsIt->first]);
435
436
    }

437

438
    // === Update number of elements, vertices, etc. ===
439

440
441
    nLeaves = 0;
    nElements = 0;
442
    nVertices = 0;
443

444
445
    if (!macroElements.empty()) {
      std::set<const DegreeOfFreedom*> allVertices;
446

447
448
449
      elInfo = stack.traverseFirst(this, -1, Mesh::CALL_EVERY_EL_PREORDER);
      while (elInfo) {
	nElements++;
450

451
452
	if (elInfo->getElement()->isLeaf()) {
	  nLeaves++;
453

454
455
456
457
458
	  for (int i = 0; i < getGeo(VERTEX); i++)
	    allVertices.insert(elInfo->getElement()->getDof(i));
	}

	elInfo = stack.traverseNext(elInfo);
459
460
      }

461
462
      nVertices = allVertices.size();
    } else {
463

464
465
466
      for (size_t i = 0; i < admin.size(); i++)
      {
	TEST_EXIT_DBG(admin[i]->getUsedSize() == admin[i]->getHoleCount())
467
	  ("All macro elements has been removed. But not all dofs are cleaned. (UsedSize = %d, HoleCount = %d)\n",
468
	    admin[i]->getUsedSize(), admin[i]->getHoleCount());
469

470
471
	admin[i]->reset();
      }
472
473
474
475
476
    }
    // === Note: Although the macro elements are removed from the mesh,   ===
    // === they are not deleted from memory. The macro elements are still ===
    // === stored in macroInfo structure. They are needed, if the mesh is ===
    // === redistributed between the ranks.                               ===
477
  }
478

479

480
  void Mesh::addDOFAdmin(DOFAdmin *localAdmin)
481
  {
482
483
484
485
    FUNCNAME("Mesh::addDOFAdmin()");

    localAdmin->setMesh(this);

486
    TEST_EXIT(find(admin.begin(), admin.end(), localAdmin) == admin.end())
Thomas Witkowski's avatar
Thomas Witkowski committed
487
488
      ("admin %s is already associated to mesh %s\n",
       localAdmin->getName().c_str(), this->getName().c_str());
489
490
491

    admin.push_back(localAdmin);

492
    nDofEl = 0;
493

494
    localAdmin->setNumberOfPreDofs(VERTEX, nDof[VERTEX]);
495
496
    nDof[VERTEX] += localAdmin->getNumberOfDofs(VERTEX);
    nDofEl += getGeo(VERTEX) * nDof[VERTEX];
497

498
    if (dim > 1) {
499
      localAdmin->setNumberOfPreDofs(EDGE, nDof[EDGE]);
500
501
      nDof[EDGE] += localAdmin->getNumberOfDofs(EDGE);
      nDofEl += getGeo(EDGE) * nDof[EDGE];
502
503
    }

504
505
    localAdmin->setNumberOfPreDofs(CENTER, nDof[CENTER]);
    nDof[CENTER] += localAdmin->getNumberOfDofs(CENTER);
506
    nDofEl += nDof[CENTER];
507

508
    TEST_EXIT_DBG(nDof[VERTEX] > 0)("no vertex dofs\n");
509

510
511
    node[VERTEX] = 0;
    nNodeEl = getGeo(VERTEX);
512

513
514
    if (dim > 1) {
      node[EDGE] = nNodeEl;
515
      if (nDof[EDGE] > 0)
516
	nNodeEl += getGeo(EDGE);
517
518
    }

519
    if (dim == 3) {
520
      localAdmin->setNumberOfPreDofs(FACE, nDof[FACE]);
521
522
      nDof[FACE] += localAdmin->getNumberOfDofs(FACE);
      nDofEl += getGeo(FACE) * nDof[FACE];
523
      node[FACE] = nNodeEl;
524
      if (nDof[FACE] > 0)
525
	nNodeEl += getGeo(FACE);
526
527
    }

528
    node[CENTER] = nNodeEl;
529
    if (nDof[CENTER] > 0)
530
      nNodeEl += 1;
531
532
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
533

534
535
  void Mesh::dofCompress()
  {
536
    FUNCNAME_DBG("Mesh::dofCompress()");
537

Thomas Witkowski's avatar
Thomas Witkowski committed
538
539
    for (unsigned int iadmin = 0; iadmin < admin.size(); iadmin++) {
      DOFAdmin* compressAdmin = admin[iadmin];
540
541

      TEST_EXIT_DBG(compressAdmin)("no admin[%d] in mesh\n", iadmin);
542

Thomas Witkowski's avatar
Thomas Witkowski committed
543
      int size = compressAdmin->getSize();
544
545
546
      if (size < 1 ||
	  compressAdmin->getUsedDofs() < 1 ||
	  compressAdmin->getHoleCount() < 1)
547
	continue;
Thomas Witkowski's avatar
Thomas Witkowski committed
548

549
      vector<DegreeOfFreedom> newDofIndex(size);
550
      compressAdmin->compress(newDofIndex);
Thomas Witkowski's avatar
Thomas Witkowski committed
551

552
      Flag fill_flag = (preserveCoarseDOFs ?
Thomas Witkowski's avatar
Thomas Witkowski committed
553
			Mesh::CALL_EVERY_EL_PREORDER | Mesh::FILL_NOTHING :
554
			Mesh::CALL_LEAF_EL | Mesh::FILL_NOTHING);
555
556
557
      TraverseStack stack;
      ElInfo *elInfo = stack.traverseFirst(this, -1, fill_flag);
      while (elInfo) {
558
	elInfo->getElement()->newDofFct1(compressAdmin, newDofIndex);
559
560
561
562
563
	elInfo = stack.traverseNext(elInfo);
      }

      elInfo = stack.traverseFirst(this, -1, fill_flag);
      while (elInfo) {
564
	elInfo->getElement()->newDofFct2(compressAdmin);
565
566
	elInfo = stack.traverseNext(elInfo);
      }
567
    }
568
569
570
  }


571
  DegreeOfFreedom *Mesh::getDof(GeoIndex position)
572
  {
573
    FUNCNAME_DBG("Mesh::getDof()");
574

575
    TEST_EXIT_DBG(position >= CENTER && position <= FACE)
576
      ("unknown position %d\n", position);
577

578
    int ndof = getNumberOfDofs(position);
579
    if (ndof <= 0)
580
      return NULL;
581

582
    DegreeOfFreedom *dof = new DegreeOfFreedom[ndof];
583

584
    for (int i = 0; i < getNumberOfDOFAdmin(); i++) {
585
      const DOFAdmin *localAdmin = &getDofAdmin(i);
586
      TEST_EXIT_DBG(localAdmin)("no admin[%d]\n", i);
587

588
      int n  = localAdmin->getNumberOfDofs(position);
589
      int n0 = localAdmin->getNumberOfPreDofs(position);
590

591
592
      TEST_EXIT_DBG(n + n0 <= ndof)
	("n = %d, n0 = %d too large: ndof = %d\n", n, n0, ndof);
593

594
      for (int j = 0; j < n; j++)
595
596
	dof[n0 + j] = const_cast<DOFAdmin*>(localAdmin)->getDOFIndex();
    }
597

598
    return dof;
599
600
601
  }


Thomas Witkowski's avatar
Thomas Witkowski committed
602
  DegreeOfFreedom **Mesh::createDofPtrs()
603
604
  {
    if (nNodeEl <= 0)
605
      return NULL;
606

607
    DegreeOfFreedom **ptrs = new DegreeOfFreedom*[nNodeEl];
608
    for (int i = 0; i < nNodeEl; i++)
609
      ptrs[i] = NULL;
610

611
    return ptrs;
612
613
  }

614

615
  void Mesh::freeDofPtrs(DegreeOfFreedom **ptrs)
616
  {
617
    FUNCNAME_DBG("Mesh::freeDofPtrs()");
618

619
    TEST_EXIT_DBG(ptrs)("ptrs is NULL!\n");
620
621
622

    if (nNodeEl <= 0)
      return;
623

624
    delete [] ptrs;
625
    ptrs = NULL;
626
627
628
  }


629
  const DOFAdmin *Mesh::createDOFAdmin(string lname, DimVec<int> lnDof)
630
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
631
    DOFAdmin *localAdmin = new DOFAdmin(this, lname);
632

Thomas Witkowski's avatar
Thomas Witkowski committed
633
    for (int i = 0; i < dim + 1; i++)
634
      localAdmin->setNumberOfDofs(i, lnDof[i]);
635
636
637

    addDOFAdmin(localAdmin);

638
    return localAdmin;
639
640
641
642
643
  }


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

Thomas Witkowski's avatar
Thomas Witkowski committed
646
    for (unsigned int i = 0; i < admin.size(); i++) {
647
      if (admin[i]->getNumberOfDofs(VERTEX)) {
648
	if (!localAdmin)
649
650
651
	  localAdmin = admin[i];
	else if (admin[i]->getSize() < localAdmin->getSize())
	  localAdmin = admin[i];
652
      }
653
654
    }

655
    return localAdmin;
656
657
  }

658

659
  void Mesh::freeDof(DegreeOfFreedom* dof, GeoIndex position)
660
  {
661
    FUNCNAME_DBG("Mesh::freeDof()");
662

663
    TEST_EXIT_DBG(position >= CENTER && position <= FACE)
664
      ("unknown position %d\n", position);
665

666
    if (nDof[position]) {
667
      TEST_EXIT_DBG(dof != NULL)("dof = NULL, but ndof = %d\n", nDof[position]);
668
    } else  {
669
      TEST_EXIT_DBG(dof == NULL)("dof != NULL, but ndof = 0\n");
670
    }
671

672
673
    TEST_EXIT_DBG(nDof[position] <= MAX_DOF)
      ("ndof too big: ndof = %d, MAX_DOF = %d\n", nDof[position], MAX_DOF);
674

675
    for (unsigned int i = 0; i < admin.size(); i++) {
676
677
      int n = admin[i]->getNumberOfDofs(position);
      int n0 = admin[i]->getNumberOfPreDofs(position);
678

679
680
      TEST_EXIT_DBG(n + n0 <= nDof[position])
	("n = %d, n0 = %d too large: ndof = %d\n", n, n0, nDof[position]);
681

682
      for (int j = 0; j < n; j++)
683
	admin[i]->freeDofIndex(dof[n0 + j]);
684
    }
685
    delete [] dof;
686
687
  }

688

689
690
  void Mesh::freeElement(Element* el)
  {
691
    freeDofPtrs(const_cast<DegreeOfFreedom**>(el->getDof()));
Thomas Witkowski's avatar
Thomas Witkowski committed
692
    delete el;
693
694
695
696
697
  }


  Element* Mesh::createNewElement(Element *parent)
  {
698
    FUNCNAME_DBG("Mesh::createNewElement()");
699
700

    TEST_EXIT_DBG(elementPrototype)("no element prototype\n");
701
702

    Element *el = parent ? parent->clone() : elementPrototype->clone();
703

Thomas Witkowski's avatar
Thomas Witkowski committed
704
    if (!parent && elementDataPrototype)
705
      el->setElementData(elementDataPrototype->clone());
Thomas Witkowski's avatar
Thomas Witkowski committed
706
    else
707
      el->setElementData(NULL); // must be done in ElementData::refineElementData()
708

709
710
711
    return el;
  }

712

713
714
  ElInfo* Mesh::createNewElInfo()
  {
715
716
    FUNCNAME("Mesh::createNewElInfo()");

717
    switch (dim) {
718
    case 1:
Thomas Witkowski's avatar
Thomas Witkowski committed
719
      return new ElInfo1d(this);
720
721
      break;
    case 2:
Thomas Witkowski's avatar
Thomas Witkowski committed
722
      return new ElInfo2d(this);
723
724
      break;
    case 3:
Thomas Witkowski's avatar
Thomas Witkowski committed
725
      return new ElInfo3d(this);
726
727
      break;
    default:
728
      ERROR_EXIT("invalid dim [%d]\n",dim);
729
      return NULL;
730
    }
731
732
  }

733

734
735
  bool Mesh::findElInfoAtPoint(const WorldVector<double>& xy,
			       ElInfo *el_info,
736
737
			       DimVec<double>& bary,
			       const MacroElement *start_mel,
738
			       const WorldVector<double> *xy0,
739
			       double *sp)
740
  {
741
    static const MacroElement *mel = NULL;
742
    DimVec<double> lambda(dim, NO_INIT);
Thomas Witkowski's avatar
Thomas Witkowski committed
743
    ElInfo *mel_info = createNewElInfo();
744

745
    if (start_mel != NULL)
746
747
      mel = start_mel;
    else
748
      if (mel == NULL || mel->getElement()->getMesh() != this)
749
750
751
	mel = *(macroElements.begin());

    mel_info->setFillFlag(Mesh::FILL_COORDS);
752
    g_xy = &xy;
753
    g_xy0 = xy0;
754
    g_sp = sp;
755
756
757

    mel_info->fillMacroInfo(mel);

758
    // We have the care about not to visite a macro element twice. In this case the
759
    // function would end up in an infinite loop. If a macro element is visited a
760
761
762
    // second time, what can happen with periodic boundary conditions, the point is
    // not within the mesh!
    std::set<int> macrosVisited;
763
    std::stack<MacroElement*> active;
764

765
//     macrosVisited.insert(mel->getIndex());
766

767
    int k;
768
    while ((k = mel_info->worldToCoord(xy, &lambda)) >= 0) {
769
      macrosVisited.insert(mel->getIndex());
770

771
      if (mel->getNeighbour(k) && !macrosVisited.count(mel->getNeighbour(k)->getIndex())) {
772
	// look for next macro-element in the direction of the coordinates xy
773
	mel = mel->getNeighbour(k);
774
775
	mel_info->fillMacroInfo(mel);
	continue;
776
      } else {
777
	// consider all neighbors of the current macro-element to visit next
778
	for (int i = 0; i < dim + 1; ++i) {
779
	  if (i !=  k && mel->getNeighbour(i) && !macrosVisited.count(mel->getNeighbour(i)->getIndex()))
780
	    active.push(mel->getNeighbour(i));
781
	}
782
	// if all neighbors are visited already
783
	if (active.empty()) {
784
785
	  if (macrosVisited.size() == static_cast<size_t>(getNumberOfMacros())) {
	    // if all macro-elements are visited -> no element found!
786
787
788
	    delete mel_info;
	    return false;
	  } else {
789
	    // go to an arbitrary macro-element to continue the search
790
791
792
	    deque<MacroElement*>::iterator it;
	    bool found = false;
	    for (it = firstMacroElement(); it != endOfMacroElements(); it++) {
793
	      if (!macrosVisited.count((*it)->getIndex())) {
794
795
796
797
		active.push(*it);
		found = true;
	      }
	    }
798
799
800
801
	    if (!found) {
	      delete mel_info;
	      return false;
	    }
802
	  }
803
	}
804
805
	mel = active.top();
	active.pop();
806
807

	mel_info->fillMacroInfo(mel);
808
809
810
811
      }
    }

    /* now, descend in tree to find leaf element at point */
812
    bool inside = findElementAtPointRecursive(mel_info, lambda, k, el_info);
813
    for (int i = 0; i <= dim; i++)
814
815
      bary[i] = final_lambda[i];

Thomas Witkowski's avatar
Thomas Witkowski committed
816
    delete mel_info;
817

818
    return inside;
819
820
  }

821

822
  bool Mesh::findElementAtPoint(const WorldVector<double>&  xy,
823
				Element **elp,
824
				DimVec<double>& bary,
825
				const MacroElement *start_mel,
826
827
				const WorldVector<double> *xy0,
				double *sp)
828
  {
829
830
    ElInfo *el_info = createNewElInfo();
    int val = findElInfoAtPoint(xy, el_info, bary, start_mel, xy0, sp);
831
832
833

    *elp = el_info->getElement();

Thomas Witkowski's avatar
Thomas Witkowski committed
834
    delete el_info;
835