MacroReader.cc 45.1 KB
Newer Older
1
2
3
#include "MacroReader.h"
#include "MacroWriter.h"
#include "MacroElement.h"
4
#include "MacroInfo.h"
5
6
7
#include "Boundary.h"
#include "FiniteElemSpace.h"
#include "Mesh.h"
8
#include <string.h>
9
10
11
12
13
14
15
16
17
18
19
20
21
#include "FixVec.h"
#include "PeriodicMap.h"
#include "ElInfo.h"
#include "Parameters.h"
#include "DOFIterator.h"
#include "LeafData.h"
#include "VertexVector.h"
#include <map>
#include <iostream>
#include <fstream>

namespace AMDiS {

22
  MacroInfo* MacroReader::readMacro(std::string filename, 
23
				    Mesh* mesh,
24
				    std::string periodicFile,
25
26
				    int check)
  {
27
    FUNCNAME("MacroReader::readMacro()");
28

29
30
    TEST_EXIT(filename != "")("no file specified; filename NULL pointer\n");  

Thomas Witkowski's avatar
Thomas Witkowski committed
31
    MacroInfo *macroInfo = new MacroInfo();
32
33
    macroInfo->readAMDiSMacro(filename, mesh);

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

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

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

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

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

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

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

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

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

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

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

127
128
129
130
131
132
133
134
	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());

135
136
137
138
	LeafDataPeriodic *ldp1 = 
	  dynamic_cast<LeafDataPeriodic*>(ld1->getElementData(PERIODIC));
	LeafDataPeriodic *ldp2 = 
	  dynamic_cast<LeafDataPeriodic*>(ld2->getElementData(PERIODIC));
139
140

	if (!ldp1) {
Thomas Witkowski's avatar
Thomas Witkowski committed
141
	  ldp1 = new LeafDataPeriodic(ld1);
142
143
144
145
	  element1->setElementData(ldp1);
	}

	if (!ldp2) {
Thomas Witkowski's avatar
Thomas Witkowski committed
146
	  ldp2 = new LeafDataPeriodic(ld2);
147
148
149
	  element2->setElementData(ldp2);
	}

150
151
	ldp1->addPeriodicInfo(mode, boundaryType, sideEl1, &periodicCoordsEl1);
	ldp2->addPeriodicInfo(mode, boundaryType, sideEl2, &periodicCoordsEl2);
152
153
154

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

156
	  if (!associated) {
Thomas Witkowski's avatar
Thomas Witkowski committed
157
	    associated = new VertexVector(mesh->getVertexAdmin(), "vertex vector");
158
159
	    mesh->periodicAssociations[boundaryType] = associated;
	    VertexVector::Iterator it(associated, ALL_DOFS);
160
	    for (it.reset2(); !it.end(); ++it)
161
162
163
	      *it = it.getDOFIndex();
	  }

164
	  for (int j = 0; j < dim; j++) {
165
166
167
168
169
170
171
172
	    (*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
173
174
      delete [] verticesEl1;
      delete [] verticesEl2;
175
176

      // change periodic vertex dofs
177
      for (int i = 0; i < mesh->getNumberOfVertices(); i++) {
178
	if (periodicMap.getEntry(i) != -1) {
179
	  mesh->freeDof(dof[i], VERTEX);
180
181
	  dof[i] = dof[periodicMap.getEntry(i)];

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

186
	  for (assoc = mesh->periodicAssociations.begin(); assoc != assocEnd; ++assoc) {
187
188
189
190
191

	    DegreeOfFreedom a = (*(assoc->second))[i];
	    if (a != i) {
	      (*(assoc->second))[i] = i;
	      (*(assoc->second))[a] = periodicMap.getEntry(i);
192
	    }
193
194
	  }

195
196
197
	}
      }

198
199
      std::map<BoundaryType, VertexVector*>::iterator assoc;
      std::map<BoundaryType, VertexVector*>::iterator assocEnd =
200
	mesh->periodicAssociations.end();
201
      for (assoc = mesh->periodicAssociations.begin(); assoc != assocEnd; ++assoc) {
Thomas Witkowski's avatar
Thomas Witkowski committed
202
	for (int i = 0; i < mesh->getNumberOfVertices(); i++)
203
204
	  if (i != (*(assoc->second))[i])
	    MSG("association %d: vertex %d -> vertex %d\n", 
Thomas Witkowski's avatar
Thomas Witkowski committed
205
		assoc->first, i, (*(assoc->second))[i]);	
206
207
      }

208
209
210
211
      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));
212
213
214

    } // periodicFile

