MacroReader.cc 48.9 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
#include <string.h>
#include <map>
#include <iostream>
#include <fstream>
26 27 28
#include "MacroReader.h"
#include "MacroWriter.h"
#include "MacroElement.h"
29
#include "MacroInfo.h"
30 31 32 33 34
#include "Boundary.h"
#include "FiniteElemSpace.h"
#include "Mesh.h"
#include "FixVec.h"
#include "ElInfo.h"
35
#include "Initfile.h"
36 37 38 39
#include "DOFIterator.h"
#include "LeafData.h"
#include "VertexVector.h"

40
namespace AMDiS { namespace io {
41

42 43 44
  void MacroReader::restoreMacroDofs(MacroInfo& macroInfo, int check)
  {
    FUNCNAME("MacroReader::restoreMacroDofs()");
45

46 47
    TEST_EXIT(macroInfo.isInitialized())
      ("Cannot set dof since macro info is not initialized.\n");
48

49 50 51 52 53
    Mesh* mesh = macroInfo.mesh;
    int nVertices = macroInfo.nVertices;
    int nElements = macroInfo.nElements;
    int **mel_vertex = macroInfo.mel_vertex;
    DegreeOfFreedom** dof = macroInfo.dof;
54 55

    TEST_EXIT(mesh->getNumberOfMacros() == nElements)
56
      ("Mesh's current macro doesn't match the macro file.\n");
57

58 59 60
    mesh->setNumberOfElements(nElements);
    mesh->setNumberOfLeaves(nElements);
    mesh->setNumberOfVertices(nVertices);
61

62 63 64
    // Vertices dofs
    for (int i = 0; i < nVertices; i++)
      dof[i] = mesh->getDof(VERTEX);
65

66
    std::deque<MacroElement*>::iterator it = mesh->firstMacroElement();
67 68 69

    for (int i = 0; i < mesh->getNumberOfMacros(); i++)
      for (int k = 0; k < mesh->getGeo(VERTEX); k++)
70 71
	const_cast<Element*>((*(it + i))->getElement())->
	  setDof(k, dof[mel_vertex[i][k]]);
72

73 74 75
    // Edge and face dofs
    if (mesh->getDim() > 1)
      boundaryDOFs(mesh);
76

77 78 79
    // Center dofs
    if (mesh->getNumberOfDofs(CENTER)) {
      it = mesh->firstMacroElement();
80

81 82 83 84
      for (int i = 0; i < mesh->getNumberOfMacros(); i++)
	const_cast<Element*>(it[i]->getElement())->
	  setDof(mesh->getNode(CENTER), mesh->getDof(CENTER));
    }
85

86 87 88
    if (check)
      checkMesh(mesh);
  }
89

90
  MacroInfo* MacroReader::readMacro(std::string filename,
91
				    Mesh* mesh,
92
				    std::string periodicFile,
93 94
				    int check)
  {
95
    FUNCNAME("MacroReader::readMacro()");
96

97
    TEST_EXIT(filename != "")("no file specified; filename NULL pointer\n");
98

Thomas Witkowski's avatar
Thomas Witkowski committed
99
    MacroInfo *macroInfo = new MacroInfo();
100 101
    macroInfo->readAMDiSMacro(filename, mesh);

102
    std::deque<MacroElement*>::iterator mel = macroInfo->mel.begin();
103 104 105
    int **melVertex = macroInfo->mel_vertex;
    WorldVector<double> *coords = macroInfo->coords;
    DegreeOfFreedom **dof = macroInfo->dof;
106 107

    // === read periodic data =================================
108
    if (periodicFile != "") {
109 110
      FILE *file = fopen(periodicFile.c_str(), "r");
      TEST_EXIT(file)("can't open file %s\n", periodicFile.c_str());
111 112 113 114 115

      int n;
      int dim = mesh->getDim();

      int el1, el2;
Thomas Witkowski's avatar
Thomas Witkowski committed
116 117
      int *verticesEl1 = new int[dim];
      int *verticesEl2 = new int[dim];
118
      int mode = -1; // 0: drop dofs, 1: associate dofs
119
      int result = 0;
120
      BoundaryType boundaryType;
121
      MacroReader::PeriodicMap periodicMap;
122

123 124
      result = fscanf(file, "%*s %d", &n);
      result = fscanf(file, "%*s %*s %*s %*s %*s %*s %*s %*s %*s %*s %*s");
125

126
      for (int i = 0; i < n; i++) {
127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 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
        std::map<int, int> vertexMapEl1;
        std::map<int, int> vertexMapEl2;

        result = fscanf(file, "%d", &mode);
        TEST_EXIT(result == 1)("mode?\n");

        result = fscanf(file, "%d", &boundaryType);
        TEST_EXIT(result == 1)("boundaryType?\n");

        result = fscanf(file, "%d", &el1);
        TEST_EXIT(result == 1)("el1?\n");

        for (int j = 0; j < dim; j++) {
          result = fscanf(file, "%d", &verticesEl1[j]);
          TEST_EXIT(result == 1)("vertEl1[%d]\n", j);
        }
        result = fscanf(file, "%d", &el2);
        TEST_EXIT(result == 1)("el2?\n");
        for (int j = 0; j < dim; j++) {
          result = fscanf(file, "%d", &verticesEl2[j]);
          TEST_EXIT(result == 1)("vertEl2[%d]\n", j);
        }
        for (int j = 0; j < dim; j++) {
          if (mode == 0)
            periodicMap.setEntry(melVertex[el1][verticesEl1[j]],
                    melVertex[el2][verticesEl2[j]]);
          vertexMapEl1[verticesEl1[j]] = verticesEl2[j];
          vertexMapEl2[verticesEl2[j]] = verticesEl1[j];
        }

        // calculate sides of periodic vertices
        int sideEl1 = 0, sideEl2 = 0;
        if (dim == 1) {
          sideEl1 = verticesEl1[0];
          sideEl2 = verticesEl2[0];
        } else {
          for (int j = 0; j < dim + 1; j++) {
            sideEl1 += j;
            sideEl2 += j;
          }
          for (int j = 0; j < dim; j++) {
            sideEl1 -= verticesEl1[j];
            sideEl2 -= verticesEl2[j];
          }
        }

        // create periodic info
        DimVec<WorldVector<double> > periodicCoordsEl1(dim - 1, NO_INIT);
        DimVec<WorldVector<double> > periodicCoordsEl2(dim - 1, NO_INIT);

        Element *element1 = const_cast<Element*>((*(mel + el1))->getElement());
        Element *element2 = const_cast<Element*>((*(mel + el2))->getElement());

        // for all vertices of this side
        for (int j = 0; j < dim; j++) {
          periodicCoordsEl1[element1->getPositionOfVertex(sideEl1, verticesEl1[j])] =
            coords[melVertex[el2][vertexMapEl1[verticesEl1[j]]]];

          periodicCoordsEl2[element2->getPositionOfVertex(sideEl2, verticesEl2[j])] =
            coords[melVertex[el1][vertexMapEl2[verticesEl2[j]]]];
        }

        // decorate leaf data
        ElementData *ld1 = element1->getElementData();
        ElementData *ld2 = element2->getElementData();

        TEST_EXIT_DBG(ld1)
          ("Should not happen: no element data pointer in macro element %d!\n",
          element1->getIndex());

        TEST_EXIT_DBG(ld2)
          ("Should not happen: no element data pointer in macro element %d!\n",
          element2->getIndex());

        LeafDataPeriodic *ldp1 =
          dynamic_cast<LeafDataPeriodic*>(ld1->getElementData(PERIODIC));
        LeafDataPeriodic *ldp2 =
          dynamic_cast<LeafDataPeriodic*>(ld2->getElementData(PERIODIC));

        if (!ldp1) {
          ldp1 = new LeafDataPeriodic(ld1);
          element1->setElementData(ldp1);
        }

        if (!ldp2) {
          ldp2 = new LeafDataPeriodic(ld2);
          element2->setElementData(ldp2);
        }

        ldp1->addPeriodicInfo(mode, boundaryType, sideEl1, &periodicCoordsEl1);
        ldp2->addPeriodicInfo(mode, boundaryType, sideEl2, &periodicCoordsEl2);

        if (mode != 0) {
          VertexVector* associated = NULL;
          DOFAdmin* localAdmin = NULL;

          if (mesh->periodicAssociations[boundaryType].empty()) {
            for (unsigned int iadmin = 0; iadmin < mesh->admin.size(); iadmin++) {
              localAdmin = mesh->admin[iadmin];
              TEST_EXIT_DBG(localAdmin->getNumberOfDofs(VERTEX) == 1)
                ("Periodic boundary supports only one dof on element vertex.\n");

              associated = new VertexVector(localAdmin, "vertex vector");
              mesh->periodicAssociations[boundaryType].push_back(associated);
              VertexVector::Iterator it(associated, ALL_DOFS);
              for (it.reset2(); !it.end(); ++it)
                *it = it.getDOFIndex();
            }
          }

          for (unsigned int iadmin = 0; iadmin < mesh->admin.size(); iadmin++) {
            associated = mesh->periodicAssociations[boundaryType][iadmin];

            for (int j = 0; j < dim; j++) {
241
#ifndef NDEBUG
242 243 244 245 246 247 248 249 250
              {
                unsigned initData(melVertex[el1][verticesEl1[j]]);
                unsigned oldData((*associated)[melVertex[el1][verticesEl1[j]]]);
                unsigned newData(melVertex[el2][vertexMapEl1[verticesEl1[j]]]);
                if( initData != oldData && newData != oldData ) {
                    MSG("warning: element %d overwrites assoc index %d: %d -> %d\n",
                            el1, initData, oldData, newData);
                }
              }
251
#endif
252 253 254 255 256 257 258
              (*associated)[melVertex[el1][verticesEl1[j]]] =
                melVertex[el2][vertexMapEl1[verticesEl1[j]]];
              (*associated)[melVertex[el2][verticesEl2[j]]] =
                melVertex[el1][vertexMapEl2[verticesEl2[j]]];
            }
          }
        }
259
      }
260

Thomas Witkowski's avatar
Thomas Witkowski committed
261 262
      delete [] verticesEl1;
      delete [] verticesEl2;
263 264

      // change periodic vertex dofs
265
      for (int i = 0; i < mesh->getNumberOfVertices(); i++) {
266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285
        if (periodicMap.getEntry(i) != -1) {
          mesh->freeDof(dof[i], VERTEX);
          dof[i] = dof[periodicMap.getEntry(i)];

          std::map<BoundaryType, std::vector<VertexVector*> >::iterator assoc;
          std::map<BoundaryType, std::vector<VertexVector*> >::iterator assocEnd =
            mesh->periodicAssociations.end();

          for (assoc = mesh->periodicAssociations.begin(); assoc != assocEnd; ++assoc) {
            for (unsigned int iadmin = 0; iadmin < mesh->admin.size(); iadmin++) {
              DegreeOfFreedom a = (*(assoc->second)[iadmin])[i];
              if (a != i) {
                (*(assoc->second)[iadmin])[i] = i;
                (*(assoc->second)[iadmin])[a] = periodicMap.getEntry(i);
              }
            } // for iadmin
          } // for assoc

        } // if
      } // for i
286

287
#ifndef NDEBUG
288 289 290
      std::map<BoundaryType, std::vector<VertexVector*> >::iterator assoc;
      std::map<BoundaryType, std::vector<VertexVector*> >::iterator assocEnd =
        mesh->periodicAssociations.end();
291
      for (assoc = mesh->periodicAssociations.begin();
292 293 294 295 296 297 298 299 300 301 302 303
           assoc != assocEnd; ++assoc)
        for (int i = 0; i < mesh->getNumberOfVertices(); i++)
          if (i != (*(assoc->second)[0])[i])
            MSG("association %d: vertex %d -> vertex %d\n",
            assoc->first, i, (*(assoc->second)[0])[i]);

      for (int i = 0; i < mesh->getNumberOfVertices(); i++) {
        if (periodicMap.getEntry(i) != -1) {
          MSG("identification : vertex %d is now vertex %d\n",
              i, periodicMap.getEntry(i));
        }
      }
304
#endif
305 306 307

    } // periodicFile

308

309
    // =========================================================
310
    // Set coordinate
311 312 313
    for (int i = 0; i < mesh->getNumberOfMacros(); i++) {
      for (int k = 0; k < mesh->getGeo(VERTEX); k++) {
	(*(mel + i))->setCoord(k, coords[melVertex[i][k]]);
314

315
	const_cast<Element*>((*(mel + i))->getElement())->
316
	  setDof(k, dof[melVertex[i][k]]);
317
      }
318
    }
319

320
    if (!macroInfo->neigh_set) {
321
      TEST_EXIT(periodicFile == "")
322
	("periodic boundary condition => element neighbours must be set\n");
323
      computeNeighbours(mesh);
324
    } else {
325 326 327
      /****************************************************************************/
      /* fill MEL oppVertex values when reading neighbour information form file  */
      /****************************************************************************/
328

329 330 331 332 333
      for (int i = 0; i < mesh->getNumberOfMacros(); i++) {
	for (int k = 0; k < mesh->getGeo(NEIGH); k++) {
	  MacroElement *neigh = const_cast<MacroElement*>(mel[i]->getNeighbour(k));

	  if (neigh) {
334 335 336
	    if (!macroInfo->neigh_inverse_set) {
	      int j = 0;
	      for (; j < mesh->getGeo(NEIGH); j++)
337
		if (neigh->getNeighbour(j) == *(mel + i))
338
		  break;
339 340

	      TEST_EXIT(j < mesh->getGeo(NEIGH))("el %d no neighbour of neighbour %d\n",
341 342 343 344 345
						mel[i]->getIndex(), neigh->getIndex());
	      mel[i]->setOppVertex(k, j);
	    } else {
	      int j = 0;
	      for (; j < mesh->getGeo(NEIGH); j++)
346
		if (neigh->getNeighbourInv(j) == *(mel + i))
347
		  break;
348 349

	      TEST_EXIT(j < mesh->getGeo(NEIGH))("el %d no neighbour-inv of neighbour %d\n",
350 351 352
						mel[i]->getIndex(), neigh->getIndex());
	      mel[i]->setOppVertex(k, j);
	    }
353 354
	  } else {
	    mel[i]->setOppVertex(k, -1);
355
	  }
356
	}
357
      }
358
    }
359

360
    if (!macroInfo->bound_set)
361 362
      macroInfo->dirichletBoundary();

363
    if (mesh->getDim() > 1)
364 365 366 367 368 369
      boundaryDOFs(mesh);

    // initial boundary projections
    int numFaces = mesh->getGeo(FACE);
    int dim = mesh->getDim();
    mel = mesh->firstMacroElement();
370
    for (int i = 0; i < mesh->getNumberOfLeaves(); i++) {
371
      MacroElement *macroEl = *(mel + i);
372
      Projection *projector = macroEl->getProjection(0);
373
      if (projector && projector->getType() == VOLUME_PROJECTION) {
Thomas Witkowski's avatar
Thomas Witkowski committed
374
	for (int j = 0; j <= dim; j++)
375
	  projector->project(macroEl->getCoord(j));
376
      } else {
377
	for (int j = 0; j < mesh->getGeo(EDGE); j++) {
378
	  projector = macroEl->getProjection(numFaces + j);
379
	  if (projector) {
380 381 382 383 384 385 386 387 388 389 390
	    int vertex0 = Global::getReferenceElement(dim)->getVertexOfEdge(j, 0);
	    int vertex1 = Global::getReferenceElement(dim)->getVertexOfEdge(j, 1);
	    projector->project(macroEl->getCoord(vertex0));
	    projector->project(macroEl->getCoord(vertex1));
	  }
	}
      }
    }

    macroInfo->fillBoundaryInfo(mesh);

391
    if (mesh->getNumberOfDofs(CENTER)) {
392
      for (int i = 0; i < mesh->getNumberOfMacros(); i++)
393
	const_cast<Element*>(mel[i]->getElement())->
394
	  setDof(mesh->getNode(CENTER), mesh->getDof(CENTER));
395 396
    }

397
    // === Domain size ===
398

399
    WorldVector<double> x_min, x_max;
400

401 402 403 404 405 406 407
    for (int j = 0; j < Global::getGeo(WORLD); j++) {
      x_min[j] =  1.E30;
      x_max[j] = -1.E30;
    }

    for (int i = 0; i < mesh->getNumberOfVertices(); i++) {
      for (int j = 0; j < Global::getGeo(WORLD); j++) {
408 409
	x_min[j] = std::min(x_min[j], coords[i][j]);
	x_max[j] = std::max(x_max[j], coords[i][j]);
410 411
      }
    }
412

413
    for (int j = 0; j < Global::getGeo(WORLD); j++)
414 415 416 417 418
      mesh->setDiameter(j, x_max[j] - x_min[j]);

    if (check) {
      checkMesh(mesh);

Thomas Witkowski's avatar
Thomas Witkowski committed
419 420
      if (mesh->getDim() > 1)
	macroTest(mesh);
421 422
    }

423
    return macroInfo;
424
  }
425

426 427
  void MacroReader::computeNeighbours(Mesh *mesh)
  {
428
    FUNCNAME("MacroReader::computeNeighbours()");
429

430
    int dim = mesh->getDim();
431
    FixVec<DegreeOfFreedom*, DIMEN> dof(dim, NO_INIT);
432

433 434 435
    for (int i = 0; i < mesh->getNumberOfLeaves(); i++) {
      for (int k = 0; k < mesh->getGeo(NEIGH); k++) {
	mesh->getMacroElement(i)->setOppVertex(k, AMDIS_UNDEFINED);
436
	mesh->getMacroElement(i)->setNeighbour(k, NULL);
437
      }
438
    }
439

440 441 442
    for (int i = 0; i < mesh->getNumberOfLeaves(); i++) {
      for (int k = 0; k < mesh->getGeo(NEIGH); k++) {
	if (mesh->getMacroElement(i)->getBoundary(k) != INTERIOR) {
443
	  mesh->getMacroElement(i)->setNeighbour(k, NULL);
444 445 446
	  mesh->getMacroElement(i)->setOppVertex(k, -1);
	  continue;
	}
447

448 449 450
	if (mesh->getMacroElement(i)->getOppVertex(k) == AMDIS_UNDEFINED) {
	  if (dim == 1) {
	    dof[0] = const_cast<DegreeOfFreedom*>(mesh->getMacroElement(i)->
451
						  getElement()->getDof(k));
452 453 454 455
	  } else {
	    for (int l = 0; l < dim; l++)
	      dof[l] = const_cast<DegreeOfFreedom*>(mesh->getMacroElement(i)->
						    getElement()->
456
						    getDof((k + l + 1) % (dim + 1)));
457
	  }
458

459 460 461 462 463 464 465 466 467 468 469 470
	  int j = 0;
	  for (j = i + 1; j < mesh->getNumberOfLeaves(); j++) {
	    int m = mesh->getMacroElement(j)->getElement()->oppVertex(dof);
	    if (m != -1) {
	      mesh->getMacroElement(i)->setNeighbour(k, mesh->getMacroElement(j));
	      mesh->getMacroElement(j)->setNeighbour(m, mesh->getMacroElement(i));
	      mesh->getMacroElement(i)->setOppVertex(k, m);
	      mesh->getMacroElement(j)->setOppVertex(m, k);
	      break;
	    }
	  }

471 472 473 474 475 476 477
	  if (j >= mesh->getNumberOfLeaves()) {
	    std::cout << "----------- ERROR ------------" << std::endl;
	    std::cout << "Cannot find neighbour " << k << " of element " << i << std::endl;
	    std::cout << "  dim = " << dim << std::endl;
	    std::cout << "  coords of element = ";
	    for (int l = 0; l <= dim; l++) {
	      std::cout << mesh->getMacroElement(i)->getCoord(l);
478
	      if (l < dim)
479
		std::cout << " / ";
480 481 482
	    }
	    std::cout << std::endl;
	    std::cout << "  dofs = ";
483
	    for (int l = 0; l < dim; l++)
484
	      std::cout << *(dof[l]) << " ";
485 486 487
	    std::cout << std::endl;

	    ERROR_EXIT("\n");
488
	  }
489
	}
490
      }
491
    }
492 493 494 495 496 497 498 499 500 501 502
  }


