Fixed error message for 2D mesh with periodic boundaries and higher order basis functions.

AMDiS/src/io/MacroReader.cc

@@ -1185,7 +1185,7 @@ namespace AMDiS {
   {
     FUNCNAME("MacroReader::checkMesh()");
     
-    int i, nused, nfree;
+    int nused, nfree;
     DOFAdmin  *localAdmin=mesh->admin[0];
     Flag fill_flag;
     int error_detected = 0;
@@ -1206,22 +1206,23 @@ namespace AMDiS {
     else 
       fill_flag = Mesh::CALL_LEAF_EL;
   
-    fill_flag |= Mesh::FILL_NEIGH;
+    fill_flag |= Mesh::FILL_NEIGH | Mesh::FILL_BOUND;
 
-    for (int iadmin = 0; iadmin < static_cast<int>(mesh->admin.size()); iadmin++) {
+    for (unsigned int iadmin = 0; iadmin < mesh->admin.size(); iadmin++) {
       localAdmin = mesh->admin[iadmin];
    
       if (localAdmin->getSize() > 0) {
 	if (static_cast<int>(mesh->dof_used.size()) < localAdmin->getSize()) 
 	  mesh->dof_used.resize(localAdmin->getSize() + 1000);
 	
-	for (i = 0; i < static_cast<int>(mesh->dof_used.size()); i++) 
+	for (unsigned int i = 0; i < mesh->dof_used.size(); i++) 
 	  mesh->dof_used[i] = 0;
 	
 	nused = nfree = 0;
 
 	TraverseStack stack;
-	ElInfo *elInfo = stack.traverseFirst(mesh, -1, fill_flag | Mesh::FILL_ADD_ALL);
+	ElInfo *elInfo = 
+	       stack.traverseFirst(mesh, -1, fill_flag | Mesh::FILL_ADD_ALL);
 	while (elInfo) {
 	  basicDOFCheckFct(elInfo, mesh, iadmin);
 	  elInfo = stack.traverseNext(elInfo);
@@ -1277,14 +1278,17 @@ namespace AMDiS {
 
     int j, k, opp_v;
     Element *el = elinfo->getElement();
-    const Element *neig;
     int error_detected = 0;
     int dim = mesh->getDim();
 
     elinfo->testFlag(Mesh::FILL_NEIGH);
 
     for (int i = 0; i < mesh->getGeo(NEIGH); i++) {
-      if ((neig = elinfo->getNeighbour(i))) {
+      const Element *neig = elinfo->getNeighbour(i);
+
+      if (neig) {
+	// If element has neighbour but also a non periodic boundary, this is
+	// an error.
 	if (elinfo->getBoundary(i) > 0) { // < 0 => periodic boundary
 	  if (!error_detected)
 	    MSG("error detected!!!\n");
@@ -1294,7 +1298,7 @@ namespace AMDiS {
 	}
 	
 	opp_v = elinfo->getOppVertex(i);
-	if (opp_v < 0  ||  opp_v >= mesh->getGeo(NEIGH)) {
+	if (opp_v < 0 || opp_v >= mesh->getGeo(NEIGH)) {
 	  if (!error_detected)
 	    MSG("error detected!!!\n");
 	  error_detected++;
@@ -1350,7 +1354,7 @@ namespace AMDiS {
     const Element *neig;
     const DegreeOfFreedom *dof;
 
-    if (0 == mesh->dof_used.size()) 
+    if (mesh->dof_used.size() == 0)
       return;
 
     int ndof = admin.getNumberOfDofs(VERTEX);
@@ -1361,13 +1365,13 @@ namespace AMDiS {
 	 j0, ndof, mesh->getNumberOfDofs(VERTEX));
       int i0 = mesh->getNode(VERTEX);
       for (int i = 0; i < mesh->getGeo(VERTEX); i++) {
-	if ((dof = el->getDof(i0+i)) == NULL) {
+	if ((dof = el->getDof(i0 + i)) == NULL) {
 	  ERROR("no vertex dof %d on element %d\n", i, el->getIndex());
 	} else {
 	  for (int j = 0; j < ndof; j++) {
 	    int jdof = dof[j0 + j];
 	    TEST(jdof >= 0 && jdof < static_cast<int>(mesh->dof_used.size()))
-	      ("vertex dof=%d invalid? size=%d\n",jdof, mesh->dof_used.size());
+	      ("vertex dof = %d invalid? size = %d\n", jdof, mesh->dof_used.size());
 	    mesh->dof_used[jdof]++;
 	  }
 	}
@@ -1378,7 +1382,10 @@ namespace AMDiS {
   
     if (mesh->getDim() > 1) {
       ndof = admin.getNumberOfDofs(EDGE);
+      
       if (ndof) {
+	// === Check for higher order DOFs on edges. ===
+
 	int j0 = admin.getNumberOfPreDofs(EDGE);
 
 	TEST_EXIT(j0 + ndof <= mesh->getNumberOfDofs(EDGE))
@@ -1404,7 +1411,12 @@ namespace AMDiS {
 	  if (el->getFirstChild() == NULL) {
 	    if (mesh->getDim() == 2) {
 	      neig = elinfo->getNeighbour(i);
-	      if (neig) {
+
+	      // Just check the edge if it is not part of a periodic boundary. In
+	      // the case of periodic boundaries the DOFs on neighbouring edges
+	      // cannot fit together!
+
+	      if (neig && elinfo->getBoundary(i) >= 0) {
 		int ov = elinfo->getOppVertex(i);
 
 		TEST(neig->getDof(i0 + ov) == dof)