MacroReader.cc 49.1 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
<<<<<<< HEAD
242
#ifndef NDEBUG
243
=======
244
#ifdef DEBUG
245
>>>>>>> 7513d640d9de009c53e95dd9e1f252958172b31d
246
247
248
249
250
251
252
253
254
              {
                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);
                }
              }
255
#endif
256
257
258
259
260
261
262
              (*associated)[melVertex[el1][verticesEl1[j]]] =
                melVertex[el2][vertexMapEl1[verticesEl1[j]]];
              (*associated)[melVertex[el2][verticesEl2[j]]] =
                melVertex[el1][vertexMapEl2[verticesEl2[j]]];
            }
          }
        }
263
      }
264

Thomas Witkowski's avatar
Thomas Witkowski committed
265
266
      delete [] verticesEl1;
      delete [] verticesEl2;
267
268

      // change periodic vertex dofs
269
      for (int i = 0; i < mesh->getNumberOfVertices(); i++) {
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
        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
290

291
<<<<<<< HEAD
292
#ifndef NDEBUG
293
=======
Thomas Witkowski's avatar
Thomas Witkowski committed
294
#if (DEBUG != 0)
295
>>>>>>> 7513d640d9de009c53e95dd9e1f252958172b31d
296
297
298
      std::map<BoundaryType, std::vector<VertexVector*> >::iterator assoc;
      std::map<BoundaryType, std::vector<VertexVector*> >::iterator assocEnd =
        mesh->periodicAssociations.end();
299
      for (assoc = mesh->periodicAssociations.begin();
300
301
302
303
304
305
306
307
308
309
310
311
           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));
        }
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
312
#endif
313
314
315

    } // periodicFile

316

317
    // =========================================================
318
    // Set coordinate
319
320
321
    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]]);
322

323
	const_cast<Element*>((*(mel + i))->getElement())->
324
	  setDof(k, dof[melVertex[i][k]]);
325
      }
326
    }
327

328
    if (!macroInfo->neigh_set) {
329
      TEST_EXIT(periodicFile == "")
330
	("periodic boundary condition => element neighbours must be set\n");
331
      computeNeighbours(mesh);
332
    } else {
333
334
335
      /****************************************************************************/
      /* fill MEL oppVertex values when reading neighbour information form file  */
      /****************************************************************************/
336

337
338
339
340
341
      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) {
342
343
344
	    if (!macroInfo->neigh_inverse_set) {
	      int j = 0;
	      for (; j < mesh->getGeo(NEIGH); j++)
345
		if (neigh->getNeighbour(j) == *(mel + i))
346
		  break;
347
348

	      TEST_EXIT(j < mesh->getGeo(NEIGH))("el %d no neighbour of neighbour %d\n",
349
350
351
352
353
						mel[i]->getIndex(), neigh->getIndex());
	      mel[i]->setOppVertex(k, j);
	    } else {
	      int j = 0;
	      for (; j < mesh->getGeo(NEIGH); j++)
354
		if (neigh->getNeighbourInv(j) == *(mel + i))
355
		  break;
356
357

	      TEST_EXIT(j < mesh->getGeo(NEIGH))("el %d no neighbour-inv of neighbour %d\n",
358
359
360
						mel[i]->getIndex(), neigh->getIndex());
	      mel[i]->setOppVertex(k, j);
	    }
361
362
	  } else {
	    mel[i]->setOppVertex(k, -1);
363
	  }
364
	}
365
      }
366
    }
367

368
    if (!macroInfo->bound_set)
369
370
      macroInfo->dirichletBoundary();

371
    if (mesh->getDim() > 1)
372
373
374
375
376
377
      boundaryDOFs(mesh);

    // initial boundary projections
    int numFaces = mesh->getGeo(FACE);
    int dim = mesh->getDim();
    mel = mesh->firstMacroElement();
