Liebe Gitlab-Nutzer, lieber Gitlab-Nutzer, es ist nun möglich sich mittels des ZIH-Logins/LDAP an unserem Dienst anzumelden. Ein Anmelden über dieses erzeugt ein neues Konto. Das alte Konto ist über den Reiter "Standard" erreichbar. Die Administratoren

Dear Gitlab user, it is now possible to log in to our service using the ZIH login/LDAP. Logging in via this will create a new account. The old account can be accessed via the "Standard" tab. The administrators

MeshStructure.cc 14.3 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


22 23 24 25 26
#include "Debug.h"
#include "DOFVector.h"
#include "Element.h"
#include "ElementDofIterator.h"
#include "ElInfo.h"
27 28 29 30
#include "MeshStructure.h"
#include "MeshStructure_ED.h"
#include "Mesh.h"
#include "RefinementManager.h"
31 32
#include "Traverse.h"

33

34 35
using namespace std;

36 37
namespace AMDiS {

38
  const int MeshStructure::structureSize = 64;
39

40

41
  void MeshStructure::insertElement(bool isLeaf)
42
  {
43
    // overflow? -> next index
44
    if (pos >= structureSize) {
45 46 47
      code.push_back(currentCode);
      pos = 0;
      currentCode = 0;
48 49 50
    }

    // insert element in binary code
51
    if (!isLeaf) {
52
      uint64_t one = 1;
53
      currentCode += (one << pos);
54
    }
55

56 57
    pos++;
    nElements++;
58 59
  }

60

61 62
  void MeshStructure::clear()
  {
63 64 65 66 67
    currentCode = 0;
    code.resize(0);
    pos = 0;
    nElements = 0;
    currentElement = 0;
68 69
  }

70

71
  void MeshStructure::init(Mesh *mesh, int macroElIndex)
72 73
  {
    clear();
74

75
    TraverseStack stack;
76

Thomas Witkowski's avatar
BLUB  
Thomas Witkowski committed
77 78 79 80
    ElInfo *elInfo;
    if (macroElIndex == -1)
      elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_EVERY_EL_PREORDER);
    else
81
      elInfo = stack.traverseFirstOneMacro(mesh, macroElIndex, -1,
Thomas Witkowski's avatar
BLUB  
Thomas Witkowski committed
82 83 84 85
					   Mesh::CALL_EVERY_EL_PREORDER);

    while (elInfo) {
      insertElement(elInfo->getElement()->isLeaf());
86 87
      elInfo = stack.traverseNext(elInfo);
    }
88

89 90 91
    commit();
  }

92

93
  void MeshStructure::init(BoundaryObject &bound, Element* element)
94 95 96
  {
    FUNCNAME("MeshStructure::init()");

97
    Element* el = (element == NULL) ? bound.el : element;
98

99 100
    TEST_EXIT_DBG(el)("No element!\n");

101 102
    clear();

103 104
    int s1 = el->getSubObjOfChild(0, bound.subObj, bound.ithObj, bound.elType);
    int s2 = el->getSubObjOfChild(1, bound.subObj, bound.ithObj, bound.elType);
105

106
    TEST_EXIT(s1 != -1 || s2 != -1)("This should not happen!\n");
107 108 109 110 111 112 113

    if (debugMode) {
      MSG("addAlondSide(%d, %d, %d, %d)\n",
	  bound.elIndex, bound.ithObj, bound.elType, bound.reverseMode);
      MSG("Element is leaf: %d\n", el->isLeaf());
      MSG("s1 = %d    s2 = %d\n", s1, s2);
    }
114

115 116 117 118 119 120
    /*    switch (bound.subObj)
    {
      case EDGE:
    */
	if (!el->isLeaf()) {
	  if (s1 == -1)
121
	    addAlongSide(el->getSecondChild(), bound.subObj, s2,
122 123
			el->getChildType(bound.elType), bound.reverseMode);
	  else if (s2 == -1)
124
	    addAlongSide(el->getFirstChild(), bound.subObj, s1,
125 126 127 128 129 130
			el->getChildType(bound.elType), bound.reverseMode);
	  else
	    addAlongSide(el, bound.subObj, bound.ithObj, bound.elType, bound.reverseMode);
	}
	  /*	break;
      case FACE:
131
	addAlongSide(el, bound.subObj, bound.ithObj, bound.elType, bound.reverseMode);
132 133 134
	break;
      default:
	ERROR_EXIT("What is this?\n");
135
    }
136
	  */
137 138

    commit();
139 140
  }

