Liebe Gitlab-Nutzer, lieber Gitlab-Nutzer,
es ist nun möglich sich mittels des ZIH-Logins/LDAP an unserem Dienst anzumelden. Die Konten der externen Nutzer:innen sind ü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. The accounts of external users 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;
283
#ifndef NDEBUG
284 285
      bool cont2;
#endif
286
      if (structure1->isLeafElement() == structure2->isLeafElement()) {
287
	cont1 = structure1->nextElement(result);
288
#ifndef NDEBUG
289
	cont2 = structure2->nextElement();
290
#endif
291
      } else {
292
	if (structure1->isLeafElement()) {
293
	  cont1 = structure1->nextElement();
294
#ifndef NDEBUG
295
	  cont2 = structure2->skipBranch(result);
296
#endif
297 298
	} else {
	  cont1 = structure1->skipBranch(result);
299
#ifndef NDEBUG
300
	  cont2 = structure2->nextElement();
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 507 508

      elInfo = stack.traverseNext(elInfo);
    }
  }


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);
522
    unsigned int counter = 0;
523 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
  }

572
}