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 39.6 KB
Newer Older
Thomas Witkowski's avatar
Thomas Witkowski committed
1 2 3 4 5 6
#include <algorithm>
#include <set>
#include <map>

#include "time.h"

7 8 9 10
#include "AdaptStationary.h"
#include "AdaptInstationary.h"
#include "FiniteElemSpace.h"
#include "ElementData.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
11
#include "ElementDofIterator.h"
12 13
#include "MacroElement.h"
#include "MacroReader.h"
14
#include "MacroInfo.h"
15 16 17 18 19 20 21 22 23 24 25
#include "Mesh.h"
#include "Traverse.h"
#include "Parameters.h"
#include "FixVec.h"
#include "DOFVector.h"
#include "CoarseningManager.h"
#include "DOFIterator.h"
#include "VertexVector.h"
#include "MacroWriter.h"
#include "PeriodicMap.h"
#include "Projection.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
26
#include "ElInfoStack.h"
27
#include "Serializer.h"
28
#include "Lagrange.h"
29 30 31 32 33 34 35 36 37

namespace AMDiS {

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

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

38 39 40 41 42 43 44 45 46 47 48 49 50 51
  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);
52 53 54 55 56 57 58 59 60 61 62

  //**************************************************************************
  //  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;
63 64
  const Flag Mesh::CALL_MG_LEVEL           = 0X4000L;  
  const Flag Mesh::CALL_REVERSE_MODE       = 0X8000L;
65 66


Thomas Witkowski's avatar
Thomas Witkowski committed
67
  std::vector<DegreeOfFreedom> Mesh::dof_used;
68
  const int Mesh::MAX_DOF = 100;
69
  std::map<std::pair<DegreeOfFreedom, int>, DegreeOfFreedom*> Mesh::serializedDOFs;
70 71 72

  struct delmem { 
    DegreeOfFreedom* ptr;
73
    int len;
74 75 76
  };


Thomas Witkowski's avatar
Thomas Witkowski committed
77
  Mesh::Mesh(std::string aName, int dimension) 
78 79 80 81 82 83 84 85
    : name(aName), 
      dim(dimension), 
      nVertices(0),
      nEdges(0),
      nLeaves(0), 
      nElements(0),
      parametric(NULL), 
      preserveCoarseDOFs(false),
86 87
      nDofEl(0),
      nDof(dimension, DEFAULT_VALUE, 0),
88 89 90 91 92 93
      nNodeEl(0),
      node(dimension, DEFAULT_VALUE, 0),
      elementPrototype(NULL),
      elementDataPrototype(NULL),
      elementIndex(-1),
      initialized(false),
94
      macroFileInfo(NULL),
95
      changeIndex(0),
96 97 98
      final_lambda(dimension, DEFAULT_VALUE, 0.0)
  {

99
    FUNCNAME("Mesh::Mesh()");
100 101 102 103

    // set default element prototype
    switch(dim) {
    case 1:
Thomas Witkowski's avatar
Thomas Witkowski committed
104
      elementPrototype = new Line(this);
105 106
      break;
    case 2:
Thomas Witkowski's avatar
Thomas Witkowski committed
107
      elementPrototype = new Triangle(this);
108 109
      break;
    case 3:
Thomas Witkowski's avatar
Thomas Witkowski committed
110
      elementPrototype = new Tetrahedron(this);
111 112 113 114 115 116 117
      break;
    default:
      ERROR_EXIT("invalid dimension\n");
    }

    elementPrototype->setIndex(-1);

118 119
    elementIndex = 0;
  }
120

121

122
  Mesh::~Mesh()
123
  {
124
    deleteMeshStructure();
125

126 127 128 129 130 131
    if (macroFileInfo != NULL)
      delete macroFileInfo;    
    if (elementPrototype)
      delete elementPrototype;    
    if (elementDataPrototype)
      delete elementDataPrototype;    
132
    
133 134
    for (int i = 0; i < static_cast<int>(admin.size()); i++)
      delete admin[i];    
135
  }
136

137

138
  Mesh& Mesh::operator=(const Mesh& m)