215

216
217
    // =========================================================

218
219
220
    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]]);
221

222
	const_cast<Element*>((*(mel + i))->getElement())->
223
	  setDOF(k, dof[melVertex[i][k]]);
224
      }
225
    }
226

227
    if (!macroInfo->neigh_set) {
228
      TEST_EXIT(periodicFile == "")
229
	("periodic boundary condition => element neighbours must be set\n");
230
      computeNeighbours(mesh);
231
    } else {
232
233
234
      /****************************************************************************/
      /* fill MEL oppVertex values when reading neighbour information form file  */
      /****************************************************************************/
235

236
237
238
239
240
      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) {
241
242
	    int j = 0;
	    for (; j < mesh->getGeo(NEIGH); j++)
243
244
	      if (neigh->getNeighbour(j) == *(mel + i))  
		break;
245
	
246
247
248
249
250
	    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);
251
	  }
252
	}
253
      }
254
    }
255

256
257
    if (!macroInfo->bound_set)
      macroInfo->dirichletBoundary();    
258
  
259
    if (mesh->getDim() > 1)
260
261
262
263
264
265
      boundaryDOFs(mesh);

    // initial boundary projections
    int numFaces = mesh->getGeo(FACE);
    int dim = mesh->getDim();
    mel = mesh->firstMacroElement();
266
    for (int i = 0; i < mesh->getNumberOfLeaves(); i++) {
267
268
      MacroElement *macroEl = *(mel+i);
      Projection *projector = macroEl->getProjection(0);
269
      if (projector && projector->getType() == VOLUME_PROJECTION) {
Thomas Witkowski's avatar
Thomas Witkowski committed
270
271
	for (int j = 0; j <= dim; j++)
	  projector->project(macroEl->getCoord(j));	
272
      } else {
273
	for (int j = 0; j < mesh->getGeo(EDGE); j++) {
274
	  projector = macroEl->getProjection(numFaces + j);
275
	  if (projector) {
276
277
278
279
280
281
282
283
284
285
286
287
	    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)) {
288
      for (int i = 0; i < mesh->getNumberOfMacros(); i++)
289
	const_cast<Element*>(mel[i]->getElement())->
290
	  setDOF(mesh->getNode(CENTER), mesh->getDof(CENTER));
291
292
293
294
295
296
    }

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

297
    WorldVector<double> x_min, x_max;
298

299
300
301
302
303
304
305
    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++) {
306
307
	x_min[j] = std::min(x_min[j], coords[i][j]);
	x_max[j] = std::max(x_max[j], coords[i][j]);
308
309
      }
    }
310

311
    for (int j = 0; j < Global::getGeo(WORLD); j++)
312
313
314
315
316
317
      mesh->setDiameter(j, x_max[j] - x_min[j]);

    if (check) {
      checkMesh(mesh);

      if (mesh->getDim() > 1) {
318
	char filenew[128];
319
	strncpy(filenew, filename.c_str(), 128); 
320
321
322
	filenew[127] = 0;
	strncat(filenew, ".new", 128);   
	filenew[127] = 0;
323
324
325
326
	macroTest(mesh, filenew);
      }
    }

327
    return macroInfo;
328
329
330
331
332
  }


  void MacroReader::computeNeighbours(Mesh *mesh)
  {
333
    FUNCNAME("MacroReader::computeNeighbours()");
334

335
    int dim = mesh->getDim();
336
    FixVec<DegreeOfFreedom*, DIMEN> dof(dim, NO_INIT);
337

338
339
340
341
    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);
342
      }
343
    }
344

345
346
347
348
349
350
351
    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;
	}
352

353
354
355
356
357
358
359
360
	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()->
361
						    getDOF((k + l + 1) % (dim + 1)));
