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

Mesh.cc 41.3 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12
//
// Software License for AMDiS
//
// Copyright (c) 2010 Dresden University of Technology 
// All rights reserved.
// Authors: Simon Vey, Thomas Witkowski et al.
//
// This file is part of AMDiS
//
// See also license.opensource.txt in the distribution.


Thomas Witkowski's avatar
Thomas Witkowski committed
13 14 15 16 17 18
#include <algorithm>
#include <set>
#include <map>

#include "time.h"

19 20 21 22
#include "io/MacroReader.h"
#include "io/MacroInfo.h"
#include "io/MacroWriter.h"

23 24 25 26
#include "AdaptStationary.h"
#include "AdaptInstationary.h"
#include "FiniteElemSpace.h"
#include "ElementData.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
27
#include "ElementDofIterator.h"
28 29 30
#include "MacroElement.h"
#include "Mesh.h"
#include "Traverse.h"
31
#include "Initfile.h"
32 33 34 35 36 37 38
#include "FixVec.h"
#include "DOFVector.h"
#include "CoarseningManager.h"
#include "DOFIterator.h"
#include "VertexVector.h"
#include "PeriodicMap.h"
#include "Projection.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
39
#include "ElInfoStack.h"
40
#include "Serializer.h"
41
#include "Lagrange.h"
42 43 44 45 46 47 48 49 50

namespace AMDiS {

#define TIME_USED(f,s) ((double)((s)-(f))/(double)CLOCKS_PER_SEC)

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

51 52 53 54 55 56 57 58 59 60 61 62 63 64
  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);
65 66 67 68 69 70 71 72 73 74 75

  //**************************************************************************
  //  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;
76 77
  const Flag Mesh::CALL_MG_LEVEL           = 0X4000L;  
  const Flag Mesh::CALL_REVERSE_MODE       = 0X8000L;
78 79


80
  vector<DegreeOfFreedom> Mesh::dof_used;
81
  const int Mesh::MAX_DOF = 100;
82
  map<pair<DegreeOfFreedom, int>, DegreeOfFreedom*> Mesh::serializedDOFs;
83

84
  Mesh::Mesh(string aName, int dimension) 
85 86 87 88 89 90 91 92
    : name(aName), 
      dim(dimension), 
      nVertices(0),
      nEdges(0),
      nLeaves(0), 
      nElements(0),
      parametric(NULL), 
      preserveCoarseDOFs(false),
93 94
      nDofEl(0),
      nDof(dimension, DEFAULT_VALUE, 0),
95 96 97 98 99 100
      nNodeEl(0),
      node(dimension, DEFAULT_VALUE, 0),
      elementPrototype(NULL),
      elementDataPrototype(NULL),
      elementIndex(-1),
      initialized(false),
101
      macroFileInfo(NULL),
102
      changeIndex(0),
103 104
      final_lambda(dimension, DEFAULT_VALUE, 0.0)
  {
105
    FUNCNAME("Mesh::Mesh()");
106

107 108 109 110
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    nParallelPreRefinements = 0;
#endif

111 112 113
    // set default element prototype
    switch(dim) {
    case 1:
Thomas Witkowski's avatar
Thomas Witkowski committed
114
      elementPrototype = new Line(this);
115 116
      break;
    case 2:
Thomas Witkowski's avatar
Thomas Witkowski committed
117
      elementPrototype = new Triangle(this);
118 119
      break;
    case 3:
Thomas Witkowski's avatar
Thomas Witkowski committed
120
      elementPrototype = new Tetrahedron(this);
121 122 123 124 125 126 127
      break;
    default:
      ERROR_EXIT("invalid dimension\n");
    }

    elementPrototype->setIndex(-1);

128 129
    elementIndex = 0;
  }
130

131

132
  Mesh::~Mesh()
133
  {
134
    deleteMeshStructure();
135

136 137 138 139 140 141
    if (macroFileInfo != NULL)
      delete macroFileInfo;    
    if (elementPrototype)
      delete elementPrototype;    
    if (elementDataPrototype)
      delete elementDataPrototype;    
142
    
143
    for (unsigned int i = 0; i < admin.size(); i++)
144
      delete admin[i];    
145
  }
146

147

148
  Mesh& Mesh::operator=(const Mesh& m)