139
  {
140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157
    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;
158 159
    nDofEl = m.nDofEl;
    nDof = m.nDof;
160 161 162 163 164 165 166 167 168 169
    nNodeEl = m.nNodeEl;
    node = m.node;
    newDOF = m.newDOF;
    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
170
      admin[i] = new DOFAdmin(this);
171
      *(admin[i]) = *(m.admin[i]);
172 173
      admin[i]->setMesh(this);
    }
174

175

176
    /* ====================== Copy macro elements =================== */
177
  
178 179 180 181 182 183 184 185 186
    // mapIndex[i] is the index of the MacroElement element in the vector
    // macroElements, for which holds: element->getIndex() = i    
    std::map<int, int> mapIndex;

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

    int insertCounter = 0;
187

188 189
    macroElements.clear();

190 191 192
    // Go through all MacroElements of mesh m, and create for every a new
    // MacroElement in this mesh.
    for (std::deque<MacroElement*>::const_iterator it = m.macroElements.begin();
193
	 it != m.macroElements.end(); ++it, insertCounter++) {
194

195
      // Create new MacroElement.
Thomas Witkowski's avatar
Thomas Witkowski committed
196
      MacroElement *el = new MacroElement(dim);
197

198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240
      // 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.
      mapIndex.insert(std::pair<int, int>(el->getIndex(), insertCounter));
    }

    // Now we have to go through all the new MacroElements, and update the neighbour
    // connections.
    insertCounter = 0;
    for (std::deque<MacroElement*>::const_iterator it = m.macroElements.begin();
	 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.
	macroElements[insertCounter]->
	  setNeighbour(i, macroElements[mapIndex[(*it)->getNeighbour(i)->getIndex()]]);
      }
    }

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

241

242 243 244 245 246 247 248 249 250 251 252 253
  void Mesh::updateNumberOfLeaves()
  {
    nLeaves = 0;

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

254

255 256 257 258 259 260
  void Mesh::addMacroElement(MacroElement* me) 
  {
    macroElements.push_back(me); 
    me->setIndex(macroElements.size());
  }

261

262
  void Mesh::removeMacroElements(std::set<MacroElement*>& macros,
Thomas Witkowski's avatar
Thomas Witkowski committed
263
				 const FiniteElemSpace *feSpace) 
264 265 266
  {
    FUNCNAME("Mesh::removeMacroElement()");

Thomas Witkowski's avatar
Thomas Witkowski committed
267 268 269 270 271 272
    typedef std::map<const DegreeOfFreedom*, std::set<MacroElement*> > DofElMap;
    typedef std::map<const DegreeOfFreedom*, GeoIndex> DofPosMap;

    TEST_EXIT(admin.size() == 1)("Not yet implemented for multiple admins!\n");
    TEST_EXIT(admin[0])("There is something wrong!\n");

273
    // === Determine to all dofs the macro elements where the dof is part of. ===
274

275
    // Map that stores for each dof pointer (which may have a list of dofs)
276
    // all macro element indices that own this dof.
Thomas Witkowski's avatar
Thomas Witkowski committed
277 278
    DofElMap dofsOwner;
    DofPosMap dofsPosIndex;
279

280 281 282 283 284
    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
285
      do {
286
	dofsOwner[elDofIter.getDofPtr()].insert(elInfo->getMacroElement());
Thomas Witkowski's avatar
Thomas Witkowski committed
287 288
	dofsPosIndex[elDofIter.getDofPtr()] = elDofIter.getPosIndex();
      } while (elDofIter.nextStrict());
289 290 291 292

      elInfo = stack.traverseNext(elInfo);
    }		   

Thomas Witkowski's avatar
Thomas Witkowski committed
293

294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312
    // === Remove macro elements from mesh macro element list. ===

    // Removing arbitrary elements from an std::deque is very slow. Therefore, we
    // 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.

    std::deque<MacroElement*> newMacroElements;
    for (std::deque<MacroElement*>::iterator elIter = macroElements.begin();
	 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())
	newMacroElements.push_back(*elIter);      
    }
    // And replace the macro element list with the new one.
    macroElements.clear();
    macroElements = newMacroElements;

Thomas Witkowski's avatar
Thomas Witkowski committed
313

314
    // === For all macro elements to be deleted, delete them also to be neighbours ===
315 316
    // === of some other macro elements. Furtheremore, delete the whole element    ===
    // === hierarchie structure of the macro element.                              ===
