Liebe Gitlab-Nutzer, lieber Gitlab-Nutzer, es ist nun möglich sich mittels des ZIH-Logins/LDAP an unserem Dienst anzumelden. Ein Anmelden über dieses erzeugt ein neues Konto. Das alte Konto ist über den Reiter "Standard" erreichbar. Die Administratoren

Dear Gitlab user, it is now possible to log in to our service using the ZIH login/LDAP. Logging in via this will create a new account. The old account can be accessed via the "Standard" tab. The administrators

MacroReader.cc 45.2 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12
//
// Software License for AMDiS
//
// Copyright (c) 2010 Dresden University of Technology 
// All rights reserved.
// Authors: Simon Vey, Thomas Witkowski et al.
//
// This file is part of AMDiS
//
// See also license.opensource.txt in the distribution.


13 14 15 16
#include <string.h>
#include <map>
#include <iostream>
#include <fstream>
17 18 19
#include "MacroReader.h"
#include "MacroWriter.h"
#include "MacroElement.h"
20
#include "MacroInfo.h"
21 22 23 24 25
#include "Boundary.h"
#include "FiniteElemSpace.h"
#include "Mesh.h"
#include "FixVec.h"
#include "ElInfo.h"
26
#include "Initfile.h"
27 28 29 30 31 32
#include "DOFIterator.h"
#include "LeafData.h"
#include "VertexVector.h"