141

142
  void MeshStructure::addAlongSide(Element *el, GeoIndex subObj, int ithObj,
143
				   int elType, bool reverseOrder)
144
  {
145 146 147 148
    FUNCNAME("MeshStructure::addAlongSide()");

    if (debugMode) {
      MSG("addAlondSide(%d, %d, %d, %d)\n",
149
	  el->getIndex(), ithObj, elType, reverseOrder);
150 151
      MSG("Element is leaf: %d\n", el->isLeaf());
    }
152

153
    insertElement(el->isLeaf());
154

155
    if (!el->isLeaf()) {
156 157
      int s1 = el->getSubObjOfChild(0, subObj, ithObj, elType);
      int s2 = el->getSubObjOfChild(1, subObj, ithObj, elType);
158

159
      if (debugMode) {
160
	MSG("Child index %d  %d\n",
161 162 163 164 165 166
	    el->getFirstChild()->getIndex(),
	    el->getSecondChild()->getIndex());
	MSG("s1 = %d    s2 = %d\n", s1, s2);
	MSG("   \n");
      }

167
      if (!reverseOrder) {
168 169
	if (s1 != -1)
	  addAlongSide(el->getFirstChild(), subObj, s1,
170
		       el->getChildType(elType), reverseOrder);
171
	if (s2 != -1)
172
	  addAlongSide(el->getSecondChild(), subObj, s2,
173
		       el->getChildType(elType), reverseOrder);
174 175
      } else {
	if (s2 != -1)
176
	  addAlongSide(el->getSecondChild(), subObj, s2,
177
		       el->getChildType(elType), reverseOrder);
178 179
	if (s1 != -1)
	  addAlongSide(el->getFirstChild(), subObj, s1,
180
		       el->getChildType(elType), reverseOrder);
181
      }
182 183 184
    }
  }

185

186
  void MeshStructure::reset()
187 188 189 190
  {
    currentIndex = 0;
    pos = 0;
    currentElement = 0;
191 192 193

    if (code.size() > 0)
      currentCode = code[0];
194
    else
195
      currentCode = 0;
196 197
  }

198

199 200
  bool MeshStructure::nextElement(MeshStructure *insert)
  {
201
    FUNCNAME_DBG("MeshStructure::nextElement()");
202 203

    if (insert)
204 205
      insert->insertElement(isLeafElement());

206 207
    pos++;
    currentElement++;
208

209
    if (currentElement >= nElements)
210
      return false;
211

212
    if (pos >= structureSize) {
213 214 215 216 217
      currentIndex++;
      TEST_EXIT_DBG(currentIndex < static_cast<int>(code.size()))
	("End of structure reached!\n");
      pos = 0;
      currentCode = code[currentIndex];
218
    } else {
219
      currentCode >>= 1;
220
    }
221

222 223 224
    return true;
  }

225

226 227 228 229 230 231 232 233 234 235 236
  int MeshStructure::lookAhead(unsigned int n)
  {
    int returnValue = 0;

    int tmp_pos = pos;
    int tmp_currentElement = currentElement;
    int tmp_currentIndex = currentIndex;
    uint64_t tmp_currentCode = currentCode;

    for (unsigned int i = 0; i < n; i++) {
      if (nextElement() == false) {
237
	returnValue = -1;
238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253
	break;
      }
    }

    if (returnValue != -1)
      returnValue = static_cast<int>(!isLeafElement());

    pos = tmp_pos;
    currentElement = tmp_currentElement;
    currentIndex = tmp_currentIndex;
    currentCode = tmp_currentCode;

    return returnValue;
  }


254 255
  bool MeshStructure::skipBranch(MeshStructure *insert)
  {
256
    FUNCNAME_DBG("MeshStructure::skipBranch()");
257

258
    if (isLeafElement()) {
259 260 261 262
      return nextElement(insert);
    } else {
      bool cont = nextElement(insert);
      cont = skipBranch(insert); // left branch
263
      TEST_EXIT_DBG(cont)("Invalid structure!\n");
264 265 266 267 268
      cont = skipBranch(insert); // righ branch
      return cont;
    }
  }

