MacroReader.cc 47.9 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
/******************************************************************************
 *
 * AMDiS - Adaptive multidimensional simulations
 *
 * Copyright (C) 2013 Dresden University of Technology. All Rights Reserved.
 * Web: https://fusionforge.zih.tu-dresden.de/projects/amdis
 *
 * Authors: 
 * 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.
 * 
 ******************************************************************************/
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
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
  
  void MacroReader::restoreMacroDofs(MacroInfo& macroInfo, int check)
  {
    FUNCNAME("MacroReader::restoreMacroDofs()");
    
    TEST_EXIT(macroInfo.isInitialized())
      ("Cannot set dof since macro info is not initialized.\n");
      
    Mesh* mesh = macroInfo.mesh;
    int nVertices = macroInfo.nVertices;
    int nElements = macroInfo.nElements;
    int **mel_vertex = macroInfo.mel_vertex;
    DegreeOfFreedom** dof = macroInfo.dof;
    
    TEST_EXIT(mesh->getNumberOfMacros() == nElements) 
      ("Mesh's current macro doesn't match the macro file.\n");
    
    mesh->setNumberOfElements(nElements);
    mesh->setNumberOfLeaves(nElements);
    mesh->setNumberOfVertices(nVertices);
    
    // Vertices dofs
    for (int i = 0; i < nVertices; i++)
      dof[i] = mesh->getDof(VERTEX);
    
    std::deque<MacroElement*>::iterator it = mesh->firstMacroElement();
    
    for (int i = 0; i < mesh->getNumberOfMacros(); i++) 
      for (int k = 0; k < mesh->getGeo(VERTEX); k++) 
	const_cast<Element*>((*(it + i))->getElement())->
	  setDof(k, dof[mel_vertex[i][k]]);
	  
    // Edge and face dofs
    if (mesh->getDim() > 1)
      boundaryDOFs(mesh);
    
    // Center dofs
    if (mesh->getNumberOfDofs(CENTER)) {
      it = mesh->firstMacroElement();
      
      for (int i = 0; i < mesh->getNumberOfMacros(); i++)
	const_cast<Element*>(it[i]->getElement())->
	  setDof(mesh->getNode(CENTER), mesh->getDof(CENTER));
    }
    
    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
	std::map<int, int> vertexMapEl1;
	std::map<int, int> vertexMapEl2;
129

130
131
	result = fscanf(file, "%d", &mode);
	TEST_EXIT(result == 1)("mode?\n");
132
      
133
134
	result = fscanf(file, "%d", &boundaryType);
	TEST_EXIT(result == 1)("boundaryType?\n");
135
      
136
137
138
	result = fscanf(file, "%d", &el1);
	TEST_EXIT(result == 1)("el1?\n");

139
	for (int j = 0; j < dim; j++) {
140
141
	  result = fscanf(file, "%d", &verticesEl1[j]);
	  TEST_EXIT(result == 1)("vertEl1[%d]\n", j);
142
	}
143
144
	result = fscanf(file, "%d", &el2);
	TEST_EXIT(result == 1)("el2?\n");
145
	for (int j = 0; j < dim; j++) {
146
147
	  result = fscanf(file, "%d", &verticesEl2[j]);
	  TEST_EXIT(result == 1)("vertEl2[%d]\n", j);
148
	}
149
	for (int j = 0; j < dim; j++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
150
	  if (mode == 0)
151
	    periodicMap.setEntry(melVertex[el1][verticesEl1[j]], 
Thomas Witkowski's avatar
Thomas Witkowski committed
152
				 melVertex[el2][verticesEl2[j]]);	  
153
154
155
156
157
158
	  vertexMapEl1[verticesEl1[j]] = verticesEl2[j];
	  vertexMapEl2[verticesEl2[j]] = verticesEl1[j];
	}

	// calculate sides of periodic vertices
	int sideEl1 = 0, sideEl2 = 0;
159
	if (dim == 1) {
160
161
162
	  sideEl1 = verticesEl1[0];
	  sideEl2 = verticesEl2[0];
	} else {
163
	  for (int j = 0; j < dim + 1; j++) {
164
165
166
	    sideEl1 += j;
	    sideEl2 += j;
	  }
167
	  for (int j = 0; j < dim; j++) {
168
169
170
171
	    sideEl1 -= verticesEl1[j];
	    sideEl2 -= verticesEl2[j];
	  }
	}
172

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

177
178
	Element *element1 = const_cast<Element*>((*(mel + el1))->getElement());
	Element *element2 = const_cast<Element*>((*(mel + el2))->getElement());
179
     
180
	// for all vertices of this side
181
	for (int j = 0; j < dim; j++) {
182
183
	  periodicCoordsEl1[element1->getPositionOfVertex(sideEl1, verticesEl1[j])] = 
	    coords[melVertex[el2][vertexMapEl1[verticesEl1[j]]]];
184

185
186
187
	  periodicCoordsEl2[element2->getPositionOfVertex(sideEl2, verticesEl2[j])] =
	    coords[melVertex[el1][vertexMapEl2[verticesEl2[j]]]];
	}
188
	     
189
190
191
192
	// decorate leaf data
	ElementData *ld1 = element1->getElementData();
	ElementData *ld2 = element2->getElementData();

193
194
195
196
197
198
199
200
	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());

