MacroReader.cc 61.9 KB
Newer Older
1
2
3
4
5
6
#include "MacroReader.h"
#include "MacroWriter.h"
#include "MacroElement.h"
#include "Boundary.h"
#include "FiniteElemSpace.h"
#include "Mesh.h"
7
#include <string.h>
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
#include "FixVec.h"
#include "FixVecConvert.h"
#include "PeriodicMap.h"
#include "ElInfo.h"
#include "Parameters.h"
#include "DOFIterator.h"
#include "SurfaceRegion_ED.h"
#include "ElementRegion_ED.h"
#include "LeafData.h"
#include "VertexVector.h"
#include <map>
#include <iostream>
#include <fstream>

namespace AMDiS {

  MacroInfo* MacroReader::readMacro(const char *filename, 
				    Mesh* mesh,
				    const char *periodicFile,
				    int check)
  {
    FUNCNAME("Mesh::readMacro()");

    TEST_EXIT(filename)("no file specified; filename NULL pointer\n");
32
   
Thomas Witkowski's avatar
Thomas Witkowski committed
33
    MacroInfo *macroInfo = new MacroInfo();
34
35
    macroInfo->readAMDiSMacro(filename, mesh);

36
    std::deque<MacroElement*>::iterator mel = macroInfo->mel.begin();
37
38
39
    int **melVertex = macroInfo->mel_vertex;
    WorldVector<double> *coords = macroInfo->coords;
    DegreeOfFreedom **dof = macroInfo->dof;
40
41

    // === read periodic data =================================
42
    if (periodicFile && (strcmp(periodicFile, "") != 0)) {
43
      WARNING("Periodic boundaries may lead to errors in small meshes if element neighbours are not set!\n");
44
45
46
47
48
49
50
51
    
      FILE *file = fopen(periodicFile, "r");
      TEST_EXIT(file)("can't open file %s\n", periodicFile);

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

      int el1, el2;
Thomas Witkowski's avatar
Thomas Witkowski committed
52
53
      int *verticesEl1 = new int[dim];
      int *verticesEl2 = new int[dim];
54
      int mode = -1; // 0: drop dofs, 1: associate dofs
55
      int result;
56
      BoundaryType boundaryType;
57
      PeriodicMap periodicMap;
58
59
60

      fscanf(file, "%*s %d", &n);
      fscanf(file, "%*s %*s %*s %*s %*s %*s %*s %*s %*s %*s %*s");
61
   
62
      for (int i = 0; i < n; i++) {
63
64
	std::map<int, int> vertexMapEl1;
	std::map<int, int> vertexMapEl2;
65

66
67
	result = fscanf(file, "%d", &mode);
	TEST_EXIT(result == 1)("mode?\n");
68
      
69
70
	result = fscanf(file, "%d", &boundaryType);
	TEST_EXIT(result == 1)("boundaryType?\n");
71
      
72
73
74
	result = fscanf(file, "%d", &el1);
	TEST_EXIT(result == 1)("el1?\n");

75
	for (int j = 0; j < dim; j++) {
76
77
	  result = fscanf(file, "%d", &verticesEl1[j]);
	  TEST_EXIT(result == 1)("vertEl1[%d]\n", j);
78
	}
79
80
	result = fscanf(file, "%d", &el2);
	TEST_EXIT(result == 1)("el2?\n");
81
	for (int j = 0; j < dim; j++) {
82
83
	  result = fscanf(file, "%d", &verticesEl2[j]);
	  TEST_EXIT(result == 1)("vertEl2[%d]\n", j);
84
	}
85
	for (int j = 0; j < dim; j++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
86
	  if (mode == 0)
87
	    periodicMap.setEntry(melVertex[el1][verticesEl1[j]], 
Thomas Witkowski's avatar
Thomas Witkowski committed
88
				 melVertex[el2][verticesEl2[j]]);	  
89
90
91
92
93
94
	  vertexMapEl1[verticesEl1[j]] = verticesEl2[j];
	  vertexMapEl2[verticesEl2[j]] = verticesEl1[j];
	}

	// calculate sides of periodic vertices
	int sideEl1 = 0, sideEl2 = 0;
95
	if (dim == 1) {
96
97
98
	  sideEl1 = verticesEl1[0];
	  sideEl2 = verticesEl2[0];
	} else {
99
	  for (int j = 0; j < dim + 1; j++) {
100
101
102
	    sideEl1 += j;
	    sideEl2 += j;
	  }
103
	  for (int j = 0; j < dim; j++) {
104
105
106
107
	    sideEl1 -= verticesEl1[j];
	    sideEl2 -= verticesEl2[j];
	  }
	}
108

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

113
114
	Element *element1 = const_cast<Element*>((*(mel + el1))->getElement());
	Element *element2 = const_cast<Element*>((*(mel + el2))->getElement());
115
     
116
	// for all vertices of this side
117
	for (int j = 0; j < dim; j++) {
118
119
	  periodicCoordsEl1[element1->getPositionOfVertex(sideEl1, verticesEl1[j])] = 
	    coords[melVertex[el2][vertexMapEl1[verticesEl1[j]]]];
120

121
122
123
	  periodicCoordsEl2[element2->getPositionOfVertex(sideEl2, verticesEl2[j])] =
	    coords[melVertex[el1][vertexMapEl2[verticesEl2[j]]]];
	}
124
	     
125
126
127
128
129
130
131
132
	// decorate leaf data
	ElementData *ld1 = element1->getElementData();
	ElementData *ld2 = element2->getElementData();

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

	if (!ldp1) {
Thomas Witkowski's avatar
Thomas Witkowski committed
133
	  ldp1 = new LeafDataPeriodic(ld1);
134
135
136
137
	  element1->setElementData(ldp1);
	}

	if (!ldp2) {
Thomas Witkowski's avatar
Thomas Witkowski committed
138
	  ldp2 = new LeafDataPeriodic(ld2);
139
140
141
	  element2->setElementData(ldp2);
	}

142
143
	ldp1->addPeriodicInfo(mode, boundaryType, sideEl1, &periodicCoordsEl1);
	ldp2->addPeriodicInfo(mode, boundaryType, sideEl2, &periodicCoordsEl2);
144
145
146

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

148
	  if (!associated) {
Thomas Witkowski's avatar
Thomas Witkowski committed
149
	    associated = new VertexVector(mesh->getVertexAdmin(), "vertex vector");
150
151
	    mesh->periodicAssociations[boundaryType] = associated;
	    VertexVector::Iterator it(associated, ALL_DOFS);
152
	    for (it.reset2(); !it.end(); ++it)
153
154
155
	      *it = it.getDOFIndex();
	  }

156
	  for (int j = 0; j < dim; j++) {
157
158
159
160
161
162
163
164
	    (*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
165
166
      delete [] verticesEl1;
      delete [] verticesEl2;
167
168

      // change periodic vertex dofs
169
      for (int i = 0; i < mesh->getNumberOfVertices(); i++) {
170
171
172
173
	if (periodicMap.getEntry(i) != -1) {
	  mesh->freeDOF(dof[i], VERTEX);
	  dof[i] = dof[periodicMap.getEntry(i)];

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

178
	  for (assoc = mesh->periodicAssociations.begin(); assoc != assocEnd; ++assoc) {
179
180
181
182
183

	    DegreeOfFreedom a = (*(assoc->second))[i];
	    if (a != i) {
	      (*(assoc->second))[i] = i;
	      (*(assoc->second))[a] = periodicMap.getEntry(i);
184
	    }
185
186
	  }

187
188
189
	}
      }

190
191
      std::map<BoundaryType, VertexVector*>::iterator assoc;
      std::map<BoundaryType, VertexVector*>::iterator assocEnd =
192
	mesh->periodicAssociations.end();
193
      for (assoc = mesh->periodicAssociations.begin(); assoc != assocEnd; ++assoc) {
Thomas Witkowski's avatar
Thomas Witkowski committed
194
	for (int i = 0; i < mesh->getNumberOfVertices(); i++)
195
196
	  if (i != (*(assoc->second))[i])
	    MSG("association %d: vertex %d -> vertex %d\n", 
Thomas Witkowski's avatar
Thomas Witkowski committed
197
		assoc->first, i, (*(assoc->second))[i]);	
198
199
      }

200
201
202
203
      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));
204
    }
205

206
207
    // =========================================================

208
209
210
    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]]);