378
    for (int i = 0; i < mesh->getNumberOfLeaves(); i++) {
379
      MacroElement *macroEl = *(mel + i);
380
      Projection *projector = macroEl->getProjection(0);
381
      if (projector && projector->getType() == VOLUME_PROJECTION) {
Thomas Witkowski's avatar
Thomas Witkowski committed
382
	for (int j = 0; j <= dim; j++)
383
	  projector->project(macroEl->getCoord(j));
384
      } else {
385
	for (int j = 0; j < mesh->getGeo(EDGE); j++) {
386
	  projector = macroEl->getProjection(numFaces + j);
387
	  if (projector) {
388
389
390
391
392
393
394
395
396
397
398
	    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);

399
    if (mesh->getNumberOfDofs(CENTER)) {
400
      for (int i = 0; i < mesh->getNumberOfMacros(); i++)
401
	const_cast<Element*>(mel[i]->getElement())->
402
	  setDof(mesh->getNode(CENTER), mesh->getDof(CENTER));
403
404
    }

405
    // === Domain size ===
406

407
    WorldVector<double> x_min, x_max;
408

409
410
411
412
413
414
415
    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++) {
416
417
	x_min[j] = std::min(x_min[j], coords[i][j]);
	x_max[j] = std::max(x_max[j], coords[i][j]);
418
419
      }
    }
420

421
    for (int j = 0; j < Global::getGeo(WORLD); j++)
422
423
424
425
426
      mesh->setDiameter(j, x_max[j] - x_min[j]);

    if (check) {
      checkMesh(mesh);

Thomas Witkowski's avatar
Thomas Witkowski committed
427
428
      if (mesh->getDim() > 1)
	macroTest(mesh);
429
430
    }

431
    return macroInfo;
432
  }
433

434
435
  void MacroReader::computeNeighbours(Mesh *mesh)
  {
436
    FUNCNAME("MacroReader::computeNeighbours()");
437

438
    int dim = mesh->getDim();
439
    FixVec<DegreeOfFreedom*, DIMEN> dof(dim, NO_INIT);
440

441
442
443
    for (int i = 0; i < mesh->getNumberOfLeaves(); i++) {
      for (int k = 0; k < mesh->getGeo(NEIGH); k++) {
	mesh->getMacroElement(i)->setOppVertex(k, AMDIS_UNDEFINED);
444
	mesh->getMacroElement(i)->setNeighbour(k, NULL);
445
      }
446
    }
447

448
449
450
    for (int i = 0; i < mesh->getNumberOfLeaves(); i++) {
      for (int k = 0; k < mesh->getGeo(NEIGH); k++) {
	if (mesh->getMacroElement(i)->getBoundary(k) != INTERIOR) {
451
	  mesh->getMacroElement(i)->setNeighbour(k, NULL);
452
453
454
	  mesh->getMacroElement(i)->setOppVertex(k, -1);
	  continue;
	}
455

456
457
458
	if (mesh->getMacroElement(i)->getOppVertex(k) == AMDIS_UNDEFINED) {
	  if (dim == 1) {
	    dof[0] = const_cast<DegreeOfFreedom*>(mesh->getMacroElement(i)->
459
						  getElement()->getDof(k));
460
461
462
463
	  } else {
	    for (int l = 0; l < dim; l++)
	      dof[l] = const_cast<DegreeOfFreedom*>(mesh->getMacroElement(i)->
						    getElement()->
464
						    getDof((k + l + 1) % (dim + 1)));
465
	  }
466

467
468
469
470
471
472
473
474
475
476
477
478
	  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;
	    }
	  }

479
480
481
482
483
484
485
	  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);
486
	      if (l < dim)
487
		std::cout << " / ";
488
489
490
	    }
	    std::cout << std::endl;
	    std::cout << "  dofs = ";
491
	    for (int l = 0; l < dim; l++)
492
	      std::cout << *(dof[l]) << " ";
493
494
495
	    std::cout << std::endl;

	    ERROR_EXIT("\n");
496
	  }
497
	}
498
      }
499
    }