201
202
203
204
	LeafDataPeriodic *ldp1 = 
	  dynamic_cast<LeafDataPeriodic*>(ld1->getElementData(PERIODIC));
	LeafDataPeriodic *ldp2 = 
	  dynamic_cast<LeafDataPeriodic*>(ld2->getElementData(PERIODIC));
205
206

	if (!ldp1) {
Thomas Witkowski's avatar
Thomas Witkowski committed
207
	  ldp1 = new LeafDataPeriodic(ld1);
208
209
210
211
	  element1->setElementData(ldp1);
	}

	if (!ldp2) {
Thomas Witkowski's avatar
Thomas Witkowski committed
212
	  ldp2 = new LeafDataPeriodic(ld2);
213
214
215
	  element2->setElementData(ldp2);
	}

216
217
	ldp1->addPeriodicInfo(mode, boundaryType, sideEl1, &periodicCoordsEl1);
	ldp2->addPeriodicInfo(mode, boundaryType, sideEl2, &periodicCoordsEl2);
218
219
220

	if (mode != 0) {
	  VertexVector *associated = mesh->periodicAssociations[boundaryType];
221

222
	  if (!associated) {
Thomas Witkowski's avatar
Thomas Witkowski committed
223
	    associated = new VertexVector(mesh->getVertexAdmin(), "vertex vector");
224
225
	    mesh->periodicAssociations[boundaryType] = associated;
	    VertexVector::Iterator it(associated, ALL_DOFS);
226
	    for (it.reset2(); !it.end(); ++it)
227
228
229
	      *it = it.getDOFIndex();
	  }

230
	  for (int j = 0; j < dim; j++) {
231
232
233
234
235
236
237
238
239
240
241
#ifdef DEBUG
	  {
		  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);
		  }
	  }
#endif
242
243
244
245
246
247
248
249
	    (*associated)[melVertex[el1][verticesEl1[j]]] =
	      melVertex[el2][vertexMapEl1[verticesEl1[j]]];
	    (*associated)[melVertex[el2][verticesEl2[j]]] =
	      melVertex[el1][vertexMapEl2[verticesEl2[j]]];
	  }
	}
      }    

Thomas Witkowski's avatar
Thomas Witkowski committed
250
251
      delete [] verticesEl1;
      delete [] verticesEl2;
252
253

      // change periodic vertex dofs
254
      for (int i = 0; i < mesh->getNumberOfVertices(); i++) {
255
	if (periodicMap.getEntry(i) != -1) {
256
	  mesh->freeDof(dof[i], VERTEX);
257
258
	  dof[i] = dof[periodicMap.getEntry(i)];

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

263
	  for (assoc = mesh->periodicAssociations.begin(); assoc != assocEnd; ++assoc) {
264
265
266
267
268

	    DegreeOfFreedom a = (*(assoc->second))[i];
	    if (a != i) {
	      (*(assoc->second))[i] = i;
	      (*(assoc->second))[a] = periodicMap.getEntry(i);
269
	    }
270
271
	  }

272
273
274
	}
      }

Thomas Witkowski's avatar
Thomas Witkowski committed
275
#if (DEBUG != 0)
276
277
      std::map<BoundaryType, VertexVector*>::iterator assoc;
      std::map<BoundaryType, VertexVector*>::iterator assocEnd =
278
	mesh->periodicAssociations.end();