149
  {
150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167
    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;
168 169
    nDofEl = m.nDofEl;
    nDof = m.nDof;
170 171 172 173 174 175 176 177 178
    nNodeEl = m.nNodeEl;
    node = m.node;
    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
    // mapIndex[i] is the index of the MacroElement element in the vector
    // macroElements, for which holds: element->getIndex() = i    
189
    map<int, int> mapIndex;
190 191 192 193 194 195

    // 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
    // Go through all MacroElements of mesh m, and create for every a new
    // MacroElement in this mesh.
201
    for (deque<MacroElement*>::const_iterator it = m.macroElements.begin();
202
	 it != m.macroElements.end(); ++it, insertCounter++) {
203

204
      // Create new MacroElement.
Thomas Witkowski's avatar
Thomas Witkowski committed
205
      MacroElement *el = new MacroElement(dim);
206

207 208 209 210 211 212 213 214 215 216
      // 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.
217
      mapIndex.insert(pair<int, int>(el->getIndex(), insertCounter));
218 219 220 221 222
    }

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

    // 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;
  }

252

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

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

265

266 267 268 269 270 271
  void Mesh::addMacroElement(MacroElement* me) 
  {
    macroElements.push_back(me); 
    me->setIndex(macroElements.size());
  }

272

273
  void Mesh::removeMacroElements(std::set<MacroElement*>& macros,
274
				 vector<const FiniteElemSpace*>& feSpaces) 
275 276 277
  {
    FUNCNAME("Mesh::removeMacroElement()");

278 279
    typedef map<const DegreeOfFreedom*, std::set<MacroElement*> > DofElMap;
    typedef map<const DegreeOfFreedom*, GeoIndex> DofPosMap;
Thomas Witkowski's avatar
Thomas Witkowski committed
280

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

283 284 285 286 287 288 289 290 291 292 293 294

    // === 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 = feSpaces[0];
    for (unsigned int i = 1; i < feSpaces.size(); i++)
      if (feSpaces[i]->getBasisFcts()->getDegree() >
	  feSpace->getBasisFcts()->getDegree())
	feSpace = feSpaces[i];

    
295 296
    // === Determine to all DOFs in mesh the macro elements where the DOF  ===
    // === is part of.                                                     ===
297

298 299
    // 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
300 301
    DofElMap dofsOwner;
    DofPosMap dofsPosIndex;
302

303 304 305 306 307
    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
308
      do {
309 310
	dofsOwner[elDofIter.getStartDof()].insert(elInfo->getMacroElement());
	dofsPosIndex[elDofIter.getStartDof()] = elDofIter.getPosIndex();
Thomas Witkowski's avatar
Thomas Witkowski committed
311
      } while (elDofIter.nextStrict());
312 313 314 315

      elInfo = stack.traverseNext(elInfo);
    }		   

Thomas Witkowski's avatar
Thomas Witkowski committed
316

317 318
    // === Remove macro elements from mesh macro element list. ===

319
    // Removing arbitrary elements from an deque is very slow. Therefore, we
320 321 322
    // 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.

323 324
    deque<MacroElement*> newMacroElements;
    for (deque<MacroElement*>::iterator elIter = macroElements.begin();
325 326 327 328 329
	 elIter != macroElements.end(); ++elIter) {
      // 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.
      if (macros.find(*elIter) == macros.end())
330
	newMacroElements.push_back(*elIter);     
331
    }
332

333 334 335 336
    // And replace the macro element list with the new one.
    macroElements.clear();
    macroElements = newMacroElements;

Thomas Witkowski's avatar
Thomas Witkowski committed
337

338
    // === For all macro elements to be deleted, delete them also to be neighbours ===
339 340
    // === of some other macro elements. Furtheremore, delete the whole element    ===
    // === hierarchie structure of the macro element.                              ===
341 342 343
    
    for (std::set<MacroElement*>::iterator macroIt = macros.begin();
	 macroIt != macros.end(); ++macroIt) {
Thomas Witkowski's avatar
Thomas Witkowski committed
344

345 346
      // Go through all neighbours of the macro element and remove this macro element
      // to be neighbour of some other macro element.
347
      for (int i = 0; i < getGeo(NEIGH); i++)
348
	if ((*macroIt)->getNeighbour(i))
349
	  for (int j = 0; j < getGeo(NEIGH); j++)
Thomas Witkowski's avatar
Thomas Witkowski committed
350
	    if ((*macroIt)->getNeighbour(i)->getNeighbour(j) == *macroIt)
351 352
	      (*macroIt)->getNeighbour(i)->setNeighbour(j, NULL);

353 354 355 356
      // Delete element hierarchie
      if (!(*macroIt)->getElement()->isLeaf()) {
	delete (*macroIt)->getElement()->getChild(0);
	delete (*macroIt)->getElement()->getChild(1);
357 358 359 360 361

	(*macroIt)->getElement()->child[0] = NULL;
	(*macroIt)->getElement()->child[1] = NULL;

	(*macroIt)->getElement()->setElementData(elementDataPrototype->clone()); 
362
      }
363 364 365

      // We will delete at least some element DOFs in the next but will keep the
      // element object still in memory. Therefore the DOF pointer must be set to be
366
      // invalid.
367
      (*macroIt)->getElement()->setDofValid(false);
368
    }     