namespace AMDiS {

33
  MacroInfo* MacroReader::readMacro(std::string filename, 
34
				    Mesh* mesh,
35
				    std::string periodicFile,
36 37
				    int check)
  {
38
    FUNCNAME("MacroReader::readMacro()");
39

40 41
    TEST_EXIT(filename != "")("no file specified; filename NULL pointer\n");  

Thomas Witkowski's avatar
Thomas Witkowski committed
42
    MacroInfo *macroInfo = new MacroInfo();
43 44
    macroInfo->readAMDiSMacro(filename, mesh);

45
    std::deque<MacroElement*>::iterator mel = macroInfo->mel.begin();
46 47 48
    int **melVertex = macroInfo->mel_vertex;
    WorldVector<double> *coords = macroInfo->coords;
    DegreeOfFreedom **dof = macroInfo->dof;
49 50

    // === read periodic data =================================
51
    if (periodicFile != "") {
52
      WARNING("Periodic boundaries may lead to errors in small meshes if element neighbours are not set!\n");
53
    
54 55
      FILE *file = fopen(periodicFile.c_str(), "r");
      TEST_EXIT(file)("can't open file %s\n", periodicFile.c_str());
56 57 58 59 60

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

      int el1, el2;
Thomas Witkowski's avatar
Thomas Witkowski committed
61 62
      int *verticesEl1 = new int[dim];
      int *verticesEl2 = new int[dim];
63
      int mode = -1; // 0: drop dofs, 1: associate dofs
64
      int result;
65
      BoundaryType boundaryType;
66
      MacroReader::PeriodicMap periodicMap;
67 68 69

      fscanf(file, "%*s %d", &n);
      fscanf(file, "%*s %*s %*s %*s %*s %*s %*s %*s %*s %*s %*s");
70
   
71
      for (int i = 0; i < n; i++) {
72 73
	std::map<int, int> vertexMapEl1;
	std::map<int, int> vertexMapEl2;
74

75 76
	result = fscanf(file, "%d", &mode);
	TEST_EXIT(result == 1)("mode?\n");
77
      
78 79
	result = fscanf(file, "%d", &boundaryType);
	TEST_EXIT(result == 1)("boundaryType?\n");
80
      
81 82 83
	result = fscanf(file, "%d", &el1);
	TEST_EXIT(result == 1)("el1?\n");

84
	for (int j = 0; j < dim; j++) {
85 86
	  result = fscanf(file, "%d", &verticesEl1[j]);
	  TEST_EXIT(result == 1)("vertEl1[%d]\n", j);
87
	}
88 89
	result = fscanf(file, "%d", &el2);
	TEST_EXIT(result == 1)("el2?\n");
90
	for (int j = 0; j < dim; j++) {
91 92
	  result = fscanf(file, "%d", &verticesEl2[j]);
	  TEST_EXIT(result == 1)("vertEl2[%d]\n", j);
93
	}
94
	for (int j = 0; j < dim; j++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
95
	  if (mode == 0)
96
	    periodicMap.setEntry(melVertex[el1][verticesEl1[j]], 
Thomas Witkowski's avatar
Thomas Witkowski committed
97
				 melVertex[el2][verticesEl2[j]]);	  
98 99 100 101 102 103
	  vertexMapEl1[verticesEl1[j]] = verticesEl2[j];
	  vertexMapEl2[verticesEl2[j]] = verticesEl1[j];
	}

	// calculate sides of periodic vertices
	int sideEl1 = 0, sideEl2 = 0;
104
	if (dim == 1) {
105 106 107
	  sideEl1 = verticesEl1[0];
	  sideEl2 = verticesEl2[0];
	} else {
108
	  for (int j = 0; j < dim + 1; j++) {
109 110 111
	    sideEl1 += j;
	    sideEl2 += j;
	  }
112
	  for (int j = 0; j < dim; j++) {
113 114 115 116
	    sideEl1 -= verticesEl1[j];
	    sideEl2 -= verticesEl2[j];
	  }
	}
117

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

122 123
	Element *element1 = const_cast<Element*>((*(mel + el1))->getElement());
	Element *element2 = const_cast<Element*>((*(mel + el2))->getElement());
124
     
125
	// for all vertices of this side
126
	for (int j = 0; j < dim; j++) {
127 128
	  periodicCoordsEl1[element1->getPositionOfVertex(sideEl1, verticesEl1[j])] = 
	    coords[melVertex[el2][vertexMapEl1[verticesEl1[j]]]];
129

130 131 132
	  periodicCoordsEl2[element2->getPositionOfVertex(sideEl2, verticesEl2[j])] =
	    coords[melVertex[el1][vertexMapEl2[verticesEl2[j]]]];
	}
133
	     
134 135 136 137
	// decorate leaf data
	ElementData *ld1 = element1->getElementData();
	ElementData *ld2 = element2->getElementData();

138 139 140 141 142 143 144 145
	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());

146 147 148 149
	LeafDataPeriodic *ldp1 = 
	  dynamic_cast<LeafDataPeriodic*>(ld1->getElementData(PERIODIC));
	LeafDataPeriodic *ldp2 = 
	  dynamic_cast<LeafDataPeriodic*>(ld2->getElementData(PERIODIC));
150 151

	if (!ldp1) {
Thomas Witkowski's avatar
Thomas Witkowski committed
152
	  ldp1 = new LeafDataPeriodic(ld1);
153 154 155 156
	  element1->setElementData(ldp1);
	}

	if (!ldp2) {
Thomas Witkowski's avatar
Thomas Witkowski committed
157
	  ldp2 = new LeafDataPeriodic(ld2);
158 159 160
	  element2->setElementData(ldp2);
	}

161 162
	ldp1->addPeriodicInfo(mode, boundaryType, sideEl1, &periodicCoordsEl1);
	ldp2->addPeriodicInfo(mode, boundaryType, sideEl2, &periodicCoordsEl2);
163 164 165

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

167
	  if (!associated) {
Thomas Witkowski's avatar
Thomas Witkowski committed
168
	    associated = new VertexVector(mesh->getVertexAdmin(), "vertex vector");
169 170
	    mesh->periodicAssociations[boundaryType] = associated;
	    VertexVector::Iterator it(associated, ALL_DOFS);
171
	    for (it.reset2(); !it.end(); ++it)
172 173 174
	      *it = it.getDOFIndex();
	  }

175
	  for (int j = 0; j < dim; j++) {
176 177 178 179 180 181 182 183
	    (*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
184 185
      delete [] verticesEl1;
      delete [] verticesEl2;
186 187

      // change periodic vertex dofs
188
      for (int i = 0; i < mesh->getNumberOfVertices(); i++) {
189
	if (periodicMap.getEntry(i) != -1) {
190
	  mesh->freeDof(dof[i], VERTEX);
191 192
	  dof[i] = dof[periodicMap.getEntry(i)];

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

197
	  for (assoc = mesh->periodicAssociations.begin(); assoc != assocEnd; ++assoc) {
198 199 200 201 202

	    DegreeOfFreedom a = (*(assoc->second))[i];
	    if (a != i) {
	      (*(assoc->second))[i] = i;
	      (*(assoc->second))[a] = periodicMap.getEntry(i);
203
	    }
204 205
	  }

206 207 208
	}
      }

Thomas Witkowski's avatar
Thomas Witkowski committed
209
#if (DEBUG != 0)
210 211
      std::map<BoundaryType, VertexVector*>::iterator assoc;
      std::map<BoundaryType, VertexVector*>::iterator assocEnd =
212
	mesh->periodicAssociations.end();
Thomas Witkowski's avatar
Thomas Witkowski committed
213 214
      for (assoc = mesh->periodicAssociations.begin(); 
	   assoc != assocEnd; ++assoc)
Thomas Witkowski's avatar
Thomas Witkowski committed
215
	for (int i = 0; i < mesh->getNumberOfVertices(); i++)
216 217
	  if (i != (*(assoc->second))[i])
	    MSG("association %d: vertex %d -> vertex %d\n", 
Thomas Witkowski's avatar
Thomas Witkowski committed
218
		assoc->first, i, (*(assoc->second))[i]);
219

220 221 222 223
      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
224
#endif
225 226 227

    } // periodicFile

228

229 230
    // =========================================================

231 232 233
    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]]);
234

235
	const_cast<Element*>((*(mel + i))->getElement())->
236
	  setDof(k, dof[melVertex[i][k]]);
237
      }