269

270 271 272 273
  void MeshStructure::merge(MeshStructure *structure1,
			    MeshStructure *structure2,
			    MeshStructure *result)
  {
274
    FUNCNAME_DBG("MeshStructure::merge()");
275

276 277 278 279 280
    result->clear();
    structure1->reset();
    structure2->reset();

    bool cont = true;
281
    while (cont) {
282
      bool cont1;
Praetorius, Simon's avatar
Praetorius, Simon committed
283
// #ifndef NDEBUG
284
      bool cont2;
Praetorius, Simon's avatar
Praetorius, Simon committed
285
// #endif
286
      if (structure1->isLeafElement() == structure2->isLeafElement()) {
287
	cont1 = structure1->nextElement(result);
Praetorius, Simon's avatar
Praetorius, Simon committed
288
// #ifndef NDEBUG
289
	cont2 = structure2->nextElement();
Praetorius, Simon's avatar
Praetorius, Simon committed
290
// #endif
291
      } else {
292
	if (structure1->isLeafElement()) {
293
	  cont1 = structure1->nextElement();
Praetorius, Simon's avatar
Praetorius, Simon committed
294
// #ifndef NDEBUG
295
	  cont2 = structure2->skipBranch(result);
Praetorius, Simon's avatar
Praetorius, Simon committed
296
// #endif
297 298
	} else {
	  cont1 = structure1->skipBranch(result);
Praetorius, Simon's avatar
Praetorius, Simon committed
299
// #ifndef NDEBUG
300
	  cont2 = structure2->nextElement();
Praetorius, Simon's avatar
Praetorius, Simon committed
301
// #endif
302 303
	}
      }
304
      TEST_EXIT_DBG(cont1 == cont2)("Structures don't match!\n");
305 306 307 308 309 310
      cont = cont1;
    }

    result->commit();
  }

311

312 313
  void MeshStructure::fitMeshToStructure(Mesh *mesh,
					 RefinementManager *manager,
314
					 bool debugMode,
Thomas Witkowski's avatar
Thomas Witkowski committed
315
					 int macroElIndex,
316
					 bool ignoreFinerMesh)
317 318 319 320 321 322 323
  {
    FUNCNAME("MeshStructure::fitMeshToStructure()");

    bool cont = true;

    // decorate leaf data
    reset();
324
    TraverseStack stack;
325
    ElInfo *elInfo = NULL;
326 327 328 329
    if (macroElIndex == -1)
      elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_EVERY_EL_PREORDER);
    else
      elInfo = stack.traverseFirstOneMacro(mesh, macroElIndex, -1, Mesh::CALL_EVERY_EL_PREORDER);
330

331

332
    while (elInfo) {
333 334
      Element *element = elInfo->getElement();

335 336
      TEST_EXIT(cont)("unexpected structure code end!\n");

337
      if (isLeafElement()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
338 339 340 341 342 343
	if (ignoreFinerMesh && !element->isLeaf()) {
	  int level = elInfo->getLevel();
	  while (elInfo && level >= elInfo->getLevel())
	    elInfo = stack.traverseNext(elInfo);
	} else {
	  TEST_EXIT(element->isLeaf())
344
	    ("Mesh is finer than strucutre code! (Element index: %d Macro element index: %d)\n",
345
	     element->getIndex(), elInfo->getMacroElement()->getIndex());
Thomas Witkowski's avatar
Thomas Witkowski committed
346
	}
347
      }
348

349 350
      TEST_EXIT_DBG(element)("Should not happen!\n");

351
      if (element->isLeaf() && !isLeafElement()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
352
	MeshStructure *structure = new MeshStructure();
353 354 355
	cont = skipBranch(structure);
	structure->commit();

356 357 358
	MeshStructure_ED *elData = new MeshStructure_ED(element->getElementData());
	elData->setStructure(structure);
	element->setElementData(elData);
359 360 361 362
      } else {
	cont = nextElement();
      }

Thomas Witkowski's avatar
Thomas Witkowski committed
363 364
      if (elInfo)
	elInfo = stack.traverseNext(elInfo);
365 366 367
    }

    // refine mesh
