~elmer-csc-ubuntu/elmercsc/devel

« back to all changes in this revision

Viewing changes to elmergrid/src/metis-5.1.0/libmetis/mesh.c

  • Committer: GitHub
  • Author(s): Peter R
  • Date: 2017-10-04 13:47:14 UTC
  • mfrom: (575.1.17)
  • Revision ID: git-v1:91780af69137b74ec539fedc4ad70c9ec52d732a
Merge pull request #115 from ElmerCSC/metis_update

Updated to new Metis API for devel branch.

Show diffs side-by-side

added added

removed removed

Lines of Context:
 
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