Thomas Witkowski's avatar
Thomas Witkowski committed
279
280
      for (assoc = mesh->periodicAssociations.begin(); 
	   assoc != assocEnd; ++assoc)
Thomas Witkowski's avatar
Thomas Witkowski committed
281
	for (int i = 0; i < mesh->getNumberOfVertices(); i++)
282
283
	  if (i != (*(assoc->second))[i])
	    MSG("association %d: vertex %d -> vertex %d\n", 
Thomas Witkowski's avatar
Thomas Witkowski committed
284
		assoc->first, i, (*(assoc->second))[i]);
285

286
287
288
289
      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
290
#endif
291
292
293

    } // periodicFile

294

295
    // =========================================================
296
    // Set coordinate
297
298
299
    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]]);
300

301
	const_cast<Element*>((*(mel + i))->getElement())->
302
	  setDof(k, dof[melVertex[i][k]]);
303
      }
304
    }
305

306
    if (!macroInfo->neigh_set) {
307
      TEST_EXIT(periodicFile == "")
308
	("periodic boundary condition => element neighbours must be set\n");
309
      computeNeighbours(mesh);
310
    } else {
311
312
313
      /****************************************************************************/
      /* fill MEL oppVertex values when reading neighbour information form file  */
      /****************************************************************************/
314

315
316
317
318
319
      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) {
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
	    if (!macroInfo->neigh_inverse_set) {
	      int j = 0;
	      for (; j < mesh->getGeo(NEIGH); j++)
		if (neigh->getNeighbour(j) == *(mel + i))  
		  break;
	  
	      TEST_EXIT(j < mesh->getGeo(NEIGH))("el %d no neighbour of neighbour %d\n", 
						mel[i]->getIndex(), neigh->getIndex());
	      mel[i]->setOppVertex(k, j);
	    } else {
	      int j = 0;
	      for (; j < mesh->getGeo(NEIGH); j++)
		if (neigh->getNeighbourInv(j) == *(mel + i))  
		  break;
	  
	      TEST_EXIT(j < mesh->getGeo(NEIGH))("el %d no neighbour-inv of neighbour %d\n", 
						mel[i]->getIndex(), neigh->getIndex());
	      mel[i]->setOppVertex(k, j);
	    }
339
340
	  } else {
	    mel[i]->setOppVertex(k, -1);
341
	  }
342
	}
343
      }
344
    }
345

346
347
    if (!macroInfo->bound_set)
      macroInfo->dirichletBoundary();    
348
  
349
    if (mesh->getDim() > 1)
350
351
352
353
354
355
      boundaryDOFs(mesh);

    // initial boundary projections
    int numFaces = mesh->getGeo(FACE);
    int dim = mesh->getDim();
    mel = mesh->firstMacroElement();
356
    for (int i = 0; i < mesh->getNumberOfLeaves(); i++) {
357
      MacroElement *macroEl = *(mel + i);
358
      Projection *projector = macroEl->getProjection(0);
359
      if (projector && projector->getType() == VOLUME_PROJECTION) {
Thomas Witkowski's avatar
Thomas Witkowski committed
360
361
	for (int j = 0; j <= dim; j++)
	  projector->project(macroEl->getCoord(j));	
362
      } else {
363
	for (int j = 0; j < mesh->getGeo(EDGE); j++) {
364
	  projector = macroEl->getProjection(numFaces + j);
365
	  if (projector) {
366
367
368
369
370
371
372
373
374
375
376
	    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);

377
    if (mesh->getNumberOfDofs(CENTER)) {
378
      for (int i = 0; i < mesh->getNumberOfMacros(); i++)
379
	const_cast<Element*>(mel[i]->getElement())->
380
	  setDof(mesh->getNode(CENTER), mesh->getDof(CENTER));
381
382
    }

383
    // === Domain size ===
384

385
    WorldVector<double> x_min, x_max;
386

387
388
389
390
391
392
393
    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++) {
394
395
	x_min[j] = std::min(x_min[j], coords[i][j]);
	x_max[j] = std::max(x_max[j], coords[i][j]);
396
397
      }
    }
398

399
    for (int j = 0; j < Global::getGeo(WORLD); j++)
400
401
402
403
404
      mesh->setDiameter(j, x_max[j] - x_min[j]);

    if (check) {
      checkMesh(mesh);

Thomas Witkowski's avatar
Thomas Witkowski committed
405
406
      if (mesh->getDim() > 1)
	macroTest(mesh);
407
408
    }

409
    return macroInfo;
410
  }