  /****************************************************************************/
  /*  boundaryDOFs:                                                           */
  /*  adds dof's at the edges of a given macro triangulation and calculates   */
  /*  the number of edges                                                     */
  /****************************************************************************/

  void MacroReader::boundaryDOFs(Mesh *mesh)
  {
503 504 505 506 507
    FUNCNAME("Mesh::boundaryDOFs()");

    int lnode = mesh->getNode(EDGE);
    int k, lne = mesh->getNumberOfLeaves();
    int max_n_neigh = 0, n_neigh, ov;
508
    std::deque<MacroElement*>::iterator mel;
509
    const MacroElement* neigh;
510
    DegreeOfFreedom *dof;
511 512 513 514 515 516

    mesh->setNumberOfEdges(0);
    mesh->setNumberOfFaces(0);

    int dim = mesh->getDim();

Thomas Witkowski's avatar
Thomas Witkowski committed
517
    switch (dim) {
518
    case 2:
519
      for (mel = mesh->firstMacroElement(); mel != mesh->endOfMacroElements(); mel++) {
520 521 522 523 524 525
	// check for periodic boundary
	Element *el = const_cast<Element*>((*mel)->getElement());
	ElementData *ed = el->getElementData(PERIODIC);

	DimVec<bool> periodic(dim, DEFAULT_VALUE, false);

526
	if (ed) {
527
	  std::list<LeafDataPeriodic::PeriodicInfo> &periodicInfos =
528
	    dynamic_cast<LeafDataPeriodic*>(ed)->getInfoList();
529 530
	  std::list<LeafDataPeriodic::PeriodicInfo>::iterator it;
	  std::list<LeafDataPeriodic::PeriodicInfo>::iterator end = periodicInfos.end();
Thomas Witkowski's avatar
Thomas Witkowski committed
531 532
	  for (it = periodicInfos.begin(); it != end; ++it)
	    if (it->type != 0)
533 534 535
	      periodic[it->elementSide] = true;
	}

536
	for (int i = 0; i < mesh->getGeo(NEIGH); i++) {
537
	  if (!(*mel)->getNeighbour(i) ||
538 539 540 541
	      ((*mel)->getNeighbour(i)->getIndex() < (*mel)->getIndex())) {

	    mesh->incrementNumberOfEdges(1);

542
	    if (mesh->getNumberOfDofs(EDGE)) {
543 544
	      dof = mesh->getDof(EDGE);
	      el->setDof(lnode + i, dof);
545

546
	      if ((*mel)->getNeighbour(i)) {
547
		Element *neigh =
Thomas Witkowski's avatar
Thomas Witkowski committed
548
		  const_cast<Element*>((*mel)->getNeighbour(i)->getElement());
549

Thomas Witkowski's avatar
Thomas Witkowski committed
550
		if (periodic[i])
551
		  neigh->setDof(lnode + (*mel)->getOppVertex(i), mesh->getDof(EDGE));
Thomas Witkowski's avatar
Thomas Witkowski committed
552
		else
553
		  neigh->setDof(lnode + (*mel)->getOppVertex(i), dof);
554
	      }
555
	    }
556
	  }
557 558 559 560 561 562
	}
      }
      break;
    case 3:
      lnode = mesh->getNode(FACE);
      mel = mesh->firstMacroElement();
563
      for (int i = 0; i < lne; i++) {
564 565 566

	// check for periodic boundary
	Element *el = const_cast<Element*>((*(mel+i))->getElement());
567

568 569 570
	ElementData *ed = el->getElementData(PERIODIC);

	DimVec<bool> periodic(dim, DEFAULT_VALUE, false);
571

572
	if (ed) {
573
	  std::list<LeafDataPeriodic::PeriodicInfo> &periodicInfos =
574
	    dynamic_cast<LeafDataPeriodic*>(ed)->getInfoList();
575 576
	  std::list<LeafDataPeriodic::PeriodicInfo>::iterator it;
	  std::list<LeafDataPeriodic::PeriodicInfo>::iterator end = periodicInfos.end();
577 578
	  for (it = periodicInfos.begin(); it != end; ++it)
	    if (it->type != 0)
579 580 581
	      periodic[it->elementSide] = true;
	}

582 583
	for (k = 0; k < mesh->getGeo(EDGE); k++) {
	  // === Check for not counted edges. ===
584 585
	  n_neigh = 1;

Thomas Witkowski's avatar
Thomas Witkowski committed
586
	  if (newEdge(mesh, (*(mel + i)), k, &n_neigh)) {
587
	    mesh->incrementNumberOfEdges(1);
588
	    max_n_neigh = std::max(max_n_neigh, n_neigh);
589 590
	  }
	}
591

592
	for (k = 0; k < mesh->getGeo(NEIGH); k++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
593
	  neigh = (*(mel + i))->getNeighbour(k);
594
	  // === Face is counted and dof is added by the element with bigger index. ===
595
	  if (neigh && (neigh->getIndex() > (*(mel + i))->getIndex()))
Thomas Witkowski's avatar
Thomas Witkowski committed
596
	    continue;
597

598
	  mesh->incrementNumberOfFaces(1);
599

600
	  if (mesh->getNumberOfDofs(FACE)) {
601
	    TEST_EXIT(!(*(mel + i))->getElement()->getDof(lnode + k))
602
	      ("dof %d on element %d already set\n",
Thomas Witkowski's avatar
Thomas Witkowski committed
603
	       lnode + k, (*(mel + i))->getIndex());
604

605
	    const_cast<Element*>((*(mel + i))->getElement())->setDof(lnode + k,
606
								     mesh->getDof(FACE));
607 608

	    if (neigh) {
Thomas Witkowski's avatar
Thomas Witkowski committed
609
	      ov = (*(mel + i))->getOppVertex(k);
610
	      TEST_EXIT(!neigh->getElement()->getDof(lnode + ov))
611
		("dof %d on neighbour %d already set\n",
Thomas Witkowski's avatar
Thomas Witkowski committed
612
		 lnode + ov, neigh->getIndex());
613 614

	      Element *neighEl =
Thomas Witkowski's avatar
Thomas Witkowski committed
615
		const_cast<Element*>((*(mel + i))->getNeighbour(k)->getElement());
616

Thomas Witkowski's avatar
Thomas Witkowski committed
617
	      if (periodic[k])
618
		neighEl->setDof(lnode+ov, mesh->getDof(FACE));
Thomas Witkowski's avatar
Thomas Witkowski committed
619
	      else
620
		neighEl->setDof(lnode+ov, const_cast<int*>((*(mel + i))->getElement()->
621
							   getDof(lnode + k)));
622 623 624 625 626
	    }
	  }
	}
      }
      break;
627
    default:
Thomas Witkowski's avatar
Thomas Witkowski committed
628
      ERROR_EXIT("invalid dim\n");
629
    }
630

631
    if (3 == dim)
Thomas Witkowski's avatar
Thomas Witkowski committed
632
      mesh->setMaxEdgeNeigh(std::max(8, 2 * max_n_neigh));
633
    else
634
      mesh->setMaxEdgeNeigh(dim - 1);
635 636
  }

637
  /*
638
     testet mesh auf auftretende Zyklen
639

640 641
     wenn Zyklus auftritt:
     ordnet Eintraege in MacroElement-Struktur um, so dass kein Zyklus auftritt
642
     erzeugt neue Macro-Datei nameneu mit umgeordnetem Netz
643
     (wenn nameneu=NULL wird keine MAcro-Datei erzeugt)
644
  */
645

Thomas Witkowski's avatar
Thomas Witkowski committed
646
  void MacroReader::macroTest(Mesh *mesh)
647
  {
648
    FUNCNAME("MacroReader::macroTest()");
649

650 651 652
    int i = macrotest(mesh);

    if (i >= 0) {
Thomas Witkowski's avatar
Thomas Witkowski committed
653 654
      WARNING("There is a cycle beginning in macro element %d\n", i);
      WARNING("Entries in MacroElement structures get reordered\n");
655
      umb(NULL, mesh, umbVkantMacro);
656
    }
657
  }
658

659 660 661 662 663
  /****************************************************************************/
  /*  macro_test():                              Author: Thomas Kastl (1998)  */
  /****************************************************************************/
  /*
    testet mesh auf auftretende Zyklen
664

665
    wenn mesh zyklenfrei -> -1
666
    sonst ->  globaler Index des Macroelementes bei dem ein Zyklus beginnt
667 668 669 670
  */