238
    }
239

240
    if (!macroInfo->neigh_set) {
241
      TEST_EXIT(periodicFile == "")
242
	("periodic boundary condition => element neighbours must be set\n");
243
      computeNeighbours(mesh);
244
    } else {
245 246 247
      /****************************************************************************/
      /* fill MEL oppVertex values when reading neighbour information form file  */
      /****************************************************************************/
248

249 250 251 252 253
      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) {
254 255
	    int j = 0;
	    for (; j < mesh->getGeo(NEIGH); j++)
256 257
	      if (neigh->getNeighbour(j) == *(mel + i))  
		break;
258
	
259 260 261 262 263
	    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);
264
	  }
265
	}
266
      }
267
    }
268

269 270
    if (!macroInfo->bound_set)
      macroInfo->dirichletBoundary();    
271
  
272
    if (mesh->getDim() > 1)
273 274 275 276 277 278
      boundaryDOFs(mesh);

    // initial boundary projections
    int numFaces = mesh->getGeo(FACE);
    int dim = mesh->getDim();
    mel = mesh->firstMacroElement();
279
    for (int i = 0; i < mesh->getNumberOfLeaves(); i++) {
280
      MacroElement *macroEl = *(mel + i);
281
      Projection *projector = macroEl->getProjection(0);
282
      if (projector && projector->getType() == VOLUME_PROJECTION) {
Thomas Witkowski's avatar
Thomas Witkowski committed
283 284
	for (int j = 0; j <= dim; j++)
	  projector->project(macroEl->getCoord(j));	
285
      } else {
286
	for (int j = 0; j < mesh->getGeo(EDGE); j++) {
287
	  projector = macroEl->getProjection(numFaces + j);
288
	  if (projector) {
289 290 291 292 293 294 295 296 297 298 299
	    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);

300
    if (mesh->getNumberOfDofs(CENTER)) {
301
      for (int i = 0; i < mesh->getNumberOfMacros(); i++)
302
	const_cast<Element*>(mel[i]->getElement())->
303
	  setDof(mesh->getNode(CENTER), mesh->getDof(CENTER));
304 305
    }

306
    // === Domain size ===
307

308
    WorldVector<double> x_min, x_max;
309

310 311 312 313 314 315 316
    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++) {
317 318
	x_min[j] = std::min(x_min[j], coords[i][j]);
	x_max[j] = std::max(x_max[j], coords[i][j]);
319 320
      }
    }
321

322
    for (int j = 0; j < Global::getGeo(WORLD); j++)
323 324 325 326 327 328
      mesh->setDiameter(j, x_max[j] - x_min[j]);

    if (check) {
      checkMesh(mesh);

      if (mesh->getDim() > 1) {
329
	char filenew[128];
330
	strncpy(filenew, filename.c_str(), 128); 
331 332 333
	filenew[127] = 0;
	strncat(filenew, ".new", 128);   
	filenew[127] = 0;
334 335 336 337
	macroTest(mesh, filenew);
      }
    }

338
    return macroInfo;
339 340 341 342 343
  }


  void MacroReader::computeNeighbours(Mesh *mesh)
  {
344
    FUNCNAME("MacroReader::computeNeighbours()");
345

346
    int dim = mesh->getDim();
347
    FixVec<DegreeOfFreedom*, DIMEN> dof(dim, NO_INIT);
348

349 350 351 352
    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);