500
501
502
503
504
505
506
507
508
509
510
  }


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

  void MacroReader::boundaryDOFs(Mesh *mesh)
  {
511
512
513
514
515
    FUNCNAME("Mesh::boundaryDOFs()");

    int lnode = mesh->getNode(EDGE);
    int k, lne = mesh->getNumberOfLeaves();
    int max_n_neigh = 0, n_neigh, ov;
516
    std::deque<MacroElement*>::iterator mel;
517
    const MacroElement* neigh;
518
    DegreeOfFreedom *dof;
519
520
521
522
523
524

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

    int dim = mesh->getDim();

Thomas Witkowski's avatar
Thomas Witkowski committed
525
    switch (dim) {
526
    case 2:
527
      for (mel = mesh->firstMacroElement(); mel != mesh->endOfMacroElements(); mel++) {
528
529
530
531
532
533
	// check for periodic boundary
	Element *el = const_cast<Element*>((*mel)->getElement());
	ElementData *ed = el->getElementData(PERIODIC);

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

534
	if (ed) {
535
	  std::list<LeafDataPeriodic::PeriodicInfo> &periodicInfos =
536
	    dynamic_cast<LeafDataPeriodic*>(ed)->getInfoList();
537
538
	  std::list<LeafDataPeriodic::PeriodicInfo>::iterator it;
	  std::list<LeafDataPeriodic::PeriodicInfo>::iterator end = periodicInfos.end();
Thomas Witkowski's avatar
Thomas Witkowski committed
539
540
	  for (it = periodicInfos.begin(); it != end; ++it)
	    if (it->type != 0)
541
542
543
	      periodic[it->elementSide] = true;
	}

544
	for (int i = 0; i < mesh->getGeo(NEIGH); i++) {
545
	  if (!(*mel)->getNeighbour(i) ||
546
547
548
549
	      ((*mel)->getNeighbour(i)->getIndex() < (*mel)->getIndex())) {

	    mesh->incrementNumberOfEdges(1);

550
	    if (mesh->getNumberOfDofs(EDGE)) {
551
552
	      dof = mesh->getDof(EDGE);
	      el->setDof(lnode + i, dof);
553

554
	      if ((*mel)->getNeighbour(i)) {
555
		Element *neigh =
Thomas Witkowski's avatar
Thomas Witkowski committed
556
		  const_cast<Element*>((*mel)->getNeighbour(i)->getElement());
557

Thomas Witkowski's avatar
Thomas Witkowski committed
558
		if (periodic[i])
559
		  neigh->setDof(lnode + (*mel)->getOppVertex(i), mesh->getDof(EDGE));
Thomas Witkowski's avatar
Thomas Witkowski committed
560
		else
561
		  neigh->setDof(lnode + (*mel)->getOppVertex(i), dof);
562
	      }
563
	    }
564
	  }
565
566
567
568
569
570
	}
      }
      break;
    case 3:
      lnode = mesh->getNode(FACE);
      mel = mesh->firstMacroElement();
571
      for (int i = 0; i < lne; i++) {
572
573
574

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

576
577
578
	ElementData *ed = el->getElementData(PERIODIC);

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

580
	if (ed) {
581
	  std::list<LeafDataPeriodic::PeriodicInfo> &periodicInfos =
582
	    dynamic_cast<LeafDataPeriodic*>(ed)->getInfoList();
583
584
	  std::list<LeafDataPeriodic::PeriodicInfo>::iterator it;
	  std::list<LeafDataPeriodic::PeriodicInfo>::iterator end = periodicInfos.end();
585
586
	  for (it = periodicInfos.begin(); it != end; ++it)
	    if (it->type != 0)
587
588
589
	      periodic[it->elementSide] = true;
	}

590
591
	for (k = 0; k < mesh->getGeo(EDGE); k++) {
	  // === Check for not counted edges. ===
592
593
	  n_neigh = 1;

Thomas Witkowski's avatar
Thomas Witkowski committed
594
	  if (newEdge(mesh, (*(mel + i)), k, &n_neigh)) {
595
	    mesh->incrementNumberOfEdges(1);
596
	    max_n_neigh = std::max(max_n_neigh, n_neigh);
597
598
	  }
	}
599

600
	for (k = 0; k < mesh->getGeo(NEIGH); k++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
601
	  neigh = (*(mel + i))->getNeighbour(k);
602
	  // === Face is counted and dof is added by the element with bigger index. ===
603
	  if (neigh && (neigh->getIndex() > (*(mel + i))->getIndex()))
Thomas Witkowski's avatar
Thomas Witkowski committed
604
	    continue;
605

606
	  mesh->incrementNumberOfFaces(1);
607

608
	  if (mesh->getNumberOfDofs(FACE)) {
609
	    TEST_EXIT(!(*(mel + i))->getElement()->getDof(lnode + k))
610
	      ("dof %d on element %d already set\n",
Thomas Witkowski's avatar
Thomas Witkowski committed
611
	       lnode + k, (*(mel + i))->getIndex());
612

613
	    const_cast<Element*>((*(mel + i))->getElement())->setDof(lnode + k,
614
								     mesh->getDof(FACE));
615
616

	    if (neigh) {
Thomas Witkowski's avatar
Thomas Witkowski committed
617
	      ov = (*(mel + i))->getOppVertex(k);
618
	      TEST_EXIT(!neigh->getElement()->getDof(lnode + ov))
619
		("dof %d on neighbour %d already set\n",
Thomas Witkowski's avatar
Thomas Witkowski committed
620
		 lnode + ov, neigh->getIndex());
621
622

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

Thomas Witkowski's avatar
Thomas Witkowski committed
625
	      if (periodic[k])
626
		neighEl->setDof(lnode+ov, mesh->getDof(FACE));
Thomas Witkowski's avatar
Thomas Witkowski committed
627
	      else
628
		neighEl->setDof(lnode+ov, const_cast<int*>((*(mel + i))->getElement()->
629
							   getDof(lnode + k)));
630
631
632
633
634
	    }
	  }
	}
      }
      break;
635
    default:
Thomas Witkowski's avatar
Thomas Witkowski committed
636
      ERROR_EXIT("invalid dim\n");
637
    }
638

639
    if (3 == dim)
Thomas Witkowski's avatar
Thomas Witkowski committed
640
      mesh->setMaxEdgeNeigh(std::max(8, 2 * max_n_neigh));
641
    else
642
      mesh->setMaxEdgeNeigh(dim - 1);
643
644
  }