317 318 319
    
    for (std::set<MacroElement*>::iterator macroIt = macros.begin();
	 macroIt != macros.end(); ++macroIt) {
Thomas Witkowski's avatar
Thomas Witkowski committed
320

321 322
      // Go through all neighbours of the macro element and remove this macro element
      // to be neighbour of some other macro element.
323 324
      for (int i = 0; i <= dim; i++)
	if ((*macroIt)->getNeighbour(i))
Thomas Witkowski's avatar
Thomas Witkowski committed
325 326
	  for (int j = 0; j <= dim; j++)
	    if ((*macroIt)->getNeighbour(i)->getNeighbour(j) == *macroIt)
327 328
	      (*macroIt)->getNeighbour(i)->setNeighbour(j, NULL);

329 330 331 332 333
      // Delete element hierarchie
      if (!(*macroIt)->getElement()->isLeaf()) {
	delete (*macroIt)->getElement()->getChild(0);
	delete (*macroIt)->getElement()->getChild(1);
      }
334
    }     
335 336


337 338
    // === Check now all the dofs, that have no owner anymore and therefore have ===
    // === to be removed.                                                        ===
339

Thomas Witkowski's avatar
Thomas Witkowski committed
340
    for (DofElMap::iterator dofsIt = dofsOwner.begin(); 
341 342 343 344 345 346 347 348 349 350 351 352 353 354
	 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)
355
	freeDof(const_cast<DegreeOfFreedom*>(dofsIt->first), 
Thomas Witkowski's avatar
Thomas Witkowski committed
356
		dofsPosIndex[dofsIt->first]);
357 358
    }

359

360
    // === Update number of elements, vertices, etc. ===
361

362 363 364 365 366 367 368 369
    nLeaves = 0;
    nElements = 0;

    std::set<const DegreeOfFreedom*> allVertices;

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

371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386
      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.                               ===
387
  }
388

389

