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,v 1.1 1998/11/27 17:59:20 karypis Exp $
13 *
14 */
15
16 #include <metis.h>
17
18 /*****************************************************************************
19 * This function creates a graph corresponding to the dual of a finite element
20 * mesh. At this point the supported elements are triangles, tetrahedrons, and
21 * bricks.
22 ******************************************************************************/
METIS_MeshToDual(int * ne,int * nn,idxtype * elmnts,int * etype,int * numflag,idxtype * dxadj,idxtype * dadjncy)23 void METIS_MeshToDual(int *ne, int *nn, idxtype *elmnts, int *etype, int *numflag,
24 idxtype *dxadj, idxtype *dadjncy)
25 {
26 int esizes[] = {-1, 3, 4, 8, 4};
27
28 if (*numflag == 1)
29 ChangeMesh2CNumbering((*ne)*esizes[*etype], elmnts);
30
31 GENDUALMETIS(*ne, *nn, *etype, elmnts, dxadj, dadjncy);
32
33 if (*numflag == 1)
34 ChangeMesh2FNumbering((*ne)*esizes[*etype], elmnts, *ne, dxadj, dadjncy);
35 }
36
37
38 /*****************************************************************************
39 * This function creates a graph corresponding to the finite element mesh.
40 * At this point the supported elements are triangles, tetrahedrons.
41 ******************************************************************************/
METIS_MeshToNodal(int * ne,int * nn,idxtype * elmnts,int * etype,int * numflag,idxtype * dxadj,idxtype * dadjncy)42 void METIS_MeshToNodal(int *ne, int *nn, idxtype *elmnts, int *etype, int *numflag,
43 idxtype *dxadj, idxtype *dadjncy)
44 {
45 int esizes[] = {-1, 3, 4, 8, 4};
46
47 if (*numflag == 1)
48 ChangeMesh2CNumbering((*ne)*esizes[*etype], elmnts);
49
50 switch (*etype) {
51 case 1:
52 TRINODALMETIS(*ne, *nn, elmnts, dxadj, dadjncy);
53 break;
54 case 2:
55 TETNODALMETIS(*ne, *nn, elmnts, dxadj, dadjncy);
56 break;
57 case 3:
58 HEXNODALMETIS(*ne, *nn, elmnts, dxadj, dadjncy);
59 break;
60 case 4:
61 QUADNODALMETIS(*ne, *nn, elmnts, dxadj, dadjncy);
62 break;
63 }
64
65 if (*numflag == 1)
66 ChangeMesh2FNumbering((*ne)*esizes[*etype], elmnts, *nn, dxadj, dadjncy);
67 }
68
69
70
71 /*****************************************************************************
72 * This function creates the dual of a finite element mesh
73 ******************************************************************************/
GENDUALMETIS(int nelmnts,int nvtxs,int etype,idxtype * elmnts,idxtype * dxadj,idxtype * dadjncy)74 void GENDUALMETIS(int nelmnts, int nvtxs, int etype, idxtype *elmnts, idxtype *dxadj, idxtype *dadjncy)
75 {
76 int i, j, jj, k, kk, kkk, l, m, n, nedges, mask;
77 idxtype *nptr, *nind;
78 idxtype *mark, ind[200], wgt[200];
79 int esize, esizes[] = {-1, 3, 4, 8, 4},
80 mgcnum, mgcnums[] = {-1, 2, 3, 4, 2};
81
82 mask = (1<<11)-1;
83 mark = idxsmalloc(mask+1, -1, "GENDUALMETIS: mark");
84
85 /* Get the element size and magic number for the particular element */
86 esize = esizes[etype];
87 mgcnum = mgcnums[etype];
88
89 /* Construct the node-element list first */
90 nptr = idxsmalloc(nvtxs+1, 0, "GENDUALMETIS: nptr");
91 for (j=esize*nelmnts, i=0; i<j; i++)
92 nptr[elmnts[i]]++;
93 MAKECSR(i, nvtxs, nptr);
94
95 nind = idxmalloc(nptr[nvtxs], "GENDUALMETIS: nind");
96 for (k=i=0; i<nelmnts; i++) {
97 for (j=0; j<esize; j++, k++)
98 nind[nptr[elmnts[k]]++] = i;
99 }
100 for (i=nvtxs; i>0; i--)
101 nptr[i] = nptr[i-1];
102 nptr[0] = 0;
103
104 for (i=0; i<nelmnts; i++)
105 dxadj[i] = esize*i;
106
107 for (i=0; i<nelmnts; i++) {
108 for (m=j=0; j<esize; j++) {
109 n = elmnts[esize*i+j];
110 for (k=nptr[n+1]-1; k>=nptr[n]; k--) {
111 if ((kk = nind[k]) <= i)
112 break;
113
114 kkk = kk&mask;
115 if ((l = mark[kkk]) == -1) {
116 ind[m] = kk;
117 wgt[m] = 1;
118 mark[kkk] = m++;
119 }
120 else if (ind[l] == kk) {
121 wgt[l]++;
122 }
123 else {
124 for (jj=0; jj<m; jj++) {
125 if (ind[jj] == kk) {
126 wgt[jj]++;
127 break;
128 }
129 }
130 if (jj == m) {
131 ind[m] = kk;
132 wgt[m++] = 1;
133 }
134 }
135 }
136 }
137 for (j=0; j<m; j++) {
138 if (wgt[j] == mgcnum) {
139 k = ind[j];
140 dadjncy[dxadj[i]++] = k;
141 dadjncy[dxadj[k]++] = i;
142 }
143 mark[ind[j]&mask] = -1;
144 }
145 }
146
147 /* Go and consolidate the dxadj and dadjncy */
148 for (j=i=0; i<nelmnts; i++) {
149 for (k=esize*i; k<dxadj[i]; k++, j++)
150 dadjncy[j] = dadjncy[k];
151 dxadj[i] = j;
152 }
153 for (i=nelmnts; i>0; i--)
154 dxadj[i] = dxadj[i-1];
155 dxadj[0] = 0;
156
157 free(mark);
158 free(nptr);
159 free(nind);
160
161 }
162
163
164
165
166 /*****************************************************************************
167 * This function creates the nodal graph of a finite element mesh
168 ******************************************************************************/
TRINODALMETIS(int nelmnts,int nvtxs,idxtype * elmnts,idxtype * dxadj,idxtype * dadjncy)169 void TRINODALMETIS(int nelmnts, int nvtxs, idxtype *elmnts, idxtype *dxadj, idxtype *dadjncy)
170 {
171 int i, j, jj, k, kk, kkk, l, m, n, nedges;
172 idxtype *nptr, *nind;
173 idxtype *mark;
174
175 /* Construct the node-element list first */
176 nptr = idxsmalloc(nvtxs+1, 0, "TRINODALMETIS: nptr");
177 for (j=3*nelmnts, i=0; i<j; i++)
178 nptr[elmnts[i]]++;
179 MAKECSR(i, nvtxs, nptr);
180
181 nind = idxmalloc(nptr[nvtxs], "TRINODALMETIS: nind");
182 for (k=i=0; i<nelmnts; i++) {
183 for (j=0; j<3; j++, k++)
184 nind[nptr[elmnts[k]]++] = i;
185 }
186 for (i=nvtxs; i>0; i--)
187 nptr[i] = nptr[i-1];
188 nptr[0] = 0;
189
190
191 mark = idxsmalloc(nvtxs, -1, "TRINODALMETIS: mark");
192
193 nedges = dxadj[0] = 0;
194 for (i=0; i<nvtxs; i++) {
195 mark[i] = i;
196 for (j=nptr[i]; j<nptr[i+1]; j++) {
197 for (jj=3*nind[j], k=0; k<3; k++, jj++) {
198 kk = elmnts[jj];
199 if (mark[kk] != i) {
200 mark[kk] = i;
201 dadjncy[nedges++] = kk;
202 }
203 }
204 }
205 dxadj[i+1] = nedges;
206 }
207
208 free(mark);
209 free(nptr);
210 free(nind);
211
212 }
213
214
215 /*****************************************************************************
216 * This function creates the nodal graph of a finite element mesh
217 ******************************************************************************/
TETNODALMETIS(int nelmnts,int nvtxs,idxtype * elmnts,idxtype * dxadj,idxtype * dadjncy)218 void TETNODALMETIS(int nelmnts, int nvtxs, idxtype *elmnts, idxtype *dxadj, idxtype *dadjncy)
219 {
220 int i, j, jj, k, kk, kkk, l, m, n, nedges;
221 idxtype *nptr, *nind;
222 idxtype *mark;
223
224 /* Construct the node-element list first */
225 nptr = idxsmalloc(nvtxs+1, 0, "TETNODALMETIS: nptr");
226 for (j=4*nelmnts, i=0; i<j; i++)
227 nptr[elmnts[i]]++;
228 MAKECSR(i, nvtxs, nptr);
229
230 nind = idxmalloc(nptr[nvtxs], "TETNODALMETIS: nind");
231 for (k=i=0; i<nelmnts; i++) {
232 for (j=0; j<4; j++, k++)
233 nind[nptr[elmnts[k]]++] = i;
234 }
235 for (i=nvtxs; i>0; i--)
236 nptr[i] = nptr[i-1];
237 nptr[0] = 0;
238
239
240 mark = idxsmalloc(nvtxs, -1, "TETNODALMETIS: mark");
241
242 nedges = dxadj[0] = 0;
243 for (i=0; i<nvtxs; i++) {
244 mark[i] = i;
245 for (j=nptr[i]; j<nptr[i+1]; j++) {
246 for (jj=4*nind[j], k=0; k<4; k++, jj++) {
247 kk = elmnts[jj];
248 if (mark[kk] != i) {
249 mark[kk] = i;
250 dadjncy[nedges++] = kk;
251 }
252 }
253 }
254 dxadj[i+1] = nedges;
255 }
256
257 free(mark);
258 free(nptr);
259 free(nind);
260
261 }
262
263
264 /*****************************************************************************
265 * This function creates the nodal graph of a finite element mesh
266 ******************************************************************************/
HEXNODALMETIS(int nelmnts,int nvtxs,idxtype * elmnts,idxtype * dxadj,idxtype * dadjncy)267 void HEXNODALMETIS(int nelmnts, int nvtxs, idxtype *elmnts, idxtype *dxadj, idxtype *dadjncy)
268 {
269 int i, j, jj, k, kk, kkk, l, m, n, nedges;
270 idxtype *nptr, *nind;
271 idxtype *mark;
272 int table[8][3] = {1, 3, 4,
273 0, 2, 5,
274 1, 3, 6,
275 0, 2, 7,
276 0, 5, 7,
277 1, 4, 6,
278 2, 5, 7,
279 3, 4, 6};
280
281 /* Construct the node-element list first */
282 nptr = idxsmalloc(nvtxs+1, 0, "HEXNODALMETIS: nptr");
283 for (j=8*nelmnts, i=0; i<j; i++)
284 nptr[elmnts[i]]++;
285 MAKECSR(i, nvtxs, nptr);
286
287 nind = idxmalloc(nptr[nvtxs], "HEXNODALMETIS: nind");
288 for (k=i=0; i<nelmnts; i++) {
289 for (j=0; j<8; j++, k++)
290 nind[nptr[elmnts[k]]++] = i;
291 }
292 for (i=nvtxs; i>0; i--)
293 nptr[i] = nptr[i-1];
294 nptr[0] = 0;
295
296
297 mark = idxsmalloc(nvtxs, -1, "HEXNODALMETIS: mark");
298
299 nedges = dxadj[0] = 0;
300 for (i=0; i<nvtxs; i++) {
301 mark[i] = i;
302 for (j=nptr[i]; j<nptr[i+1]; j++) {
303 jj=8*nind[j];
304 for (k=0; k<8; k++) {
305 if (elmnts[jj+k] == i)
306 break;
307 }
308 ASSERT(k != 8);
309
310 /* You found the index, now go and put the 3 neighbors */
311 kk = elmnts[jj+table[k][0]];
312 if (mark[kk] != i) {
313 mark[kk] = i;
314 dadjncy[nedges++] = kk;
315 }
316 kk = elmnts[jj+table[k][1]];
317 if (mark[kk] != i) {
318 mark[kk] = i;
319 dadjncy[nedges++] = kk;
320 }
321 kk = elmnts[jj+table[k][2]];
322 if (mark[kk] != i) {
323 mark[kk] = i;
324 dadjncy[nedges++] = kk;
325 }
326 }
327 dxadj[i+1] = nedges;
328 }
329
330 free(mark);
331 free(nptr);
332 free(nind);
333
334 }
335
336
337 /*****************************************************************************
338 * This function creates the nodal graph of a finite element mesh
339 ******************************************************************************/
QUADNODALMETIS(int nelmnts,int nvtxs,idxtype * elmnts,idxtype * dxadj,idxtype * dadjncy)340 void QUADNODALMETIS(int nelmnts, int nvtxs, idxtype *elmnts, idxtype *dxadj, idxtype *dadjncy)
341 {
342 int i, j, jj, k, kk, kkk, l, m, n, nedges;
343 idxtype *nptr, *nind;
344 idxtype *mark;
345 int table[4][2] = {1, 3,
346 0, 2,
347 1, 3,
348 0, 2};
349
350 /* Construct the node-element list first */
351 nptr = idxsmalloc(nvtxs+1, 0, "QUADNODALMETIS: nptr");
352 for (j=4*nelmnts, i=0; i<j; i++)
353 nptr[elmnts[i]]++;
354 MAKECSR(i, nvtxs, nptr);
355
356 nind = idxmalloc(nptr[nvtxs], "QUADNODALMETIS: nind");
357 for (k=i=0; i<nelmnts; i++) {
358 for (j=0; j<4; j++, k++)
359 nind[nptr[elmnts[k]]++] = i;
360 }
361 for (i=nvtxs; i>0; i--)
362 nptr[i] = nptr[i-1];
363 nptr[0] = 0;
364
365
366 mark = idxsmalloc(nvtxs, -1, "QUADNODALMETIS: mark");
367
368 nedges = dxadj[0] = 0;
369 for (i=0; i<nvtxs; i++) {
370 mark[i] = i;
371 for (j=nptr[i]; j<nptr[i+1]; j++) {
372 jj=4*nind[j];
373 for (k=0; k<4; k++) {
374 if (elmnts[jj+k] == i)
375 break;
376 }
377 ASSERT(k != 4);
378
379 /* You found the index, now go and put the 2 neighbors */
380 kk = elmnts[jj+table[k][0]];
381 if (mark[kk] != i) {
382 mark[kk] = i;
383 dadjncy[nedges++] = kk;
384 }
385 kk = elmnts[jj+table[k][1]];
386 if (mark[kk] != i) {
387 mark[kk] = i;
388 dadjncy[nedges++] = kk;
389 }
390 }
391 dxadj[i+1] = nedges;
392 }
393
394 free(mark);
395 free(nptr);
396 free(nind);
397
398 }
399