645
  /*
646
     testet mesh auf auftretende Zyklen
647

648
649
     wenn Zyklus auftritt:
     ordnet Eintraege in MacroElement-Struktur um, so dass kein Zyklus auftritt
650
     erzeugt neue Macro-Datei nameneu mit umgeordnetem Netz
651
     (wenn nameneu=NULL wird keine MAcro-Datei erzeugt)
652
  */
653

Thomas Witkowski's avatar
Thomas Witkowski committed
654
  void MacroReader::macroTest(Mesh *mesh)
655
  {
656
    FUNCNAME("MacroReader::macroTest()");
657

658
659
660
    int i = macrotest(mesh);

    if (i >= 0) {
Thomas Witkowski's avatar
Thomas Witkowski committed
661
662
      WARNING("There is a cycle beginning in macro element %d\n", i);
      WARNING("Entries in MacroElement structures get reordered\n");
663
      umb(NULL, mesh, umbVkantMacro);
664
    }
665
  }
666

667
668
669
670
671
  /****************************************************************************/
  /*  macro_test():                              Author: Thomas Kastl (1998)  */
  /****************************************************************************/
  /*
    testet mesh auf auftretende Zyklen
672

673
    wenn mesh zyklenfrei -> -1
674
    sonst ->  globaler Index des Macroelementes bei dem ein Zyklus beginnt
675
676
677
678
  */

  int MacroReader::macrotest(Mesh *mesh)
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
679
680
    std::deque<MacroElement*>::const_iterator macro, mac;
    int flg = 0;
681
    int dim = mesh->getDim();
Thomas Witkowski's avatar
Thomas Witkowski committed
682
683
    int *test = new int[mesh->getNumberOfMacros()];
    int *zykl = new int[mesh->getNumberOfMacros()];
684

Thomas Witkowski's avatar
Thomas Witkowski committed
685
    for (int i = 0; i < mesh->getNumberOfMacros(); i++)
686
      test[i] = 0;
687

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

691
692
693
694
    while (macrolfd != mesh->endOfMacroElements()) {
      if (test[(*macrolfd)->getIndex()] == 1) {
	macrolfd++;
      } else {
695
	for (int i = 0; i < mesh->getNumberOfMacros(); i++)
696
697
	  zykl[i] = 0;

698
699
700
701
702
703
704
705
	macro = macrolfd;
	flg = 2;
	do {
	  if (zykl[(*macro)->getIndex()] == 1) {
	    flg = 0;
	    zykstart = (*macro)->getIndex();
	  } else {
	    zykl[(*macro)->getIndex()] = 1;
706

707
708
	    if (test[(*macro)->getIndex()] == 1) {
	      flg = 1;
709
	    } else if ((*macro)->getNeighbour(dim) == NULL) {
710
711
	      flg = 1;
	      test[(*macro)->getIndex()] = 1;
Thomas Witkowski's avatar
Thomas Witkowski committed
712
	    } else if ((*macro) == (*macro)->getNeighbour(dim)->getNeighbour(dim)) {
713
714
715
716
	      flg = 1;
	      test[(*macro)->getIndex()] = 1;
	      test[(*macro)->getNeighbour(dim)->getIndex()] = 1;
	    } else {
717
	      for (mac = mesh->firstMacroElement();
Thomas Witkowski's avatar
Thomas Witkowski committed
718
		   (*mac) != (*macro)->getNeighbour(dim); mac++);
719
	      macro = mac;
720
721
	    }
	  }
722
	} while(flg == 2);
723

724
	if (flg == 1)
725
	  macrolfd++;
726
727
	else
	  macrolfd=mesh->endOfMacroElements();
728
      }
729
    }