369 370


371 372
    // === Check now all the DOFs that have no owner anymore and therefore have ===
    // === to be removed.                                                       ===
373

Thomas Witkowski's avatar
Thomas Witkowski committed
374
    for (DofElMap::iterator dofsIt = dofsOwner.begin(); 
375 376 377 378 379 380 381 382 383 384 385 386 387 388
	 dofsIt != dofsOwner.end(); ++dofsIt) {
      
      bool deleteDof = true;

      for (std::set<MacroElement*>::iterator elIter = dofsIt->second.begin();
	   elIter != dofsIt->second.end(); ++elIter) {
	std::set<MacroElement*>::iterator mIt = macros.find(*elIter);
	if (mIt == macros.end()) {
	  deleteDof = false;
	  break;
	}
      }

      if (deleteDof)
389
	freeDof(const_cast<DegreeOfFreedom*>(dofsIt->first), 
390
		dofsPosIndex[dofsIt->first]);      
391 392
    }

393

394
    // === Update number of elements, vertices, etc. ===
395

396 397 398 399 400 401 402 403
    nLeaves = 0;
    nElements = 0;

    std::set<const DegreeOfFreedom*> allVertices;

    elInfo = stack.traverseFirst(this, -1, Mesh::CALL_EVERY_EL_PREORDER);
    while (elInfo) {
      nElements++;
404

405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420
      if (elInfo->getElement()->isLeaf()) {
	nLeaves++;

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

      elInfo = stack.traverseNext(elInfo);
    }

    nVertices = allVertices.size();

    // === 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.                               ===
421
  }
422

423