353
      }
354
    }
355

356 357 358 359 360 361 362
    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;
	}
363

364 365 366
	if (mesh->getMacroElement(i)->getOppVertex(k) == AMDIS_UNDEFINED) {
	  if (dim == 1) {
	    dof[0] = const_cast<DegreeOfFreedom*>(mesh->getMacroElement(i)->
367
						  getElement()->getDof(k));
368 369 370 371
	  } else {
	    for (int l = 0; l < dim; l++)
	      dof[l] = const_cast<DegreeOfFreedom*>(mesh->getMacroElement(i)->
						    getElement()->
372
						    getDof((k + l + 1) % (dim + 1)));
373
	  }
374 375 376 377 378 379 380 381 382 383 384 385 386
	  
	  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;
	    }
	  }

387 388 389 390 391 392 393
	  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);
394 395
	      if (l < dim)
		std::cout << " / ";	      
396 397 398
	    }
	    std::cout << std::endl;
	    std::cout << "  dofs = ";
399 400
	    for (int l = 0; l < dim; l++)
	      std::cout << *(dof[l]) << " ";	    
401 402 403 404
	    std::cout << std::endl;

	    ERROR_EXIT("\n");
	  }    
405
	}
406
      }
407
    }
408 409 410 411 412 413 414 415 416 417 418
  }


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

  void MacroReader::boundaryDOFs(Mesh *mesh)
  {
419 420 421 422 423
    FUNCNAME("Mesh::boundaryDOFs()");

    int lnode = mesh->getNode(EDGE);
    int k, lne = mesh->getNumberOfLeaves();
    int max_n_neigh = 0, n_neigh, ov;
424
    std::deque<MacroElement*>::iterator mel;
425
    const MacroElement* neigh;
426
    DegreeOfFreedom *dof;
427 428 429 430 431 432

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

    int dim = mesh->getDim();

Thomas Witkowski's avatar
Thomas Witkowski committed
433
    switch (dim) {
434
    case 2:
435
      for (mel = mesh->firstMacroElement(); mel != mesh->endOfMacroElements(); mel++) {
436 437 438 439 440 441
	// check for periodic boundary
	Element *el = const_cast<Element*>((*mel)->getElement());
	ElementData *ed = el->getElementData(PERIODIC);

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

442
	if (ed) {
443
	  std::list<LeafDataPeriodic::PeriodicInfo> &periodicInfos = 
444
	    dynamic_cast<LeafDataPeriodic*>(ed)->getInfoList();
445 446
	  std::list<LeafDataPeriodic::PeriodicInfo>::iterator it;
	  std::list<LeafDataPeriodic::PeriodicInfo>::iterator end = periodicInfos.end();
Thomas Witkowski's avatar
Thomas Witkowski committed
447 448
	  for (it = periodicInfos.begin(); it != end; ++it)
	    if (it->type != 0)
449 450 451
	      periodic[it->elementSide] = true;
	}

452
	for (int i = 0; i < mesh->getGeo(NEIGH); i++) {
453
	  if (!(*mel)->getNeighbour(i) || 
454 455 456 457
	      ((*mel)->getNeighbour(i)->getIndex() < (*mel)->getIndex())) {

	    mesh->incrementNumberOfEdges(1);

458
	    if (mesh->getNumberOfDofs(EDGE)) {
459 460
	      dof = mesh->getDof(EDGE);
	      el->setDof(lnode + i, dof);
461 462
      
	      if ((*mel)->getNeighbour(i)) {
Thomas Witkowski's avatar
Thomas Witkowski committed
463 464
		Element *neigh = 
		  const_cast<Element*>((*mel)->getNeighbour(i)->getElement());
465

Thomas Witkowski's avatar
Thomas Witkowski committed
466
		if (periodic[i])
467
		  neigh->setDof(lnode + (*mel)->getOppVertex(i), mesh->getDof(EDGE));
Thomas Witkowski's avatar
Thomas Witkowski committed
468
		else
469
		  neigh->setDof(lnode + (*mel)->getOppVertex(i), dof);		
470
	      }
471 472
	    }
	  }  
473 474 475 476 477 478
	}
      }
      break;
    case 3:
      lnode = mesh->getNode(FACE);
      mel = mesh->firstMacroElement();
479
      for (int i = 0; i < lne; i++) {
480 481 482 483 484 485 486

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

	DimVec<bool> periodic(dim, DEFAULT_VALUE, false);
      
487
	if (ed) {
488
	  std::list<LeafDataPeriodic::PeriodicInfo> &periodicInfos = 
489
	    dynamic_cast<LeafDataPeriodic*>(ed)->getInfoList();
490 491
	  std::list<LeafDataPeriodic::PeriodicInfo>::iterator it;
	  std::list<LeafDataPeriodic::PeriodicInfo>::iterator end = periodicInfos.end();
492 493
	  for (it = periodicInfos.begin(); it != end; ++it)
	    if (it->type != 0)
494 495 496
	      periodic[it->elementSide] = true;
	}

497 498
	for (k = 0; k < mesh->getGeo(EDGE); k++) {
	  // === Check for not counted edges. ===
499 500
	  n_neigh = 1;

Thomas Witkowski's avatar
Thomas Witkowski committed
501
	  if (newEdge(mesh, (*(mel + i)), k, &n_neigh)) {
502
	    mesh->incrementNumberOfEdges(1);
503
	    max_n_neigh = std::max(max_n_neigh, n_neigh);
504 505 506 507
	  }
	}
      
	for (k = 0; k < mesh->getGeo(NEIGH); k++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
508
	  neigh = (*(mel + i))->getNeighbour(k);
509
	  // === Face is counted and dof is added by the element with bigger index. ===
Thomas Witkowski's avatar
Thomas Witkowski committed
510 511
	  if (neigh && (neigh->getIndex() > (*(mel + i))->getIndex()))  
	    continue;
512 513 514
	
	  mesh->incrementNumberOfFaces(1);
	
515
	  if (mesh->getNumberOfDofs(FACE)) {
516
	    TEST_EXIT(!(*(mel + i))->getElement()->getDof(lnode + k))
517
	      ("dof %d on element %d already set\n", 
Thomas Witkowski's avatar
Thomas Witkowski committed
518
	       lnode + k, (*(mel + i))->getIndex());
519
	  
520
	    const_cast<Element*>((*(mel + i))->getElement())->setDof(lnode + k,
521
								     mesh->getDof(FACE));
522 523

	    if (neigh) {
Thomas Witkowski's avatar
Thomas Witkowski committed
524
	      ov = (*(mel + i))->getOppVertex(k);
525
	      TEST_EXIT(!neigh->getElement()->getDof(lnode + ov))
526
		("dof %d on neighbour %d already set\n", 
Thomas Witkowski's avatar
Thomas Witkowski committed
527
		 lnode + ov, neigh->getIndex());
528 529
	    
	      Element *neighEl = 
Thomas Witkowski's avatar
Thomas Witkowski committed
530
		const_cast<Element*>((*(mel + i))->getNeighbour(k)->getElement());
531

Thomas Witkowski's avatar
Thomas Witkowski committed
532
	      if (periodic[k])
533
		neighEl->setDof(lnode+ov, mesh->getDof(FACE));
Thomas Witkowski's avatar
Thomas Witkowski committed
534
	      else
535 536
		neighEl->setDof(lnode+ov, const_cast<int*>((*(mel + i))->getElement()->
							   getDof(lnode + k)));	      
537 538 539 540 541
	    }
	  }
	}
      }
      break;
Thomas Witkowski's avatar
Thomas Witkowski committed
542 543
    default: 
      ERROR_EXIT("invalid dim\n");
544 545
    }
    
546
    if (3 == dim)
Thomas Witkowski's avatar
Thomas Witkowski committed
547
      mesh->setMaxEdgeNeigh(std::max(8, 2 * max_n_neigh));
548
    else
Thomas Witkowski's avatar
Thomas Witkowski committed
549
      mesh->setMaxEdgeNeigh(dim - 1);    
550 551 552 553 554 555 556 557 558 559 560 561 562
  }

  /* 
     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)
  {
563
    FUNCNAME("MacroReader::macroTest()");
564
   
565 566 567 568 569 570 571
    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
572
      if (nameneu)
573 574
	ERROR_EXIT("mesh->feSpace\n");
    }
575 576 577 578 579 580 581 582 583 584 585 586 587 588
  }
  
  /****************************************************************************/
  /*  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)
  {
589
    FUNCNAME("MacroReader::macrotest()");
590

Thomas Witkowski's avatar
Thomas Witkowski committed
591 592
    std::deque<MacroElement*>::const_iterator macro, mac;
    int flg = 0;
593
    int dim = mesh->getDim();
Thomas Witkowski's avatar
Thomas Witkowski committed
594 595
    int *test = new int[mesh->getNumberOfMacros()];
    int *zykl = new int[mesh->getNumberOfMacros()];
596
 
Thomas Witkowski's avatar
Thomas Witkowski committed
597
    for (int i = 0; i < mesh->getNumberOfMacros(); i++)
598
      test[i] = 0;
599

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
673
    int *test = new int[mesh->getNumberOfMacros()];
674
  
Thomas Witkowski's avatar
Thomas Witkowski committed
675 676
    for (int i = 0; i < static_cast<int>(mesh->getNumberOfMacros()); i++)
      test[i] = 0;
677

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

Thomas Witkowski's avatar
Thomas Witkowski committed
680
    delete [] test;
681 682 683 684 685
  }

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

Thomas Witkowski's avatar
Thomas Witkowski committed
701
    int mel_index = mel->getIndex();
702 703 704 705

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

706
    if (mesh->getNumberOfDofs(EDGE)) {
707
      node = mesh->getNode(EDGE);
708
      if (el->getDof(node+edge_no)) {
709 710 711 712 713 714
	/****************************************************************************/
	/*  edge was counted by another macro element and dof was added on the      */
	/*  complete patch                                                          */
	/****************************************************************************/
	return false;
      } else {
715 716
	edge_dof = mesh->getDof(EDGE);
	el->setDof(node+edge_no, edge_dof);
717 718 719
      }
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
720
    for (j = 0; j < 2; j++)
