mesh.c

Go to the documentation of this file.

/*
 * Copyright 1997, Regents of the University of Minnesota
 *
 * mesh.c
 *
 * This file contains routines for converting 3D and 4D finite element
 * meshes into dual or nodal graphs
 *
 * Started 8/18/97
 * George
 *
 * $Id: mesh.c 13804 2013-03-04 23:49:08Z karypis $
 *
 */

#include "metislib.h"


/*****************************************************************************/
/*****************************************************************************/
int METIS_MeshToDual(idx_t *ne, idx_t *nn, idx_t *eptr, idx_t *eind, 
          idx_t *ncommon, idx_t *numflag,  idx_t **r_xadj, idx_t **r_adjncy)
{
  int sigrval=0, renumber=0;

  /* set up malloc cleaning code and signal catchers */
  if (!gk_malloc_init()) 
    return METIS_ERROR_MEMORY;

  gk_sigtrap();

  if ((sigrval = gk_sigcatch()) != 0) 
    goto SIGTHROW;


  /* renumber the mesh */
  if (*numflag == 1) {
    ChangeMesh2CNumbering(*ne, eptr, eind);
    renumber = 1;
  }

  /* create dual graph */
  *r_xadj = *r_adjncy = NULL;
  CreateGraphDual(*ne, *nn, eptr, eind, *ncommon, r_xadj, r_adjncy);


SIGTHROW:
  if (renumber)
    ChangeMesh2FNumbering(*ne, eptr, eind, *ne, *r_xadj, *r_adjncy);

  gk_siguntrap();
  gk_malloc_cleanup(0);

  if (sigrval != 0) {
    if (*r_xadj != NULL)
      free(*r_xadj);
    if (*r_adjncy != NULL)
      free(*r_adjncy);
    *r_xadj = *r_adjncy = NULL;
  }

  return metis_rcode(sigrval);
}


/*****************************************************************************/
/*****************************************************************************/
int METIS_MeshToNodal(idx_t *ne, idx_t *nn, idx_t *eptr, idx_t *eind, 
          idx_t *numflag,  idx_t **r_xadj, idx_t **r_adjncy)
{
  int sigrval=0, renumber=0;

  /* set up malloc cleaning code and signal catchers */
  if (!gk_malloc_init()) 
    return METIS_ERROR_MEMORY;

  gk_sigtrap();

  if ((sigrval = gk_sigcatch()) != 0) 
    goto SIGTHROW;


  /* renumber the mesh */
  if (*numflag == 1) {
    ChangeMesh2CNumbering(*ne, eptr, eind);
    renumber = 1;
  }

  /* create nodal graph */
  *r_xadj = *r_adjncy = NULL;
  CreateGraphNodal(*ne, *nn, eptr, eind, r_xadj, r_adjncy);


SIGTHROW:
  if (renumber)
    ChangeMesh2FNumbering(*ne, eptr, eind, *nn, *r_xadj, *r_adjncy);

  gk_siguntrap();
  gk_malloc_cleanup(0);

  if (sigrval != 0) {
    if (*r_xadj != NULL)
      free(*r_xadj);
    if (*r_adjncy != NULL)
      free(*r_adjncy);
    *r_xadj = *r_adjncy = NULL;
  }

  return metis_rcode(sigrval);
}