211

212
	const_cast<Element*>((*(mel + i))->getElement())->
213
	  setDOF(k, dof[melVertex[i][k]]);
214
      }
215
    }
216

217
218
219
    if (!macroInfo->neigh_set) {
      TEST_EXIT(!periodicFile)
	("periodic boundary condition => element neighbours must be set\n");
220
      computeNeighbours(mesh);
221
    } else {
222
223
224
      /****************************************************************************/
      /* fill MEL oppVertex values when reading neighbour information form file  */
      /****************************************************************************/
225

226
227
228
229
230
      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) {
231
232
	    int j = 0;
	    for (; j < mesh->getGeo(NEIGH); j++)
233
234
	      if (neigh->getNeighbour(j) == *(mel + i))  
		break;
235
	
236
237
238
239
240
	    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 {
	    mel[i]->setOppVertex(k, -1);
241
	  }
242
	}
243
      }
244
    }
245

246
    if (!macroInfo->bound_set) {
247
      macroInfo->dirichletBoundary();
248
    }
249
  
250
    if (mesh->getDim() > 1)
251
252
253
254
255
256
      boundaryDOFs(mesh);

    // initial boundary projections
    int numFaces = mesh->getGeo(FACE);
    int dim = mesh->getDim();
    mel = mesh->firstMacroElement();