411
  
412
413
  void MacroReader::computeNeighbours(Mesh *mesh)
  {
414
    FUNCNAME("MacroReader::computeNeighbours()");
415

416
    int dim = mesh->getDim();
417
    FixVec<DegreeOfFreedom*, DIMEN> dof(dim, NO_INIT);
418

419
420
421
    for (int i = 0; i < mesh->getNumberOfLeaves(); i++) {
      for (int k = 0; k < mesh->getGeo(NEIGH); k++) {
	mesh->getMacroElement(i)->setOppVertex(k, AMDIS_UNDEFINED);
422
	mesh->getMacroElement(i)->setNeighbour(k, NULL);
423
      }
424
    }
425

426
427
428
    for (int i = 0; i < mesh->getNumberOfLeaves(); i++) {
      for (int k = 0; k < mesh->getGeo(NEIGH); k++) {
	if (mesh->getMacroElement(i)->getBoundary(k) != INTERIOR) {
429
	  mesh->getMacroElement(i)->setNeighbour(k, NULL);
430
431
432
	  mesh->getMacroElement(i)->setOppVertex(k, -1);
	  continue;
	}
433

434
435
436
	if (mesh->getMacroElement(i)->getOppVertex(k) == AMDIS_UNDEFINED) {
	  if (dim == 1) {
	    dof[0] = const_cast<DegreeOfFreedom*>(mesh->getMacroElement(i)->
437
						  getElement()->getDof(k));
438
439
440
441
	  } else {
	    for (int l = 0; l < dim; l++)
	      dof[l] = const_cast<DegreeOfFreedom*>(mesh->getMacroElement(i)->
						    getElement()->
442
						    getDof((k + l + 1) % (dim + 1)));
443
	  }
444
445
446
447
448
449
450
451
452
453
454
455
456
	  
	  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;
	    }
	  }

457
458
459
460
461
462
463
	  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);
464
465
	      if (l < dim)
		std::cout << " / ";	      
466
467
468
	    }
	    std::cout << std::endl;
	    std::cout << "  dofs = ";
469
470
	    for (int l = 0; l < dim; l++)
	      std::cout << *(dof[l]) << " ";	    
471
472
473
474
	    std::cout << std::endl;

	    ERROR_EXIT("\n");
	  }    
475
	}
476
      }
477
    }