424 425 426 427 428 429
  void Mesh::addDOFAdmin(DOFAdmin *localAdmin)
  {    
    FUNCNAME("Mesh::addDOFAdmin()");

    localAdmin->setMesh(this);

430
    TEST_EXIT(find(admin.begin(), admin.end(), localAdmin) == admin.end())
Thomas Witkowski's avatar
Thomas Witkowski committed
431 432
      ("admin %s is already associated to mesh %s\n",
       localAdmin->getName().c_str(), this->getName().c_str());
433 434 435

    admin.push_back(localAdmin);

436
    nDofEl = 0;
437

438
    localAdmin->setNumberOfPreDofs(VERTEX, nDof[VERTEX]);
439 440
    nDof[VERTEX] += localAdmin->getNumberOfDofs(VERTEX);
    nDofEl += getGeo(VERTEX) * nDof[VERTEX];
441

442
    if (dim > 1) {
443
      localAdmin->setNumberOfPreDofs(EDGE, nDof[EDGE]);
444 445
      nDof[EDGE] += localAdmin->getNumberOfDofs(EDGE);
      nDofEl += getGeo(EDGE) * nDof[EDGE];
446 447
    }

448 449
    localAdmin->setNumberOfPreDofs(CENTER, nDof[CENTER]);
    nDof[CENTER] += localAdmin->getNumberOfDofs(CENTER);
450
    nDofEl += nDof[CENTER];
451

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

454 455
    node[VERTEX] = 0;
    nNodeEl = getGeo(VERTEX);
456

457 458
    if (dim > 1) {
      node[EDGE] = nNodeEl;
459
      if (nDof[EDGE] > 0) 
460
	nNodeEl += getGeo(EDGE);
461 462
    }

463
    if (dim == 3) {
464
      localAdmin->setNumberOfPreDofs(FACE, nDof[FACE]);
465 466
      nDof[FACE] += localAdmin->getNumberOfDofs(FACE);
      nDofEl += getGeo(FACE) * nDof[FACE];
467
      node[FACE] = nNodeEl;
468
      if (nDof[FACE] > 0) 
469
	nNodeEl += getGeo(FACE);
470 471
    }

472
    node[CENTER] = nNodeEl;
473
    if (nDof[CENTER] > 0)
474
      nNodeEl += 1;
475 476
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
477

478 479
  void Mesh::dofCompress()
  {
480
    FUNCNAME("Mesh::dofCompress()");
481

Thomas Witkowski's avatar
Thomas Witkowski committed
482 483
    for (unsigned int iadmin = 0; iadmin < admin.size(); iadmin++) {
      DOFAdmin* compressAdmin = admin[iadmin];
484 485

      TEST_EXIT_DBG(compressAdmin)("no admin[%d] in mesh\n", iadmin);
486
      
Thomas Witkowski's avatar
Thomas Witkowski committed
487
      int size = compressAdmin->getSize();
Thomas Witkowski's avatar
Thomas Witkowski committed
488
      if (size < 1 || 
489
	  compressAdmin->getUsedDofs() < 1 || 
Thomas Witkowski's avatar
Thomas Witkowski committed
490
	  compressAdmin->getHoleCount() < 1)    
491
	continue;
Thomas Witkowski's avatar
Thomas Witkowski committed
492

493
      vector<int> newDofIndex(size);     
494
      compressAdmin->compress(newDofIndex);
Thomas Witkowski's avatar
Thomas Witkowski committed
495 496 497 498

      Flag fill_flag = (preserveCoarseDOFs ?  
			Mesh::CALL_EVERY_EL_PREORDER | Mesh::FILL_NOTHING :
			Mesh::CALL_LEAF_EL | Mesh::FILL_NOTHING);          
499 500 501
      TraverseStack stack;
      ElInfo *elInfo = stack.traverseFirst(this, -1, fill_flag);
      while (elInfo) {
502
	elInfo->getElement()->newDofFct1(compressAdmin, newDofIndex);
503 504 505 506 507
	elInfo = stack.traverseNext(elInfo);
      }

      elInfo = stack.traverseFirst(this, -1, fill_flag);
      while (elInfo) {
508
	elInfo->getElement()->newDofFct2(compressAdmin);
509 510
	elInfo = stack.traverseNext(elInfo);
      }
511
    }       
512 513 514
  }


515
  DegreeOfFreedom *Mesh::getDof(GeoIndex position)
516
  {
517
    FUNCNAME("Mesh::getDof()");
518

519
    TEST_EXIT_DBG(position >= CENTER && position <= FACE)
520
      ("unknown position %d\n", position);
521

522
    int ndof = getNumberOfDofs(position);
523
    if (ndof <= 0) 
524
      return NULL;
525

526
    DegreeOfFreedom *dof = new DegreeOfFreedom[ndof];
527

528
    for (int i = 0; i < getNumberOfDOFAdmin(); i++) {
529
      const DOFAdmin *localAdmin = &getDofAdmin(i);
530
      TEST_EXIT_DBG(localAdmin)("no admin[%d]\n", i);
531
      
532
      int n  = localAdmin->getNumberOfDofs(position);
533
      int n0 = localAdmin->getNumberOfPreDofs(position);
534
      
535
      TEST_EXIT_DBG(n + n0 <= ndof)("n=%d, n0=%d too large: ndof=%d\n", n, n0, ndof);
536
      
537
      for (int j = 0; j < n; j++)
538 539
	dof[n0 + j] = const_cast<DOFAdmin*>(localAdmin)->getDOFIndex();
    }
540
 
541
    return dof;
542 543 544
  }


Thomas Witkowski's avatar
Thomas Witkowski committed
545
  DegreeOfFreedom **Mesh::createDofPtrs()
546
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
547
    FUNCNAME("Mesh::createDofPtrs()");
548 549

    if (nNodeEl <= 0)
550
      return NULL;
551

552
    DegreeOfFreedom **ptrs = new DegreeOfFreedom*[nNodeEl];
553
    for (int i = 0; i < nNodeEl; i++)
554 555
      ptrs[i] = NULL;

556
    return ptrs;
557 558
  }

559

560
  void Mesh::freeDofPtrs(DegreeOfFreedom **ptrs)
561
  {
562
    FUNCNAME("Mesh::freeDofPtrs()");
563

564
    TEST_EXIT_DBG(ptrs)("ptrs is NULL!\n");
565 566 567 568

    if (nNodeEl <= 0)
      return;
  
569
    delete [] ptrs;
570 571 572
  }


573
  const DOFAdmin *Mesh::createDOFAdmin(string lname, DimVec<int> lnDof)