390 391 392 393 394 395
  void Mesh::addDOFAdmin(DOFAdmin *localAdmin)
  {    
    FUNCNAME("Mesh::addDOFAdmin()");

    localAdmin->setMesh(this);

396
    std::vector<DOFAdmin*>::iterator dai = 
397
      std::find(admin.begin(), admin.end(), localAdmin);
398

Thomas Witkowski's avatar
Thomas Witkowski committed
399 400 401
    TEST_EXIT(dai == admin.end())
      ("admin %s is already associated to mesh %s\n",
       localAdmin->getName().c_str(), this->getName().c_str());
402

403
    // if this will be required, see the untested code in revision < 224
404
    //    TEST_EXIT(!initialized)("Adding DOFAdmins to initilized meshes does not work yet!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
405

406 407
    admin.push_back(localAdmin);

408
    nDofEl = 0;
409

410 411 412
    localAdmin->setNumberOfPreDOFs(VERTEX,nDof[VERTEX]);
    nDof[VERTEX] += localAdmin->getNumberOfDofs(VERTEX);
    nDofEl += getGeo(VERTEX) * nDof[VERTEX];
413

414
    if (dim > 1) {
415 416 417
      localAdmin->setNumberOfPreDOFs(EDGE,nDof[EDGE]);
      nDof[EDGE] += localAdmin->getNumberOfDofs(EDGE);
      nDofEl += getGeo(EDGE) * nDof[EDGE];
418 419
    }

420 421 422
    localAdmin->setNumberOfPreDOFs(CENTER,nDof[CENTER]);
    nDof[CENTER]  += localAdmin->getNumberOfDofs(CENTER);
    nDofEl += nDof[CENTER];
423

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

426 427
    node[VERTEX] = 0;
    nNodeEl = getGeo(VERTEX);
428

429 430
    if (dim > 1) {
      node[EDGE] = nNodeEl;
431
      if (nDof[EDGE] > 0) 
432
	nNodeEl += getGeo(EDGE);
433 434
    }

435
    if (dim == 3) {
436 437 438
      localAdmin->setNumberOfPreDOFs(FACE,nDof[FACE]);
      nDof[FACE] += localAdmin->getNumberOfDofs(FACE);
      nDofEl += getGeo(FACE) * nDof[FACE];
439
      node[FACE] = nNodeEl;
440
      if (nDof[FACE] > 0) 
441
	nNodeEl += getGeo(FACE);
442 443
    }

444
    node[CENTER] = nNodeEl;
445
    if (nDof[CENTER] > 0)
446
      nNodeEl += 1;
447 448
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
449

450 451
  void Mesh::dofCompress()
  {
452
    FUNCNAME("Mesh::dofCompress()");
453

Thomas Witkowski's avatar
Thomas Witkowski committed
454 455
    for (unsigned int iadmin = 0; iadmin < admin.size(); iadmin++) {
      DOFAdmin* compressAdmin = admin[iadmin];
456 457

      TEST_EXIT_DBG(compressAdmin)("no admin[%d] in mesh\n", iadmin);
458
      
Thomas Witkowski's avatar
Thomas Witkowski committed
459
      int size = compressAdmin->getSize();
Thomas Witkowski's avatar
Thomas Witkowski committed
460
      if (size < 1 || 
461
	  compressAdmin->getUsedDofs() < 1 || 
Thomas Witkowski's avatar
Thomas Witkowski committed
462
	  compressAdmin->getHoleCount() < 1)    
463
	continue;
Thomas Witkowski's avatar
Thomas Witkowski committed
464

465 466 467
      newDOF.resize(size);
      
      compressAdmin->compress(newDOF);
Thomas Witkowski's avatar
Thomas Witkowski committed
468 469 470 471

      Flag fill_flag = (preserveCoarseDOFs ?  
			Mesh::CALL_EVERY_EL_PREORDER | Mesh::FILL_NOTHING :
			Mesh::CALL_LEAF_EL | Mesh::FILL_NOTHING);          
472 473 474
      TraverseStack stack;
      ElInfo *elInfo = stack.traverseFirst(this, -1, fill_flag);
      while (elInfo) {
475
	elInfo->getElement()->newDofFct1(compressAdmin);
476 477 478 479 480
	elInfo = stack.traverseNext(elInfo);
      }

      elInfo = stack.traverseFirst(this, -1, fill_flag);
      while (elInfo) {
481
	elInfo->getElement()->newDofFct2(compressAdmin);
482 483 484
	elInfo = stack.traverseNext(elInfo);
      }

485
      newDOF.resize(0);
486
    }       
487 488 489
  }


490
  DegreeOfFreedom *Mesh::getDof(GeoIndex position)
491
  {
492
    FUNCNAME("Mesh::getDof()");
493

494
    TEST_EXIT_DBG(position >= CENTER && position <= FACE)
495
      ("unknown position %d\n", position);
496

497
    int ndof = getNumberOfDofs(position);
498
    if (ndof <= 0) 
499
      return NULL;
500

501
    DegreeOfFreedom *dof = new DegreeOfFreedom[ndof];
502

503
    for (int i = 0; i < getNumberOfDOFAdmin(); i++) {
504
      const DOFAdmin *localAdmin = &getDofAdmin(i);
505
      TEST_EXIT_DBG(localAdmin)("no admin[%d]\n", i);
506
      
507
      int n  = localAdmin->getNumberOfDofs(position);
508 509
      int n0 = localAdmin->getNumberOfPreDOFs(position);
      
510
      TEST_EXIT_DBG(n + n0 <= ndof)("n=%d, n0=%d too large: ndof=%d\n", n, n0, ndof);
511
      
512
      for (int j = 0; j < n; j++)
513 514
	dof[n0 + j] = const_cast<DOFAdmin*>(localAdmin)->getDOFIndex();
    }
515
 
516
    return dof;
517 518 519
  }


Thomas Witkowski's avatar
Thomas Witkowski committed
520
  DegreeOfFreedom **Mesh::createDofPtrs()
521
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
522
    FUNCNAME("Mesh::createDofPtrs()");
523 524

    if (nNodeEl <= 0)
525
      return NULL;
526

527
    DegreeOfFreedom **ptrs = new DegreeOfFreedom*[nNodeEl];
528
    for (int i = 0; i < nNodeEl; i++)
529 530
      ptrs[i] = NULL;

531
    return ptrs;
532 533
  }

534

535
  void Mesh::freeDofPtrs(DegreeOfFreedom **ptrs)
