Liebe Gitlab-Nutzerin, 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 10.3 KB
Newer Older
1 2 3 4 5 6 7
#include "MeshStructure.h"
#include "MeshStructure_ED.h"
#include "Mesh.h"
#include "Element.h"
#include "Traverse.h"
#include "ElInfo.h"
#include "RefinementManager.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
8
#include "Debug.h"
9
#include "DOFVector.h"
10 11 12

namespace AMDiS {

13
  const int MeshStructure::structureSize = 64;
14

15

16 17
  void MeshStructure::insertElement(bool isLeaf) 
  {
18
    // overflow? -> next index
19
    if (pos >= structureSize) {
20 21 22
      code.push_back(currentCode);
      pos = 0;
      currentCode = 0;
23 24 25
    }

    // insert element in binary code
26
    if (!isLeaf) {
27
      uint64_t one = 1;
28
      currentCode += (one << pos);
29 30
    } 

31 32
    pos++;
    nElements++;
33 34
  }

35

36 37
  void MeshStructure::clear()
  {
38 39 40 41 42
    currentCode = 0;
    code.resize(0);
    pos = 0;
    nElements = 0;
    currentElement = 0;
43 44
  }

45

46
  void MeshStructure::init(Mesh *mesh, int macroElIndex) 
47 48
  {
    clear();
49

50 51
    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_EVERY_EL_PREORDER);
52
    while (elInfo) {
53 54 55 56
      if (macroElIndex == -1 ||
	  elInfo->getMacroElement()->getIndex() == macroElIndex)
	insertElement(elInfo->getElement()->isLeaf());

57 58 59 60 61 62
      elInfo = stack.traverseNext(elInfo);
    }
  
    commit();
  }

63 64 65 66 67 68 69

  void MeshStructure::init(BoundaryObject &bound)
  {
    FUNCNAME("MeshStructure::init()");

    Element *el = bound.el;

70 71
    TEST_EXIT_DBG(el)("No element!\n");

72 73
    clear();

74 75
    int s1 = el->getSubObjOfChild(0, bound.subObj, bound.ithObj, bound.elType);
    int s2 = el->getSubObjOfChild(1, bound.subObj, bound.ithObj, bound.elType);
76 77
    
    TEST_EXIT(s1 != -1 || s2 != -1)("This should not happen!\n");
78 79 80 81 82 83 84

    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);
    }
85 86 87
    
    if (!el->isLeaf()) {
      if (s1 == -1)
88
	addAlongSide(el->getSecondChild(), bound.subObj, s2, 
89 90
		     el->getChildType(bound.elType), bound.reverseMode);
      else if (s2 == -1)
91
	addAlongSide(el->getFirstChild(), bound.subObj, s1, 
92 93
		     el->getChildType(bound.elType), bound.reverseMode);
      else
94
	addAlongSide(el, bound.subObj, bound.ithObj, bound.elType, bound.reverseMode);
95
    }
96 97 98 99

    commit();    
  }

100 101 102

  void MeshStructure::addAlongSide(Element *el, GeoIndex subObj, int ithObj, 
				   int elType, bool reverseOrder)
103
  {
104 105 106 107
    FUNCNAME("MeshStructure::addAlongSide()");

    if (debugMode) {
      MSG("addAlondSide(%d, %d, %d, %d)\n",
108
	  el->getIndex(), ithObj, elType, reverseOrder);
109 110 111
      MSG("Element is leaf: %d\n", el->isLeaf());
    }
    
112 113 114
    insertElement(el->isLeaf());
    
    if (!el->isLeaf()) {
115 116
      int s1 = el->getSubObjOfChild(0, subObj, ithObj, elType);
      int s2 = el->getSubObjOfChild(1, subObj, ithObj, elType);
117

118 119 120 121 122 123 124 125
      if (debugMode) {
	MSG("Child index %d  %d\n", 
	    el->getFirstChild()->getIndex(),
	    el->getSecondChild()->getIndex());
	MSG("s1 = %d    s2 = %d\n", s1, s2);
	MSG("   \n");
      }

126 127
      if (!reverseOrder) {
	if (s1 != -1) 
128
	  addAlongSide(el->getFirstChild(), subObj, s1, el->getChildType(elType), reverseOrder);
129
	if (s2 != -1)
130
	  addAlongSide(el->getSecondChild(), subObj, s2, el->getChildType(elType), reverseOrder);
131 132
      } else {
	if (s2 != -1)
133
	  addAlongSide(el->getSecondChild(), subObj, s2, el->getChildType(elType), reverseOrder);
134
	if (s1 != -1) 
135
	  addAlongSide(el->getFirstChild(), subObj, s1, el->getChildType(elType), reverseOrder);
136
      }
137 138 139
    }
  }

140