  int MacroReader::macrotest(Mesh *mesh)
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
671 672
    std::deque<MacroElement*>::const_iterator macro, mac;
    int flg = 0;
673
    int dim = mesh->getDim();
Thomas Witkowski's avatar
Thomas Witkowski committed
674 675
    int *test = new int[mesh->getNumberOfMacros()];
    int *zykl = new int[mesh->getNumberOfMacros()];
676

Thomas Witkowski's avatar
Thomas Witkowski committed
677
    for (int i = 0; i < mesh->getNumberOfMacros(); i++)
678
      test[i] = 0;
679

Thomas Witkowski's avatar
Thomas Witkowski committed
680 681
    int zykstart = -1;
    std::deque<MacroElement*>::const_iterator macrolfd = mesh->firstMacroElement();
682

683 684 685 686
    while (macrolfd != mesh->endOfMacroElements()) {
      if (test[(*macrolfd)->getIndex()] == 1) {
	macrolfd++;
      } else {
687
	for (int i = 0; i < mesh->getNumberOfMacros(); i++)
688 689
	  zykl[i] = 0;

690 691 692 693 694 695 696 697
	macro = macrolfd;
	flg = 2;
	do {
	  if (zykl[(*macro)->getIndex()] == 1) {
	    flg = 0;
	    zykstart = (*macro)->getIndex();
	  } else {
	    zykl[(*macro)->getIndex()] = 1;
698

699 700
	    if (test[(*macro)->getIndex()] == 1) {
	      flg = 1;
701
	    } else if ((*macro)->getNeighbour(dim) == NULL) {
702 703
	      flg = 1;
	      test[(*macro)->getIndex()] = 1;
Thomas Witkowski's avatar
Thomas Witkowski committed
704
	    } else if ((*macro) == (*macro)->getNeighbour(dim)->getNeighbour(dim)) {
705 706 707 708
	      flg = 1;
	      test[(*macro)->getIndex()] = 1;
	      test[(*macro)->getNeighbour(dim)->getIndex()] = 1;
	    } else {
709
	      for (mac = mesh->firstMacroElement();
Thomas Witkowski's avatar
Thomas Witkowski committed
710
		   (*mac) != (*macro)->getNeighbour(dim); mac++);
711
	      macro = mac;
712 713
	    }
	  }
714
	} while(flg == 2);
715

716
	if (flg == 1)
717
	  macrolfd++;
718 719
	else
	  macrolfd=mesh->endOfMacroElements();
720
      }
721
    }
722

Thomas Witkowski's avatar
Thomas Witkowski committed
723 724
    delete [] zykl;
    delete [] test;
725

726 727 728 729 730 731 732
    return zykstart;
  }