536
  {
537
    FUNCNAME("Mesh::freeDofPtrs()");
538

539
    TEST_EXIT_DBG(ptrs)("ptrs is NULL!\n");
540 541 542 543

    if (nNodeEl <= 0)
      return;
  
544
    delete [] ptrs;
545 546 547
  }


548
  const DOFAdmin *Mesh::createDOFAdmin(std::string lname, DimVec<int> lnDof)
549
  {
550
    FUNCNAME("Mesh::createDOFAdmin()");
551

Thomas Witkowski's avatar
Thomas Witkowski committed
552
    DOFAdmin *localAdmin = new DOFAdmin(this, lname);
553

Thomas Witkowski's avatar
Thomas Witkowski committed
554
    for (int i = 0; i < dim + 1; i++)
555
      localAdmin->setNumberOfDOFs(i, lnDof[i]);
556 557 558

    addDOFAdmin(localAdmin);

559
    return localAdmin;
560 561 562 563 564 565 566
  }


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

567
    for (int i = 0; i < static_cast<int>(admin.size()); i++) {
568
      if (admin[i]->getNumberOfDofs(VERTEX)) {
569 570 571 572
	if (!localAdmin)  
	  localAdmin = admin[i];
	else if (admin[i]->getSize() < localAdmin->getSize())
	  localAdmin = admin[i];
573
      }
574 575
    }

576
    return localAdmin;
577 578
  }

579

580
  void Mesh::freeDof(DegreeOfFreedom* dof, GeoIndex position)
581
  {
582
    FUNCNAME("Mesh::freeDof()");
583

584
    TEST_EXIT_DBG(position >= CENTER && position <= FACE)
585
      ("unknown position %d\n", position);
586

587
    int ndof = nDof[position];
588 589
    if (ndof) {
      if (!dof) {
590
	MSG("dof = NULL, but ndof = %d\n", ndof);
591 592
	return;
      }
593
    } else  {
594 595 596
      if (dof)
	MSG("dof != NULL, but ndof = 0\n");

597 598
      return;
    }
599

600
    TEST_EXIT_DBG(ndof <= MAX_DOF)
601
      ("ndof too big: ndof = %d, MAX_DOF = %d\n", ndof, MAX_DOF);
602

603
    for (unsigned int i = 0; i < admin.size(); i++) {
604
      DOFAdmin *localAdmin = admin[i];
605
      int n = localAdmin->getNumberOfDofs(position);
606 607
      int n0 = localAdmin->getNumberOfPreDOFs(position);
      
608 609
      TEST_EXIT_DBG(n + n0 <= ndof)
	("n = %d, n0 = %d too large: ndof = %d\n", n, n0, ndof);
610 611
      
      for (int j = 0; j < n; j++)
612
	localAdmin->freeDofIndex(dof[n0 + j]);      
613
    }
614

615
    delete [] dof;
616
    dof = NULL;
617 618
  }

619

620 621
  void Mesh::freeElement(Element* el)
  {
622
    freeDofPtrs(const_cast<DegreeOfFreedom**>(el->getDof()));
Thomas Witkowski's avatar
Thomas Witkowski committed
623
    delete el;
624 625 626 627 628 629
  }


  Element* Mesh::createNewElement(Element *parent)
  {
    FUNCNAME("Mesh::createNewElement()");
630 631

    TEST_EXIT_DBG(elementPrototype)("no element prototype\n");
632 633

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

Thomas Witkowski's avatar
Thomas Witkowski committed
635
    if (!parent && elementDataPrototype)
636
      el->setElementData(elementDataPrototype->clone()); 
Thomas Witkowski's avatar
Thomas Witkowski committed
637
    else
638
      el->setElementData(NULL); // must be done in ElementData::refineElementData()
639

640 641 642
    return el;
  }

643

644 645
  ElInfo* Mesh::createNewElInfo()
  {
646 647
    FUNCNAME("Mesh::createNewElInfo()");

648 649
    switch(dim) {
    case 1:
Thomas Witkowski's avatar
Thomas Witkowski committed
650
      return new ElInfo1d(this);
651 652
      break;
    case 2:
Thomas Witkowski's avatar
Thomas Witkowski committed
653
      return new ElInfo2d(this);
654 655
      break;
    case 3:
Thomas Witkowski's avatar
Thomas Witkowski committed
656
      return new ElInfo3d(this);
657 658 659 660
      break;
    default:
      ERROR_EXIT("invalid dim\n");
      return NULL;
661
    }
662 663
  }