362
	  }
363
364
365
366
367
368
369
370
371
372
373
374
375
	  
	  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;
	    }
	  }

376
377
378
379
380
381
382
	  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);
383
384
	      if (l < dim)
		std::cout << " / ";	      
385
386
387
	    }
	    std::cout << std::endl;
	    std::cout << "  dofs = ";
388
389
	    for (int l = 0; l < dim; l++)
	      std::cout << *(dof[l]) << " ";	    
390
391
392
393
	    std::cout << std::endl;

	    ERROR_EXIT("\n");
	  }    
394
	}
395
      }
396
    }
397
398
399
400
401
402
403
404
405
406
407
  }


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

  void MacroReader::boundaryDOFs(Mesh *mesh)
  {
408
409
410
411
412
    FUNCNAME("Mesh::boundaryDOFs()");

    int lnode = mesh->getNode(EDGE);
    int k, lne = mesh->getNumberOfLeaves();
    int max_n_neigh = 0, n_neigh, ov;
413
    std::deque<MacroElement*>::iterator mel;
414
    const MacroElement* neigh;
415
    DegreeOfFreedom *dof;
416
417
418
419
420
421

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

    int dim = mesh->getDim();

Thomas Witkowski's avatar
Thomas Witkowski committed
422
    switch (dim) {
423
    case 2:
424
      for (mel = mesh->firstMacroElement(); mel != mesh->endOfMacroElements(); mel++) {
425
426
427
428
429
430
	// check for periodic boundary
	Element *el = const_cast<Element*>((*mel)->getElement());
	ElementData *ed = el->getElementData(PERIODIC);

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

431
	if (ed) {
432
	  std::list<LeafDataPeriodic::PeriodicInfo> &periodicInfos = 
433
	    dynamic_cast<LeafDataPeriodic*>(ed)->getInfoList();
434
435
	  std::list<LeafDataPeriodic::PeriodicInfo>::iterator it;
	  std::list<LeafDataPeriodic::PeriodicInfo>::iterator end = periodicInfos.end();
Thomas Witkowski's avatar
Thomas Witkowski committed
436
437
	  for (it = periodicInfos.begin(); it != end; ++it)
	    if (it->type != 0)
438
439
440
	      periodic[it->elementSide] = true;
	}

441
	for (int i = 0; i < mesh->getGeo(NEIGH); i++) {
442
	  if (!(*mel)->getNeighbour(i) || 
443
444
445
446
447
	      ((*mel)->getNeighbour(i)->getIndex() < (*mel)->getIndex())) {

	    mesh->incrementNumberOfEdges(1);

	    if (mesh->getNumberOfDOFs(EDGE)) {
448
	      dof = el->setDOF(lnode + i, mesh->getDof(EDGE));
449
450
      
	      if ((*mel)->getNeighbour(i)) {
Thomas Witkowski's avatar
Thomas Witkowski committed
451
452
		Element *neigh = 
		  const_cast<Element*>((*mel)->getNeighbour(i)->getElement());
453

Thomas Witkowski's avatar
Thomas Witkowski committed
454
		if (periodic[i])
455
		  neigh->setDOF(lnode + (*mel)->getOppVertex(i), mesh->getDof(EDGE));
Thomas Witkowski's avatar
Thomas Witkowski committed
456
457
		else
		  neigh->setDOF(lnode + (*mel)->getOppVertex(i), dof);		
458
	      }
459
460
	    }
	  }  
461
462
463
464
465
466
	}
      }
      break;
    case 3:
      lnode = mesh->getNode(FACE);
      mel = mesh->firstMacroElement();
467
      for (int i = 0; i < lne; i++) {
468
469
470
471
472
473
474

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

	DimVec<bool> periodic(dim, DEFAULT_VALUE, false);
      
475
	if (ed) {
476
	  std::list<LeafDataPeriodic::PeriodicInfo> &periodicInfos = 
477
	    dynamic_cast<LeafDataPeriodic*>(ed)->getInfoList();
478
479
	  std::list<LeafDataPeriodic::PeriodicInfo>::iterator it;
	  std::list<LeafDataPeriodic::PeriodicInfo>::iterator end = periodicInfos.end();
480
481
	  for (it = periodicInfos.begin(); it != end; ++it)
	    if (it->type != 0)
482
483
484
485
486
487
488
489
490
	      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
491
	  if (newEdge(mesh, (*(mel + i)), k, &n_neigh)) {
492
493
494
495
496
497
	    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
498
	  neigh = (*(mel + i))->getNeighbour(k);
499
500
501
	  /*********************************************************************/
	  /* face is counted and dof is added by the element with bigger index */
	  /*********************************************************************/
Thomas Witkowski's avatar
Thomas Witkowski committed
502
503
	  if (neigh && (neigh->getIndex() > (*(mel + i))->getIndex()))  
	    continue;
504
505
506
507
	
	  mesh->incrementNumberOfFaces(1);
	
	  if (mesh->getNumberOfDOFs(FACE)) {
Thomas Witkowski's avatar
Thomas Witkowski committed
508
	    TEST_EXIT(!(*(mel + i))->getElement()->getDOF(lnode + k))
509
	      ("dof %d on element %d already set\n", 
Thomas Witkowski's avatar
Thomas Witkowski committed
510
	       lnode + k, (*(mel + i))->getIndex());
511
	  
Thomas Witkowski's avatar
Thomas Witkowski committed
512
	    const_cast<Element*>((*(mel + i))->getElement())->setDOF(lnode + k,
513
								     mesh->getDof(FACE));
514
515

	    if (neigh) {
Thomas Witkowski's avatar
Thomas Witkowski committed
516
517
	      ov = (*(mel + i))->getOppVertex(k);
	      TEST_EXIT(!neigh->getElement()->getDOF(lnode + ov))
518
		("dof %d on neighbour %d already set\n", 
Thomas Witkowski's avatar
Thomas Witkowski committed
519
		 lnode + ov, neigh->getIndex());
520
521
	    
	      Element *neighEl = 
Thomas Witkowski's avatar
Thomas Witkowski committed
522
		const_cast<Element*>((*(mel + i))->getNeighbour(k)->getElement());
523

Thomas Witkowski's avatar
Thomas Witkowski committed
524
	      if (periodic[k])
525
		neighEl->setDOF(lnode+ov, mesh->getDof(FACE));
Thomas Witkowski's avatar
Thomas Witkowski committed
526
527
528
	      else
		neighEl->setDOF(lnode+ov, const_cast<int*>((*(mel + i))->getElement()->
							   getDOF(lnode + k)));	      
529
530
531
532
533
	    }
	  }
	}
      }
      break;
Thomas Witkowski's avatar
Thomas Witkowski committed
534
535
    default: 
      ERROR_EXIT("invalid dim\n");
536
537
    }
    
538
    if (3 == dim)
Thomas Witkowski's avatar
Thomas Witkowski committed
539
      mesh->setMaxEdgeNeigh(std::max(8, 2 * max_n_neigh));
540
    else
Thomas Witkowski's avatar
Thomas Witkowski committed
541
      mesh->setMaxEdgeNeigh(dim - 1);    
542
543
544
545
546
547
548
549
550
551
552
553
554
  }

  /* 
     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)
  {
555
    FUNCNAME("MacroReader::macroTest()");
556
   
557
558
559
560
561
562
563
    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
564
      if (nameneu)
565
566
	ERROR_EXIT("mesh->feSpace\n");
    }
567
568
569
570
571
572
573
574
575
576
577
578
579
580
  }
  
  /****************************************************************************/
  /*  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)
  {
581
    FUNCNAME("MacroReader::macrotest()");
582

Thomas Witkowski's avatar
Thomas Witkowski committed
583
584
    std::deque<MacroElement*>::const_iterator macro, mac;
    int flg = 0;
585
    int dim = mesh->getDim();
Thomas Witkowski's avatar
Thomas Witkowski committed
586
587
    int *test = new int[mesh->getNumberOfMacros()];
    int *zykl = new int[mesh->getNumberOfMacros()];
588
 
Thomas Witkowski's avatar
Thomas Witkowski committed
589
    for (int i = 0; i < mesh->getNumberOfMacros(); i++)
590
      test[i] = 0;
591

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

595
596
597
598
    while (macrolfd != mesh->endOfMacroElements()) {
      if (test[(*macrolfd)->getIndex()] == 1) {
	macrolfd++;
      } else {
599
600
	for (int i = 0; i < mesh->getNumberOfMacros(); i++)
	  zykl[i] = 0;	
601
    
602
603
604
605
606
607
608
609
	macro = macrolfd;
	flg = 2;
	do {
	  if (zykl[(*macro)->getIndex()] == 1) {
	    flg = 0;
	    zykstart = (*macro)->getIndex();
	  } else {
	    zykl[(*macro)->getIndex()] = 1;
610
      
611
612
613
614
615
	    if (test[(*macro)->getIndex()] == 1) {
	      flg = 1;
	    } else if ((*macro)->getNeighbour(dim) == NULL) {
	      flg = 1;
	      test[(*macro)->getIndex()] = 1;
Thomas Witkowski's avatar
Thomas Witkowski committed
616
	    } else if ((*macro) == (*macro)->getNeighbour(dim)->getNeighbour(dim)) {
617
618
619
620
	      flg = 1;
	      test[(*macro)->getIndex()] = 1;
	      test[(*macro)->getNeighbour(dim)->getIndex()] = 1;
	    } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
621
622
	      for (mac = mesh->firstMacroElement(); 
		   (*mac) != (*macro)->getNeighbour(dim); mac++);
623
624
625
626
	      macro = mac;
	    } 
	  }	  
	} while(flg == 2);
627
 
628
	if (flg == 1)
629
	  macrolfd++;
630
631
	else 
	  macrolfd=mesh->endOfMacroElements();	
632
      }
633
    }
634
  
Thomas Witkowski's avatar
Thomas Witkowski committed
635
636
    delete [] zykl;
    delete [] test;
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
 
    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
661
			void (*umbvk)(Mesh*, MacroElement*, int, int*))
662
663
664
  {
    FUNCNAME("MacroReader::umb");

Thomas Witkowski's avatar
Thomas Witkowski committed
665
    int *test = new int[mesh->getNumberOfMacros()];
666
  
Thomas Witkowski's avatar
Thomas Witkowski committed
667
668
    for (int i = 0; i < static_cast<int>(mesh->getNumberOfMacros()); i++)
      test[i] = 0;
669

Thomas Witkowski's avatar
Thomas Witkowski committed
670
    recumb(mesh, (*mesh->firstMacroElement()), NULL, test, 0, 0, ele, umbvk);
671

Thomas Witkowski's avatar
Thomas Witkowski committed
672
    delete [] test;
673
674
675
676
677
  }

  bool MacroReader::newEdge(Mesh *mesh, MacroElement *mel,
			    int mel_edge_no, int *n_neigh)
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
678
679
    FUNCNAME("MacroElement::newEdge()"); 
    MacroElement *nei;
680
681
    const DegreeOfFreedom *dof[2];
    DegreeOfFreedom *edge_dof = NULL;
Thomas Witkowski's avatar
Thomas Witkowski committed
682
683
    int j, k, opp_v, node = 0;
    BoundaryType lbound = INTERIOR;
684
    Projection *lproject = NULL;
Thomas Witkowski's avatar
Thomas Witkowski committed
685
    const int max_no_list_el = 100;
686
687
688
689
    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
690
    static int next_el[6][2] = {{3,2},{1,3},{1,2},{0,3},{0,2},{0,1}};
691
692
    int vertices = mesh->getGeo(VERTEX);

Thomas Witkowski's avatar
Thomas Witkowski committed
693
    int mel_index = mel->getIndex();
694
695
696
697
698
699
700
701
702
703
704
705
706

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

    if (mesh->getNumberOfDOFs(EDGE)) {
      node = mesh->getNode(EDGE);
      if (el->getDOF(node+edge_no)) {
	/****************************************************************************/
	/*  edge was counted by another macro element and dof was added on the      */
	/*  complete patch                                                          */
	/****************************************************************************/
	return false;
      } else {
707
	edge_dof = el->setDOF(node+edge_no,mesh->getDof(EDGE));
708
709
710
      }
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
711
    for (j = 0; j < 2; j++)