478
479
480
481
482
483
484
485
486
487
488
  }


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

  void MacroReader::boundaryDOFs(Mesh *mesh)
  {
489
490
491
492
493
    FUNCNAME("Mesh::boundaryDOFs()");

    int lnode = mesh->getNode(EDGE);
    int k, lne = mesh->getNumberOfLeaves();
    int max_n_neigh = 0, n_neigh, ov;
494
    std::deque<MacroElement*>::iterator mel;
495
    const MacroElement* neigh;
496
    DegreeOfFreedom *dof;
497
498
499
500
501
502

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

    int dim = mesh->getDim();

Thomas Witkowski's avatar
Thomas Witkowski committed
503
    switch (dim) {
504
    case 2:
505
      for (mel = mesh->firstMacroElement(); mel != mesh->endOfMacroElements(); mel++) {
506
507
508
509
510
511
	// check for periodic boundary
	Element *el = const_cast<Element*>((*mel)->getElement());
	ElementData *ed = el->getElementData(PERIODIC);

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

512
	if (ed) {
513
	  std::list<LeafDataPeriodic::PeriodicInfo> &periodicInfos = 
514
	    dynamic_cast<LeafDataPeriodic*>(ed)->getInfoList();
515
516
	  std::list<LeafDataPeriodic::PeriodicInfo>::iterator it;
	  std::list<LeafDataPeriodic::PeriodicInfo>::iterator end = periodicInfos.end();
Thomas Witkowski's avatar
Thomas Witkowski committed
517
518
	  for (it = periodicInfos.begin(); it != end; ++it)
	    if (it->type != 0)
519
520
521
	      periodic[it->elementSide] = true;
	}

522
	for (int i = 0; i < mesh->getGeo(NEIGH); i++) {
523
	  if (!(*mel)->getNeighbour(i) || 
524
525
526
527
	      ((*mel)->getNeighbour(i)->getIndex() < (*mel)->getIndex())) {

	    mesh->incrementNumberOfEdges(1);

528
	    if (mesh->getNumberOfDofs(EDGE)) {
529
530
	      dof = mesh->getDof(EDGE);
	      el->setDof(lnode + i, dof);
531
532
      
	      if ((*mel)->getNeighbour(i)) {
Thomas Witkowski's avatar
Thomas Witkowski committed
533
534
		Element *neigh = 
		  const_cast<Element*>((*mel)->getNeighbour(i)->getElement());
535

Thomas Witkowski's avatar
Thomas Witkowski committed
536
		if (periodic[i])
537
		  neigh->setDof(lnode + (*mel)->getOppVertex(i), mesh->getDof(EDGE));
Thomas Witkowski's avatar
Thomas Witkowski committed
538
		else
539
		  neigh->setDof(lnode + (*mel)->getOppVertex(i), dof);		
540
	      }
541
542
	    }
	  }  
543
544
545
546
547
548
	}
      }
      break;
    case 3:
      lnode = mesh->getNode(FACE);
      mel = mesh->firstMacroElement();
549
      for (int i = 0; i < lne; i++) {
550
551
552

	// check for periodic boundary
	Element *el = const_cast<Element*>((*(mel+i))->getElement());
553
	
554
555
556
557
	ElementData *ed = el->getElementData(PERIODIC);

	DimVec<bool> periodic(dim, DEFAULT_VALUE, false);
      
558
	if (ed) {
559
	  std::list<LeafDataPeriodic::PeriodicInfo> &periodicInfos = 
560
	    dynamic_cast<LeafDataPeriodic*>(ed)->getInfoList();
561
562
	  std::list<LeafDataPeriodic::PeriodicInfo>::iterator it;
	  std::list<LeafDataPeriodic::PeriodicInfo>::iterator end = periodicInfos.end();
563
564
	  for (it = periodicInfos.begin(); it != end; ++it)
	    if (it->type != 0)
565
566
567
	      periodic[it->elementSide] = true;
	}

568
569
	for (k = 0; k < mesh->getGeo(EDGE); k++) {
	  // === Check for not counted edges. ===
570
571
	  n_neigh = 1;

Thomas Witkowski's avatar
Thomas Witkowski committed
572
	  if (newEdge(mesh, (*(mel + i)), k, &n_neigh)) {
573
	    mesh->incrementNumberOfEdges(1);
574
	    max_n_neigh = std::max(max_n_neigh, n_neigh);
575
576
577
578
	  }
	}
      
	for (k = 0; k < mesh->getGeo(NEIGH); k++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
579
	  neigh = (*(mel + i))->getNeighbour(k);
580
	  // === Face is counted and dof is added by the element with bigger index. ===
Thomas Witkowski's avatar
Thomas Witkowski committed
581
582
	  if (neigh && (neigh->getIndex() > (*(mel + i))->getIndex()))  
	    continue;
583
584
585
	
	  mesh->incrementNumberOfFaces(1);
	
586
	  if (mesh->getNumberOfDofs(FACE)) {
587
	    TEST_EXIT(!(*(mel + i))->getElement()->getDof(lnode + k))
588
	      ("dof %d on element %d already set\n", 
Thomas Witkowski's avatar
Thomas Witkowski committed
589
	       lnode + k, (*(mel + i))->getIndex());
590
	  
591
	    const_cast<Element*>((*(mel + i))->getElement())->setDof(lnode + k,
592
								     mesh->getDof(FACE));
593
594

	    if (neigh) {
Thomas Witkowski's avatar
Thomas Witkowski committed
595
	      ov = (*(mel + i))->getOppVertex(k);
596
	      TEST_EXIT(!neigh->getElement()->getDof(lnode + ov))
597
		("dof %d on neighbour %d already set\n", 
Thomas Witkowski's avatar
Thomas Witkowski committed
598
		 lnode + ov, neigh->getIndex());
599
600
	    
	      Element *neighEl = 
Thomas Witkowski's avatar
Thomas Witkowski committed
601
		const_cast<Element*>((*(mel + i))->getNeighbour(k)->getElement());
602

Thomas Witkowski's avatar
Thomas Witkowski committed
603
	      if (periodic[k])
604
		neighEl->setDof(lnode+ov, mesh->getDof(FACE));
Thomas Witkowski's avatar
Thomas Witkowski committed
605
	      else
606
		neighEl->setDof(lnode+ov, const_cast<int*>((*(mel + i))->getElement()->
607
							   getDof(lnode + k)));	      
608
609
610
611
612
	    }
	  }
	}
      }
      break;
Thomas Witkowski's avatar
Thomas Witkowski committed
613
614
    default: 
      ERROR_EXIT("invalid dim\n");
615
616
    }
    
617
    if (3 == dim)
Thomas Witkowski's avatar
Thomas Witkowski committed
618
      mesh->setMaxEdgeNeigh(std::max(8, 2 * max_n_neigh));
619
    else
Thomas Witkowski's avatar
Thomas Witkowski committed
620
      mesh->setMaxEdgeNeigh(dim - 1);    
621
622
623
624
625
626
627
628
  }

  /* 
     testet mesh auf auftretende Zyklen
  
     wenn Zyklus auftritt:
     ordnet Eintraege in MacroElement-Struktur um, so dass kein Zyklus auftritt
     erzeugt neue Macro-Datei nameneu mit umgeordnetem Netz 
629
     (wenn nameneu=NULL wird keine MAcro-Datei erzeugt)
630
631
  */      