  //   waehlt geeignete Verfeinerungskanten, so dass kein Zyklus auftritt (recumb)

  //   ele     Integer-Vektor der Dimension Anzahl der Macro-Elemente
  //           zur Speicherung der neuen Verfeinerungskanten
733 734
  //           (wird nur benoetigt, wenn umbvk=umb_vkant_macrodat)

735 736 737 738 739
  //   umbvk   Fkt. zur Neuordnung der Verfeinerungskanten
  //           = umb_vkant_macro :
  //               Eintraege in MacroElement-Struktur und Eintraege in macro->el
  //               werden tatsaechlich umgeordnet
  //               -> ALBERT-Routine write_macro kann zum erzeugen einer
740
  //                  neuen Macro-Datei angewendet werden
741 742 743 744 745 746 747 748
  //           = umb_vkant_macrodat :
  //               Eintraege werden nicht veraendert, es werden nur die lokalen
  //               Indices der Kanten, die zu Verfeinerungskanten werden im
  //               Integer-Vektor ele abgespeichert
  //               -> print_Macrodat zur Erzeugung einer zyklenfreien neuen
  //                  Macro-Datei kann angewendet werden

  void MacroReader::umb(int *ele, Mesh *mesh,
Thomas Witkowski's avatar
Thomas Witkowski committed
749
			void (*umbvk)(Mesh*, MacroElement*, int, int*))
750
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
751
    int *test = new int[mesh->getNumberOfMacros()];
752

Thomas Witkowski's avatar
Thomas Witkowski committed
753 754
    for (int i = 0; i < static_cast<int>(mesh->getNumberOfMacros()); i++)
      test[i] = 0;
755

756
    recumb(mesh, (*mesh->firstMacroElement()), NULL, test, 0, 0, ele, umbvk);
757

Thomas Witkowski's avatar
Thomas Witkowski committed
758
    delete [] test;
759 760 761 762 763
  }

  bool MacroReader::newEdge(Mesh *mesh, MacroElement *mel,
			    int mel_edge_no, int *n_neigh)
  {
764
    FUNCNAME("MacroElement::newEdge()");
Thomas Witkowski's avatar
Thomas Witkowski committed
765
    MacroElement *nei;
766
    const DegreeOfFreedom *dof[2];
767
    DegreeOfFreedom *edge_dof = NULL;
Thomas Witkowski's avatar
Thomas Witkowski committed
768 769
    int j, k, opp_v, node = 0;
    BoundaryType lbound = INTERIOR;
770
    Projection *lproject = NULL;
Thomas Witkowski's avatar
Thomas Witkowski committed
771
    const int max_no_list_el = 100;
772 773 774 775
    BoundaryType *list_bound[100];
    Projection **list_project[100];
    Element *el = const_cast<Element*>(mel->getElement());
    int edge_no = mel_edge_no;
Thomas Witkowski's avatar
Thomas Witkowski committed
776
    static int next_el[6][2] = {{3,2},{1,3},{1,2},{0,3},{0,2},{0,1}};
777 778
    int vertices = mesh->getGeo(VERTEX);

Thomas Witkowski's avatar
Thomas Witkowski committed
779
    int mel_index = mel->getIndex();
780 781 782 783

    list_bound[0] = &(mel->boundary[mesh->getGeo(FACE)+edge_no]);
    list_project[0] = &(mel->projection[mesh->getGeo(FACE)+edge_no]);

784
    if (mesh->getNumberOfDofs(EDGE)) {
785
      node = mesh->getNode(EDGE);
786
      if (el->getDof(node+edge_no)) {
787 788 789 790 791 792
	/****************************************************************************/
	/*  edge was counted by another macro element and dof was added on the      */
	/*  complete patch                                                          */
	/****************************************************************************/
	return false;
      } else {
793 794
	edge_dof = mesh->getDof(EDGE);
	el->setDof(node+edge_no, edge_dof);
795 796 797
      }
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
798
    for (j = 0; j < 2; j++)
799
      dof[j] = el->getDof(el->getVertexOfEdge(edge_no, j));
800

801 802 803 804 805 806 807 808 809 810 811 812
    /****************************************************************************/
    /*  first look for all neighbours in one direction until a boundary is      */
    /*  reached :-( or we are back to mel :-)                                   */
    /*  if the index of a neighbour is smaller than the element index, the edge */
    /*  is counted by this neighbour, return 0.                                 */
    /*  If we are back to element, return 1, to count the edge                  */
    /****************************************************************************/

    nei = mel->getNeighbour(next_el[edge_no][0]);
    opp_v = mel->getOppVertex(next_el[edge_no][0]);


813
    if (mel->getBoundary(next_el[edge_no][0])) {
814 815 816 817 818 819
      lbound = newBound(mel->getBoundary(next_el[edge_no][0]), lbound);
      lproject = mel->getProjection(next_el[edge_no][0]);
    }

    while (nei  &&  nei != mel) {
      for (j = 0; j < vertices; j++)
820
	if (nei->getElement()->getDof(j) == dof[0])  break;
821
      for (k = 0; k < vertices; k++)
822
	if (nei->getElement()->getDof(k) == dof[1])  break;
823 824

      // check for periodic boundary
825
      if (j == 4 || k == 4) {
826
	nei = NULL;
827 828 829
	break;
      }

830
      if (mesh->getNumberOfDofs(EDGE)) {
831 832
	TEST_EXIT(nei->index > mel_index)
	  ("neighbour index %d < element index %d\n", nei->getIndex(), mel_index);
833
      }
834

835
      if (!mesh->getNumberOfDofs(EDGE) && nei->getIndex() < mel_index)
836
	return false;
837

838
      edge_no = Tetrahedron::edgeOfDofs[j][k];
839 840 841 842

      TEST_EXIT(*n_neigh < max_no_list_el)
	("too many neigbours for local list\n");

843
      list_bound[(*n_neigh)] =
844 845
	&(nei->boundary[mesh->getGeo(FACE)+edge_no]);

846
      list_project[(*n_neigh)++] =
847 848
	&(nei->projection[mesh->getGeo(FACE)+edge_no]);

849
      if (mesh->getNumberOfDofs(EDGE))
850
	nei->element->setDof(node+edge_no,edge_dof);
851

852
      if (next_el[edge_no][0] != opp_v) {
853
	if (nei->getBoundary(next_el[edge_no][0])) {
854 855 856 857 858 859
	  lbound = newBound(nei->getBoundary(next_el[edge_no][0]), lbound);
	  Projection *neiProject = nei->getProjection(next_el[edge_no][0]);
	  if (!lproject) {
	    lproject = neiProject;
	  } else {
	    if (neiProject && (lproject->getID() < neiProject->getID()))
860 861 862
	      lproject = neiProject;
	  }
	}
863 864 865 866 867
	opp_v = nei->getOppVertex(next_el[edge_no][0]);
	nei = nei->getNeighbour(next_el[edge_no][0]);
      } else {
	if (nei->getBoundary(next_el[edge_no][1])) {
	  lbound = newBound(nei->getBoundary(next_el[edge_no][1]), lbound);
868
	  Projection *neiProject = nei->getProjection(next_el[edge_no][1]);
869 870 871 872
	  if (!lproject) {
	    lproject = neiProject;
	  } else {
	    if (neiProject && (lproject->getID() < neiProject->getID()))
873 874 875
	      lproject = neiProject;
	  }
	}
876 877 878
	opp_v = nei->getOppVertex(next_el[edge_no][1]);
	nei = nei->getNeighbour(next_el[edge_no][1]);
      }
879 880
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
881 882 883 884 885 886
    if (!nei) {
      /****************************************************************************/
      /*  while looping around the edge the domain's boundary was reached. Now,   */
      /*  loop in the other direction to the domain's boundary		    */
      /****************************************************************************/
      edge_no = mel_edge_no;
887

Thomas Witkowski's avatar
Thomas Witkowski committed
888 889 890
      nei = mel->getNeighbour(next_el[edge_no][1]);
      opp_v = mel->getOppVertex(next_el[edge_no][1]);
      if (mel->getBoundary(next_el[edge_no][1])) {
891
	lbound = newBound(mel->getBoundary(next_el[edge_no][1]), lbound);
Thomas Witkowski's avatar
Thomas Witkowski committed
892 893 894 895 896
	Projection *neiProject =  mel->getProjection(next_el[edge_no][1]);
	if (!lproject) {
	  lproject = neiProject;
	} else {
	  if (neiProject && (lproject->getID() < neiProject->getID()))
897
	    lproject = neiProject;
Thomas Witkowski's avatar
Thomas Witkowski committed
898 899
	}
      }
900

Thomas Witkowski's avatar
Thomas Witkowski committed
901 902
      while (nei) {
	for (j = 0; j < vertices; j++)
903
	  if (nei->getElement()->getDof(j) == dof[0])  break;
Thomas Witkowski's avatar
Thomas Witkowski committed
904
	for (k = 0; k < vertices; k++)
905
	  if (nei->getElement()->getDof(k) == dof[1])  break;
906

Thomas Witkowski's avatar
Thomas Witkowski committed
907 908 909
	// check for periodic boundary
	if (j == 4 || k == 4)
	  return false;
910

911
	if (mesh->getNumberOfDofs(EDGE)) {
Thomas Witkowski's avatar
Thomas Witkowski committed
912 913 914
	  TEST_EXIT(nei->getIndex() > mel_index)
	    ("neighbour index %d < element index %d\n", nei->getIndex(),
	     mel_index);
915
	}
916 917

	if (nei->getIndex() < mel_index)
Thomas Witkowski's avatar
Thomas Witkowski committed
918
	  return false;
919

920
	edge_no = Tetrahedron::edgeOfDofs[j][k];
921

Thomas Witkowski's avatar
Thomas Witkowski committed
922
	TEST_EXIT(*n_neigh < max_no_list_el)("too many neigbours for local list\n");
923

Thomas Witkowski's avatar
Thomas Witkowski committed
924 925
	list_bound[(*n_neigh)] = &(nei->boundary[mesh->getGeo(FACE) + edge_no]);
	list_project[(*n_neigh)++] = &(nei->projection[mesh->getGeo(FACE) + edge_no]);
926

927
	if (mesh->getNumberOfDofs(EDGE)) {
928
	  TEST_EXIT(!nei->getElement()->getDof(node+edge_no))
Thomas Witkowski's avatar
Thomas Witkowski committed
929 930
	    ("dof %d on element %d is already set, but not on element %d\n",
	     node + edge_no, nei->getIndex(), mel_index);
931

932
	  nei->element->setDof(node+edge_no, edge_dof);
Thomas Witkowski's avatar
Thomas Witkowski committed
933
	}
934

Thomas Witkowski's avatar
Thomas Witkowski committed
935 936 937 938 939 940 941 942
	if (next_el[edge_no][0] != opp_v) {
	  if (nei->getBoundary(next_el[edge_no][0])) {
	    lbound = newBound(nei->getBoundary(next_el[edge_no][0]), lbound);
	    Projection *neiProject = nei->getProjection(next_el[edge_no][0]);
	    if (!lproject) {
	      lproject = neiProject;
	    } else {
	      if (neiProject &&( lproject->getID() < neiProject->getID()))
943
		lproject = neiProject;
944 945
	    }
	  }
946

Thomas Witkowski's avatar
Thomas Witkowski committed
947 948 949 950
	  opp_v = nei->getOppVertex(next_el[edge_no][0]);
	  nei = nei->getNeighbour(next_el[edge_no][0]);
	} else {
	  if (nei->getBoundary(next_el[edge_no][1])) {
951
	    lbound = newBound(nei->getBoundary(next_el[edge_no][1]), lbound);
Thomas Witkowski's avatar
Thomas Witkowski committed
952 953 954 955 956
	    Projection *neiProject = nei->getProjection(next_el[edge_no][1]);
	    if (!lproject) {
	      lproject = neiProject;
	    } else {
	      if (neiProject && (lproject->getID() < neiProject->getID()))
957
		lproject = neiProject;
Thomas Witkowski's avatar
Thomas Witkowski committed
958 959
	    }
	  }
960

Thomas Witkowski's avatar
Thomas Witkowski committed
961 962 963
	  opp_v = nei->getOppVertex(next_el[edge_no][1]);
	  nei = nei->getNeighbour(next_el[edge_no][1]);
	}
964
      }