730

Thomas Witkowski's avatar
Thomas Witkowski committed
731
732
    delete [] zykl;
    delete [] test;
733

734
735
736
737
738
739
740
    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
741
742
  //           (wird nur benoetigt, wenn umbvk=umb_vkant_macrodat)

743
744
745
746
747
  //   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
748
  //                  neuen Macro-Datei angewendet werden
749
750
751
752
753
754
755
756
  //           = 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
757
			void (*umbvk)(Mesh*, MacroElement*, int, int*))
758
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
759
    int *test = new int[mesh->getNumberOfMacros()];
760

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

764
    recumb(mesh, (*mesh->firstMacroElement()), NULL, test, 0, 0, ele, umbvk);
765

Thomas Witkowski's avatar
Thomas Witkowski committed
766
    delete [] test;
767
768
769
770
771
  }

  bool MacroReader::newEdge(Mesh *mesh, MacroElement *mel,
			    int mel_edge_no, int *n_neigh)
  {
772
    FUNCNAME("MacroElement::newEdge()");
Thomas Witkowski's avatar
Thomas Witkowski committed
773
    MacroElement *nei;
774
    const DegreeOfFreedom *dof[2];
775
    DegreeOfFreedom *edge_dof = NULL;
Thomas Witkowski's avatar
Thomas Witkowski committed
776
777
    int j, k, opp_v, node = 0;
    BoundaryType lbound = INTERIOR;
778
    Projection *lproject = NULL;
Thomas Witkowski's avatar
Thomas Witkowski committed
779
    const int max_no_list_el = 100;
780
781
782
783
    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
784
    static int next_el[6][2] = {{3,2},{1,3},{1,2},{0,3},{0,2},{0,1}};
785
786
    int vertices = mesh->getGeo(VERTEX);

Thomas Witkowski's avatar
Thomas Witkowski committed
787
    int mel_index = mel->getIndex();
788
789
790
791

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

792
    if (mesh->getNumberOfDofs(EDGE)) {
793
      node = mesh->getNode(EDGE);
794
      if (el->getDof(node+edge_no)) {
795
796
797
798
799
800
	/****************************************************************************/
	/*  edge was counted by another macro element and dof was added on the      */
	/*  complete patch                                                          */
	/****************************************************************************/
	return false;
      } else {
801
802
	edge_dof = mesh->getDof(EDGE);
	el->setDof(node+edge_no, edge_dof);
803
804
805
      }
    }

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

809
810
811
812
813
814
815
816
817
818
819
820
    /****************************************************************************/
    /*  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]);


821
    if (mel->getBoundary(next_el[edge_no][0])) {
822
823
824
825
826
827
      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++)
828
	if (nei->getElement()->getDof(j) == dof[0])  break;
829
      for (k = 0; k < vertices; k++)
830
	if (nei->getElement()->getDof(k) == dof[1])  break;
831
832

      // check for periodic boundary
833
      if (j == 4 || k == 4) {
834
	nei = NULL;
835
836
837
	break;
      }

838
      if (mesh->getNumberOfDofs(EDGE)) {
839
840
	TEST_EXIT(nei->index > mel_index)
	  ("neighbour index %d < element index %d\n", nei->getIndex(), mel_index);
841
      }
842

843
      if (!mesh->getNumberOfDofs(EDGE) && nei->getIndex() < mel_index)
844
	return false;
845

846
      edge_no = Tetrahedron::edgeOfDofs[j][k];
847
848
849
850

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

851
      list_bound[(*n_neigh)] =
852
853
	&(nei->boundary[mesh->getGeo(FACE)+edge_no]);

854
      list_project[(*n_neigh)++] =
855
856
	&(nei->projection[mesh->getGeo(FACE)+edge_no]);

857
      if (mesh->getNumberOfDofs(EDGE))
858
	nei->element->setDof(node+edge_no,edge_dof);
859

860
      if (next_el[edge_no][0] != opp_v) {
861
	if (nei->getBoundary(next_el[edge_no][0])) {
862
863
864
865
866
867
	  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()))
868
869
870