257
    for (int i = 0; i < mesh->getNumberOfLeaves(); i++) {
258
259
      MacroElement *macroEl = *(mel+i);
      Projection *projector = macroEl->getProjection(0);
260
      if (projector && projector->getType() == VOLUME_PROJECTION) {
Thomas Witkowski's avatar
Thomas Witkowski committed
261
262
	for (int j = 0; j <= dim; j++)
	  projector->project(macroEl->getCoord(j));	
263
      } else {
264
	for (int j = 0; j < mesh->getGeo(EDGE); j++) {
265
	  projector = macroEl->getProjection(numFaces + j);
266
	  if (projector) {
267
268
269
270
271
272
273
274
275
276
277
278
	    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);

    if (mesh->getNumberOfDOFs(CENTER)) {
279
      for (int i = 0; i < mesh->getNumberOfMacros(); i++)
280
	const_cast<Element*>(mel[i]->getElement())->
281
	  setDOF(mesh->getNode(CENTER), mesh->getDOF(CENTER));
282
283
284
285
286
287
    }

    /****************************************************************************/
    /* domain size                                                              */
    /****************************************************************************/

288
    WorldVector<double> x_min, x_max;
289

290
291
292
293
294
295
296
    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++) {
297
298
	x_min[j] = std::min(x_min[j], coords[i][j]);
	x_max[j] = std::max(x_max[j], coords[i][j]);
299
300
      }
    }
301

302
    for (int j = 0; j < Global::getGeo(WORLD); j++)
303
304
305
306
307
308
      mesh->setDiameter(j, x_max[j] - x_min[j]);

    if (check) {
      checkMesh(mesh);

      if (mesh->getDim() > 1) {
309
310
311
312
313
	char filenew[128];
	strncpy(filenew, filename, 128); 
	filenew[127] = 0;
	strncat(filenew, ".new", 128);   
	filenew[127] = 0;
314
315
316
317
	macroTest(mesh, filenew);
      }
    }

318
    return macroInfo;
319
320
321
322
323
324
325
326
  }

  /****************************************************************************/
  /*  fill macro info structure and some pointers in mesh ...                 */
  /****************************************************************************/

  void MacroInfo::fill(Mesh *pmesh, int ne, int nv)
  {
327
328
    FUNCNAME("MacroInfo::fill()");

329
330
    TEST_EXIT(pmesh)("no mesh\n");

331
    int dim = pmesh->getDim(); 
332
    mesh = pmesh;
333
334
335
336
337

    mesh->setNumberOfElements(ne);
    mesh->setNumberOfLeaves(ne);
    mesh->setNumberOfVertices(nv);

338
    for (int i = 0; i < ne; i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
339
      MacroElement *newMacro = new MacroElement(mesh->getDim());
340
341
342
343
      mel.push_back(newMacro);
      mesh->addMacroElement(mel[i]);
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
344
345
346
    dof = new DegreeOfFreedom*[nv];
    coords = new WorldVector<double>[nv];
    mel_vertex = new int*[ne];
347

Thomas Witkowski's avatar
Thomas Witkowski committed
348
349
    for (int i = 0; i < ne; i++)
      mel_vertex[i] = new int[mesh->getGeo(VERTEX)];
350

Thomas Witkowski's avatar
Thomas Witkowski committed
351
    for (int i = 0; i < nv; i++)
352
353
      dof[i] = mesh->getDOF(VERTEX);

354
    for (int i = 0; i < ne; i++) {
355
      mel[i]->element = mesh->createNewElement();
Thomas Witkowski's avatar
Thomas Witkowski committed
356
      mel[i]->index = i;
357

Thomas Witkowski's avatar
Thomas Witkowski committed
358
      if (dim == 3)
Thomas Witkowski's avatar
Thomas Witkowski committed
359
	mel[i]->elType = 0;
360
361
362
363
364
    }
    neigh_set = false;
    bound_set = false;
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
365
  void MacroInfo::clear()
366
367
  {
    for (int i = 0; i < mesh->getNumberOfMacros(); i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
368
      delete [] mel_vertex[i];
369

Thomas Witkowski's avatar
Thomas Witkowski committed
370
371
    delete [] mel_vertex;
    delete [] coords;
372
    coords = NULL;  
Thomas Witkowski's avatar
Thomas Witkowski committed
373
    delete [] dof;
374
    dof = NULL;
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395

    mesh = NULL;
    neigh_set = false;
  }

  /****************************************************************************/
  /****************************************************************************/
  /*  tool for reading macro triangulations in ALBERT-format                  */
  /****************************************************************************/
  /****************************************************************************/

  /****************************************************************************/
  /*  read_indices()  reads dim+1 indices from  file  into  id[0-dim],        */
  /*    returns true if dim+1 inputs arguments could be read successfully by  */
  /*    fscanf(), else false                                                  */
  /****************************************************************************/

  int  MacroInfo::read_indices(FILE *file, DimVec<int> &id)
  {
    int dim = mesh->getDim();

Thomas Witkowski's avatar
Thomas Witkowski committed
396
    for (int i = 0; i <= dim; i++)
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
      if (fscanf(file, "%d", &id[i]) != 1)
	return(false);

    return(true);
  }

#define N_KEYS      14
#define N_MIN_KEYS  7
  static const char *keys[N_KEYS] = {
    "DIM",                   //  0 
    "DIM_OF_WORLD",          //  1
    "number of vertices",    //  2
    "number of elements",    //  3
    "vertex coordinates",    //  4
    "element vertices",      //  5
    "element boundaries",    //  6
    "element neighbours",    //  7
    "element type",          //  8
    "projections",           //  9
    "element region",        // 10
    "surface region",        // 11
    "mesh name",             // 12
    "time"                   // 13
  };

  static int get_key_no(const char *key)
  {
424
425
    for (int i = 0; i < N_KEYS; i++)
      if (!strcmp(keys[i], key))  
426
	return i;
427

428
    return -1;
429
430
431
432
433
434
  }

#include <ctype.h>

  static const char *read_key(const char *line)
  {
435
436
    static char key[100];
    char *k = key;
437
438

    while (isspace(*line)) 
439
      line++;
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
    while ((*k++ = *line++) != ':');
    *--k = '\0';
  
    return(const_cast<const char *>( key));
  }

  /****************************************************************************/
  /*  read_albert_macro():                                                    */
  /*    read macro triangulation from ascii file in ALBERT format             */
  /*    fills macro_info structure                                            */
  /*    called by read_macro(), fills missing information                     */
  /****************************************************************************/


  void MacroInfo::readAMDiSMacro(const char *filename, Mesh* mesh)
  {
456
457
458
459
460
    FUNCNAME("MacroInfo::readAMDiSMacro()");

    int dim;
    int dow, nv, ne, j, k;
    double dbl;
461
462
    char line[256];
    int line_no, n_keys, sort_key[N_KEYS], nv_key, ne_key;
463
    int key_def[N_KEYS] = {0,0,0,0,0,0,0,0,0,0,0,0};
464
465
466
467
468
469
470
    const char *key;
    DimVec<int> *ind = NULL;

    TEST_EXIT(filename)("no file specified; filename NULL pointer\n");
    TEST_EXIT(strlen(filename) < static_cast<unsigned int>(127))
      ("can only handle filenames up to 127 characters\n");

Thomas Witkowski's avatar
Thomas Witkowski committed
471
    FILE *file = fopen(filename, "r");
472
    TEST_EXIT(file)("cannot open file %s\n", filename);
473
474
475
476
477
478
479

    /****************************************************************************/
    /*  looking for all keys in the macro file ...                              */
    /****************************************************************************/

    line_no = n_keys = 0;
    while (fgets(line, 255, file)) {
480
      line_no++;
Thomas Witkowski's avatar
Thomas Witkowski committed
481
482
      if (!strchr(line, ':'))  
	continue;
483
484
485
486
487
488
489
490
491
      key = read_key(line);
      int i_key = get_key_no(key);
      TEST_EXIT(i_key >= 0)
	("macro file %s must not contain key %s on line %d\n",
	 filename, key, line_no);
      TEST_EXIT(!key_def[i_key])("key %s defined second time on line %d in file %s\n");

      sort_key[n_keys++] = i_key;
      key_def[i_key] = true;
492
493
    }
    fclose(file);
494

495

496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
    /*******************************************************************************/
    /*  Test, if there is data for every key and if all is defined in right order. */
    /*******************************************************************************/

    for (int i_key = 0; i_key < N_MIN_KEYS; i_key++) {
      for (j = 0; j < n_keys; j++)
	if (sort_key[j] == i_key)  break;
      TEST_EXIT(j < n_keys)("You do not have specified data for %s in %s\n",
			    keys[i_key], filename);

      for (j = 0; j < n_keys; j++)
	if (sort_key[j] == 2)  break;
      nv_key = j;
      for (j = 0; j < n_keys; j++)
	if (sort_key[j] == 3)  break;
      ne_key = j;
512
    
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
      switch (i_key) {
      case 0:
      case 1:
	TEST_EXIT(sort_key[i_key] < 2)
	  ("You have to specify DIM or mesh->getGeo(WORLD) before all other data\n");
	break;
      case 4: 
	TEST_EXIT(nv_key < i_key)
	  ("Before reading data for %s, you have to specify the %s in file\n",
	   keys[4], keys[2], filename);
	break;
      case 5: 
	TEST_EXIT(nv_key < i_key  &&  ne_key < i_key)
	  ("Before reading data for %s, you have to specify the %s and %s in file %s\n",
	   keys[5], keys[3], keys[2], filename);
      case 6:
      case 7:
      case 8:
	TEST_EXIT(ne_key < i_key)
	  ("Before reading data for %s, you have to specify the %s in file %s\n",
	   keys[i_key], keys[3], filename);
      }
535
536
    }

537
    for (int i_key = 0; i_key < N_KEYS; i_key++)
538
539
540
541
542
543
      key_def[i_key] = false;

    /****************************************************************************/
    /*  and now, reading data ...                                               */
    /****************************************************************************/
	
544
545
    file = fopen(filename, "r");
    TEST_EXIT(file)("cannot open file %s\n", filename);
546
547

    int result;
548

549
    for (int i_key = 0; i_key < n_keys; i_key++) {
550

551
552
553
554
555
556
      switch (sort_key[i_key]) {
	
      case 0:
	// line "DIM"
	result = fscanf(file, "%*s %d", &dim);
	TEST_EXIT(result == 1)("cannot read DIM correctly in file %s\n", filename);
557

558
	ind = new DimVec<int>(dim, NO_INIT);
559

560
561
	key_def[0] = true;
	break;
562

563
564
565
566
567
568
569
570
      case 1:
	// line "DIM_OF_WORLD"
	result = fscanf(file, "%*s %d", &dow);
	TEST_EXIT(result == 1)
	  ("cannot read Global::getGeo(WORLD) correctly in file %s\n", filename);
	TEST_EXIT(dow == Global::getGeo(WORLD))
	  ("dimension of world = %d != Global::getGeo(WORLD) = %d\n", 
	   dow, Global::getGeo(WORLD));
571

572
573
574
575
576
577
578
579
580
581
582
583
	key_def[1] = true;
	break;

      case 2:
	// line "number of vertices"
	result = fscanf(file, "%*s %*s %*s %d", &nv);
	TEST_EXIT(result == 1)
	  ("cannot read number of vertices correctly in file %s\n", filename);
	TEST_EXIT(nv > 0)
	  ("number of vertices = %d must be bigger than 0\n", nv);

	key_def[2] = true;
584
585
	if (key_def[3])
	  fill(mesh, ne, nv);
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
	break;

      case 3:
	// line "number of elements"
	result = fscanf(file, "%*s %*s %*s %d", &ne);
	TEST_EXIT(result == 1)
	  ("cannot read number of elements correctly in file %s\n", filename);
	TEST_EXIT(ne > 0)
	  ("number of elements = %d must be bigger than 0\n", ne);

	key_def[3] = true;
	if (key_def[2])
	  fill(mesh, ne, nv);
	break;

      case 4:
	// block "vertex coordinates"
	fscanf(file, "%*s %*s");
	for (int i = 0; i < nv; i++) {
	  for (j = 0; j <Global::getGeo(WORLD) ; j++) {
	    result = fscanf(file, "%lf", &dbl);
607
	    TEST_EXIT(result == 1)
608
609
610
611
612
613
	      ("error while reading coordinates, check file %s\n", filename);
	    coords[i][j] = dbl;
	  }
	}
	key_def[4] = true;
	break;
614

615
616
617
      case 5:
	// block "element vertices"
	fscanf(file, "%*s %*s");
618

619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
	/****************************************************************************/
	/* global number of vertices on a single element                            */
	/****************************************************************************/

	for (int i = 0; i < ne; i++) {
	  result = read_indices(file, *ind);
	  TEST_EXIT(result)
	    ("cannot read vertex indices of element %d in file %s\n",  i, filename);

	  for (k = 0; k < mesh->getGeo(VERTEX); k++)
	    mel_vertex[i][k] = (*ind)[k];
	}

	key_def[5] = true;
	break;

      case 6:
	// block "element boundaries"
	fscanf(file, "%*s %*s");

	/****************************************************************************/
	/* MEL boundary pointers                                                    */
	/****************************************************************************/
	for (int i = 0; i < ne; i++) {  
	  // boundary information of ith element 
644

645
646
647
648
649
650
651
652
653
654
655
	  result = read_indices(file, *ind);
	  TEST_EXIT(result)
	    ("cannot read boundary type of element %d in file %s\n", i, filename);

	  // fill boundary of macro-element
	  MacroReader::fillMelBoundary(mesh, 
				       mel[i], 
				       VecConv<int,NEIGH,PARTS>::convertVec((*ind), mesh));
	}

	this->fillBoundaryInfo(mesh);
656
                   
657
658
659
	bound_set = true;
	key_def[6] = true;
	break;
660

661
      case 7:
662
663
	// block "element neighbours"
	fscanf(file, "%*s %*s");
664

665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
	/****************************************************************************/
	/* fill MEL neighbour pointers:                                             */
	/*   if they are specified in the file: read them from file,                */
	/*   else init them by a call of fill_mel_neighbour()                       */
	/****************************************************************************/
	neigh_set = true;
	for (int i = 0; i < ne; i++) {
	  //  neighbour information about ith element

	  if (read_indices(file, *ind)) {
	    MacroReader::fillMelNeigh(mel[i], mel, 
				      VecConv<int,NEIGH,PARTS>::convertVec((*ind), 
									   mesh));
	  } else {
	    neigh_set = false; /* setting of neighbours fails :-( */
680
	    break;
681
682
	  }
	}
683

684
685
	key_def[7] = true;
	break;
686

687
688
689
690
691
692
      case 8:
	// block "element type"
	fscanf(file, "%*s %*s");
	/****************************************************************************/
	/* MEL elType                                                               */
	/****************************************************************************/
693

694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
	if (dim == 2 || dim == 1)
	  ERROR("there is no element type in 2d and 2d; ignoring data for elType\n");

	for (int i = 0; i < ne; i++) {
	  result = fscanf(file, "%d", &j);
	  TEST_EXIT(result == 1)
	    ("cannot read elType of element %d in file %s\n", i, filename);
	  if (dim == 3)
	    (mel)[i]->elType = j;
	}

	key_def[8] = true;
	break;

      case 9:
	// block "projections"
	{
	  fscanf(file, "%*s");

	  int numFaces = mesh->getGeo(FACE);
	  int numEdgesAtBoundary = 0;

	  for (k = 1; k < dim; k++)
	    numEdgesAtBoundary += k;

	  for (int i = 0; i < ne; i++) {
	    result = read_indices(file, *ind);
	    TEST_EXIT(result)
	      ("cannot read boundary projector of element %d in file %s\n", i, filename);
723
	
724
725
726
727
728
	    Projection *projector = Projection::getProjection((*ind)[0]);

	    if (projector && projector->getType() == VOLUME_PROJECTION) {
	      mel[i]->setProjection(0, projector);
	    } else { // boundary projection
729
	      for (j = 0; j < mesh->getGeo(NEIGH); j++) {
730
		projector = Projection::getProjection((*ind)[j]);
731
		if (projector) {
732
		  mel[i]->setProjection(j, projector);
733
734
		  if (dim > 2) {
		    for (k = 0; k < numEdgesAtBoundary; k++) {
735
736
		      int edgeNr = Global::getReferenceElement(dim)->getEdgeOfFace(j, k);
		      mel[i]->setProjection(numFaces + edgeNr, projector);
737
738
739
740
741
		    }
		  }
		}
	      }
	    }
742
743
744
745
	  }
	}
	key_def[9] = true;
	break;
746

747
748
749
750
751
752
      case 10:
	// block "element region"
	fscanf(file, "%*s %*s");
	/****************************************************************************/
	/* MEL regions                                                              */
	/****************************************************************************/
753

754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
	for (int i = 0; i < ne; i++) {
	  result = fscanf(file, "%d", &j);
	  TEST_EXIT(result == 1)
	    ("cannot read region of element %d in file %s\n", i, filename);
	  if (j >= 0) {
	    Element *el = mel[i]->getElement();
	    ElementRegion_ED *elementRegion = 
	      new ElementRegion_ED(el->getElementData());
	    elementRegion->setRegion(j);
	    el->setElementData(elementRegion);
	  }
	}
	key_def[10] = true;
	break;

      case 11:
	// block "surface region"
	fscanf(file, "%*s %*s");
	for (int i = 0; i < ne; i++) {
	  result = read_indices(file, *ind);
	  TEST_EXIT(result)
	    ("cannot read surface regions of element %d in file %s\n", i, filename);

	  Element *el = mel[i]->getElement();

	  for (j = 0; j < mesh->getGeo(NEIGH); j++) {
	    if ((*ind)[j] >= 0) {
	      SurfaceRegion_ED *surfaceRegion = 
		new SurfaceRegion_ED(el->getElementData());
	      surfaceRegion->setSide(j);
	      surfaceRegion->setRegion((*ind)[j]);
	      el->setElementData(surfaceRegion);
786
787
	    }
	  }
788
789
790
791
792
793
794
795
796
797
798
799
800
	}
	key_def[11] = true;
	break;

      case 12:
	// line "mesh name"
	fscanf(file, "%*s %*s %*s");
	break;

      case 13:
	// line "time"
	fscanf(file, "%*s %*s");
	break;
801
      }
802
    }
803

804
    if (ind)
Thomas Witkowski's avatar
Thomas Witkowski committed
805
      delete ind;
806

807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
    fclose(file);
  }


  int macro_type(const char *filename, const char *type)
  {
    const char *fn, *t;
  
    if (strlen(filename) <= strlen(type))
      return(false);
  
    fn = filename;
    while (*fn) fn++;
    t = type;
    while (*t) t++;
  
    while (t != type  &&  *t == *fn) t--;
  
    return(t == type);
  }


  /****************************************************************************/
  /*  sets the boundary of all edges/faces with no neigbour to a straight     */
  /*  line/face with Dirichlet boundary type                                  */
  /****************************************************************************/

  void MacroInfo::dirichletBoundary()
  {
836
837
    for (int i = 0; i < static_cast<int>( mel.size()); i++) {
      for (int k = 0; k < mesh->getGeo(NEIGH); k++) {
838
839
840
841
842
843
844
845
846
847
848
	if (mel[i]->neighbour[k])
	  mel[i]->boundary[k] = INTERIOR;
	else
	  mel[i]->boundary[k] = DIRICHLET;
      }
    }
  }


  void MacroInfo::fillBoundaryInfo(Mesh *mesh)
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
849
    int i, j, k, nv = mesh->getNumberOfVertices();
850
    std::deque<MacroElement*>::iterator melIt;
Thomas Witkowski's avatar
Thomas Witkowski committed
851
    BoundaryType *bound = new BoundaryType[nv];
852
853
854
855
856
857
858
859
860
    int dim = mesh->getDim();

    switch(dim) {
    case 1:
      break;
    case 2:
      for (i = 0; i < nv; i++)
	bound[i] = INTERIOR;

Thomas Witkowski's avatar
Thomas Witkowski committed
861
      for (i = 0, melIt = mesh->firstMacroElement(); 
862
	   melIt != mesh->endOfMacroElements(); 
Thomas Witkowski's avatar
Thomas Witkowski committed
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
	   ++melIt, i++) {
	for (j = 0; j < mesh->getGeo(NEIGH); j++) {
	  if ((*melIt)->getBoundary(j) != INTERIOR) {
	    if ((*melIt)->getBoundary(j) >= DIRICHLET) {
	      int j1 = mel_vertex[i][(j+1)%3];
	      int j2 = mel_vertex[i][(j+2)%3];
	      
	      bound[j1] = max(bound[j1], (*melIt)->getBoundary(j));
	      bound[j2] = max(bound[j2], (*melIt)->getBoundary(j));
	    } else if ((*melIt)->getBoundary(j) <= NEUMANN) {
	      int j1 = mel_vertex[i][(j+1)%3];
	      int j2 = mel_vertex[i][(j+2)%3];
	      
	      if (bound[j1] != INTERIOR)
		bound[j1] = max(bound[j1], (*melIt)->getBoundary(j));
	      else
		bound[j1] = (*melIt)->getBoundary(j);
	      
	      if (bound[j2] != INTERIOR)
		bound[j2] = max(bound[j2], (*melIt)->getBoundary(j));
	      else
		bound[j2] = (*melIt)->getBoundary(j);
885
886
887
	    }
	  }
	}
Thomas Witkowski's avatar
Thomas Witkowski committed
888
      }
889

Thomas Witkowski's avatar
Thomas Witkowski committed
890
      for (i = 0, melIt = mesh->firstMacroElement(); 
891
892
	   melIt != mesh->endOfMacroElements(); 
	   ++melIt, i++) 
Thomas Witkowski's avatar
Thomas Witkowski committed
893
894
895
	for (j = 0; j < mesh->getGeo(VERTEX); j++)
	  (*melIt)->setBoundary(3 + j, bound[mel_vertex[i][j]]);
	
896
897
898
899
900
      break;
    case 3:
      for (i = 0; i < nv; i++)
	bound[i] = INTERIOR;

Thomas Witkowski's avatar
Thomas Witkowski committed
901
      for (i = 0, melIt = mesh->firstMacroElement(); 
902
	   melIt != mesh->endOfMacroElements(); 
Thomas Witkowski's avatar
Thomas Witkowski committed
903
904
905
906
907
908
909
910
	   ++melIt, i++) {
	for (j = 0; j < mesh->getGeo(NEIGH); j++) {
	  for (k = 1; k < 4; k++)
	    bound[mel_vertex[i][(j+k)%4]] =
	      ((*melIt)->getBoundary(j) != INTERIOR) ?
	      newBound((*melIt)->getBoundary(j),
		       bound[mel_vertex[i][(j+k)%4]]) :
	      bound[mel_vertex[i][(j+k)%4]];
911
	}
Thomas Witkowski's avatar
Thomas Witkowski committed
912
      }
913
914
915

      for (i = 0, melIt = mesh->firstMacroElement(); 
	   melIt != mesh->endOfMacroElements(); 
Thomas Witkowski's avatar
Thomas Witkowski committed
916
917
918
919
	   ++melIt, i++) 	
	for (j = 0; j < mesh->getGeo(VERTEX); j++)
	  (*melIt)->setBoundary(10 + j, bound[mel_vertex[i][j]]);
      
920
      break;
Thomas Witkowski's avatar
Thomas Witkowski committed
921
922
    default: 
      ERROR_EXIT("invalid dim\n");
923
924
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
925
    delete [] bound;
926
927
928
929
  }

  void MacroReader::computeNeighbours(Mesh *mesh)
  {
930
    FUNCNAME("MacroReader::computeNeighbours()");
931

932
    int dim = mesh->getDim();
933
    FixVec<DegreeOfFreedom*, DIMEN> dof(dim, NO_INIT);
934

935
936
937
938
    for (int i = 0; i < mesh->getNumberOfLeaves(); i++) {
      for (int k = 0; k < mesh->getGeo(NEIGH); k++) {
	mesh->getMacroElement(i)->setOppVertex(k, AMDIS_UNDEFINED);
	mesh->getMacroElement(i)->setNeighbour(k, NULL);
939
      }
940
    }
941

942
943
944
945
946
947
948
    for (int i = 0; i < mesh->getNumberOfLeaves(); i++) {
      for (int k = 0; k < mesh->getGeo(NEIGH); k++) {
	if (mesh->getMacroElement(i)->getBoundary(k) != INTERIOR) {
	  mesh->getMacroElement(i)->setNeighbour(k, NULL);
	  mesh->getMacroElement(i)->setOppVertex(k, -1);
	  continue;
	}
949

950
951
952
953
954
955
956
957
	if (mesh->getMacroElement(i)->getOppVertex(k) == AMDIS_UNDEFINED) {
	  if (dim == 1) {
	    dof[0] = const_cast<DegreeOfFreedom*>(mesh->getMacroElement(i)->
						  getElement()->getDOF(k));
	  } else {
	    for (int l = 0; l < dim; l++)
	      dof[l] = const_cast<DegreeOfFreedom*>(mesh->getMacroElement(i)->
						    getElement()->
958
						    getDOF((k + l + 1) % (dim + 1)));
959
	  }
960
961
962
963
964
965
966
967
968
969
970
971
972
	  
	  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;
	    }
	  }

973
974
975
976
977
978
979
	  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);
980
981
	      if (l < dim)
		std::cout << " / ";	      
982
983
984
	    }
	    std::cout << std::endl;
	    std::cout << "  dofs = ";
985
986
	    for (int l = 0; l < dim; l++)
	      std::cout << *(dof[l]) << " ";	    
987
988
989
990
	    std::cout << std::endl;

	    ERROR_EXIT("\n");
	  }    
991
	}
992
      }