141 142 143 144 145
  void MeshStructure::reset() 
  {
    currentIndex = 0;
    pos = 0;
    currentElement = 0;
146 147 148 149 150

    if (code.size() > 0)
      currentCode = code[0];
    else 
      currentCode = 0;
151 152
  }

153

154 155
  bool MeshStructure::nextElement(MeshStructure *insert)
  {
156 157 158
    FUNCNAME("MeshStructure::nextElement()");

    if (insert)
159 160
      insert->insertElement(isLeafElement());

161 162
    pos++;
    currentElement++;
163

164 165
    if (currentElement >= nElements) 
      return false;
166

167
    if (pos >= structureSize) {
168 169 170 171 172
      currentIndex++;
      TEST_EXIT_DBG(currentIndex < static_cast<int>(code.size()))
	("End of structure reached!\n");
      pos = 0;
      currentCode = code[currentIndex];
173
    } else {
174
      currentCode >>= 1;
175
    }
176

177 178 179
    return true;
  }

180

181 182
  bool MeshStructure::skipBranch(MeshStructure *insert)
  {
183 184
    FUNCNAME("MeshStructure::skipBranch()");

185
    if (isLeafElement()) {
186 187 188 189
      return nextElement(insert);
    } else {
      bool cont = nextElement(insert);
      cont = skipBranch(insert); // left branch
190
      TEST_EXIT_DBG(cont)("Invalid structure!\n");
191 192 193 194 195
      cont = skipBranch(insert); // righ branch
      return cont;
    }
  }

196

197 198 199 200
  void MeshStructure::merge(MeshStructure *structure1,
			    MeshStructure *structure2,
			    MeshStructure *result)
  {
201 202
    FUNCNAME("MeshStructure::merge()");

203 204 205 206 207
    result->clear();
    structure1->reset();
    structure2->reset();

    bool cont = true;
208
    while (cont) {
209
      bool cont1, cont2;
210
      if (structure1->isLeafElement() == structure2->isLeafElement()) {
211 212 213
	cont1 = structure1->nextElement(result);
	cont2 = structure2->nextElement();
      } else {
214
	if (structure1->isLeafElement()) {
215 216 217 218 219 220 221
	  cont1 = structure1->nextElement();
	  cont2 = structure2->skipBranch(result);
	} else {
	  cont1 = structure1->skipBranch(result);
	  cont2 = structure2->nextElement();
	}
      }
222
      TEST_EXIT_DBG(cont1 == cont2)("Structures don't match!\n");
223 224 225 226 227 228
      cont = cont1;
    }

    result->commit();
  }

229

230 231
  void MeshStructure::fitMeshToStructure(Mesh *mesh,
					 RefinementManager *manager,
232
					 bool debugMode,
Thomas Witkowski's avatar
Thomas Witkowski committed
233 234
					 int macroElIndex,
					 bool ignoreFinerMesh) 
235 236 237 238 239 240 241
  {
    FUNCNAME("MeshStructure::fitMeshToStructure()");

    bool cont = true;

    // decorate leaf data
    reset();
242
    TraverseStack stack;
243 244 245 246 247
    ElInfo *elInfo = NULL;
    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);
248

249

250
    while (elInfo) {
251 252
      Element *element = elInfo->getElement();

253 254
      TEST_EXIT(cont)("unexpected structure code end!\n");

255
      if (isLeafElement()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
256 257 258 259 260 261
	if (ignoreFinerMesh && !element->isLeaf()) {
	  int level = elInfo->getLevel();
	  while (elInfo && level >= elInfo->getLevel())
	    elInfo = stack.traverseNext(elInfo);
	} else {
	  TEST_EXIT(element->isLeaf())
262 263
	    ("Mesh is finer than strucutre code! (Element index: %d Macro element index: %d)\n", 
	     element->getIndex(), elInfo->getMacroElement()->getIndex());
Thomas Witkowski's avatar
Thomas Witkowski committed
264 265
	}
      } 
266

267 268
      TEST_EXIT_DBG(element)("Should not happen!\n");

269
      if (element->isLeaf() && !isLeafElement()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
270
	MeshStructure *structure = new MeshStructure();
271 272 273
	cont = skipBranch(structure);
	structure->commit();

274 275 276
	MeshStructure_ED *elData = new MeshStructure_ED(element->getElementData());
	elData->setStructure(structure);
	element->setElementData(elData);
277 278 279 280
      } else {
	cont = nextElement();
      }

Thomas Witkowski's avatar
Thomas Witkowski committed
281 282
      if (elInfo)
	elInfo = stack.traverseNext(elInfo);
283 284 285
    }

    // refine mesh
286
    bool finished = true;
287