/*****************************************************************************/
/*****************************************************************************/
void CreateGraphDual(idx_t ne, idx_t nn, idx_t *eptr, idx_t *eind, idx_t ncommon, 
          idx_t **r_xadj, idx_t **r_adjncy)
{
  idx_t i, j, nnbrs;
  idx_t *nptr, *nind;
  idx_t *xadj, *adjncy;
  idx_t *marker, *nbrs;

  if (ncommon < 1) {
    printf("  Increased ncommon to 1, as it was initially %"PRIDX"\n", ncommon);
    ncommon = 1;
  }

  /* construct the node-element list first */
  nptr = ismalloc(nn+1, 0, "CreateGraphDual: nptr");
  nind = imalloc(eptr[ne], "CreateGraphDual: nind");

  for (i=0; i<ne; i++) {
    for (j=eptr[i]; j<eptr[i+1]; j++)
      nptr[eind[j]]++;
  }
  MAKECSR(i, nn, nptr);

  for (i=0; i<ne; i++) {
    for (j=eptr[i]; j<eptr[i+1]; j++)
      nind[nptr[eind[j]]++] = i;
  }
  SHIFTCSR(i, nn, nptr);


  /* Allocate memory for xadj, since you know its size.
     These are done using standard malloc as they are returned
     to the calling function */
  if ((xadj = (idx_t *)malloc((ne+1)*sizeof(idx_t))) == NULL) 
    gk_errexit(SIGMEM, "***Failed to allocate memory for xadj.\n");
  *r_xadj = xadj;
  iset(ne+1, 0, xadj);

  /* allocate memory for working arrays used by FindCommonElements */
  marker = ismalloc(ne, 0, "CreateGraphDual: marker");
  nbrs   = imalloc(ne, "CreateGraphDual: nbrs");

  for (i=0; i<ne; i++) {
    xadj[i] = FindCommonElements(i, eptr[i+1]-eptr[i], eind+eptr[i], nptr, 
                  nind, eptr, ncommon, marker, nbrs);
  }
  MAKECSR(i, ne, xadj);

  /* Allocate memory for adjncy, since you now know its size.
     These are done using standard malloc as they are returned
     to the calling function */
  if ((adjncy = (idx_t *)malloc(xadj[ne]*sizeof(idx_t))) == NULL) {
    free(xadj);
    *r_xadj = NULL;
    gk_errexit(SIGMEM, "***Failed to allocate memory for adjncy.\n");
  }
  *r_adjncy = adjncy;

  for (i=0; i<ne; i++) {
    nnbrs = FindCommonElements(i, eptr[i+1]-eptr[i], eind+eptr[i], nptr, 
                nind, eptr, ncommon, marker, nbrs);
    for (j=0; j<nnbrs; j++)
      adjncy[xadj[i]++] = nbrs[j];
  }
  SHIFTCSR(i, ne, xadj);
  
  gk_free((void **)&nptr, &nind, &marker, &nbrs, LTERM);
}


/*****************************************************************************/
/*****************************************************************************/
idx_t FindCommonElements(idx_t qid, idx_t elen, idx_t *eind, idx_t *nptr, 
          idx_t *nind, idx_t *eptr, idx_t ncommon, idx_t *marker, idx_t *nbrs)
{
  idx_t i, ii, j, jj, k, l, overlap;

  /* find all elements that share at least one node with qid */
  for (k=0, i=0; i<elen; i++) {
    j = eind[i];
    for (ii=nptr[j]; ii<nptr[j+1]; ii++) {
      jj = nind[ii];

      if (marker[jj] == 0) 
        nbrs[k++] = jj;
      marker[jj]++;
    }
  }

  /* put qid into the neighbor list (in case it is not there) so that it
     will be removed in the next step */
  if (marker[qid] == 0)
    nbrs[k++] = qid;
  marker[qid] = 0;

  /* compact the list to contain only those with at least ncommon nodes */
  for (j=0, i=0; i<k; i++) {
    overlap = marker[l = nbrs[i]];
    if (overlap >= ncommon || 
        overlap >= elen-1 || 
        overlap >= eptr[l+1]-eptr[l]-1)
      nbrs[j++] = l;
    marker[l] = 0;
  }

  return j;
}