993
    }
994
995
996
997
998
999
1000
1001
1002
1003
1004
  }


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

  void MacroReader::boundaryDOFs(Mesh *mesh)
  {
1005
1006
1007
1008
1009
    FUNCNAME("Mesh::boundaryDOFs()");

    int lnode = mesh->getNode(EDGE);
    int k, lne = mesh->getNumberOfLeaves();
    int max_n_neigh = 0, n_neigh, ov;
1010
    std::deque<MacroElement*>::iterator mel;
1011
    const MacroElement* neigh;
1012
    DegreeOfFreedom *dof;
1013
1014
1015
1016
1017
1018

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

    int dim = mesh->getDim();

Thomas Witkowski's avatar
Thomas Witkowski committed
1019
    switch (dim) {
1020
    case 2:
1021
      for (mel = mesh->firstMacroElement(); mel != mesh->endOfMacroElements(); mel++) {
1022
1023
1024
1025
1026
1027
	// check for periodic boundary
	Element *el = const_cast<Element*>((*mel)->getElement());
	ElementData *ed = el->getElementData(PERIODIC);

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

1028
	if (ed) {
1029
	  std::list<LeafDataPeriodic::PeriodicInfo> &periodicInfos = 
1030
	    dynamic_cast<LeafDataPeriodic*>(ed)->getInfoList();
1031
1032
	  std::list<LeafDataPeriodic::PeriodicInfo>::iterator it;
	  std::list<LeafDataPeriodic::PeriodicInfo>::iterator end = periodicInfos.end();
Thomas Witkowski's avatar
Thomas Witkowski committed
1033
1034
	  for (it = periodicInfos.begin(); it != end; ++it)
	    if (it->type != 0)
1035
1036
1037
	      periodic[it->elementSide] = true;
	}

1038
	for (int i = 0; i < mesh->getGeo(NEIGH); i++) {
1039
	  if (!(*mel)->getNeighbour(i) || 
1040
1041
1042
1043
1044
1045
1046
1047
	      ((*mel)->getNeighbour(i)->getIndex() < (*mel)->getIndex())) {

	    mesh->incrementNumberOfEdges(1);

	    if (mesh->getNumberOfDOFs(EDGE)) {
	      dof = el->setDOF(lnode + i, mesh->getDOF(EDGE));
      
	      if ((*mel)->getNeighbour(i)) {
Thomas Witkowski's avatar
Thomas Witkowski committed
1048
1049
		Element *neigh = 
		  const_cast<Element*>((*mel)->getNeighbour(i)->getElement());
1050

Thomas Witkowski's avatar
Thomas Witkowski committed
1051
		if (periodic[i])
1052
		  neigh->setDOF(lnode + (*mel)->getOppVertex(i), mesh->getDOF(EDGE));
Thomas Witkowski's avatar
Thomas Witkowski committed
1053
1054
		else
		  neigh->setDOF(lnode + (*mel)->getOppVertex(i), dof);		
1055
	      }
1056
1057
	    }
	  }  
1058
1059
1060
1061
1062
1063
	}
      }
      break;
    case 3:
      lnode = mesh->getNode(FACE);
      mel = mesh->firstMacroElement();
1064
      for (int i = 0; i < lne; i++) {
1065
1066
1067
1068
1069
1070
1071

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

	DimVec<bool> periodic(dim, DEFAULT_VALUE, false);
      
1072
	if (ed) {
1073
	  std::list<LeafDataPeriodic::PeriodicInfo> &periodicInfos = 
1074
	    dynamic_cast<LeafDataPeriodic*>(ed)->getInfoList();
1075
1076
	  std::list<LeafDataPeriodic::PeriodicInfo>::iterator it;
	  std::list<LeafDataPeriodic::PeriodicInfo>::iterator end = periodicInfos.end();
1077
1078
	  for (it = periodicInfos.begin(); it != end; ++it)
	    if (it->type != 0)
1079
1080
1081
1082
1083
1084
1085
1086
1087
	      periodic[it->elementSide] = true;
	}

	for (k = 0; k < mesh->getGeo(EDGE); k++) {      
	  /*********************************************************************/
	  /* check for not counted edges                                       */
	  /*********************************************************************/
	  n_neigh = 1;

Thomas Witkowski's avatar
Thomas Witkowski committed
1088
	  if (newEdge(mesh, (*(mel + i)), k, &n_neigh)) {
1089
1090
1091
1092
1093
1094
	    mesh->incrementNumberOfEdges(1);
	    max_n_neigh = max(max_n_neigh, n_neigh);
	  }
	}
      
	for (k = 0; k < mesh->getGeo(NEIGH); k++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
1095
	  neigh = (*(mel + i))->getNeighbour(k);
1096
1097
1098
	  /*********************************************************************/
	  /* face is counted and dof is added by the element with bigger index */
	  /*********************************************************************/
Thomas Witkowski's avatar
Thomas Witkowski committed
1099
1100
	  if (neigh && (neigh->getIndex() > (*(mel + i))->getIndex()))  
	    continue;
1101
1102
1103
1104
	
	  mesh->incrementNumberOfFaces(1);
	
	  if (mesh->getNumberOfDOFs(FACE)) {
Thomas Witkowski's avatar
Thomas Witkowski committed
1105
	    TEST_EXIT(!(*(mel + i))->getElement()->getDOF(lnode + k))
1106
	      ("dof %d on element %d already set\n", 
Thomas Witkowski's avatar
Thomas Witkowski committed
1107
	       lnode + k, (*(mel + i))->getIndex());
1108
	  
Thomas Witkowski's avatar
Thomas Witkowski committed
1109
1110
	    const_cast<Element*>((*(mel + i))->getElement())->setDOF(lnode + k,
								     mesh->getDOF(FACE));
1111
1112

	    if (neigh) {
Thomas Witkowski's avatar
Thomas Witkowski committed
1113
1114
	      ov = (*(mel + i))->getOppVertex(k);
	      TEST_EXIT(!neigh->getElement()->getDOF(lnode + ov))
1115
		("dof %d on neighbour %d already set\n", 
Thomas Witkowski's avatar
Thomas Witkowski committed
1116
		 lnode + ov, neigh->getIndex());
1117
1118
	    
	      Element *neighEl = 
Thomas Witkowski's avatar
Thomas Witkowski committed
1119
		const_cast<Element*>((*(mel + i))->getNeighbour(k)->getElement());
1120

Thomas Witkowski's avatar
Thomas Witkowski committed
1121
	      if (periodic[k])
1122
		neighEl->setDOF(lnode+ov, mesh->getDOF(FACE));
Thomas Witkowski's avatar
Thomas Witkowski committed
1123
1124
1125
	      else
		neighEl->setDOF(lnode+ov, const_cast<int*>((*(mel + i))->getElement()->
							   getDOF(lnode + k)));	      
1126
1127
1128
1129
1130
	    }
	  }
	}
      }
      break;
Thomas Witkowski's avatar
Thomas Witkowski committed
1131
1132
    default: 
      ERROR_EXIT("invalid dim\n");
1133
1134
    }
    
1135
    if (3 == dim)
Thomas Witkowski's avatar
Thomas Witkowski committed
1136
      mesh->setMaxEdgeNeigh(std::max(8, 2 * max_n_neigh));
1137
    else
Thomas Witkowski's avatar
Thomas Witkowski committed
1138
      mesh->setMaxEdgeNeigh(dim - 1);    
1139
1140
1141
1142
1143
1144
1145
1146
1147
1148
1149
1150
1151
  }

  /* 
     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 
     (wenn nameneu=NULL wird keine MAcro-Datei erzeugt)
  */      

  void MacroReader::macroTest(Mesh *mesh, const char *nameneu)
  {
1152
    FUNCNAME("MacroReader::macroTest()");
1153
   
1154
1155
1156
1157
1158
1159
1160
    int i = macrotest(mesh);

    if (i >= 0) {
      ERROR("There is a cycle beginning in macro element %d\n", i);
      ERROR("Entries in MacroElement structures get reordered\n");
      umb(NULL, mesh, umbVkantMacro);

Thomas Witkowski's avatar
Thomas Witkowski committed
1161
      if (nameneu)