Thomas Witkowski's avatar
Thomas Witkowski committed
632
  void MacroReader::macroTest(Mesh *mesh)
633
  {
634
    FUNCNAME("MacroReader::macroTest()");
635
   
636
637
638
    int i = macrotest(mesh);

    if (i >= 0) {
Thomas Witkowski's avatar
Thomas Witkowski committed
639
640
      WARNING("There is a cycle beginning in macro element %d\n", i);
      WARNING("Entries in MacroElement structures get reordered\n");
641
      umb(NULL, mesh, umbVkantMacro);
642
    }
643
644
645
646
647
648
649
650
651
652
653
654
655
656
  }
  
  /****************************************************************************/
  /*  macro_test():                              Author: Thomas Kastl (1998)  */
  /****************************************************************************/
  /*
    testet mesh auf auftretende Zyklen
  
    wenn mesh zyklenfrei -> -1
    sonst ->  globaler Index des Macroelementes bei dem ein Zyklus beginnt 
  */

  int MacroReader::macrotest(Mesh *mesh)
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
657
658
    std::deque<MacroElement*>::const_iterator macro, mac;
    int flg = 0;
659
    int dim = mesh->getDim();
Thomas Witkowski's avatar
Thomas Witkowski committed
660
661
    int *test = new int[mesh->getNumberOfMacros()];
    int *zykl = new int[mesh->getNumberOfMacros()];
662
 
Thomas Witkowski's avatar
Thomas Witkowski committed
663
    for (int i = 0; i < mesh->getNumberOfMacros(); i++)
664
      test[i] = 0;
665

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

669
670
671
672
    while (macrolfd != mesh->endOfMacroElements()) {
      if (test[(*macrolfd)->getIndex()] == 1) {
	macrolfd++;
      } else {
673
674
	for (int i = 0; i < mesh->getNumberOfMacros(); i++)
	  zykl[i] = 0;	
675
    
676
677
678
679
680
681
682
683
	macro = macrolfd;
	flg = 2;
	do {
	  if (zykl[(*macro)->getIndex()] == 1) {
	    flg = 0;
	    zykstart = (*macro)->getIndex();
	  } else {
	    zykl[(*macro)->getIndex()] = 1;
684
      
685
686
	    if (test[(*macro)->getIndex()] == 1) {
	      flg = 1;
687
	    } else if ((*macro)->getNeighbour(dim) == NULL) {
688
689
	      flg = 1;
	      test[(*macro)->getIndex()] = 1;
Thomas Witkowski's avatar
Thomas Witkowski committed
690
	    } else if ((*macro) == (*macro)->getNeighbour(dim)->getNeighbour(dim)) {
691
692
693
694
	      flg = 1;
	      test[(*macro)->getIndex()] = 1;
	      test[(*macro)->getNeighbour(dim)->getIndex()] = 1;
	    } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
695
696
	      for (mac = mesh->firstMacroElement(); 
		   (*mac) != (*macro)->getNeighbour(dim); mac++);
697
698
699
700
	      macro = mac;
	    } 
	  }	  
	} while(flg == 2);
701
 
702
	if (flg == 1)
703
	  macrolfd++;
704
705
	else 
	  macrolfd=mesh->endOfMacroElements();	
706
      }