664

665 666
  bool Mesh::findElInfoAtPoint(const WorldVector<double>& xy,
			       ElInfo *el_info,
667 668
			       DimVec<double>& bary,
			       const MacroElement *start_mel,
669
			       const WorldVector<double> *xy0,
670
			       double *sp)
671 672 673
  {
    static const MacroElement *mel = NULL;
    DimVec<double> lambda(dim, NO_INIT);
Thomas Witkowski's avatar
Thomas Witkowski committed
674
    ElInfo *mel_info = createNewElInfo();
675 676 677 678

    if (start_mel != NULL)
      mel = start_mel;
    else
Thomas Witkowski's avatar
Thomas Witkowski committed
679
      if (mel == NULL || mel->getElement()->getMesh() != this)
680 681 682
	mel = *(macroElements.begin());

    mel_info->setFillFlag(Mesh::FILL_COORDS);
683
    g_xy = &xy;
684
    g_xy0 = xy0;
685
    g_sp = sp;
686 687 688

    mel_info->fillMacroInfo(mel);

689 690 691 692 693 694 695
    // 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;
    macrosVisited.insert(mel->getIndex());

696
    int k;
697 698
    while ((k = mel_info->worldToCoord(xy, &lambda)) >= 0) {
      if (mel->getNeighbour(k)) {
699 700 701 702 703
	mel = mel->getNeighbour(k);
	if (macrosVisited.count(mel->getIndex()))
	  return false;

	macrosVisited.insert(mel->getIndex());
704 705 706 707 708 709 710
	mel_info->fillMacroInfo(mel);
	continue;
      }
      break;
    }

    /* now, descend in tree to find leaf element at point */
711
    bool inside = findElementAtPointRecursive(mel_info, lambda, k, el_info);
712 713
    for (int i = 0; i <= dim; i++)
      bary[i] = final_lambda[i];   
714
  
Thomas Witkowski's avatar
Thomas Witkowski committed
715
    delete mel_info;
716

717
    return inside;
718 719 720
  }

  bool Mesh::findElementAtPoint(const WorldVector<double>&  xy,
721 722
				Element **elp, 
				DimVec<double>& bary,
723
				const MacroElement *start_mel,
724 725
				const WorldVector<double> *xy0,
				double *sp)
726
  {
727 728
    ElInfo *el_info = createNewElInfo();
    int val = findElInfoAtPoint(xy, el_info, bary, start_mel, xy0, sp);
729 730 731

    *elp = el_info->getElement();

Thomas Witkowski's avatar
Thomas Witkowski committed
732
    delete el_info;
733

734
    return val;
735 736
  }

737
  bool Mesh::findElementAtPointRecursive(ElInfo *el_info,
738
					 const DimVec<double>& lambda,
739
					 int outside,
740 741
					 ElInfo* final_el_info)
  {
742
    FUNCNAME("Mesh::findElementAtPointRecursive()");
Thomas Witkowski's avatar
Thomas Witkowski committed
743

744 745
    Element *el = el_info->getElement();
    DimVec<double> c_lambda(dim, NO_INIT);
746 747
    int inside;
    int ichild, c_outside;
748 749 750 751

    if (el->isLeaf()) {
      *final_el_info = *el_info;
      if (outside < 0) {
Thomas Witkowski's avatar
Thomas Witkowski committed
752
	for (int i = 0; i <= dim; i++)
753
	  final_lambda[i] = lambda[i];
Thomas Witkowski's avatar
Thomas Witkowski committed
754
	
755
	return true;
756 757 758 759
      }  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
760 761 762 763
	  for (int i = 0; i <= dim; i++) 
	    final_lambda[i] = s * c_lambda[i] + (1.0 - s) * lambda[i];
	  
	  if (g_sp)
764
	    *(g_sp) = s;
Thomas Witkowski's avatar
Thomas Witkowski committed
765
	  
766
	  if (dim == 3) 
767
	    MSG("Outside finest level on el %d: s = %.3e\n", el->getIndex(), s);
768

769
	  return false;  /* ??? */
770
	} else {
771
	  return false;
772
	}
773
      }