574
  {
575
    FUNCNAME("Mesh::createDOFAdmin()");
576
    
Thomas Witkowski's avatar
Thomas Witkowski committed
577
    DOFAdmin *localAdmin = new DOFAdmin(this, lname);
578

Thomas Witkowski's avatar
Thomas Witkowski committed
579
    for (int i = 0; i < dim + 1; i++)
580
      localAdmin->setNumberOfDofs(i, lnDof[i]);
581 582 583

    addDOFAdmin(localAdmin);

584
    return localAdmin;
585 586 587 588 589 590 591
  }


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

592
    for (int i = 0; i < static_cast<int>(admin.size()); i++) {
593
      if (admin[i]->getNumberOfDofs(VERTEX)) {
594 595 596 597
	if (!localAdmin)  
	  localAdmin = admin[i];
	else if (admin[i]->getSize() < localAdmin->getSize())
	  localAdmin = admin[i];
598
      }
599 600
    }

601
    return localAdmin;
602 603
  }

604

605
  void Mesh::freeDof(DegreeOfFreedom* dof, GeoIndex position)
606
  {
607
    FUNCNAME("Mesh::freeDof()");
608

609
    TEST_EXIT_DBG(position >= CENTER && position <= FACE)
610
      ("unknown position %d\n", position);
611

612
    int ndof = nDof[position];
613 614
    if (ndof) {
      if (!dof) {
615
	MSG("dof = NULL, but ndof = %d\n", ndof);
616 617
	return;
      }
618
    } else  {
619 620 621
      if (dof)
	MSG("dof != NULL, but ndof = 0\n");

622 623
      return;
    }
624

625
    TEST_EXIT_DBG(ndof <= MAX_DOF)
626
      ("ndof too big: ndof = %d, MAX_DOF = %d\n", ndof, MAX_DOF);
627

628
    for (unsigned int i = 0; i < admin.size(); i++) {
629 630
      int n = admin[i]->getNumberOfDofs(position);
      int n0 = admin[i]->getNumberOfPreDofs(position);
631
      
632 633
      TEST_EXIT_DBG(n + n0 <= ndof)
	("n = %d, n0 = %d too large: ndof = %d\n", n, n0, ndof);
634 635
      
      for (int j = 0; j < n; j++)
636
	admin[i]->freeDofIndex(dof[n0 + j]);      
637
    }
638

639
    delete [] dof;
640
    dof = NULL;
641 642
  }

643

644 645
  void Mesh::freeElement(Element* el)
  {
646
    freeDofPtrs(const_cast<DegreeOfFreedom**>(el->getDof()));
Thomas Witkowski's avatar
Thomas Witkowski committed
647
    delete el;
648 649 650 651 652 653
  }


  Element* Mesh::createNewElement(Element *parent)
  {
    FUNCNAME("Mesh::createNewElement()");
654 655

    TEST_EXIT_DBG(elementPrototype)("no element prototype\n");
656 657

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

Thomas Witkowski's avatar
Thomas Witkowski committed
659
    if (!parent && elementDataPrototype)
660
      el->setElementData(elementDataPrototype->clone()); 
Thomas Witkowski's avatar
Thomas Witkowski committed
661
    else
662
      el->setElementData(NULL); // must be done in ElementData::refineElementData()
663

664 665 666
    return el;
  }

667

668 669
  ElInfo* Mesh::createNewElInfo()
  {
670 671
    FUNCNAME("Mesh::createNewElInfo()");

672
    switch (dim) {
673
    case 1:
Thomas Witkowski's avatar
Thomas Witkowski committed
674
      return new ElInfo1d(this);
675 676
      break;
    case 2:
Thomas Witkowski's avatar
Thomas Witkowski committed
677
      return new ElInfo2d(this);
678 679
      break;
    case 3:
Thomas Witkowski's avatar
Thomas Witkowski committed
680
      return new ElInfo3d(this);
681 682 683 684
      break;
    default:
      ERROR_EXIT("invalid dim\n");
      return NULL;
685
    }
686 687
  }

688

689 690
  bool Mesh::findElInfoAtPoint(const WorldVector<double>& xy,
			       ElInfo *el_info,
691 692
			       DimVec<double>& bary,
			       const MacroElement *start_mel,
693
			       const WorldVector<double> *xy0,
694
			       double *sp)