707
    }
708
  
Thomas Witkowski's avatar
Thomas Witkowski committed
709
710
    delete [] zykl;
    delete [] test;
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
 
    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
  //           (wird nur benoetigt, wenn umbvk=umb_vkant_macrodat) 
  
  //   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
  //                  neuen Macro-Datei angewendet werden 
  //           = 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
735
			void (*umbvk)(Mesh*, MacroElement*, int, int*))
736
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
737
    int *test = new int[mesh->getNumberOfMacros()];
738
  
Thomas Witkowski's avatar
Thomas Witkowski committed
739
740
    for (int i = 0; i < static_cast<int>(mesh->getNumberOfMacros()); i++)
      test[i] = 0;
741

742
    recumb(mesh, (*mesh->firstMacroElement()), NULL, test, 0, 0, ele, umbvk);
743

Thomas Witkowski's avatar
Thomas Witkowski committed
744
    delete [] test;
745
746
747
748
749
  }

  bool MacroReader::newEdge(Mesh *mesh, MacroElement *mel,
			    int mel_edge_no, int *n_neigh)
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
750
751
    FUNCNAME("MacroElement::newEdge()"); 
    MacroElement *nei;
752
    const DegreeOfFreedom *dof[2];
753
    DegreeOfFreedom *edge_dof = NULL;
Thomas Witkowski's avatar
Thomas Witkowski committed
754
755
    int j, k, opp_v, node = 0;
    BoundaryType lbound = INTERIOR;
756
    Projection *lproject = NULL;
Thomas Witkowski's avatar
Thomas Witkowski committed
757
    const int max_no_list_el = 100;
758
759
760
761
    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
762
    static int next_el[6][2] = {{3,2},{1,3},{1,2},{0,3},{0,2},{0,1}};
763
764
    int vertices = mesh->getGeo(VERTEX);

Thomas Witkowski's avatar
Thomas Witkowski committed
765
    int mel_index = mel->getIndex();
766
767
768
769

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