368
    bool finished = true;
369

370
    do {
371
      finished = true;
372 373 374
      if (macroElIndex == -1)
	elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
      else
375
	elInfo = stack.traverseFirstOneMacro(mesh, macroElIndex, -1, Mesh::CALL_LEAF_EL);
376
      while (elInfo) {
377
	Element *element = elInfo->getElement();
378
	if (element->getElementData(MESH_STRUCTURE) != NULL) {
379 380 381 382 383 384 385
	  element->setMark(1);
	  finished = false;
	} else {
	  element->setMark(0);
	}
	elInfo = stack.traverseNext(elInfo);
      }
386 387

      if (!finished) {
388
#ifndef NDEBUG
389 390
	int oldMeshIndex = mesh->getChangeIndex();
#endif
391

392 393 394 395
	if (macroElIndex == -1)
	  manager->refineMesh(mesh);
	else
	  manager->refineMacroElement(mesh, macroElIndex);
396

397
#ifndef NDEBUG
398 399 400 401
	TEST_EXIT(oldMeshIndex != mesh->getChangeIndex())
	  ("Mesh has not been changed by refinement procedure!\n");
#endif
      }
402
    } while (!finished);
403 404
  }

405

Thomas Witkowski's avatar
Thomas Witkowski committed
406
  string MeshStructure::toStr(bool resetCode)
