Robust/src/Benchmarks/SingleTM/Yada/mesh.java

   1 /* =============================================================================
   2  *
   3  * mesh.c
   4  *
   5  * =============================================================================
   6  *
   7  * Copyright (C) Stanford University, 2006.  All Rights Reserved.
   8  * Author: Chi Cao Minh
   9  *
  10  * =============================================================================
  11  *
  12  * For the license of bayes/sort.h and bayes/sort.c, please see the header
  13  * of the files.
  14  *
  15  * ------------------------------------------------------------------------
  16  *
  17  * For the license of kmeans, please see kmeans/LICENSE.kmeans
  18  *
  19  * ------------------------------------------------------------------------
  20  *
  21  * For the license of ssca2, please see ssca2/COPYRIGHT
  22  *
  23  * ------------------------------------------------------------------------
  24  *
  25  * For the license of lib/mt19937ar.c and lib/mt19937ar.h, please see the
  26  * header of the files.
  27  *
  28  * ------------------------------------------------------------------------
  29  *
  30  * For the license of lib/rbtree.h and lib/rbtree.c, please see
  31  * lib/LEGALNOTICE.rbtree and lib/LICENSE.rbtree
  32  *
  33  * ------------------------------------------------------------------------
  34  *
  35  * Unless otherwise noted, the following license applies to STAMP files:
  36  *
  37  * Copyright (c) 2007, Stanford University
  38  * All rights reserved.
  39  *
  40  * Redistribution and use in source and binary forms, with or without
  41  * modification, are permitted provided that the following conditions are
  42  * met:
  43  *
  44  *     * Redistributions of source code must retain the above copyright
  45  *       notice, this list of conditions and the following disclaimer.
  46  *
  47  *     * Redistributions in binary form must reproduce the above copyright
  48  *       notice, this list of conditions and the following disclaimer in
  49  *       the documentation and/or other materials provided with the
  50  *       distribution.
  51  *
  52  *     * Neither the name of Stanford University nor the names of its
  53  *       contributors may be used to endorse or promote products derived
  54  *       from this software without specific prior written permission.
  55  *
  56  * THIS SOFTWARE IS PROVIDED BY STANFORD UNIVERSITY ``AS IS'' AND ANY
  57  * EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
  58  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
  59  * PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL STANFORD UNIVERSITY BE LIABLE
  60  * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
  61  * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
  62  * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
  63  * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
  64  * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
  65  * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
  66  * THE POSSIBILITY OF SUCH DAMAGE.
  67  *
  68  * =============================================================================
  69  */
  70
  71 public class mesh {
  72     element rootElementPtr;
  73     Queue_t initBadQueuePtr;
  74     int size;
  75     SET_T* boundarySetPtr;
  76
  77 /* =============================================================================
  78  * mesh_alloc
  79  * =============================================================================
  80  */
  81   public mesh() {
  82     rootElementPtr = null;
  83     initBadQueuePtr = new Queue_t(-1);
  84     size = 0;
  85     boundarySetPtr = SET_ALLOC(null, &element_listCompareEdge);
  86   }
  87
  88
  89   /* =============================================================================
  90    * TMmesh_insert
  91    * =============================================================================
  92    */
  93   void TMmesh_insert (element elementPtr, MAP_T edgeMapPtr) {
  94     /*
  95      * Assuming fully connected graph, we just need to record one element.
  96      * The root element is not needed for the actual refining, but rather
  97      * for checking the validity of the final mesh.
  98      */
  99     if (rootElementPtr==null) {
 100       rootElementPtr=elementPtr);
 101   }
 102
 103     /*
 104      * Record existence of each of this element's edges
 105      */
 106   int numEdge = elementPtr.element_getNumEdge();
 107   for (int i = 0; i < numEdge; i++) {
 108     edge edgePtr = elementPtr.element_getEdge(i);
 109     if (!MAP_CONTAINS(edgeMapPtr, (void*)edgePtr)) {
 110       /* Record existance of this edge */
 111       boolean isSuccess =
 112         PMAP_INSERT(edgeMapPtr, (void*)edgePtr, (void*)elementPtr);
 113       assert(isSuccess);
 114     } else {
 115       /*
 116        * Shared edge; update each element's neighborList
 117        */
 118       boolean isSuccess;
 119       element sharerPtr = (element_t*)MAP_FIND(edgeMapPtr, edgePtr);
 120       assert(sharerPtr); /* cannot be shared by >2 elements */
 121       elementPtr.element_addNeighbor(sharerPtr);
 122       sharerPtr.element_addNeighbor(elementPtr);
 123       isSuccess = PMAP_REMOVE(edgeMapPtr, edgePtr);
 124       assert(isSuccess);
 125       isSuccess = PMAP_INSERT(edgeMapPtr,
 126                               edgePtr,
 127                               NULL); /* marker to check >2 sharers */
 128       assert(isSuccess);
 129     }
 130   }
 131
 132   /*
 133    * Check if really encroached
 134    */
 135
 136   edge encroachedPtr = elementPtr.element_getEncroachedPtr();
 137   if (encroachedPtr!=null) {
 138     if (!TMSET_CONTAINS(meshPtr.boundarySetPtr, encroachedPtr)) {
 139       element_clearEncroached(elementPtr);
 140     }
 141   }
 142 }
 143
 144
 145 /* =============================================================================
 146  * TMmesh_remove
 147  * =============================================================================
 148  */
 149 public void TMmesh_remove(element elementPtr) {
 150   assert(!elementPtr.element_isGarbage());
 151
 152   /*
 153    * If removing root, a new root is selected on the next mesh_insert, which
 154    * always follows a call a mesh_remove.
 155    */
 156   if (rootElementPtr == elementPtr) {
 157     rootElementPtr=null;
 158   }
 159
 160   /*
 161    * Remove from neighbors
 162    */
 163   list_iter_t it;
 164   List neighborListPtr = elementPtr.element_getNeighborListPtr();
 165   TMLIST_ITER_RESET(&it, neighborListPtr);
 166   while (TMLIST_ITER_HASNEXT(&it, neighborListPtr)) {
 167       element neighborPtr = (element)TMLIST_ITER_NEXT(&it, neighborListPtr);
 168       List_t neighborNeighborListPtr = neighborPtr.element_getNeighborListPtr();
 169       boolean status = neighborNeighborListPtr.remove(elementPtr);
 170       assert(status);
 171   }
 172
 173   elementPtr.element_isGarbage(true);
 174 }
 175
 176
 177 /* =============================================================================
 178  * TMmesh_insertBoundary
 179  * =============================================================================
 180  */
 181 boolean TMmesh_insertBoundary (meshPtr, edge boundaryPtr) {
 182   return TMSET_INSERT(boundarySetPtr, boundaryPtr);
 183 }
 184
 185
 186 /* =============================================================================
 187  * TMmesh_removeBoundary
 188  * =============================================================================
 189  */
 190 boolean TMmesh_removeBoundary (meshPtr, edge boundaryPtr) {
 191   return TMSET_REMOVE(boundarySetPtr, boundaryPtr);
 192 }
 193
 194
 195 /* =============================================================================
 196  * createElement
 197  * =============================================================================
 198  */
 199 static void createElement (coordinate coordinates,
 200                int numCoordinate,
 201                            MAP_T edgeMapPtr) {
 202     element elementPtr = new element(coordinates, numCoordinate);
 203     assert(elementPtr);
 204
 205     if (numCoordinate == 2) {
 206         edge boundaryPtr = elementPtr.element_getEdge(0);
 207         boolean status = SET_INSERT(boundarySetPtr, boundaryPtr);
 208         assert(status);
 209     }
 210
 211     mesh_insert(elementPtr, edgeMapPtr);
 212
 213     if (elementPtr.element_isBad()) {
 214         boolean status = initBadQueuePtr.queue_push(elementPtr);
 215         assert(status);
 216     }
 217 }
 218
 219
 220 /* =============================================================================
 221  * mesh_read
 222  *
 223  * Returns number of elements read from file
 224  *
 225  * Refer to http://www.cs.cmu.edu/~quake/triangle.html for file formats.
 226  * =============================================================================
 227  */
 228 int mesh_read(String fileNamePrefix) {
 229     FILE* inputFile;
 230     coordinate_t* coordinates;
 231     char fileName[256];
 232     int fileNameSize = sizeof(fileName) / sizeof(fileName[0]);
 233     char inputBuff[256];
 234     int inputBuffSize = sizeof(inputBuff) / sizeof(inputBuff[0]);
 235     int numEntry;
 236     int numDimension;
 237     int numCoordinate;
 238     int i;
 239     int numElement = 0;
 240
 241     MAP_T* edgeMapPtr = MAP_ALLOC(NULL, &element_mapCompareEdge);
 242     assert(edgeMapPtr);
 243
 244     /*
 245      * Read .node file
 246      */
 247     snprintf(fileName, fileNameSize, "%s.node", fileNamePrefix);
 248     inputFile = fopen(fileName, "r");
 249     assert(inputFile);
 250     fgets(inputBuff, inputBuffSize, inputFile);
 251     sscanf(inputBuff, "%li %li", &numEntry, &numDimension);
 252     assert(numDimension == 2); /* must be 2-D */
 253     numCoordinate = numEntry + 1; /* numbering can start from 1 */
 254     coordinates = (coordinate_t*)malloc(numCoordinate * sizeof(coordinate_t));
 255     assert(coordinates);
 256     for (i = 0; i < numEntry; i++) {
 257         int id;
 258         double x;
 259         double y;
 260         if (!fgets(inputBuff, inputBuffSize, inputFile)) {
 261             break;
 262         }
 263         if (inputBuff[0] == '#') {
 264             continue; /* TODO: handle comments correctly */
 265         }
 266         sscanf(inputBuff, "%li %lf %lf", &id, &x, &y);
 267         coordinates[id].x = x;
 268         coordinates[id].y = y;
 269     }
 270     assert(i == numEntry);
 271     fclose(inputFile);
 272
 273     /*
 274      * Read .poly file, which contains boundary segments
 275      */
 276     snprintf(fileName, fileNameSize, "%s.poly", fileNamePrefix);
 277     inputFile = fopen(fileName, "r");
 278     assert(inputFile);
 279     fgets(inputBuff, inputBuffSize, inputFile);
 280     sscanf(inputBuff, "%li %li", &numEntry, &numDimension);
 281     assert(numEntry == 0); /* .node file used for vertices */
 282     assert(numDimension == 2); /* must be edge */
 283     fgets(inputBuff, inputBuffSize, inputFile);
 284     sscanf(inputBuff, "%li", &numEntry);
 285     for (i = 0; i < numEntry; i++) {
 286         int id;
 287         int a;
 288         int b;
 289         coordinate_t insertCoordinates[2];
 290         if (!fgets(inputBuff, inputBuffSize, inputFile)) {
 291             break;
 292         }
 293         if (inputBuff[0] == '#') {
 294             continue; /* TODO: handle comments correctly */
 295         }
 296         sscanf(inputBuff, "%li %li %li", &id, &a, &b);
 297         assert(a >= 0 && a < numCoordinate);
 298         assert(b >= 0 && b < numCoordinate);
 299         insertCoordinates[0] = coordinates[a];
 300         insertCoordinates[1] = coordinates[b];
 301         createElement(meshPtr, insertCoordinates, 2, edgeMapPtr);
 302     }
 303     assert(i == numEntry);
 304     numElement += numEntry;
 305     fclose(inputFile);
 306
 307     /*
 308      * Read .ele file, which contains triangles
 309      */
 310     snprintf(fileName, fileNameSize, "%s.ele", fileNamePrefix);
 311     inputFile = fopen(fileName, "r");
 312     assert(inputFile);
 313     fgets(inputBuff, inputBuffSize, inputFile);
 314     sscanf(inputBuff, "%li %li", &numEntry, &numDimension);
 315     assert(numDimension == 3); /* must be triangle */
 316     for (i = 0; i < numEntry; i++) {
 317         int id;
 318         int a;
 319         int b;
 320         int c;
 321         coordinate_t insertCoordinates[3];
 322         if (!fgets(inputBuff, inputBuffSize, inputFile)) {
 323             break;
 324         }
 325         if (inputBuff[0] == '#') {
 326             continue; /* TODO: handle comments correctly */
 327         }
 328         sscanf(inputBuff, "%li %li %li %li", &id, &a, &b, &c);
 329         assert(a >= 0 && a < numCoordinate);
 330         assert(b >= 0 && b < numCoordinate);
 331         assert(c >= 0 && c < numCoordinate);
 332         insertCoordinates[0] = coordinates[a];
 333         insertCoordinates[1] = coordinates[b];
 334         insertCoordinates[2] = coordinates[c];
 335         createElement(meshPtr, insertCoordinates, 3, edgeMapPtr);
 336     }
 337     assert(i == numEntry);
 338     numElement += numEntry;
 339     fclose(inputFile);
 340
 341     free(coordinates);
 342     MAP_FREE(edgeMapPtr);
 343
 344     return numElement;
 345 }
 346
 347
 348 /* =============================================================================
 349  * mesh_getBad
 350  * -- Returns NULL if none
 351  * =============================================================================
 352  */
 353 element mesh_getBad() {
 354   return (element)queue_pop(initBadQueuePtr);
 355 }
 356
 357
 358 /* =============================================================================
 359  * mesh_shuffleBad
 360  * =============================================================================
 361  */
 362 void mesh_shuffleBad (Random randomPtr) {
 363   queue_shuffle(initBadQueuePtr, randomPtr);
 364 }
 365
 366
 367 /* =============================================================================
 368  * mesh_check
 369  * =============================================================================
 370  */
 371 boolean mesh_check(int expectedNumElement) {
 372     int numBadTriangle = 0;
 373     int numFalseNeighbor = 0;
 374     int numElement = 0;
 375
 376     System.out.println("Checking final mesh:");
 377
 378     Queue_t searchQueuePtr = new Queue_t(-1);
 379     assert(searchQueuePtr);
 380     MAP_T visitedMapPtr = MAP_ALLOC(NULL, &element_mapCompare);
 381     assert(visitedMapPtr);
 382
 383     /*
 384      * Do breadth-first search starting from rootElementPtr
 385      */
 386     assert(rootElementPtr!=null);
 387     searchQueuePtr.queue_push(rootElementPtr);
 388     while (!searchQueuePtr.queue_isEmpty()) {
 389         list_iter_t it;
 390         List_t neighborListPtr;
 391
 392         element currentElementPtr = (element)queue_pop(searchQueuePtr);
 393         if (MAP_CONTAINS(visitedMapPtr, (void*)currentElementPtr)) {
 394             continue;
 395         }
 396         boolean isSuccess = MAP_INSERT(visitedMapPtr, (void*)currentElementPtr, NULL);
 397         assert(isSuccess);
 398         if (!currentElementPtr.checkAngles()) {
 399             numBadTriangle++;
 400         }
 401         neighborListPtr = currentElementPtr.element_getNeighborListPtr();
 402
 403         list_iter_reset(&it, neighborListPtr);
 404         while (list_iter_hasNext(&it, neighborListPtr)) {
 405             element neighborElementPtr =
 406                 (element)list_iter_next(&it, neighborListPtr);
 407             /*
 408              * Continue breadth-first search
 409              */
 410             if (!MAP_CONTAINS(visitedMapPtr, (void*)neighborElementPtr)) {
 411                 boolean isSuccess = searchQueuePtr.queue_push(neighborElementPtr);
 412                 assert(isSuccess);
 413             }
 414         } /* for each neighbor */
 415
 416         numElement++;
 417
 418     } /* breadth-first search */
 419
 420     System.out.println("Number of elements      = "+ numElement);
 421     System.out.println("Number of bad triangles = "+ numBadTriangle);
 422
 423     MAP_FREE(visitedMapPtr);
 424
 425     return (!(numBadTriangle > 0 ||
 426               numFalseNeighbor > 0 ||
 427               numElement != expectedNumElement));
 428 }