695 696 697
  {
    static const MacroElement *mel = NULL;
    DimVec<double> lambda(dim, NO_INIT);
Thomas Witkowski's avatar
Thomas Witkowski committed
698
    ElInfo *mel_info = createNewElInfo();
699 700 701 702

    if (start_mel != NULL)
      mel = start_mel;
    else
Thomas Witkowski's avatar
Thomas Witkowski committed
703
      if (mel == NULL || mel->getElement()->getMesh() != this)
704 705 706
	mel = *(macroElements.begin());

    mel_info->setFillFlag(Mesh::FILL_COORDS);
707
    g_xy = &xy;
708
    g_xy0 = xy0;
709
    g_sp = sp;
710 711 712

    mel_info->fillMacroInfo(mel);

713 714 715 716 717
    // We have the care about not to visite a macro element twice. In this case the
    // function would end up in an infinite loop. If a macro element is visited a 
    // second time, what can happen with periodic boundary conditions, the point is
    // not within the mesh!
    std::set<int> macrosVisited;
718 719 720
    std::stack<MacroElement*> active;
    
//     macrosVisited.insert(mel->getIndex());
721

722
    int k;
723
    while ((k = mel_info->worldToCoord(xy, &lambda)) >= 0) {
724
      macrosVisited.insert(mel->getIndex());
725
      if (mel->getNeighbour(k) && !macrosVisited.count(mel->getNeighbour(k)->getIndex())) {
726
	mel = mel->getNeighbour(k);
727 728
// 	if (macrosVisited.count(mel->getIndex())) // Nur für periodische RB interessant, muss noch implementiert werden
// 	  return false;
729

730 731
	mel_info->fillMacroInfo(mel);
	continue;
732 733
      } else {
	for (int i = 0; i < dim + 1; ++i) {
734
	  if (i !=  k && mel->getNeighbour(i) && !macrosVisited.count(mel->getNeighbour(i)->getIndex()))
735 736 737 738
	    active.push(mel->getNeighbour(i));    
	}
	if (active.empty())
	  return false;
739 740
	mel = active.top();
	active.pop();
741 742

	mel_info->fillMacroInfo(mel);
743 744 745 746
      }
    }

    /* now, descend in tree to find leaf element at point */
747
    bool inside = findElementAtPointRecursive(mel_info, lambda, k, el_info);
748 749
    for (int i = 0; i <= dim; i++)
      bary[i] = final_lambda[i];   
750
  
Thomas Witkowski's avatar
Thomas Witkowski committed
751
    delete mel_info;
752

753
    return inside;
754 755
  }

756

757
  bool Mesh::findElementAtPoint(const WorldVector<double>&  xy,
758 759
				Element **elp, 
				DimVec<double>& bary,
760
				const MacroElement *start_mel,
761 762
				const WorldVector<double> *xy0,
				double *sp)
763
  {
764 765
    ElInfo *el_info = createNewElInfo();
    int val = findElInfoAtPoint(xy, el_info, bary, start_mel, xy0, sp);
766 767 768

    *elp = el_info->getElement();

Thomas Witkowski's avatar
Thomas Witkowski committed
769
    delete el_info;
770

771
    return val;
772 773
  }

774

775
  bool Mesh::findElementAtPointRecursive(ElInfo *el_info,
776
					 const DimVec<double>& lambda,
777
					 int outside,
778 779
					 ElInfo* final_el_info)
  {
780
    FUNCNAME("Mesh::findElementAtPointRecursive()");
Thomas Witkowski's avatar
Thomas Witkowski committed
781

782 783
    Element *el = el_info->getElement();
    DimVec<double> c_lambda(dim, NO_INIT);
784 785
    int inside;
    int ichild, c_outside;
786 787 788 789

    if (el->isLeaf()) {
      *final_el_info = *el_info;
      if (outside < 0) {
Thomas Witkowski's avatar
Thomas Witkowski committed
790
	for (int i = 0; i <= dim; i++)
791
	  final_lambda[i] = lambda[i];
Thomas Witkowski's avatar
Thomas Witkowski committed
792
	
793
	return true;
794 795 796 797
      }  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]);
Thomas Witkowski's avatar
Thomas Witkowski committed
798 799 800 801
	  for (int i = 0; i <= dim; i++) 
	    final_lambda[i] = s * c_lambda[i] + (1.0 - s) * lambda[i];
	  
	  if (g_sp)
802
	    *(g_sp) = s;
Thomas Witkowski's avatar
Thomas Witkowski committed
803
	  
804
	  if (dim == 3)