770
    if (mesh->getNumberOfDofs(EDGE)) {
771
      node = mesh->getNode(EDGE);
772
      if (el->getDof(node+edge_no)) {
773
774
775
776
777
778
	/****************************************************************************/
	/*  edge was counted by another macro element and dof was added on the      */
	/*  complete patch                                                          */
	/****************************************************************************/
	return false;
      } else {
779
780
	edge_dof = mesh->getDof(EDGE);
	el->setDof(node+edge_no, edge_dof);
781
782
783
      }
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
784
    for (j = 0; j < 2; j++)
785
      dof[j] = el->getDof(el->getVertexOfEdge(edge_no, j));
Thomas Witkowski's avatar
Thomas Witkowski committed
786
    
787
788
789
790
791
792
793
794
795
796
797
798
    /****************************************************************************/
    /*  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]);


799
    if (mel->getBoundary(next_el[edge_no][0])) {
800
801
802
803
804
805
      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++)
806
	if (nei->getElement()->getDof(j) == dof[0])  break;
807
      for (k = 0; k < vertices; k++)
808
	if (nei->getElement()->getDof(k) == dof[1])  break;
809
810

      // check for periodic boundary
811
      if (j == 4 || k == 4) {
812
	nei = NULL;
813
814
815
	break;
      }

816
      if (mesh->getNumberOfDofs(EDGE)) {
817
818
	TEST_EXIT(nei->index > mel_index)
	  ("neighbour index %d < element index %d\n", nei->getIndex(), mel_index);
819
      }
820

821
      if (!mesh->getNumberOfDofs(EDGE) && nei->getIndex() < mel_index) 
822
	return false;
823
    
824
      edge_no = Tetrahedron::edgeOfDofs[j][k];
825
826
827
828
829
830
831
832
833
834

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

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

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

835
      if (mesh->getNumberOfDofs(EDGE))
836
	nei->element->setDof(node+edge_no,edge_dof);      
837

838
      if (next_el[edge_no][0] != opp_v) {
839
	if (nei->getBoundary(next_el[edge_no][0])) {
840
841
842
843
844
845
	  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()))
846
847
848
	      lproject = neiProject;
	  }
	}
849
850
851
852
853
854
855
856
857
858
	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);
	  Projection *neiProject = nei->getProjection(next_el[edge_no][1]);	
	  if (!lproject) {
	    lproject = neiProject;
	  } else {
	    if (neiProject && (lproject->getID() < neiProject->getID()))
859
860
861
	      lproject = neiProject;
	  }
	}
862
863
864
	opp_v = nei->getOppVertex(next_el[edge_no][1]);
	nei = nei->getNeighbour(next_el[edge_no][1]);
      }
865
866
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
    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;
      
      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])) {
	lbound = newBound(mel->getBoundary(next_el[edge_no][1]), lbound); 
	Projection *neiProject =  mel->getProjection(next_el[edge_no][1]);
	if (!lproject) {
	  lproject = neiProject;
	} else {
	  if (neiProject && (lproject->getID() < neiProject->getID()))
	    lproject = neiProject;	  
	}
      }
      
      while (nei) {
	for (j = 0; j < vertices; j++)
889
	  if (nei->getElement()->getDof(j) == dof[0])  break;
Thomas Witkowski's avatar
Thomas Witkowski committed
890
	for (k = 0; k < vertices; k++)
891
	  if (nei->getElement()->getDof(k) == dof[1])  break;
892
	
Thomas Witkowski's avatar
Thomas Witkowski committed
893
894
895
896
	// check for periodic boundary
	if (j == 4 || k == 4)
	  return false;
	
897
	if (mesh->getNumberOfDofs(EDGE)) {
Thomas Witkowski's avatar
Thomas Witkowski committed
898
899
900
	  TEST_EXIT(nei->getIndex() > mel_index)
	    ("neighbour index %d < element index %d\n", nei->getIndex(),
	     mel_index);
901
	}
Thomas Witkowski's avatar
Thomas Witkowski committed
902
903
904
905
	
	if (nei->getIndex() < mel_index)  
	  return false;
	
906
	edge_no = Tetrahedron::edgeOfDofs[j][k];
Thomas Witkowski's avatar
Thomas Witkowski committed
907
908
909
910
911
	
	TEST_EXIT(*n_neigh < max_no_list_el)("too many neigbours for local list\n");
	
	list_bound[(*n_neigh)] = &(nei->boundary[mesh->getGeo(FACE) + edge_no]);
	list_project[(*n_neigh)++] = &(nei->projection[mesh->getGeo(FACE) + edge_no]);
912

913
	if (mesh->getNumberOfDofs(EDGE)) {
914
	  TEST_EXIT(!nei->getElement()->getDof(node+edge_no))
Thomas Witkowski's avatar
Thomas Witkowski committed
915
916
917
	    ("dof %d on element %d is already set, but not on element %d\n",
	     node + edge_no, nei->getIndex(), mel_index);
	  
918
	  nei->element->setDof(node+edge_no, edge_dof);
Thomas Witkowski's avatar
Thomas Witkowski committed
919
	}
920

Thomas Witkowski's avatar
Thomas Witkowski committed
921
922
923
924
925
926
927
928
929
	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()))
		lproject = neiProject;	      
930
931
	    }
	  }
Thomas Witkowski's avatar
Thomas Witkowski committed
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
	  
	  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); 
	    Projection *neiProject = nei->getProjection(next_el[edge_no][1]);
	    if (!lproject) {
	      lproject = neiProject;
	    } else {
	      if (neiProject && (lproject->getID() < neiProject->getID()))
		lproject = neiProject;	      
	    }
	  }
	  
	  opp_v = nei->getOppVertex(next_el[edge_no][1]);
	  nei = nei->getNeighbour(next_el[edge_no][1]);
	}
950
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
951
952
    }
    
953
954
955
956
957
958
959
960
961
    for (j = 0; j < *n_neigh; j++) {
      *(list_bound[j]) = lbound;
      *(list_project[j]) = lproject;
    }
  
    return true;
  }

  void MacroReader::fillMelBoundary(Mesh *mesh, MacroElement *mel, 
962
				    FixVec<BoundaryType,NEIGH> const& ind)
963
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
964
965
    for (int i = 0; i < mesh->getGeo(NEIGH); i++)
      mel->boundary[i] = ind[i];        
966
967
968
969
  }


  void MacroReader::fillMelNeigh(MacroElement *mel,
970
				 std::deque<MacroElement*>& macro_elements, 
971
				 FixVec<int,NEIGH> const& ind)
972
973
974
  {
    int dim = mel->element->getMesh()->getDim();