774 775
    }

776
    ElInfo *c_el_info = createNewElInfo();
777

778
    if (dim == 1) {
779 780 781 782
      if (lambda[0] >= lambda[1]) {
	c_el_info->fillElInfo(0, el_info);
	if (outside >= 0) {
	  outside = el_info->worldToCoord(*(g_xy), &c_lambda);
Thomas Witkowski's avatar
Thomas Witkowski committed
783
	  TEST_EXIT(outside == 0)("point outside domain\n");
784 785 786 787 788 789 790 791
	} else {
	  c_lambda[0] = lambda[0] - lambda[1];
	  c_lambda[1] = 2.0 * lambda[1];
	}
      } else {
	c_el_info->fillElInfo(1, el_info);
	if (outside >= 0)  {
	  outside = el_info->worldToCoord(*(g_xy), &c_lambda);
Thomas Witkowski's avatar
Thomas Witkowski committed
792
	  TEST_EXIT(outside == 0)("point outside domain\n");
793 794 795 796 797 798 799
	} else {
	  c_lambda[1] = lambda[1] - lambda[0];
	  c_lambda[0] = 2.0 * lambda[0];
	}
      }
    } /* DIM == 1 */

800
    if (dim == 2) {
801 802 803 804
      if (lambda[0] >= lambda[1]) {
	c_el_info->fillElInfo(0, el_info);
	if (el->isNewCoordSet()) {
	  outside = c_el_info->worldToCoord(*(g_xy), &c_lambda);
Thomas Witkowski's avatar
Thomas Witkowski committed
805
	  TEST_EXIT(outside == 0)("outside curved boundary child 0\n");	  
806 807 808 809 810 811 812 813 814
	} else {
	  c_lambda[0] = lambda[2];
	  c_lambda[1] = lambda[0] - lambda[1];
	  c_lambda[2] = 2.0 * lambda[1];
	}
      } else {
	c_el_info->fillElInfo(1, el_info);
	if (el->isNewCoordSet()) {
	  outside = c_el_info->worldToCoord(*(g_xy), &c_lambda);
Thomas Witkowski's avatar
Thomas Witkowski committed
815
	  TEST_EXIT(outside == 0)("outside curved boundary child 1\n");	  
816 817 818 819 820 821 822 823
	} else {
	  c_lambda[0] = lambda[1] - lambda[0];
	  c_lambda[1] = lambda[2];
	  c_lambda[2] = 2.0 * lambda[0];
	}
      }
    } /* DIM == 2 */

824
    if (dim == 3) {
825 826 827 828 829 830 831 832 833 834 835 836 837
      if (el->isNewCoordSet()) {
	if (lambda[0] >= lambda[1])
	  ichild = 0;
	else
	  ichild = 1;
	c_el_info->fillElInfo(ichild, el_info);
	c_outside = c_el_info->worldToCoord(*(g_xy), &c_lambda);

	if (c_outside>=0) {  /* test is other child is better... */
	  DimVec<double> c_lambda2(dim, NO_INIT);
	  ElInfo *c_el_info2 = createNewElInfo();

	  c_el_info2->fillElInfo(1-ichild, el_info);
Thomas Witkowski's avatar
Thomas Witkowski committed
838
	  int c_outside2 = c_el_info2->worldToCoord(*(g_xy), &c_lambda2);
839 840

	  MSG("new_coord CHILD %d: outside=%d, lambda=(%.2f %.2f %.2f %.2f)\n",
841
	      ichild, c_outside, c_lambda[0], c_lambda[1], c_lambda[2], c_lambda[3]);
842
	  MSG("new_coord CHILD %d: outside=%d, lambda=(%.2f %.2f %.2f %.2f)\n",
843
	      1 - ichild, c_outside2, c_lambda2[0], c_lambda2[1], c_lambda2[2],
844 845 846
	      c_lambda2[3]);

	  if ((c_outside2 < 0) || (c_lambda2[c_outside2] > c_lambda[c_outside])) {
847 848
	    for (int i = 0; i <= dim; i++) 
	      c_lambda[i]