src/libs/ck-libs/metis/libmetis/mesh.c

   1 /*
   2  * Copyright 1997, Regents of the University of Minnesota
   3  *
   4  * mesh.c
   5  *
   6  * This file contains routines for converting 3D and 4D finite element
   7  * meshes into dual or nodal graphs
   8  *
   9  * Started 8/18/97
  10  * George
  11  *
  12  * $Id: mesh.c 13804 2013-03-04 23:49:08Z karypis $
  13  *
  14  */
  15
  16 #include "metislib.h"
  17
  18
  19 /*****************************************************************************/
  20 /*! This function creates a graph corresponding to the dual of a finite element
  21     mesh.
  22
  23     \param ne is the number of elements in the mesh.
  24     \param nn is the number of nodes in the mesh.
  25     \param eptr is an array of size ne+1 used to mark the start and end
  26            locations in the nind array.
  27     \param eind is an array that stores for each element the set of node IDs
  28            (indices) that it is made off. The length of this array is equal
  29            to the total number of nodes over all the mesh elements.
  30     \param ncommon is the minimum number of nodes that two elements must share
  31            in order to be connected via an edge in the dual graph.
  32     \param numflag is either 0 or 1 indicating if the numbering of the nodes
  33            starts from 0 or 1, respectively. The same numbering is used for the
  34            returned graph as well.
  35     \param r_xadj indicates where the adjacency list of each vertex is stored
  36            in r_adjncy. The memory for this array is allocated by this routine.
  37            It can be freed by calling METIS_free().
  38     \param r_adjncy stores the adjacency list of each vertex in the generated
  39            dual graph. The memory for this array is allocated by this routine.
  40            It can be freed by calling METIS_free().
  41
  42 */
  43 /*****************************************************************************/
  44 int METIS_MeshToDual(idx_t *ne, idx_t *nn, idx_t *eptr, idx_t *eind,
  45           idx_t *ncommon, idx_t *numflag,  idx_t **r_xadj, idx_t **r_adjncy)
  46 {
  47   int sigrval=0, renumber=0;
  48
  49   /* set up malloc cleaning code and signal catchers */
  50   if (!gk_malloc_init())
  51     return METIS_ERROR_MEMORY;
  52
  53   gk_sigtrap();
  54
  55   if ((sigrval = gk_sigcatch()) != 0)
  56     goto SIGTHROW;
  57
  58
  59   /* renumber the mesh */
  60   if (*numflag == 1) {
  61     ChangeMesh2CNumbering(*ne, eptr, eind);
  62     renumber = 1;
  63   }
  64
  65   /* create dual graph */
  66   *r_xadj = *r_adjncy = NULL;
  67   CreateGraphDual(*ne, *nn, eptr, eind, *ncommon, r_xadj, r_adjncy);
  68
  69
  70 SIGTHROW:
  71   if (renumber)
  72     ChangeMesh2FNumbering(*ne, eptr, eind, *ne, *r_xadj, *r_adjncy);
  73
  74   gk_siguntrap();
  75   gk_malloc_cleanup(0);
  76
  77   if (sigrval != 0) {
  78     if (*r_xadj != NULL)
  79       free(*r_xadj);
  80     if (*r_adjncy != NULL)
  81       free(*r_adjncy);
  82     *r_xadj = *r_adjncy = NULL;
  83   }
  84
  85   return metis_rcode(sigrval);
  86 }
  87
  88
  89 /*****************************************************************************/
  90 /*! This function creates a graph corresponding to (almost) the nodal of a
  91     finite element mesh. In the nodal graph, each node is connected to the
  92     nodes corresponding to the union of nodes present in all the elements
  93     in which that node belongs.
  94
  95     \param ne is the number of elements in the mesh.
  96     \param nn is the number of nodes in the mesh.
  97     \param eptr is an array of size ne+1 used to mark the start and end
  98            locations in the nind array.
  99     \param eind is an array that stores for each element the set of node IDs
 100            (indices) that it is made off. The length of this array is equal
 101            to the total number of nodes over all the mesh elements.
 102     \param numflag is either 0 or 1 indicating if the numbering of the nodes
 103            starts from 0 or 1, respectively. The same numbering is used for the
 104            returned graph as well.
 105     \param r_xadj indicates where the adjacency list of each vertex is stored
 106            in r_adjncy. The memory for this array is allocated by this routine.
 107            It can be freed by calling METIS_free().
 108     \param r_adjncy stores the adjacency list of each vertex in the generated
 109            dual graph. The memory for this array is allocated by this routine.
 110            It can be freed by calling METIS_free().
 111
 112 */
 113 /*****************************************************************************/
 114 int METIS_MeshToNodal(idx_t *ne, idx_t *nn, idx_t *eptr, idx_t *eind,
 115           idx_t *numflag,  idx_t **r_xadj, idx_t **r_adjncy)
 116 {
 117   int sigrval=0, renumber=0;
 118
 119   /* set up malloc cleaning code and signal catchers */
 120   if (!gk_malloc_init())
 121     return METIS_ERROR_MEMORY;
 122
 123   gk_sigtrap();
 124
 125   if ((sigrval = gk_sigcatch()) != 0)
 126     goto SIGTHROW;
 127
 128
 129   /* renumber the mesh */
 130   if (*numflag == 1) {
 131     ChangeMesh2CNumbering(*ne, eptr, eind);
 132     renumber = 1;
 133   }
 134
 135   /* create nodal graph */
 136   *r_xadj = *r_adjncy = NULL;
 137   CreateGraphNodal(*ne, *nn, eptr, eind, r_xadj, r_adjncy);
 138
 139
 140 SIGTHROW:
 141   if (renumber)
 142     ChangeMesh2FNumbering(*ne, eptr, eind, *nn, *r_xadj, *r_adjncy);
 143
 144   gk_siguntrap();
 145   gk_malloc_cleanup(0);
 146
 147   if (sigrval != 0) {
 148     if (*r_xadj != NULL)
 149       free(*r_xadj);
 150     if (*r_adjncy != NULL)
 151       free(*r_adjncy);
 152     *r_xadj = *r_adjncy = NULL;
 153   }
 154
 155   return metis_rcode(sigrval);
 156 }
 157
 158
 159 /*****************************************************************************/
 160 /*! This function creates the dual of a finite element mesh */
 161 /*****************************************************************************/
 162 void CreateGraphDual(idx_t ne, idx_t nn, idx_t *eptr, idx_t *eind, idx_t ncommon,
 163           idx_t **r_xadj, idx_t **r_adjncy)
 164 {
 165   idx_t i, j, nnbrs;
 166   idx_t *nptr, *nind;
 167   idx_t *xadj, *adjncy;
 168   idx_t *marker, *nbrs;
 169
 170   if (ncommon < 1) {
 171     printf("  Increased ncommon to 1, as it was initially %"PRIDX"\n", ncommon);
 172     ncommon = 1;
 173   }
 174
 175   /* construct the node-element list first */
 176   nptr = ismalloc(nn+1, 0, "CreateGraphDual: nptr");
 177   nind = imalloc(eptr[ne], "CreateGraphDual: nind");
 178
 179   for (i=0; i<ne; i++) {
 180     for (j=eptr[i]; j<eptr[i+1]; j++)
 181       nptr[eind[j]]++;
 182   }
 183   MAKECSR(i, nn, nptr);
 184
 185   for (i=0; i<ne; i++) {
 186     for (j=eptr[i]; j<eptr[i+1]; j++)
 187       nind[nptr[eind[j]]++] = i;
 188   }
 189   SHIFTCSR(i, nn, nptr);
 190
 191
 192   /* Allocate memory for xadj, since you know its size.
 193      These are done using standard malloc as they are returned
 194      to the calling function */
 195   if ((xadj = (idx_t *)malloc((ne+1)*sizeof(idx_t))) == NULL)
 196     gk_errexit(SIGMEM, "***Failed to allocate memory for xadj.\n");
 197   *r_xadj = xadj;
 198   iset(ne+1, 0, xadj);
 199
 200   /* allocate memory for working arrays used by FindCommonElements */
 201   marker = ismalloc(ne, 0, "CreateGraphDual: marker");
 202   nbrs   = imalloc(ne, "CreateGraphDual: nbrs");
 203
 204   for (i=0; i<ne; i++) {
 205     xadj[i] = FindCommonElements(i, eptr[i+1]-eptr[i], eind+eptr[i], nptr,
 206                   nind, eptr, ncommon, marker, nbrs);
 207   }
 208   MAKECSR(i, ne, xadj);
 209
 210   /* Allocate memory for adjncy, since you now know its size.
 211      These are done using standard malloc as they are returned
 212      to the calling function */
 213   if ((adjncy = (idx_t *)malloc(xadj[ne]*sizeof(idx_t))) == NULL) {
 214     free(xadj);
 215     *r_xadj = NULL;
 216     gk_errexit(SIGMEM, "***Failed to allocate memory for adjncy.\n");
 217   }
 218   *r_adjncy = adjncy;
 219
 220   for (i=0; i<ne; i++) {
 221     nnbrs = FindCommonElements(i, eptr[i+1]-eptr[i], eind+eptr[i], nptr,
 222                 nind, eptr, ncommon, marker, nbrs);
 223     for (j=0; j<nnbrs; j++)
 224       adjncy[xadj[i]++] = nbrs[j];
 225   }
 226   SHIFTCSR(i, ne, xadj);
 227
 228   gk_free((void **)&nptr, &nind, &marker, &nbrs, LTERM);
 229 }
 230
 231
 232 /*****************************************************************************/
 233 /*! This function finds all elements that share at least ncommon nodes with
 234     the ``query'' element.
 235 */
 236 /*****************************************************************************/
 237 idx_t FindCommonElements(idx_t qid, idx_t elen, idx_t *eind, idx_t *nptr,
 238           idx_t *nind, idx_t *eptr, idx_t ncommon, idx_t *marker, idx_t *nbrs)
 239 {
 240   idx_t i, ii, j, jj, k, l, overlap;
 241
 242   /* find all elements that share at least one node with qid */
 243   for (k=0, i=0; i<elen; i++) {
 244     j = eind[i];
 245     for (ii=nptr[j]; ii<nptr[j+1]; ii++) {
 246       jj = nind[ii];
 247
 248       if (marker[jj] == 0)
 249         nbrs[k++] = jj;
 250       marker[jj]++;
 251     }
 252   }
 253
 254   /* put qid into the neighbor list (in case it is not there) so that it
 255      will be removed in the next step */
 256   if (marker[qid] == 0)
 257     nbrs[k++] = qid;
 258   marker[qid] = 0;
 259
 260   /* compact the list to contain only those with at least ncommon nodes */
 261   for (j=0, i=0; i<k; i++) {
 262     overlap = marker[l = nbrs[i]];
 263     if (overlap >= ncommon ||
 264         overlap >= elen-1 ||
 265         overlap >= eptr[l+1]-eptr[l]-1)
 266       nbrs[j++] = l;
 267     marker[l] = 0;
 268   }
 269
 270   return j;
 271 }
 272
 273
 274 /*****************************************************************************/
 275 /*! This function creates the (almost) nodal of a finite element mesh */
 276 /*****************************************************************************/
 277 void CreateGraphNodal(idx_t ne, idx_t nn, idx_t *eptr, idx_t *eind,
 278           idx_t **r_xadj, idx_t **r_adjncy)
 279 {
 280   idx_t i, j, nnbrs;
 281   idx_t *nptr, *nind;
 282   idx_t *xadj, *adjncy;
 283   idx_t *marker, *nbrs;
 284
 285
 286   /* construct the node-element list first */
 287   nptr = ismalloc(nn+1, 0, "CreateGraphNodal: nptr");
 288   nind = imalloc(eptr[ne], "CreateGraphNodal: nind");
 289
 290   for (i=0; i<ne; i++) {
 291     for (j=eptr[i]; j<eptr[i+1]; j++)
 292       nptr[eind[j]]++;
 293   }
 294   MAKECSR(i, nn, nptr);
 295
 296   for (i=0; i<ne; i++) {
 297     for (j=eptr[i]; j<eptr[i+1]; j++)
 298       nind[nptr[eind[j]]++] = i;
 299   }
 300   SHIFTCSR(i, nn, nptr);
 301
 302
 303   /* Allocate memory for xadj, since you know its size.
 304      These are done using standard malloc as they are returned
 305      to the calling function */
 306   if ((xadj = (idx_t *)malloc((nn+1)*sizeof(idx_t))) == NULL)
 307     gk_errexit(SIGMEM, "***Failed to allocate memory for xadj.\n");
 308   *r_xadj = xadj;
 309   iset(nn+1, 0, xadj);
 310
 311   /* allocate memory for working arrays used by FindCommonElements */
 312   marker = ismalloc(nn, 0, "CreateGraphNodal: marker");
 313   nbrs   = imalloc(nn, "CreateGraphNodal: nbrs");
 314
 315   for (i=0; i<nn; i++) {
 316     xadj[i] = FindCommonNodes(i, nptr[i+1]-nptr[i], nind+nptr[i], eptr,
 317                   eind, marker, nbrs);
 318   }
 319   MAKECSR(i, nn, xadj);
 320
 321   /* Allocate memory for adjncy, since you now know its size.
 322      These are done using standard malloc as they are returned
 323      to the calling function */
 324   if ((adjncy = (idx_t *)malloc(xadj[nn]*sizeof(idx_t))) == NULL) {
 325     free(xadj);
 326     *r_xadj = NULL;
 327     gk_errexit(SIGMEM, "***Failed to allocate memory for adjncy.\n");
 328   }
 329   *r_adjncy = adjncy;
 330
 331   for (i=0; i<nn; i++) {
 332     nnbrs = FindCommonNodes(i, nptr[i+1]-nptr[i], nind+nptr[i], eptr,
 333                 eind, marker, nbrs);
 334     for (j=0; j<nnbrs; j++)
 335       adjncy[xadj[i]++] = nbrs[j];
 336   }
 337   SHIFTCSR(i, nn, xadj);
 338
 339   gk_free((void **)&nptr, &nind, &marker, &nbrs, LTERM);
 340 }
 341
 342
 343 /*****************************************************************************/
 344 /*! This function finds the union of nodes that are in the same elements with
 345     the ``query'' node.
 346 */
 347 /*****************************************************************************/
 348 idx_t FindCommonNodes(idx_t qid, idx_t nelmnts, idx_t *elmntids, idx_t *eptr,
 349           idx_t *eind, idx_t *marker, idx_t *nbrs)
 350 {
 351   idx_t i, ii, j, jj, k;
 352
 353   /* find all nodes that share at least one element with qid */
 354   marker[qid] = 1;  /* this is to prevent self-loops */
 355   for (k=0, i=0; i<nelmnts; i++) {
 356     j = elmntids[i];
 357     for (ii=eptr[j]; ii<eptr[j+1]; ii++) {
 358       jj = eind[ii];
 359       if (marker[jj] == 0) {
 360         nbrs[k++] = jj;
 361         marker[jj] = 1;
 362       }
 363     }
 364   }
 365
 366   /* reset the marker */
 367   marker[qid] = 0;
 368   for (i=0; i<k; i++) {
 369     marker[nbrs[i]] = 0;
 370   }
 371
 372   return k;
 373 }
 374
 375
 376
 377 /*************************************************************************/
 378 /*! This function creates and initializes a mesh_t structure */
 379 /*************************************************************************/
 380 mesh_t *CreateMesh(void)
 381 {
 382   mesh_t *mesh;
 383
 384   mesh = (mesh_t *)gk_malloc(sizeof(mesh_t), "CreateMesh: mesh");
 385
 386   InitMesh(mesh);
 387
 388   return mesh;
 389 }
 390
 391
 392 /*************************************************************************/
 393 /*! This function initializes a mesh_t data structure */
 394 /*************************************************************************/
 395 void InitMesh(mesh_t *mesh)
 396 {
 397   memset((void *)mesh, 0, sizeof(mesh_t));
 398 }
 399
 400
 401 /*************************************************************************/
 402 /*! This function deallocates any memory stored in a mesh */
 403 /*************************************************************************/
 404 void FreeMesh(mesh_t **r_mesh)
 405 {
 406   mesh_t *mesh = *r_mesh;
 407
 408   gk_free((void **)&mesh->eptr, &mesh->eind, &mesh->ewgt, &mesh, LTERM);
 409
 410   *r_mesh = NULL;
 411 }
 412