712
      dof[j] = el->getDOF(el->getVertexOfEdge(edge_no, j));
Thomas Witkowski's avatar
Thomas Witkowski committed
713
    
714
715
716
717
718
719
720
721
722
723
724
725
    /****************************************************************************/
    /*  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]);


726
    if (mel->getBoundary(next_el[edge_no][0])) {
727
728
729
730
731
732
733
734
735
736
737
      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++)
	if (nei->getElement()->getDOF(j) == dof[0])  break;
      for (k = 0; k < vertices; k++)
	if (nei->getElement()->getDOF(k) == dof[1])  break;

      // check for periodic boundary
738
      if (j == 4 || k == 4) {
739
740
741
742
	nei = NULL;
	break;
      }

743
      if (mesh->getNumberOfDOFs(EDGE)) {
744
745
	TEST_EXIT(nei->index > mel_index)
	  ("neighbour index %d < element index %d\n", nei->getIndex(), mel_index);
746
      }
747

748
749
      if (!mesh->getNumberOfDOFs(EDGE) && nei->getIndex() < mel_index) 
	return false;
750
    
751
      edge_no = Tetrahedron::edgeOfDofs[j][k];
752
753
754
755
756
757
758
759
760
761
762

      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]);

      if (mesh->getNumberOfDOFs(EDGE))
763
	nei->element->setDOF(node+edge_no,edge_dof);      
764

765
      if (next_el[edge_no][0] != opp_v) {
766
	if (nei->getBoundary(next_el[edge_no][0])) {
767
768
769
770
771
772
	  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()))
773
774
775
	      lproject = neiProject;
	  }
	}
776
777
778
779
780
781
782
783
784
785
	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()))
786
787
788
	      lproject = neiProject;
	  }
	}
789
790
791
	opp_v = nei->getOppVertex(next_el[edge_no][1]);
	nei = nei->getNeighbour(next_el[edge_no][1]);
      }
792
793
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
    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++)
	  if (nei->getElement()->getDOF(j) == dof[0])  break;
	for (k = 0; k < vertices; k++)
	  if (nei->getElement()->getDOF(k) == dof[1])  break;
819
	
Thomas Witkowski's avatar
Thomas Witkowski committed
820
821
822
823
	// check for periodic boundary
	if (j == 4 || k == 4)
	  return false;
	
824
	if (mesh->getNumberOfDOFs(EDGE)) {
Thomas Witkowski's avatar
Thomas Witkowski committed
825
826
827
	  TEST_EXIT(nei->getIndex() > mel_index)
	    ("neighbour index %d < element index %d\n", nei->getIndex(),
	     mel_index);
828
	}
Thomas Witkowski's avatar
Thomas Witkowski committed
829
830
831
832
	
	if (nei->getIndex() < mel_index)  
	  return false;
	
833
	edge_no = Tetrahedron::edgeOfDofs[j][k];
Thomas Witkowski's avatar
Thomas Witkowski committed
834
835
836
837
838
	
	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]);
839

Thomas Witkowski's avatar
Thomas Witkowski committed
840
841
842
843
844
	if (mesh->getNumberOfDOFs(EDGE)) {
	  TEST_EXIT(!nei->getElement()->getDOF(node+edge_no))
	    ("dof %d on element %d is already set, but not on element %d\n",
	     node + edge_no, nei->getIndex(), mel_index);
	  
845
	  nei->element->setDOF(node+edge_no, edge_dof);
Thomas Witkowski's avatar
Thomas Witkowski committed
846
	}
847

Thomas Witkowski's avatar
Thomas Witkowski committed
848
849
850
851
852
853
854
855
856
	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;	      
857
858
	    }
	  }
Thomas Witkowski's avatar
Thomas Witkowski committed
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
	  
	  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]);
	}
877
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
878
879
    }
    
880
881
882
883
884
885
886
887
888
889
890
    for (j = 0; j < *n_neigh; j++) {
      *(list_bound[j]) = lbound;
      *(list_project[j]) = lproject;
    }
  
    return true;
  }

  void MacroReader::fillMelBoundary(Mesh *mesh, MacroElement *mel, 
				    FixVec<BoundaryType ,NEIGH> ind)
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
891
892
    for (int i = 0; i < mesh->getGeo(NEIGH); i++)
      mel->boundary[i] = ind[i];        
893
894
895
896
  }


  void MacroReader::fillMelNeigh(MacroElement *mel,
897
				 std::deque<MacroElement*>& macro_elements, 
898
899
900
901
				 FixVec<int,NEIGH> ind)
  {
    int dim = mel->element->getMesh()->getDim();

902
903
904
905
906
907
    for (int k = 0; k < Global::getGeo(NEIGH, dim); k++) {
      if (ind[k] >= 0) 
	mel->neighbour[k] = macro_elements[ind[k]];
      else
	mel->neighbour[k] = NULL;
    }
908
909
910
911
912
913
914
915
916
917
918
919
920
921
  }


  //   ordnet Eintraege in Macro-Element macro bzgl. Verfeinerungskante ka um
  //   (coord, bound, boundary, neigh, oppVertex)

  //   ordnet Eintraege in macro->el bzgl. Verfeinerungskante ka um
  //   (Element-DOF's (macro->el->dof) in Ecken und auf Kanten,
  //    wenn NEIGH_IN_EL macro->el->neigh, macro->el->oppVertex)
  //   (wird fuer ALBERT-Routine write_macro benoetigt)

  //   ele wird nicht benoetigt (es kann NULL uebergeben werden)    
  void MacroReader::umbVkantMacro(Mesh *mesh, MacroElement* me, int ka, int *)
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
922
    MacroElement* macr=new MacroElement(mesh->getDim());
923
924
925
926
927
928
929
930
931
    int i;
    int n0;
    DegreeOfFreedom *d[7];
  
    int vertices = mesh->getGeo(VERTEX);
    int facesPlusEdges = mesh->getGeo(EDGE) + mesh->getGeo(FACE);

    if (ka == 2);
    else { 
Thomas Witkowski's avatar
Thomas Witkowski committed
932
      for (i = 0; i < 3; i++) {
933
934
935
936
937
938
939
	macr->coord[i]=me->coord[i];
	macr->setBoundary(facesPlusEdges + i, me->getBoundary(facesPlusEdges + i));
	macr->setBoundary(i, me->getBoundary(i));
	macr->setNeighbour(i, me->getNeighbour(i));
	macr->setOppVertex(i,me->getOppVertex(i));
      }    
  
Thomas Witkowski's avatar
Thomas Witkowski committed
940
941
      for (i = 0; i < 7; i++)
	d[i] = const_cast<int*>(me->getElement()->getDOF(i));      
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964

      if (ka == 1) { 
	me->coord[0] = macr->coord[2];
	me->coord[1] = macr->coord[0];
	me->coord[2] = macr->coord[1];

	me->setBoundary(facesPlusEdges + 0,macr->getBoundary(facesPlusEdges + 2));
	me->setBoundary(facesPlusEdges + 1,macr->getBoundary(facesPlusEdges + 0));
	me->setBoundary(facesPlusEdges + 2,macr->getBoundary(facesPlusEdges + 1));

	me->setBoundary(0, macr->getBoundary(2));
	me->setBoundary(1, macr->getBoundary(0));
	me->setBoundary(2, macr->getBoundary(1));

	me->setNeighbour(0,const_cast<MacroElement*>(macr->getNeighbour(2)));
	me->setNeighbour(1,const_cast<MacroElement*>(macr->getNeighbour(0)));
	me->setNeighbour(2,const_cast<MacroElement*>(macr->getNeighbour(1)));

	me->setOppVertex(0,macr->getOppVertex(2));
	me->setOppVertex(1,macr->getOppVertex(0));
	me->setOppVertex(2,macr->getOppVertex(1));


Thomas Witkowski's avatar
Thomas Witkowski committed
965
966
	if (mesh->getNumberOfDOFs(VERTEX)) {                /* Ecken */
	  n0 = mesh->getNode(VERTEX);              
967