721
      dof[j] = el->getDof(el->getVertexOfEdge(edge_no, j));
Thomas Witkowski's avatar
Thomas Witkowski committed
722
    
723 724 725 726 727 728 729 730 731 732 733 734
    /****************************************************************************/
    /*  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]);


735
    if (mel->getBoundary(next_el[edge_no][0])) {
736 737 738 739 740 741
      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++)
742
	if (nei->getElement()->getDof(j) == dof[0])  break;
743
      for (k = 0; k < vertices; k++)
744
	if (nei->getElement()->getDof(k) == dof[1])  break;
745 746

      // check for periodic boundary
747
      if (j == 4 || k == 4) {
748 749 750 751
	nei = NULL;
	break;
      }

752
      if (mesh->getNumberOfDofs(EDGE)) {
753 754
	TEST_EXIT(nei->index > mel_index)
	  ("neighbour index %d < element index %d\n", nei->getIndex(), mel_index);
755
      }
756

757
      if (!mesh->getNumberOfDofs(EDGE) && nei->getIndex() < mel_index) 
758
	return false;
759
    
760
      edge_no = Tetrahedron::edgeOfDofs[j][k];
761 762 763 764 765 766 767 768 769 770

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

771
      if (mesh->getNumberOfDofs(EDGE))
772
	nei->element->setDof(node+edge_no,edge_dof);      
773

774
      if (next_el[edge_no][0] != opp_v) {
775
	if (nei->getBoundary(next_el[edge_no][0])) {
776 777 778 779 780 781
	  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()))
782 783 784
	      lproject = neiProject;
	  }
	}
785 786 787 788 789 790 791 792 793 794
	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()))
795 796 797
	      lproject = neiProject;
	  }
	}
798 799 800
	opp_v = nei->getOppVertex(next_el[edge_no][1]);
	nei = nei->getNeighbour(next_el[edge_no][1]);
      }
801 802
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
803 804 805 806 807 808 809 810 811 812 813 814 815 816 817 818 819 820 821 822 823 824
    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<