288
    do {
289
      finished = true;
290 291 292 293
      if (macroElIndex == -1)
	elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
      else
	elInfo = stack.traverseFirstOneMacro(mesh, macroElIndex, -1, Mesh::CALL_LEAF_EL);      
294
      while (elInfo) {
295
	Element *element = elInfo->getElement();
296
	if (element->getElementData(MESH_STRUCTURE) != NULL) {
297 298 299 300 301 302 303
	  element->setMark(1);
	  finished = false;
	} else {
	  element->setMark(0);
	}
	elInfo = stack.traverseNext(elInfo);
      }
304 305 306 307 308

      if (!finished) {
#if (DEBUG != 0)
	int oldMeshIndex = mesh->getChangeIndex();
#endif
309

310 311 312 313
	if (macroElIndex == -1)
	  manager->refineMesh(mesh);
	else
	  manager->refineMacroElement(mesh, macroElIndex);
314

315 316 317 318 319
#if (DEBUG != 0)
	TEST_EXIT(oldMeshIndex != mesh->getChangeIndex())
	  ("Mesh has not been changed by refinement procedure!\n");
#endif
      }
320
    } while (!finished);
321 322
  }

323

324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345
  void MeshStructure::print(bool resetCode)
  {
    FUNCNAME("MeshStructure::print()");
    
    std::stringstream oss;
    
    if (empty()) {
      oss << "-" << std::endl;
    }	else {	
      if (resetCode)
	reset();

      bool cont = true;
      while (cont) {
	if (isLeafElement())
	  oss << "0";
	else
	  oss << "1";
	
	cont = nextElement();
      }
    }
346 347 348 349 350 351 352 353 354 355
   
    if (oss.str().length() < 255) {
      MSG("Mesh structure code: %s\n", oss.str().c_str());
    } else {
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
      std::cout << "[" << MPI::COMM_WORLD.Get_rank() << "]                Mesh structure code: " << oss.str() << "\n";
#else
      std::cout << "                Mesh structure code: " << oss.str() << "\n";
#endif
    }
356 357 358
  }


359 360 361 362
  bool MeshStructure::compare(MeshStructure &other)
  {
    return (other.getCode() == code);
  }
363 364 365 366 367 368 369 370 371


  void MeshStructure::getMeshStructureValues(Mesh *mesh,
					     int macroElIndex,
					     const DOFVector<double>* vec,
					     std::vector<double>& values)
  {
    FUNCNAME("MeshStructure::getMeshStructureValues()");

372 373 374
    TEST_EXIT_DBG(mesh)("No mesh defined!\n");
    TEST_EXIT_DBG(vec)("No DOFVector defined!\n");

375 376 377 378 379 380
    values.clear();

    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirstOneMacro(mesh, macroElIndex, -1, 
						 Mesh::CALL_EVERY_EL_PREORDER);
    while (elInfo) {
381
      if (elInfo->getLevel() == 0)
382 383
	for (int i = 0; i < mesh->getGeo(VERTEX); i++)
	  values.push_back((*vec)[elInfo->getElement()->getDof(i, 0)]);
384 385 386

      if (!elInfo->getElement()->isLeaf())
	values.push_back((*vec)[elInfo->getElement()->getChild(0)->getDof(mesh->getDim(), 0)]);      
387 388 389 390 391 392 393 394 395 396 397 398 399

      elInfo = stack.traverseNext(elInfo);
    }
  }


  void MeshStructure::setMeshStructureValues(Mesh *mesh,
					     int macroElIndex,
					     DOFVector<double>* vec,
					     const std::vector<double>& values)
  {
    FUNCNAME("MeshStructure::setMeshStructureValues()");

400 401 402 403
    TEST_EXIT_DBG(mesh)("No mesh defined!\n");
    TEST_EXIT_DBG(vec)("No DOFVector defined!\n");
    TEST_EXIT_DBG(static_cast<int>(values.size()) >= mesh->getGeo(VERTEX))
      ("Should not happen!\n");
404 405 406 407 408 409 410

    unsigned int counter = 0;

    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirstOneMacro(mesh, macroElIndex, -1, 
						 Mesh::CALL_EVERY_EL_PREORDER);
    while (elInfo) {
411
      if (elInfo->getLevel() == 0)
412
	for (int i = 0; i < mesh->getGeo(VERTEX); i++)
413
	  (*vec)[elInfo->getElement()->getDof(i, 0)] = values[counter++];      
414

415 416 417 418 419
      if (!elInfo->getElement()->isLeaf()) {
	TEST_EXIT_DBG(counter < values.size())("Should not happen!\n");
	
	(*vec)[elInfo->getElement()->getChild(0)->getDof(mesh->getDim(), 0)] =
	  values[counter++];      
420 421 422 423 424 425 426
      }

      elInfo = stack.traverseNext(elInfo);
    }
      
  }

427
}