407 408
  {
    std::stringstream oss;
409

410 411
    if (empty()) {
      oss << "-" << std::endl;
412
    }	else {
413 414 415 416 417 418 419 420 421
      if (resetCode)
	reset();

      bool cont = true;
      while (cont) {
	if (isLeafElement())
	  oss << "0";
	else
	  oss << "1";
422

423 424 425
	cont = nextElement();
      }
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
426 427 428 429 430 431 432 433

    return oss.str();
  }


  void MeshStructure::print(bool resetCode)
  {
    FUNCNAME("MeshStructure::print()");
434

Thomas Witkowski's avatar
Thomas Witkowski committed
435
    string str = toStr(resetCode);
436

Thomas Witkowski's avatar
Thomas Witkowski committed
437 438
    if (str.length() < 255) {
      MSG("Mesh structure code: %s\n", str.c_str());
439 440
    } else {
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
Thomas Witkowski's avatar
Thomas Witkowski committed
441
      std::cout << "[" << MPI::COMM_WORLD.Get_rank() << "]                Mesh structure code: " << str << "\n";
442
#else
Thomas Witkowski's avatar
Thomas Witkowski committed
443
      std::cout << "                Mesh structure code: " << str << "\n";
444 445
#endif
    }
446 447 448
  }


449 450 451 452
  bool MeshStructure::compare(MeshStructure &other)
  {
    return (other.getCode() == code);
  }
453 454


455
  void MeshStructure::getMeshStructureValues(int macroElIndex,
456
					     const DOFVector<double>* vec,
457 458
					     std::vector<double>& values,
					     bool withElIndex)
459
  {
460
    FUNCNAME_DBG("MeshStructure::getMeshStructureValues()");
461

462
    TEST_EXIT_DBG(vec)("No DOFVector defined!\n");
463

464 465
    const FiniteElemSpace *feSpace = vec->getFeSpace();
    Mesh *mesh = feSpace->getMesh();
466
    int nVertexPreDofs = feSpace->getAdmin()->getNumberOfPreDofs(VERTEX);
467
    bool feSpaceHasNonVertexDofs = (feSpace->getBasisFcts()->getDegree() > 1);
468
    values.clear();
469

470
    // In debug mode we add the macro element index to the value code.
471 472
    if (withElIndex)
      values.push_back(static_cast<double>(macroElIndex));
473

474
    ElementDofIterator elDofIter(feSpace);
475
    TraverseStack stack;
476
    ElInfo *elInfo = stack.traverseFirstOneMacro(mesh, macroElIndex, -1,
477 478
						 Mesh::CALL_EVERY_EL_PREORDER);
    while (elInfo) {
479
      // For the macro element the mesh structure code stores all vertex values.
480
      if (elInfo->getLevel() == 0)
481
	for (int i = 0; i < mesh->getGeo(VERTEX); i++)
482
	  values.push_back((*vec)[elInfo->getElement()->getDof(i, nVertexPreDofs)]);
483

484 485 486 487
      if (!elInfo->getElement()->isLeaf()) {
	// If no leaf element store the value of the "new" DOF that is created
	// by bisectioning of this element.

488
	DegreeOfFreedom dof0 =
489
	  elInfo->getElement()->getChild(0)->getDof(mesh->getDim(), nVertexPreDofs);
490 491 492 493 494 495 496 497
	values.push_back((*vec)[dof0]);
      } else {
	// If leaf element store all non vertex values of this element, thus
	// only relevant for higher order basis functions.

	if (feSpaceHasNonVertexDofs) {
	  elDofIter.reset(elInfo->getElement());
	  do {
498
	    if (elDofIter.getPosIndex() != VERTEX)
499 500 501 502
	      values.push_back((*vec)[elDofIter.getDof()]);
	  } while (elDofIter.next());
	}
      }
503 504 505 506

      elInfo = stack.traverseNext(elInfo);
    }
  }
Praetorius, Simon's avatar
Praetorius, Simon committed
507
  
508

509
  void MeshStructure::setMeshStructureValues(int macroElIndex,
510
					     DOFVector<double>* vec,
511 512
					     const std::vector<double>& values,
					     bool withElIndex)
513
  {
514
    FUNCNAME_DBG("MeshStructure::setMeshStructureValues()");
515

516
    TEST_EXIT_DBG(vec)("No DOFVector defined!\n");
517

518 519 520
    const FiniteElemSpace *feSpace = vec->getFeSpace();
    Mesh *mesh = feSpace->getMesh();
    bool feSpaceHasNonVertexDofs = (feSpace->getBasisFcts()->getDegree() > 1);
521
    int nVertexPreDofs = feSpace->getAdmin()->getNumberOfPreDofs(VERTEX);
Praetorius, Simon's avatar
Praetorius, Simon committed
522 523
    
    std::size_t counter = 0;
524 525 526 527 528 529
    if (withElIndex) {
      TEST_EXIT(static_cast<int>(values[0]) == macroElIndex)
	("Value structure code was created for macro element index %d, but should be set to macro element index %d\n",
	static_cast<int>(values[0]), macroElIndex);
      counter++;
    }
530

531 532 533 534
    TEST_EXIT_DBG(static_cast<int>(values.size()) >= mesh->getGeo(VERTEX))
      ("Should not happen!\n");

    ElementDofIterator elDofIter(feSpace);
535
    TraverseStack stack;
536
    ElInfo *elInfo = stack.traverseFirstOneMacro(mesh, macroElIndex, -1,
537 538
						 Mesh::CALL_EVERY_EL_PREORDER);
    while (elInfo) {
539
      // For the macro element all vertex nodes are set first.
540
      if (elInfo->getLevel() == 0)
541
	for (int i = 0; i < mesh->getGeo(VERTEX); i++)
542
	  (*vec)[elInfo->getElement()->getDof(i, nVertexPreDofs)] =
543
	    values[counter++];
544

545
      if (!elInfo->getElement()->isLeaf()) {
546 547
	// If no leaf element set the value of the "new" DOF that is created
	// by bisectioning of this element.
548
	TEST_EXIT_DBG(counter < values.size())("Should not happen!\n");
549

550
	(*vec)[elInfo->getElement()->getChild(0)->getDof(mesh->getDim(), nVertexPreDofs)] =
551
	  values[counter++];
552 553 554 555 556 557 558
      } else {
	// On leaf elements set all non vertex values (thus DOFs of higher order
	// basis functions).

	if (feSpaceHasNonVertexDofs) {
	  elDofIter.reset(elInfo->getElement());
	  do {
559
	    if (elDofIter.getPosIndex() != VERTEX)
560 561 562
	      (*vec)[elDofIter.getDof()] = values[counter++];
	  } while (elDofIter.next());
	}
563 564 565 566
      }

      elInfo = stack.traverseNext(elInfo);
    }
567

568 569
    TEST_EXIT_DBG(values.size() == counter)
      ("Should not happen! values size %d, counter %d\n", values.size(), counter);
570
  }
571
}