/*****************************************************************************/
/*****************************************************************************/
void CreateGraphNodal(idx_t ne, idx_t nn, idx_t *eptr, idx_t *eind, 
          idx_t **r_xadj, idx_t **r_adjncy)
{
  idx_t i, j, nnbrs;
  idx_t *nptr, *nind;
  idx_t *xadj, *adjncy;
  idx_t *marker, *nbrs;


  /* construct the node-element list first */
  nptr = ismalloc(nn+1, 0, "CreateGraphNodal: nptr");
  nind = imalloc(eptr[ne], "CreateGraphNodal: nind");

  for (i=0; i<ne; i++) {
    for (j=eptr[i]; j<eptr[i+1]; j++)
      nptr[eind[j]]++;
  }
  MAKECSR(i, nn, nptr);

  for (i=0; i<ne; i++) {
    for (j=eptr[i]; j<eptr[i+1]; j++)
      nind[nptr[eind[j]]++] = i;
  }
  SHIFTCSR(i, nn, nptr);


  /* Allocate memory for xadj, since you know its size.
     These are done using standard malloc as they are returned
     to the calling function */
  if ((xadj = (idx_t *)malloc((nn+1)*sizeof(idx_t))) == NULL)
    gk_errexit(SIGMEM, "***Failed to allocate memory for xadj.\n");
  *r_xadj = xadj;
  iset(nn+1, 0, xadj);

  /* allocate memory for working arrays used by FindCommonElements */
  marker = ismalloc(nn, 0, "CreateGraphNodal: marker");
  nbrs   = imalloc(nn, "CreateGraphNodal: nbrs");

  for (i=0; i<nn; i++) {
    xadj[i] = FindCommonNodes(i, nptr[i+1]-nptr[i], nind+nptr[i], eptr, 
                  eind, marker, nbrs);
  }
  MAKECSR(i, nn, xadj);

  /* Allocate memory for adjncy, since you now know its size.
     These are done using standard malloc as they are returned
     to the calling function */
  if ((adjncy = (idx_t *)malloc(xadj[nn]*sizeof(idx_t))) == NULL) {
    free(xadj);
    *r_xadj = NULL;
    gk_errexit(SIGMEM, "***Failed to allocate memory for adjncy.\n");
  }
  *r_adjncy = adjncy;

  for (i=0; i<nn; i++) {
    nnbrs = FindCommonNodes(i, nptr[i+1]-nptr[i], nind+nptr[i], eptr, 
                eind, marker, nbrs);
    for (j=0; j<nnbrs; j++)
      adjncy[xadj[i]++] = nbrs[j];
  }
  SHIFTCSR(i, nn, xadj);
  
  gk_free((void **)&nptr, &nind, &marker, &nbrs, LTERM);
}


/*****************************************************************************/
/*****************************************************************************/
idx_t FindCommonNodes(idx_t qid, idx_t nelmnts, idx_t *elmntids, idx_t *eptr, 
          idx_t *eind, idx_t *marker, idx_t *nbrs)
{
  idx_t i, ii, j, jj, k;

  /* find all nodes that share at least one element with qid */
  marker[qid] = 1;  /* this is to prevent self-loops */
  for (k=0, i=0; i<nelmnts; i++) {
    j = elmntids[i];
    for (ii=eptr[j]; ii<eptr[j+1]; ii++) {
      jj = eind[ii];
      if (marker[jj] == 0) {
        nbrs[k++] = jj;
        marker[jj] = 1;
      }
    }
  }

  /* reset the marker */
  marker[qid] = 0;
  for (i=0; i<k; i++) {
    marker[nbrs[i]] = 0;
  }

  return k;
}



/*************************************************************************/
/*************************************************************************/
mesh_t *CreateMesh(void)
{
  mesh_t *mesh;

  mesh = (mesh_t *)gk_malloc(sizeof(mesh_t), "CreateMesh: mesh");

  InitMesh(mesh);

  return mesh;
}


/*************************************************************************/
/*************************************************************************/
void InitMesh(mesh_t *mesh) 
{
  memset((void *)mesh, 0, sizeof(mesh_t));
}


/*************************************************************************/
/*************************************************************************/
void FreeMesh(mesh_t **r_mesh) 
{
  mesh_t *mesh = *r_mesh;
  
  gk_free((void **)&mesh->eptr, &mesh->eind, &mesh->ewgt, &mesh, LTERM);

  *r_mesh = NULL;
}