1 /*
2 ElmerGrid - A simple mesh generation and manipulation utility
3 Copyright (C) 1995- , CSC - IT Center for Science Ltd.
4
5 Author: Peter Raback
6 Email: elmeradm@csc.fi
7 Address: CSC - IT Center for Science Ltd.
8 Keilaranta 14
9 02101 Espoo, Finland
10
11 This program is free software; you can redistribute it and/or
12 modify it under the terms of the GNU General Public License
13 as published by the Free Software Foundation; either version 2
14 of the License, or (at your option) any later version.
15
16 This program is distributed in the hope that it will be useful,
17 but WITHOUT ANY WARRANTY; without even the implied warranty of
18 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
19 GNU General Public License for more details.
20
21 You should have received a copy of the GNU General Public License
22 along with this program; if not, write to the Free Software
23 Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
24 */
25
26 /* -------------------------------: egparallel.c :----------------------------
27 This module includes routines related to parallel file format of Elmer suite.
28 */
29
30 #include <stdio.h>
31 #include <unistd.h>
32 #include <stddef.h>
33 #include <stdlib.h>
34 #include <string.h>
35 #include <math.h>
36
37 #if HAVE_UNISTD_H
38 #include <unistd.h>
39 #endif
40
41 #include <stdarg.h>
42 #include <sys/stat.h>
43 #include <sys/types.h>
44
45 #include "egutils.h"
46 #include "egdef.h"
47 #include "egtypes.h"
48 #include "egmesh.h"
49 #include "egparallel.h"
50 #include "../config.h"
51
52 #if USE_METIS
53 #include <metis.h>
54 #endif
55
56
FuseSolutionElmerPartitioned(char * prefix,char * outfile,int decimals,int parts,int minstep,int maxstep,int dstep,int info)57 int FuseSolutionElmerPartitioned(char *prefix,char *outfile,int decimals,int parts,
58 int minstep, int maxstep, int dstep, int info)
59 {
60 #define LONGLINE 2048
61 int *noknots,*noelements,novctrs,elemcode;
62 int totknots,totelements,sumknots,sumelements;
63 int timesteps,i,j,k,l,step;
64 int ind[MAXNODESD3];
65 int nofiles,activestep;
66 Real *res, x, y, z;
67 FILE *in[MAXPARTITIONS+1],*intest,*out;
68 char line[LONGLINE],filename[MAXFILESIZE],text[MAXNAMESIZE],outstyle[MAXFILESIZE];
69 char *cp, *charend;
70
71 if(minstep || maxstep || dstep) {
72 if(info) printf("Saving results in the interval from %d to %d with step %d\n",minstep,maxstep,dstep);
73 }
74
75 for(i=0;;i++) {
76 sprintf(filename,"%s.ep.%d",prefix,i);
77 if ((intest = fopen(filename,"r")) == NULL) break;
78 if(i<=MAXPARTITIONS) in[i] = intest;
79 }
80 if(i > MAXPARTITIONS) {
81 printf("**********************************************************\n");
82 printf("Only data for %d partitions is fused (%d)\n",MAXPARTITIONS,i);
83 printf("**********************************************************\n");
84 i = MAXPARTITIONS;
85 }
86 nofiles = i;
87
88 if(nofiles < 2) {
89 printf("Opening of partitioned data from file %s wasn't successful!\n",
90 filename);
91 return(2);
92 } else {
93 if(parts > 0) nofiles = MIN(parts,nofiles);
94 if(info) printf("Loading Elmer results from %d partitions.\n",nofiles);
95 }
96
97 if(minstep || maxstep || dstep) {
98 if(info) printf("Saving results in the interval from %d to %d with step %d\n",minstep,maxstep,dstep);
99 }
100
101 noknots = Ivector(0,nofiles-1);
102 noelements = Ivector(0,nofiles-1);
103
104 sumknots = 0;
105 sumelements = 0;
106
107 for(i=0;i<nofiles;i++) {
108 charend = fgets(line,LONGLINE,in[i]);
109 if(i==0) {
110 cp = line;
111 noknots[i] = next_int(&cp);
112 noelements[i] = next_int(&cp);
113 novctrs = next_int(&cp);
114 timesteps = next_int(&cp);
115 }
116 else {
117 sscanf(line,"%d %d",&noknots[i],&noelements[i]);
118 }
119 sumknots += noknots[i];
120 sumelements += noelements[i];
121 }
122 totknots = sumknots;
123 totelements = sumelements;
124 res = Rvector(1,novctrs);
125
126 if(info) printf("There are altogether %d nodes and %d elements.\n",totknots,sumelements);
127
128
129 AddExtension(outfile,filename,"ep");
130 if(info) printf("Saving ElmerPost data to %s.\n",filename);
131 out = fopen(filename,"w");
132 if(out == NULL) {
133 printf("opening of file was not successful\n");
134 return(3);
135 }
136
137 i = timesteps;
138 if(minstep || maxstep || dstep) {
139 if ( dstep>1 ) {
140 i = 0;
141 for(step = minstep; step <= maxstep; step++)
142 if((step-minstep)%dstep==0) i++;
143 } else i=maxstep-minstep+1;
144 }
145 fprintf(out,"%d %d %d %d %s %s",totknots,totelements,novctrs+1,i,"scalar: Partition",cp);
146
147 if(info) printf("Reading and writing %d coordinates.\n",totknots);
148 sprintf(outstyle,"%%.%dg %%.%dg %%.%dg\n",decimals,decimals,decimals);
149
150 for(j=0; j < nofiles; j++) {
151 for(i=1; i <= noknots[j]; i++) {
152 do {
153 charend = fgets(line,LONGLINE,in[j]);
154 } while(line[0] == '#');
155
156 sscanf(line,"%le %le %le",&x,&y,&z);
157 fprintf(out,outstyle,x,y,z);
158 }
159 }
160
161 if(info) printf("Reading and writing %d element topologies.\n",totelements);
162 sumknots = 0;
163
164 for(j=0; j < nofiles; j++) {
165 for(i=1; i <= noelements[j]; i++) {
166 do {
167 charend = fgets(line,LONGLINE,in[j]);
168 } while (line[0] == '#');
169
170 sscanf(line,"%s",text);
171 cp = strstr(line," ");
172
173 elemcode = next_int(&cp);
174
175 for(k=0;k< elemcode%100 ;k++) {
176 /* Dirty trick for long lines */
177 l = strspn(cp," ");
178 if( l == 0) {
179 charend = fgets(line,LONGLINE,in[j]);
180 cp = line;
181 }
182 ind[k] = next_int(&cp);
183 }
184 if(elemcode == 102) elemcode = 101;
185
186 fprintf(out,"%s %d",text,elemcode);
187 for(k=0;k < elemcode%100 ;k++)
188 fprintf(out," %d",ind[k]+sumknots);
189 fprintf(out,"\n");
190 }
191 sumknots += noknots[j];
192 }
193
194 if(info) printf("Reading and writing %d degrees of freedom.\n",novctrs);
195 sprintf(outstyle,"%%.%dg ",decimals);
196
197 activestep = FALSE;
198 if(maxstep) timesteps = MIN(timesteps, maxstep);
199
200 for(step = 1; step <= timesteps; step++) {
201
202 if (step>=minstep) {
203 if ( dstep>0 ) {
204 activestep=((step-minstep)%dstep==0);
205 } else activestep=TRUE;
206 }
207
208 for(k=0;k<nofiles;k++)
209 for(i=1; i <= noknots[k]; i++) {
210 do {
211 charend = fgets(line,LONGLINE,in[k]);
212 if (activestep) {
213 if(k==0 && strstr(line,"#time")) {
214 fprintf(out,"%s",line);
215 fprintf(stderr,"%s",line);
216 }
217 }
218 }
219 while (line[0] == '#');
220
221 if(activestep) {
222 cp = line;
223 for(j=1;j <= novctrs;j++)
224 res[j] = next_real(&cp);
225
226 fprintf(out,"%d ",k+1);
227 for(j=1;j <= novctrs;j++)
228 fprintf(out,outstyle,res[j]);
229 fprintf(out,"\n");
230 }
231
232 }
233 }
234
235
236 for(i=0;i<nofiles;i++)
237 fclose(in[i]);
238 fclose(out);
239
240 if(info) printf("Successfully fused partitioned Elmer results\n");
241
242 return(0);
243 }
244
245
246
247
CreatePartitionTable(struct FemType * data,int info)248 static int CreatePartitionTable(struct FemType *data,int info)
249 {
250 int i,j,k,m,noelements,noknots,partitions,nonodes,periodic;
251 int maxneededtimes,sharings,part,ind,hit,notinany,debug;
252 int *indxper;
253
254 printf("Creating a table showing all parenting partitions of nodes.\n");
255
256 if(data->maxpartitiontable) {
257 printf("The partition table already exists!\n");
258 smallerror("Partition table not done");
259 return(0);
260 }
261
262 partitions = data->nopartitions;
263 if( partitions <= 1 ) {
264 bigerror("Cannot do partition table for less than two partitions");
265 }
266
267
268 noelements = data->noelements;
269 periodic = data->periodicexist;
270 noknots = data->noknots;
271 if(periodic) indxper = data->periodic;
272
273 maxneededtimes = 0;
274 sharings = 0;
275
276 /* Make the partition list taking into account periodicity */
277 for(i=1;i<=noelements;i++) {
278 part = data->elempart[i];
279 nonodes = data->elementtypes[i] % 100;
280
281 for(j=0;j < nonodes;j++) {
282 ind = data->topology[i][j];
283 if(periodic) ind = indxper[ind];
284
285 debug = FALSE;
286 if(debug) printf("ind=%d i=%d j=%d part=%d\n",ind,i,j,part);
287
288 hit = 0;
289 for(k=1;k<=maxneededtimes;k++) {
290 if(data->partitiontable[k][ind] == 0) hit = k;
291 if(data->partitiontable[k][ind] == part) hit = -k;
292 if(hit) break;
293 }
294 if( hit > 0) {
295 data->partitiontable[hit][ind] = part;
296 if(hit == 2) sharings++;
297 }
298 else if(hit == 0) {
299 maxneededtimes++;
300 data->partitiontable[maxneededtimes] = Ivector(1,noknots);
301 for(m=1;m<=noknots;m++)
302 data->partitiontable[maxneededtimes][m] = 0;
303 data->partitiontable[maxneededtimes][ind] = part;
304 if(maxneededtimes == 2) sharings++;
305 }
306 if(debug) printf("hit = %d\n",hit);
307 }
308 }
309
310
311 /* Make the partitiontable such that the owner node is the first one in the list */
312 notinany = 0;
313 for(i=1;i<=noknots;i++) {
314
315 /* Skip the periodic nodes and take care of them later */
316 if(periodic)
317 if(i != indxper[i]) continue;
318
319 hit = FALSE;
320 for(k=1;k<=maxneededtimes;k++) {
321 if(!data->partitiontable[k][i]) break;
322 if(data->partitiontable[k][i] == data->nodepart[i]) {
323 hit = k;
324 break;
325 }
326 }
327 if( hit > 1 ) {
328 data->partitiontable[hit][i] = data->partitiontable[1][i];
329 data->partitiontable[1][i] = data->nodepart[i];
330 }
331 else if(!hit) {
332 if(0) {
333 printf("Node %d in partition %d not in the table!\n",i,data->nodepart[i]);
334 if(periodic) printf("indexper: %d\n",indxper[i]);
335 }
336
337 notinany++;
338 data->nodepart[i] = data->partitiontable[1][i];
339 }
340 }
341
342
343 /* For periodic counterparts copy the table and ownership */
344 if(periodic) {
345 for(i=1;i<=noknots;i++) {
346 ind = indxper[i];
347 if(ind == i) continue;
348 for(k=1;k<=maxneededtimes;k++) {
349 if( !data->partitiontable[k][ind] ) break;
350 data->partitiontable[k][i] = data->partitiontable[k][ind];
351 }
352 data->nodepart[i] = data->nodepart[ind];
353 }
354 }
355
356
357 if(info) {
358 printf("Nodes belong to %d partitions in maximum\n",maxneededtimes);
359 printf("There are %d shared nodes which is %.2f %% of all nodes.\n",
360 sharings,(100.*sharings)/noknots);
361 printf("The initial owner was not any of the elements for %d nodes\n",notinany);
362 }
363
364
365 if(0) for(i=1;i<=noknots;i++) {
366 if(data->partitiontable[maxneededtimes][i]) {
367 printf("node %d parts: ",i);
368 for(k=1;k<=maxneededtimes;k++)
369 printf("%d ",data->partitiontable[k][i]);
370 printf("\n");
371 }
372 }
373
374 data->maxpartitiontable = maxneededtimes;
375 data->partitiontableexists = TRUE;
376 return(0);
377 }
378
379
380
PartitionElementsByNodes(struct FemType * data,int info)381 static int PartitionElementsByNodes(struct FemType *data,int info)
382 /* Given nodal partitioning determine the elemental partitioning.
383 This is usually suboptimal, it is preferable to do the other way around. */
384 {
385 int i,j,noelements,nopartitions,part,maxpart,maxpart2,minpart;
386 int *elempart,*nodepart,*nodesinpart,*cuminpart,**knows,**cumknows,set;
387
388 if(!data->partitionexist) return(1);
389
390 noelements = data->noelements;
391 nopartitions = data->nopartitions;
392 elempart = data->elempart;
393 nodepart = data->nodepart;
394
395 nodesinpart = Ivector(1,nopartitions);
396 cuminpart = Ivector(1,nopartitions);
397 for(j=1;j<=nopartitions;j++)
398 cuminpart[j] = 0;
399
400 knows = Imatrix(1,nopartitions,1,nopartitions);
401 cumknows = Imatrix(1,nopartitions,1,nopartitions);
402 for(i=1;i<=nopartitions;i++)
403 for(j=1;j<=nopartitions;j++)
404 knows[i][j] = cumknows[i][j] = 0;
405
406 set = FALSE;
407
408 omstart:
409
410 /* In the first round count the equally joined elements and
411 on the second round split them equally using cumulative numbering */
412
413 for(i=1;i<=noelements;i++) {
414 for(j=1;j<=nopartitions;j++)
415 nodesinpart[j] = 0;
416 for(j=0;j<data->elementtypes[i] % 100;j++) {
417 part = nodepart[data->topology[i][j]];
418 nodesinpart[part] += 1;
419 }
420 maxpart = maxpart2 = 1;
421 for(j=1;j<=nopartitions;j++)
422 if(nodesinpart[j] > nodesinpart[maxpart]) maxpart = j;
423 if(maxpart == 1) maxpart2 = 2;
424 for(j=1;j<=nopartitions;j++) {
425 if(j == maxpart) continue;
426 if(nodesinpart[j] > nodesinpart[maxpart2]) maxpart2 = j;
427 }
428
429 if(nodesinpart[maxpart] > nodesinpart[maxpart2]) {
430 if(set)
431 elempart[i] = maxpart;
432 else
433 cuminpart[maxpart] += 1;
434 }
435 else {
436 if(set) {
437 cumknows[maxpart][maxpart2] += 1;
438 if( cumknows[maxpart][maxpart2] > knows[maxpart][maxpart2] / 2) {
439 elempart[i] = maxpart2;
440 cuminpart[maxpart2] += 1;
441 }
442 else {
443 elempart[i] = maxpart;
444 cuminpart[maxpart] += 1;
445 }
446 }
447 else
448 knows[maxpart][maxpart2] += 1;
449 }
450 }
451
452 if(!set) {
453 set = TRUE;
454 goto omstart;
455 }
456
457 minpart = maxpart = cuminpart[1];
458 for(j=1;j<=nopartitions;j++) {
459 minpart = MIN( minpart, cuminpart[j]);
460 maxpart = MAX( maxpart, cuminpart[j]);
461 }
462
463 if(info) {
464 printf("Set the element partitions by the dominating nodal partition\n");
465 printf("There are from %d to %d elements in the %d partitions.\n",minpart,maxpart,nopartitions);
466 }
467
468 free_Ivector(nodesinpart,1,nopartitions);
469 free_Ivector(cuminpart,1,nopartitions);
470 free_Imatrix(knows,1,nopartitions,1,nopartitions);
471 free_Imatrix(cumknows,1,nopartitions,1,nopartitions);
472
473 return(0);
474 }
475
476
PartitionNodesByElements(struct FemType * data,int info)477 static int PartitionNodesByElements(struct FemType *data,int info)
478 /* Given the elemental partitioning set the nodal ownership.
479 This is optimal for Elmer since the elemental partitioning is primary. */
480 {
481 int i,j,k,noknots,nopartitions,part,minpart,maxpart;
482 int maxpart2,*cuminpart,**knows,**cumknows,set;
483 int *elempart,*nodepart,*nodesinpart;
484 int *invrow,*invcol;
485
486 if(!data->partitionexist) return(1);
487
488 CreateInverseTopology(data, info);
489
490 invrow = data->invtopo.rows;
491 invcol = data->invtopo.cols;
492
493 noknots = data->noknots;
494 nopartitions = data->nopartitions;
495 elempart = data->elempart;
496 nodepart = data->nodepart;
497
498 if(info) printf("Number of nodal partitions: %d\n",nopartitions);
499
500 nodesinpart = Ivector(1,nopartitions);
501 cuminpart = Ivector(1,nopartitions);
502 for(j=1;j<=nopartitions;j++)
503 cuminpart[j] = 0;
504
505 knows = Imatrix(1,nopartitions,1,nopartitions);
506 cumknows = Imatrix(1,nopartitions,1,nopartitions);
507 for(i=1;i<=nopartitions;i++)
508 for(j=1;j<=nopartitions;j++)
509 knows[i][j] = cumknows[i][j] = 0;
510
511 set = FALSE;
512
513 omstart:
514
515 for(i=1;i<=noknots;i++) {
516
517 for(j=1;j<=nopartitions;j++)
518 nodesinpart[j] = 0;
519
520 /* Tag the number of owner partitions */
521 for(j=invrow[i-1];j<invrow[i];j++) {
522 k = invcol[j]+1;
523 part = elempart[k];
524 if( part > 0 ) nodesinpart[part] += 1;
525 }
526
527 /* Find the partition with maximum number of hits */
528 maxpart = maxpart2 = 0;
529 for(j=1;j<=nopartitions;j++)
530 if(nodesinpart[j] > 0 ) {
531 if(!maxpart) {
532 maxpart = j;
533 }
534 else if(nodesinpart[j] > nodesinpart[maxpart]) {
535 maxpart = j;
536 }
537 }
538
539 /* Find the partition with 2nd largest number of hits */
540 for(j=1;j<=nopartitions;j++)
541 if(nodesinpart[j] > 0 ) {
542 if(j == maxpart) continue;
543 if(!maxpart2) {
544 maxpart2 = j;
545 }
546 else if(nodesinpart[j] > nodesinpart[maxpart2]) {
547 maxpart2 = j;
548 }
549 }
550
551 /* If there is a clear dominator use that */
552 if(!maxpart && !maxpart2)
553 printf("Node is not included in any partition: %d\n",i);
554 else if(!maxpart2) {
555 nodepart[i] = maxpart;
556 cuminpart[maxpart] += 1;
557 }
558 else
559 if(nodesinpart[maxpart] > nodesinpart[maxpart2]) {
560 if(set)
561 nodepart[i] = maxpart;
562 else
563 cuminpart[maxpart] += 1;
564 }
565
566 /* Otherwise make a half and half split between the major owners */
567 else {
568 if(set) {
569 cumknows[maxpart][maxpart2] += 1;
570 if( cumknows[maxpart][maxpart2] > knows[maxpart][maxpart2] / 2) {
571 nodepart[i] = maxpart2;
572 cuminpart[maxpart2] += 1;
573 }
574 else {
575 nodepart[i] = maxpart;
576 cuminpart[maxpart] += 1;
577 }
578 }
579 else {
580 knows[maxpart][maxpart2] += 1;
581 }
582 }
583 }
584
585 if(!set) {
586 set = TRUE;
587 goto omstart;
588 }
589
590 minpart = maxpart = cuminpart[1];
591 for(j=1;j<=nopartitions;j++) {
592 minpart = MIN( minpart, cuminpart[j]);
593 maxpart = MAX( maxpart, cuminpart[j]);
594 }
595
596 if(info) {
597 printf("Set the node partitions by the dominating element partition.\n");
598 printf("There are from %d to %d nodes in the %d partitions.\n",minpart,maxpart,nopartitions);
599 }
600
601 free_Ivector(nodesinpart,1,nopartitions);
602 free_Ivector(cuminpart,1,nopartitions);
603 free_Imatrix(knows,1,nopartitions,1,nopartitions);
604 free_Imatrix(cumknows,1,nopartitions,1,nopartitions);
605
606 return(0);
607 }
608
609
610
PartitionSimpleElements(struct FemType * data,struct ElmergridType * eg,struct BoundaryType * bound,int dimpart[],int dimper[],int partorder,Real corder[],Real parttol,int info)611 int PartitionSimpleElements(struct FemType *data,struct ElmergridType *eg,struct BoundaryType *bound,
612 int dimpart[],int dimper[],int partorder, Real corder[],
613 Real parttol, int info)
614 /* Partition elements recursively in major directions.
615 This may be the optimal method of partitioning for simple geometries. */
616 {
617 int i,j,k,l,kprev,ind,minpart,maxpart;
618 int noknots,noelements,nonodes,elemsinpart,periodic;
619 int partitions1,partitions2,partitions3,partitions;
620 int vpartitions1,vpartitions2,vpartitions3,vpartitions;
621 int noelements0,noelements1,noparts0;
622 int *indx,*nopart,*inpart,*elemconnect;
623 Real *arrange;
624 Real x,y,z,cx,cy,cz;
625
626 if(info) printf("PartitionSimpleElements\n");
627
628 noelements = data->noelements;
629 noknots = data->noknots;
630
631 partitions1 = dimpart[0];
632 partitions2 = dimpart[1];
633 partitions3 = dimpart[2];
634 if(data->dim < 3) partitions3 = 1;
635 partitions = partitions1 * partitions2 * partitions3;
636
637 if(partitions1 < 2 && partitions2 < 2 && partitions3 < 2 && !eg->connect) {
638 printf("Some of the divisions must be larger than one: %d %d %d\n",
639 partitions1, partitions2, partitions3 );
640 bigerror("Partitioning not performed");
641 }
642
643 if(partitions >= noelements) {
644 printf("There must be fever partitions than elements (%d vs %d)!\n",
645 partitions,noelements);
646 bigerror("Partitioning not performed");
647 }
648
649 if( eg->partbcz > 1 || eg->partbcr )
650 PartitionConnectedElements1D(data,bound,eg,info);
651 else if( eg->partbcmetis > 1 )
652 PartitionConnectedElementsMetis(data,bound,eg->partbcmetis,3,info);
653 else if( eg->connect )
654 PartitionConnectedElementsStraight(data,bound,eg,info);
655
656
657 if( data->nodeconnectexist )
658 ExtendBoundaryPartitioning(data,bound,eg->partbclayers,info);
659
660 if(!data->partitionexist) {
661 data->partitionexist = TRUE;
662 data->elempart = Ivector(1,noelements);
663 data->nodepart = Ivector(1,noknots);
664 data->nopartitions = partitions;
665 }
666 inpart = data->elempart;
667
668 vpartitions1 = partitions1;
669 vpartitions2 = partitions2;
670 vpartitions3 = partitions3;
671
672 if(0) printf("connect: %d %d\n",data->elemconnectexist,data->nodeconnectexist);
673
674 noparts0 = 0;
675 noelements0 = 0;
676 if( data->elemconnectexist || data->nodeconnectexist) {
677 elemconnect = data->elemconnect;
678 noparts0 = 0;
679 for(i=1;i<=data->noelements;i++) {
680 if( elemconnect[i] ) noelements0 = noelements0 + 1;
681 noparts0 = MAX( noparts0, elemconnect[i] );
682 }
683 if(info) printf("There are %d initial partitions in the connected mesh\n",noparts0);
684 if(info) printf("There are %d initial elements in the connected mesh\n",noelements0);
685 }
686 noelements1 = noelements - noelements0;
687
688 vpartitions1 = partitions1;
689 vpartitions2 = partitions2;
690 vpartitions3 = partitions3;
691
692 periodic = dimper[0] || dimper[1] || dimper[2];
693 if(periodic) {
694 if(dimper[0] && partitions1 > 1) vpartitions1 *= 2;
695 if(dimper[1] && partitions2 > 1) vpartitions2 *= 2;
696 if(dimper[2] && partitions3 > 1) vpartitions3 *= 2;
697 }
698 vpartitions = vpartitions1 * vpartitions2 * vpartitions3;
699 nopart = Ivector(1,vpartitions+noparts0);
700
701 if( vpartitions == 1 && data->elemconnectexist ) {
702 if(info) printf("Only one regular partitions requested, skipping 2nd part of hybrid partitioning\n");
703 for(j=1;j<=noelements;j++) {
704 if( elemconnect[j] > 0 )
705 inpart[j] = elemconnect[j];
706 else
707 inpart[j] = noparts0 + 1;
708 }
709 partitions = noparts0 + 1;
710 data->nopartitions = partitions;
711 goto skippart;
712 }
713
714
715
716 if(info) printf("Making a simple partitioning for %d elements in %d-dimensions.\n",
717 noelements1,data->dim);
718
719 arrange = Rvector(1,noelements);
720 indx = Ivector(1,noelements);
721
722 if(partorder) {
723 cx = corder[0];
724 cy = corder[1];
725 cz = corder[2];
726 }
727 else {
728 cx = 1.0;
729 cy = 0.0001;
730 cz = cy*cy;
731 }
732 z = 0.0;
733
734 for(i=1;i<=noelements;i++)
735 inpart[i] = 1;
736
737 if(vpartitions1 > 1) {
738
739 if(info) printf("Ordering 1st direction with (%.3lg*x + %.3lg*y + %.3lg*z)\n",cx,cy,cz);
740
741 for(j=1;j<=noelements;j++) {
742 if( data->elemconnectexist ) {
743 if( elemconnect[j] > 0 ) {
744 arrange[j] = 1.0e9;
745 continue;
746 }
747 }
748 nonodes = data->elementtypes[j]%100;
749 x = y = z = 0.0;
750 for(i=0;i<nonodes;i++) {
751 k = data->topology[j][i];
752 x += data->x[k];
753 y += data->y[k];
754 z += data->z[k];
755 }
756 arrange[j] = (cx*x + cy*y + cz*z) / nonodes;
757 }
758
759 SortIndex(noelements,arrange,indx);
760
761 for(i=1;i<=noelements1;i++) {
762 ind = indx[i];
763 k = (i*vpartitions1-1)/noelements1+1;
764 inpart[ind] = k;
765 }
766 }
767
768
769 /* Partition in the 2nd direction taking into account the 1st direction */
770 if(vpartitions2 > 1) {
771 if(info) printf("Ordering in the 2nd direction.\n");
772
773 for(j=1;j<=noelements;j++) {
774 if( data->elemconnectexist ) {
775 if( elemconnect[j] > 0 ) {
776 arrange[j] = 1.0e9;
777 continue;
778 }
779 }
780 nonodes = data->elementtypes[j]%100;
781 x = y = z = 0.0;
782 for(i=0;i<nonodes;i++) {
783 k = data->topology[j][i];
784 x += data->x[k];
785 y += data->y[k];
786 z += data->z[k];
787 }
788 arrange[j] = (-cy*x + cx*y + cz*z) / nonodes;
789 }
790 SortIndex(noelements,arrange,indx);
791
792 for(i=1;i<=vpartitions;i++)
793 nopart[i] = 0;
794
795 elemsinpart = noelements1 / (vpartitions1*vpartitions2);
796 for(i=1;i<=noelements1;i++) {
797 j = 0;
798 ind = indx[i];
799 do {
800 j++;
801 k = (inpart[ind]-1) * vpartitions2 + j;
802 }
803 while(nopart[k] >= elemsinpart && j < vpartitions2);
804
805 nopart[k] += 1;
806 inpart[ind] = (inpart[ind]-1)*vpartitions2 + j;
807 }
808 }
809
810 /* Partition in the 3rd direction taking into account the 1st and 2nd direction */
811 if(vpartitions3 > 1) {
812 if(info) printf("Ordering in the 3rd direction.\n");
813
814 for(j=1;j<=noelements;j++) {
815 if( data->elemconnectexist ) {
816 if( elemconnect[j] > 0 ) {
817 arrange[j] = 1.0e9;
818 continue;
819 }
820 }
821
822 nonodes = data->elementtypes[j]%100;
823 x = y = z = 0.0;
824 for(i=0;i<nonodes;i++) {
825 k = data->topology[j][i];
826 x += data->x[k];
827 y += data->y[k];
828 z += data->z[k];
829 }
830 arrange[j] = (-cz*x - cy*y + cx*z) / nonodes;
831 }
832
833 SortIndex(noelements,arrange,indx);
834
835 for(i=1;i<=vpartitions;i++)
836 nopart[i] = 0;
837
838 elemsinpart = noelements1 / (vpartitions1*vpartitions2*vpartitions3);
839 for(i=1;i<=noelements1;i++) {
840 j = 0;
841 ind = indx[i];
842 do {
843 j++;
844 k = (inpart[ind]-1)*vpartitions3 + j;
845 }
846 while(nopart[k] >= elemsinpart && j < vpartitions3);
847
848 nopart[k] += 1;
849 inpart[ind] = (inpart[ind]-1)*vpartitions3 + j;
850 }
851 }
852
853 /* This piece of code may be used to assign "strides" back to same partition.
854 After performing the geometric partitioning we revisit whether elements
855 are too close to each other and yet in different partition in the chosen measure.
856 For example, if chosen measure has the smallest coefficient in z-direction then
857 a suitable value will revert elements to being on the same stride. */
858 if(parttol > 1.0e-20 ) {
859 if(0) {
860 printf("Original partition counts:\n");
861 for(i=1;i<=partitions;i++)
862 printf("nopart[%d] = %d\n",i,nopart[i]);
863 }
864
865 for(j=1;j<=noelements;j++) {
866 if( data->elemconnectexist ) {
867 if( elemconnect[j] > 0 ) {
868 arrange[j] = 1.0e9;
869 continue;
870 }
871 }
872 nonodes = data->elementtypes[j]%100;
873 x = y = z = 0.0;
874 for(i=0;i<nonodes;i++) {
875 k = data->topology[j][i];
876 x += data->x[k];
877 y += data->y[k];
878 z += data->z[k];
879 }
880 arrange[j] = (cx*x + cy*y + cz*z) / nonodes;
881 }
882 SortIndex(noelements,arrange,indx);
883
884 kprev = inpart[indx[1]];
885 l = 0;
886 for(i=2;i<=noelements1;i++) {
887 ind = indx[i];
888 k = inpart[ind];
889 if(k != kprev ) {
890 if(fabs(arrange[i]-arrange[i-1]) < parttol ) {
891 nopart[k] -= 1;
892 nopart[kprev] += 1;
893 inpart[ind] = kprev;
894 k = kprev;
895 if(0) printf("arrange %d %d %d %d %.6lg %.6lg\n",i,ind,k,kprev,arrange[i],arrange[i-1]);
896 l++;
897 }
898 }
899 kprev = k;
900 }
901 if(0) {
902 printf("Updated partition counts:\n");
903 for(i=1;i<=partitions;i++)
904 printf("nopart[%d] = %d\n",i,nopart[i]);
905 }
906 if(info) printf("Number of partition changes due to tolerance: %d\n",l);
907 }
908
909
910 free_Rvector(arrange,1,noelements);
911 free_Ivector(indx,1,noelements);
912
913
914 if( data->elemconnectexist ) {
915 for(j=1;j<=noelements;j++) {
916 if( elemconnect[j] > 0 )
917 inpart[j] = elemconnect[j];
918 else
919 inpart[j] = inpart[j] + noparts0;
920 }
921 partitions = partitions + noparts0;
922 data->nopartitions = partitions;
923 }
924
925
926 /* For periodic systems the number of virtual partitions is larger. Now map the mesh so that the
927 1st and last partition for each direction will be joined */
928 if(periodic) {
929 int *partmap;
930 int p1,p2,p3,q1,q2,q3;
931 int P,Q;
932
933 if(data->elemconnectexist ) {
934 bigerror("Cannot use connect flag with periodic systems\n");
935 }
936
937 p1=p2=p3=1;
938 partmap = Ivector(1,vpartitions);
939 for(i=1;i<=vpartitions;i++)
940 partmap[i] = 0;
941 for(p1=1;p1<=vpartitions1;p1++) {
942 q1 = p1;
943 if(dimper[0] && vpartitions1 > 1) {
944 if(q1==vpartitions1) q1 = 0;
945 q1 = q1/2 + 1;
946 }
947 for(p2=1;p2<=vpartitions2;p2++) {
948 q2 = p2;
949 if(dimper[1] && vpartitions2 > 1) {
950 if(q2==vpartitions2) q2 = 0;
951 q2 = q2/2 + 1;
952 }
953 for(p3=1;p3<=vpartitions3;p3++) {
954 q3 = p3;
955 if(dimper[2] && vpartitions3 > 1) {
956 if(q3==vpartitions3) q3 = 0;
957 q3 = q3/2 + 1;
958 }
959
960 P = vpartitions3 * vpartitions2 * (p1 - 1) + vpartitions3 * (p2-1) + p3;
961 Q = partitions3 * partitions2 * (q1 - 1) + partitions3 * (q2-1) + q3;
962
963 partmap[P] = Q;
964 }
965 }
966 }
967 for(i=1;i<=noelements;i++)
968 inpart[i] = partmap[inpart[i]];
969 free_Ivector(partmap,1,vpartitions);
970 }
971
972 skippart:
973
974 for(i=1;i<=partitions;i++)
975 nopart[i] = 0;
976 for(i=1;i<=noelements;i++)
977 nopart[inpart[i]] += 1;
978
979 minpart = maxpart = nopart[1];
980 for(i=1;i<=partitions;i++) {
981 minpart = MIN( nopart[i], minpart );
982 maxpart = MAX( nopart[i], maxpart );
983 }
984
985 free_Ivector(nopart,1,vpartitions);
986 PartitionNodesByElements(data,info);
987
988 if(info) printf("Successfully made a partitioning with %d to %d elements.\n",minpart,maxpart);
989
990 return(0);
991 }
992
993
PartitionSimpleElementsNonRecursive(struct FemType * data,int dimpart[],int dimper[],int info)994 int PartitionSimpleElementsNonRecursive(struct FemType *data,int dimpart[],int dimper[],
995 int info)
996 /* This is similar to the previous except the partitioning is not recursive.
997 This kind of partitioning is optimal for uniform grids with blockwise structure.
998 An example would be L-type geometry where 3/4 of potential partitions would be active. */
999 {
1000 int i,j,k,minpart,maxpart;
1001 int noknots, noelements,nonodes,periodic;
1002 int partitions1,partitions2,partitions3,partitions;
1003 int IndX,IndY,IndZ;
1004 int *nopart,*inpart;
1005 Real x,y,z,MaxX,MinX,MaxY,MinY,MaxZ,MinZ;
1006
1007 partitions1 = dimpart[0];
1008 partitions2 = dimpart[1];
1009 partitions3 = dimpart[2];
1010 if(data->dim < 3) partitions3 = 1;
1011 partitions = partitions1 * partitions2 * partitions3;
1012
1013 if(partitions1 < 2 && partitions2 < 2 && partitions3 < 2) {
1014 printf("Some of the divisions must be larger than one: %d %d %d\n",
1015 partitions1, partitions2, partitions3 );
1016 bigerror("Partitioning not performed");
1017 }
1018
1019 if(!data->partitionexist) {
1020 data->partitionexist = TRUE;
1021 data->elempart = Ivector(1,data->noelements);
1022 data->nodepart = Ivector(1,data->noknots);
1023 data->nopartitions = partitions;
1024 }
1025 inpart = data->elempart;
1026 noelements = data->noelements;
1027 noknots = data->noknots;
1028
1029 periodic = dimper[0] || dimper[1] || dimper[2];
1030 if(periodic) {
1031 printf("Implement periodicity for this partitioning routine!\n");
1032 }
1033
1034 if(info) {
1035 printf("Making a simple partitioning for %d elements in %d-dimensions.\n",
1036 noelements,data->dim);
1037 printf("There can be at maximum %d partitions\n",partitions);
1038 }
1039
1040 nopart = Ivector(1,partitions);
1041 for(i=1;i<=partitions;i++)
1042 nopart[i] = 0;
1043
1044 z = 0.0;
1045 IndZ = 1;
1046 for(i=1;i<=noelements;i++)
1047 inpart[i] = 0;
1048
1049 MaxX = MinX = data->x[1];
1050 MaxY = MinY = data->y[1];
1051 MaxZ = MinZ = data->z[1];
1052
1053 for(i=1;i<=noknots;i++) {
1054 x = data->x[i];
1055 y = data->y[i];
1056 z = data->z[i];
1057
1058 MaxX = MAX( MaxX, x);
1059 MinX = MIN( MinX, x);
1060 MaxY = MAX( MaxY, y);
1061 MinY = MIN( MinY, y);
1062 MaxZ = MAX( MaxZ, z);
1063 MinZ = MIN( MinZ, z);
1064 }
1065 if( info ) {
1066 printf("Range in x-direction: %12.5le %12.5le\n",MinX,MaxX);
1067 printf("Range in y-direction: %12.5le %12.5le\n",MinY,MaxY);
1068 printf("Range in z-direction: %12.5le %12.5le\n",MinZ,MaxZ);
1069 }
1070
1071 for(j=1;j<=noelements;j++) {
1072 nonodes = data->elementtypes[j]%100;
1073 x = y = z = 0.0;
1074 for(i=0;i<nonodes;i++) {
1075 k = data->topology[j][i];
1076 x += data->x[k];
1077 y += data->y[k];
1078 z += data->z[k];
1079 }
1080 x = x / nonodes;
1081 y = y / nonodes;
1082 z = z / nonodes;
1083
1084 IndX = ceil( partitions1 * ( x - MinX ) / ( MaxX - MinX ) );
1085 IndY = ceil( partitions2 * ( y - MinY ) / ( MaxY - MinY ) );
1086 if( data->dim == 3 ) {
1087 IndZ = ceil( partitions3 * ( z - MinZ ) / ( MaxZ - MinZ ) );
1088 }
1089 k = (IndZ-1) * partitions1 * partitions2 +
1090 (IndY-1) * partitions1 + IndX;
1091
1092 inpart[j] = k;
1093 nopart[k] += 1;
1094 }
1095
1096 maxpart = 0;
1097 minpart = noelements;
1098 j = 0;
1099 if( info ) printf("Renumbering the partitions (showing only 64):\n");
1100 for(i=1;i<=partitions;i++) {
1101 k = nopart[i];
1102 if( k ) {
1103 maxpart = MAX( k, maxpart );
1104 minpart = MIN( k, minpart );
1105 j += 1;
1106 nopart[i] = j;
1107 if( info && j<=64) printf("%d -> %d\n",i,j);
1108 }
1109 }
1110 if(info) printf("There are %d active partitions out of %d possible ones\n",j,partitions);
1111
1112 data->nopartitions = j;
1113
1114 for(i=1;i<=noelements;i++) {
1115 j = inpart[i];
1116 inpart[i] = nopart[j];
1117 }
1118
1119 free_Ivector(nopart,1,partitions);
1120
1121 PartitionNodesByElements(data,info);
1122
1123 if(info) printf("Successfully made a partitioning with %d to %d elements.\n",minpart,maxpart);
1124
1125 return(0);
1126 }
1127
1128
1129 #define MAXCATEGORY 500
PartitionSimpleElementsRotational(struct FemType * data,int dimpart[],int dimper[],int info)1130 int PartitionSimpleElementsRotational(struct FemType *data,int dimpart[],int dimper[],
1131 int info) {
1132 int i,j,k,minpart,maxpart,dim,hit;
1133 int noknots, noelements,nonodes,periodic;
1134 int partitions1,partitions2,partitions3,partitions;
1135 int IndR,IndF,IndZ,connect;
1136 int *nopart,*inpart,*cumf,*cumr,*cumz = NULL,*elemconnect = NULL;
1137 Real x,y,z,r,f,MaxR,MinR,MaxF,MinF,MaxZ,MinZ;
1138
1139 dim = data->dim;
1140 partitions1 = dimpart[0];
1141 partitions2 = dimpart[1];
1142 partitions3 = dimpart[2];
1143 if(dim < 3) partitions3 = 1;
1144 partitions = partitions1 * partitions2 * partitions3;
1145
1146 /* See if there are connected elements */
1147 connect = FALSE;
1148 if(data->nodeconnectexist) {
1149 SetConnectedElements(data,info);
1150 connect = TRUE;
1151 elemconnect = data->elemconnect;
1152 partitions = partitions + partitions3;
1153 }
1154
1155 if(partitions1 < 2 && partitions2 < 2 && partitions3 < 2) {
1156 printf("Some of the divisions must be larger than one: %d %d %d\n",
1157 partitions1, partitions2, partitions3 );
1158 bigerror("Partitioning not performed");
1159 }
1160
1161 if(!data->partitionexist) {
1162 data->partitionexist = TRUE;
1163 data->elempart = Ivector(1,data->noelements);
1164 data->nodepart = Ivector(1,data->noknots);
1165 data->nopartitions = partitions;
1166 }
1167 inpart = data->elempart;
1168 noelements = data->noelements;
1169 noknots = data->noknots;
1170
1171 periodic = dimper[0] || dimper[1] || dimper[2];
1172 if(periodic) {
1173 printf("Implement periodicity for this partitioning routine!\n");
1174 }
1175
1176 if(info) {
1177 printf("Making a simple rotational partitioning for %d elements in %d-dimensions.\n",
1178 noelements,data->dim);
1179 printf("There can be at maximum %d partitions\n",partitions);
1180 }
1181
1182 nopart = Ivector(1,partitions);
1183 for(i=1;i<=partitions;i++)
1184 nopart[i] = 0;
1185
1186 z = 0.0;
1187 IndZ = 1;
1188 for(i=1;i<=noelements;i++)
1189 inpart[i] = 0;
1190
1191 x = data->x[1];
1192 y = data->y[1];
1193 z = data->z[1];
1194
1195 r = sqrt(x*x+y*y);
1196 f = 180 * atan2(y,x)/FM_PI;
1197 if( f < 0.0 ) f = f + 360.0;
1198 MaxR = MinR = r;
1199 MaxF = MinF = f;
1200 MaxZ = MinZ = z;
1201
1202 for(i=1;i<=noknots;i++) {
1203 x = data->x[i];
1204 y = data->y[i];
1205 z = data->z[i];
1206
1207 r = sqrt(x*x+y*y);
1208 f = 180 * atan2(y,x)/FM_PI;
1209 if( f < 0.0 ) f = f + 360.0;
1210
1211 MaxR = MAX( MaxR, r);
1212 MinR = MIN( MinR, r);
1213 MaxF = MAX( MaxF, f);
1214 MinF = MIN( MinF, f);
1215 MaxZ = MAX( MaxZ, z);
1216 MinZ = MIN( MinZ, z);
1217 }
1218
1219 if( info ) {
1220 printf("Range in r-direction: %12.5le %12.5le\n",MinR,MaxR);
1221 printf("Range in f-direction: %12.5le %12.5le\n",MinF,MaxF);
1222 printf("Range in z-direction: %12.5le %12.5le\n",MinZ,MaxZ);
1223 }
1224 if( MaxF - MinF > 180.0 ) {
1225 MaxF = 360.0;
1226 MinF = 0.0;
1227 }
1228
1229 /* Zero is the 1st value so that recursive algos can be used. */
1230 cumr = Ivector(0,MAXCATEGORY);
1231 cumf = Ivector(0,MAXCATEGORY);
1232 if(dim==3) cumz = Ivector(0,MAXCATEGORY);
1233
1234 for(i=0;i<=MAXCATEGORY;i++) {
1235 cumf[i] = cumr[i] = 0;
1236 if( dim == 3) cumz[i] = 0;
1237 }
1238
1239 for(j=1;j<=noelements;j++) {
1240 nonodes = data->elementtypes[j]%100;
1241 x = y = z = 0.0;
1242 for(i=0;i<nonodes;i++) {
1243 k = data->topology[j][i];
1244 x += data->x[k];
1245 y += data->y[k];
1246 z += data->z[k];
1247 }
1248 x = x / nonodes;
1249 y = y / nonodes;
1250 z = z / nonodes;
1251
1252 r = sqrt(x*x+y*y);
1253 f = 180 * atan2(y,x)/FM_PI;
1254 if( f < 0.0 ) f = f + 360.0;
1255
1256 IndR = ceil( MAXCATEGORY * ( r - MinR ) / ( MaxR - MinR ) );
1257 IndF = ceil( MAXCATEGORY * ( f - MinF ) / ( MaxF - MinF ) );
1258 if(dim==3) IndZ = ceil( MAXCATEGORY * ( z - MinZ ) / ( MaxZ - MinZ ) );
1259
1260 if( IndR < 1 || IndR > MAXCATEGORY ) {
1261 printf("IndR out of bounds : %d\n",IndR );
1262 IndR = MIN( MAX( IndR, 1 ), MAXCATEGORY );
1263 }
1264 if( IndF < 1 || IndF > MAXCATEGORY ) {
1265 printf("IndF out of bounds : %d\n",IndF );
1266 IndF = MIN( MAX( IndF, 1 ), MAXCATEGORY );
1267 }
1268 if( dim == 3 ) {
1269 if( IndZ < 1 || IndZ > MAXCATEGORY ) {
1270 printf("IndZ out of bounds : %d\n",IndZ );
1271 IndZ = MIN( MAX( IndZ, 1 ), MAXCATEGORY );
1272 }
1273 cumz[IndZ] += 1;
1274 }
1275 if( connect && elemconnect[j] < 0 ) continue;
1276
1277 cumr[IndR] += 1;
1278 cumf[IndF] += 1;
1279 }
1280
1281 /* Count the cumulative numbers */
1282 for(i=1;i<=MAXCATEGORY;i++) {
1283 cumr[i] = cumr[i] + cumr[i-1];
1284 cumf[i] = cumf[i] + cumf[i-1];
1285 if(dim==3) cumz[i] = cumz[i] + cumz[i-1];
1286 }
1287
1288 j = noelements - data->elemconnectexist;
1289 for(i=1;i<=MAXCATEGORY;i++) {
1290 cumr[i] = ceil( 1.0 * partitions1 * cumr[i] / noelements );
1291 cumf[i] = ceil( 1.0 * partitions2 * cumf[i] / j );
1292 if(dim==3) cumz[i] = ceil( 1.0 * partitions3 * cumz[i] / j );
1293 }
1294
1295
1296 for(j=1;j<=noelements;j++) {
1297 nonodes = data->elementtypes[j]%100;
1298 x = y = z = 0.0;
1299 for(i=0;i<nonodes;i++) {
1300 k = data->topology[j][i];
1301 x += data->x[k];
1302 y += data->y[k];
1303 z += data->z[k];
1304 }
1305 x = x / nonodes;
1306 y = y / nonodes;
1307 z = z / nonodes;
1308
1309 r = sqrt(x*x+y*y);
1310 f = 180 * atan2(y,x)/FM_PI;
1311 if( f < 0.0 ) f = f + 360.0;
1312
1313 IndR = ceil( MAXCATEGORY * ( r - MinR ) / ( MaxR - MinR ) );
1314 IndF = ceil( MAXCATEGORY * ( f - MinF ) / ( MaxF - MinF ) );
1315 if(dim==3) IndZ = ceil( MAXCATEGORY * ( z - MinZ ) / ( MaxZ - MinZ ) );
1316
1317 IndR = MIN( MAX( IndR, 1 ), MAXCATEGORY );
1318 IndF = MIN( MAX( IndF, 1 ), MAXCATEGORY );
1319 if(dim==3) IndZ = MIN( MAX( IndZ, 1 ), MAXCATEGORY );
1320
1321 IndR = cumr[IndR];
1322 IndF = cumf[IndF];
1323 if(dim==3) {
1324 IndZ = cumz[IndZ];
1325 }
1326 else {
1327 IndZ = 1;
1328 }
1329
1330 k = (IndZ-1) * partitions1 * partitions2 +
1331 (IndF-1) * partitions1 + IndR;
1332 if( connect ) {
1333 if( elemconnect[j] < 0 )
1334 k = IndZ;
1335 else
1336 k += 1;
1337 }
1338 inpart[j] = k;
1339 nopart[k] += 1;
1340 }
1341
1342 maxpart = 0;
1343 minpart = noelements;
1344 j = 0;
1345
1346 /* Check whether some partition was not used */
1347 hit = FALSE;
1348 for(i=1;i<=partitions;i++)
1349 if(!nopart[i]) hit = TRUE;
1350
1351 if( hit ) {
1352 if( info ) printf("Renumbering the partitions (showing only 64):\n");
1353 for(i=1;i<=partitions;i++) {
1354 k = nopart[i];
1355 if( k ) {
1356 maxpart = MAX( k, maxpart );
1357 minpart = MIN( k, minpart );
1358 j += 1;
1359 nopart[i] = j;
1360 if( info && i!=j && j<=64) printf("%d -> %d %d\n",i,j,k);
1361 }
1362 }
1363 if(info) printf("There are %d active partitions out of %d possible ones\n",j,partitions);
1364 data->nopartitions = j;
1365
1366 for(i=1;i<=noelements;i++) {
1367 j = inpart[i];
1368 inpart[i] = nopart[j];
1369 }
1370 }
1371 else {
1372 if(info) printf("All possible %d out of %d partitions are active\n",j,partitions);
1373 data->nopartitions = partitions;
1374 }
1375
1376
1377 free_Ivector(nopart,1,partitions);
1378 free_Ivector( cumr,0,MAXCATEGORY);
1379 free_Ivector( cumf,0,MAXCATEGORY);
1380 if(dim==3) free_Ivector( cumz,0,MAXCATEGORY);
1381
1382 PartitionNodesByElements(data,info);
1383
1384 if(info) printf("Successfully made a partitioning with %d to %d elements.\n",minpart,maxpart);
1385
1386 return(0);
1387 }
1388
1389
1390
PartitionConnectedElementsStraight(struct FemType * data,struct BoundaryType * bound,struct ElmergridType * eg,int info)1391 int PartitionConnectedElementsStraight(struct FemType *data,struct BoundaryType *bound,
1392 struct ElmergridType *eg, int info) {
1393 int i,j,k,l,dim,allocated,debug,partz,hit,bctype;
1394 int noknots, noelements,bcelem,bc,maxbcelem;
1395 int IndZ,noconnect,totpartelems,sideelemtype,sidenodes,sidehits,nohits;
1396 int *cumz,*elemconnect,*partelems,*nodeconnect;
1397 int sideind[MAXNODESD2];
1398 Real z,MaxZ,MinZ;
1399
1400
1401 debug = FALSE;
1402
1403 if(info) {
1404 printf("Link the connected set directly to the partition\n");
1405 }
1406
1407 if(!data->nodeconnectexist) {
1408 printf("There are no connected boundary nodes?\n");
1409 return(1);
1410 }
1411 nodeconnect = data->nodeconnect;
1412 noknots = data->noknots;
1413 noelements = data->noelements;
1414 totpartelems = 0;
1415
1416 bcelem = 0;
1417 data->elemconnect = Ivector(1,noelements);
1418 elemconnect = data->elemconnect;
1419 for(i=1;i<=noelements;i++) elemconnect[i] = 0;
1420
1421 for(bc=0;bc<MAXBOUNDARIES;bc++) {
1422 if(bound[bc].created == FALSE) continue;
1423 if(bound[bc].nosides == 0) continue;
1424
1425 for(i=1;i<=bound[bc].nosides;i++) {
1426
1427 GetElementSide(bound[bc].parent[i],bound[bc].side[i],bound[bc].normal[i],
1428 data,sideind,&sideelemtype);
1429
1430 sidenodes = sideelemtype % 100;
1431
1432 hit = FALSE;
1433
1434 for(k=1;k<=eg->connect;k++) {
1435 bctype = eg->connectbounds[k-1];
1436 hit = (bound[bc].types[i] == bctype);
1437 if(hit) break;
1438 }
1439 if(!hit) continue;
1440 bcelem += 1;
1441
1442 k = bound[bc].parent[i];
1443 l = bound[bc].parent2[i];
1444 if(k) {
1445 if(!elemconnect[k]) {
1446 elemconnect[k] = 1;
1447 totpartelems += 1;
1448 }
1449 }
1450 if(l) {
1451 if(!elemconnect[l]) {
1452 elemconnect[l] = 1;
1453 totpartelems += 1;
1454 }
1455 }
1456 }
1457 }
1458
1459 data->elemconnectexist = totpartelems;
1460
1461 if(info) printf("Number of constrained boundary elements = %d\n",bcelem);
1462 if(info) printf("Number of constrained bulk elements = %d\n",totpartelems);
1463 if(info) printf("Successfully made connected elements to a partiotion\n");
1464
1465 return(0);
1466 }
1467
1468
1469
1470
PartitionConnectedElements1D(struct FemType * data,struct BoundaryType * bound,struct ElmergridType * eg,int info)1471 int PartitionConnectedElements1D(struct FemType *data,struct BoundaryType *bound,
1472 struct ElmergridType *eg, int info) {
1473 int i,j,k,l,dim,allocated,debug,partz,partr,parts,hit,bctype;
1474 int noknots, noelements,bcelem,bc,maxbcelem;
1475 int IndZ,noconnect,totpartelems,sideelemtype,sidenodes,sidehits,nohits;
1476 int *cumz,*elemconnect,*partelems,*nodeconnect;
1477 int sideind[MAXNODESD2];
1478 Real val,z,MaxZ,MinZ;
1479
1480
1481 debug = FALSE;
1482
1483 partz = eg->partbcz;
1484 partr = eg->partbcr;
1485
1486 if( partz == 0 && partr == 0) return(0);
1487
1488 parts = MAX( partz, partr );
1489
1490 if(info) {
1491 if( partz )
1492 printf("Making a simple 1D partitioning in z for the connected elements only\n");
1493 else
1494 printf("Making a simple 1D partitioning in r for the connected elements only\n");
1495 }
1496
1497 if(!data->nodeconnectexist) {
1498 printf("There are no connected boundary nodes?\n");
1499 return(1);
1500 }
1501 nodeconnect = data->nodeconnect;
1502
1503 dim = data->dim;
1504 if( dim < 3 ) {
1505 printf("PartitionConnectedElements1D is only applicable in 3D\n");
1506 return(2);
1507 }
1508
1509 noknots = data->noknots;
1510 noelements = data->noelements;
1511 totpartelems = 0;
1512
1513 /* Because we don't want to change the code too much use 'z' for the
1514 coordinate also in the radial case. */
1515 if( partr )
1516 z = sqrt( data->x[1] * data->z[1] + data->y[1] * data->y[1]);
1517 else
1518 z = data->z[1];
1519 MaxZ = MinZ = z;
1520
1521 for(i=1;i<=noknots;i++) {
1522 if( partr )
1523 z = sqrt( data->x[i] * data->x[i] + data->y[i] * data->y[i]);
1524 else
1525 z = data->z[i];
1526 MaxZ = MAX( MaxZ, z);
1527 MinZ = MIN( MinZ, z);
1528 }
1529
1530 if( info ) {
1531 printf("Range in coordinate extent: %12.5le %12.5le\n",MinZ,MaxZ);
1532 }
1533
1534 /* Zero is the 1st value so that recursive algos can be used. */
1535 cumz = Ivector(0,MAXCATEGORY);
1536 for(i=0;i<=MAXCATEGORY;i++)
1537 cumz[i] = 0;
1538
1539 allocated = FALSE;
1540 bcelem = 0;
1541
1542 omstart:
1543
1544 for(bc=0;bc<MAXBOUNDARIES;bc++) {
1545 if(bound[bc].created == FALSE) continue;
1546 if(bound[bc].nosides == 0) continue;
1547
1548 for(i=1;i<=bound[bc].nosides;i++) {
1549
1550 GetElementSide(bound[bc].parent[i],bound[bc].side[i],bound[bc].normal[i],
1551 data,sideind,&sideelemtype);
1552
1553 sidenodes = sideelemtype % 100;
1554 hit = FALSE;
1555
1556 for(k=1;k<=eg->connect;k++) {
1557 bctype = eg->connectbounds[k-1];
1558 if( bctype > 0 ) {
1559 if(bound[bc].types[i] == bctype) hit = TRUE;
1560 }
1561 else if( bctype == -1 ) {
1562 if( bound[bc].parent[i] ) hit = TRUE;
1563 }
1564 else if( bctype == -2 ) {
1565 if( bound[bc].parent[i] && bound[bc].parent2[i] ) hit = TRUE;
1566 }
1567 else if( bctype == -3 ) {
1568 if( bound[bc].parent[i] && !bound[bc].parent2[i] ) hit = TRUE;
1569 }
1570 if(hit) break;
1571 }
1572 if(!hit) continue;
1573
1574 z = 0.0;
1575 for(j=0;j<sidenodes;j++) {
1576 k = sideind[j];
1577 if( partr )
1578 val = sqrt( data->x[k]*data->x[k] + data->y[k]*data->y[k]);
1579 else
1580 val = data->z[k];
1581 z += val;
1582 }
1583
1584 z = z / sidenodes;
1585 IndZ = ceil( MAXCATEGORY * ( z - MinZ ) / ( MaxZ - MinZ ) );
1586
1587 /* To be on the safe side */
1588 IndZ = MIN( MAX( IndZ, 1 ), MAXCATEGORY );
1589
1590
1591 if(allocated) {
1592 IndZ = cumz[IndZ];
1593 if(0) partelems[IndZ] += 1;
1594
1595 /* Also set the connected elements if available */
1596 k = bound[bc].parent[i];
1597 l = bound[bc].parent2[i];
1598
1599 if(k) {
1600 elemconnect[k] = IndZ;
1601 totpartelems += 1;
1602 }
1603 if(l) {
1604 elemconnect[l] = IndZ;
1605 totpartelems += 1;
1606 }
1607 }
1608 else {
1609 bcelem += 1;
1610 cumz[IndZ] += 1;
1611 }
1612 }
1613 }
1614
1615 if( !allocated ) {
1616 maxbcelem = bcelem;
1617 if( debug ) {
1618 printf("Number of constrained boundary elements = %d\n",bcelem);
1619 printf("Differential categories (showing only 20 active ones from %d)\n",MAXCATEGORY);
1620 k = 0;
1621 for(j=0;j<=MAXCATEGORY;j++) {
1622 if( cumz[j] > 0 ) {
1623 k++;
1624 printf("%d : %d\n",j,cumz[j]);
1625 if( k == 20 ) break;
1626 }
1627 }
1628 }
1629
1630 /* Count the cumulative numbers and map them to number of partitions */
1631 for(i=1;i<=MAXCATEGORY;i++)
1632 cumz[i] = cumz[i] + cumz[i-1];
1633
1634 if( debug ) {
1635 printf("Cumulative categories\n");
1636 for(j=0;j<=MAXCATEGORY;j++) {
1637 printf("%d : %d\n",j,cumz[j]);
1638 }
1639 }
1640
1641 noconnect = bcelem;
1642 for(i=1;i<=MAXCATEGORY;i++)
1643 cumz[i] = ceil( 1.0 * parts * cumz[i] / noconnect );
1644
1645 if( debug ) {
1646 printf("Partition categories\n");
1647 for(j=0;j<=MAXCATEGORY;j++) {
1648 printf("%d : %d\n",j,cumz[j]);
1649 }
1650 }
1651
1652 data->elemconnect = Ivector(1,noelements);
1653 elemconnect = data->elemconnect;
1654
1655 for(i=1;i<=noelements;i++) elemconnect[i] = 0;
1656
1657 allocated = TRUE;
1658 if(info) printf("Allocated and now setting the entries\n");
1659 goto omstart;
1660 }
1661
1662 data->elemconnectexist = totpartelems;
1663 free_Ivector( cumz,0,MAXCATEGORY);
1664
1665 if(info) printf("Successfully made a partitioning 1D with %d elements\n",totpartelems);
1666
1667 return(0);
1668 }
1669
1670
1671 #if USE_METIS
PartitionConnectedElementsMetis(struct FemType * data,struct BoundaryType * bound,int nparts,int metisopt,int info)1672 int PartitionConnectedElementsMetis(struct FemType *data,struct BoundaryType *bound,
1673 int nparts,int metisopt,int info) {
1674 int i,j,k,l,n,m,dim;
1675 int noknots,noelements,sidenodes,ind,nohits,totpartelems;
1676 int dualmaxcon,invmaxcon,totcon,step,bc,bcelem,bcelem2,hit,set;
1677 int noconnect,minpartelems,maxpartelems,sideelemtype,con,maxbcelem;
1678 int *elemconnect,*partelems,*nodeperm,*neededby;
1679 int *bcdualgraph[MAXCONNECTIONS],*bcinvtopo[MAXCONNECTIONS];
1680 int sideind[MAXNODESD1];
1681
1682 int nn;
1683 idx_t options[METIS_NOPTIONS];
1684 int *xadj,*adjncy,*vwgt,*adjwgt,wgtflag,*npart;
1685 int numflag,edgecut,ncon;
1686 int *nodepart;
1687
1688
1689 if(info) {
1690 printf("Making a Metis partitioning for the connected BC elements only\n");
1691 }
1692
1693 if(!data->nodeconnectexist) {
1694 printf("There are no connected elements\n");
1695 return(1);
1696 }
1697
1698 dim = data->dim;
1699 noknots = data->noknots;
1700 noelements = data->noelements;
1701
1702 /* Count the connected nodes and set the permutation */
1703 noconnect = 0;
1704 for(i=1;i<=noknots;i++)
1705 if( data->nodeconnect[i] ) noconnect++;
1706 if( noconnect == 0 ) {
1707 printf("There are really no connected nodes?\n");
1708 return(2);
1709 }
1710 if(info) printf("Number of connected nodes is %d\n",noconnect);
1711
1712 /* Permute the nodes so that they only include the connected nodes. */
1713 nodeperm = Ivector(1,noknots);
1714 for(i=1;i<=noknots;i++) nodeperm[i] = 0;
1715 noconnect = 0;
1716 for(i=1;i<=noknots;i++)
1717 if( data->nodeconnect[i] ) {
1718 noconnect++;
1719 nodeperm[i] = noconnect;
1720 }
1721 if(info) printf("Permuted the connected nodes\n");
1722
1723 /* Cycle over the boundary elements
1724 1) First cycle create a table showing the elements sharing a node
1725 2) Second cycle create a table showing the element-to-element connections */
1726
1727 invmaxcon = 0;
1728 dualmaxcon = 0;
1729 totcon = 0;
1730
1731 neededby = Ivector(1,noconnect);
1732 for(i=1;i<=noconnect;i++) neededby[i] = 0;
1733
1734 for(step=1;step<=2;step++) {
1735
1736 bcelem = 0;
1737 for(bc=0;bc<MAXBOUNDARIES;bc++) {
1738 if(bound[bc].created == FALSE) continue;
1739 if(bound[bc].nosides == 0) continue;
1740
1741 for(i=1;i<=bound[bc].nosides;i++) {
1742
1743 GetBoundaryElement(i,&bound[bc],data,sideind,&sideelemtype);
1744 /* GetElementSide(bound[bc].parent[i],bound[bc].side[i],bound[bc].normal[i],
1745 data,sideind,&sideelemtype); */
1746 sidenodes = sideelemtype % 100;
1747 nohits = 0;
1748
1749 for(j=0;j<sidenodes;j++)
1750 if( nodeperm[sideind[j]] ) nohits++;
1751 if( nohits < sidenodes ) continue;
1752 bcelem++;
1753
1754 for(j=0;j<sidenodes;j++) {
1755 ind = nodeperm[sideind[j]];
1756
1757 /* Create table that shows all the elements sharing a node. */
1758 if( step == 1 ) {
1759 neededby[ind] += 1;
1760 k = neededby[ind];
1761
1762 if(k > invmaxcon) {
1763 invmaxcon++;
1764 bcinvtopo[invmaxcon] = Ivector(1,noconnect);
1765 if(0) printf("allocating invtopo %d %d\n",invmaxcon,noconnect);
1766 for(m=1;m<=noconnect;m++)
1767 bcinvtopo[invmaxcon][m] = 0;
1768 }
1769 bcinvtopo[k][ind] = bcelem;
1770 }
1771 else if( step == 2 ) {
1772 /* Now create a table that shows how elements are connected to other elements,
1773 this is the dual graph. */
1774 for(k=1;k<=invmaxcon;k++) {
1775 bcelem2 = bcinvtopo[k][ind];
1776
1777 if( bcelem2 == 0 ) break;
1778 if( bcelem2 == bcelem ) continue;
1779
1780 hit = FALSE;
1781 for(l=1;l<=dualmaxcon;l++) {
1782 if(bcdualgraph[l][bcelem] == bcelem2) hit = TRUE;
1783 if(bcdualgraph[l][bcelem] == 0) break;
1784 }
1785 if(!hit) {
1786 if(l > dualmaxcon) {
1787 dualmaxcon++;
1788 if( l >= MAXCONNECTIONS ) {
1789 printf("Number of connections %d vs. static limit %d\n",l,MAXCONNECTIONS-1);
1790 bigerror("Maximum of connections in dual graph larger than the static limit!");
1791 }
1792 if(0) printf("allocating dual topo %d %d\n",dualmaxcon,maxbcelem);
1793 bcdualgraph[dualmaxcon] = Ivector(1,maxbcelem);
1794 for(m=1;m<=maxbcelem;m++)
1795 bcdualgraph[dualmaxcon][m] = 0;
1796 }
1797 bcdualgraph[l][bcelem] = bcelem2;
1798 totcon++;
1799 }
1800 }
1801 }
1802 }
1803 }
1804 }
1805 maxbcelem = bcelem;
1806 }
1807
1808 if(info) {
1809 printf("There are %d connected boundary elements.\n",maxbcelem);
1810 printf("There are %d connections altogether in the dual graph.\n",totcon);
1811 }
1812
1813 /* Create the sparse matrix format graph for Metis */
1814 xadj = Ivector(0,maxbcelem);
1815 adjncy = Ivector(0,totcon-1);
1816 for(i=0;i<totcon;i++)
1817 adjncy[i] = 0;
1818
1819 totcon = 0;
1820 for(i=1;i<=maxbcelem;i++) {
1821 xadj[i-1] = totcon;
1822 for(j=1;j<=dualmaxcon;j++) {
1823 con = bcdualgraph[j][i];
1824 if(con) {
1825 adjncy[totcon] = con-1;
1826 totcon++;
1827 if(0) printf("con: %d %d %d\n",i,totcon,con);
1828 }
1829 }
1830 }
1831 xadj[maxbcelem] = totcon;
1832
1833 /* Deallocate the temporal graphs */
1834 free_Ivector( neededby, 1, noconnect );
1835 for(i=1;i<=invmaxcon;i++)
1836 free_Ivector(bcinvtopo[i],1,noconnect);
1837 for(i=1;i<=dualmaxcon;i++)
1838 free_Ivector(bcdualgraph[i],1,maxbcelem);
1839
1840 /* Parameters for Metis */
1841 numflag = 0;
1842 nn = maxbcelem;
1843 npart = Ivector(0,nn-1);
1844 wgtflag = 0;
1845
1846 METIS_SetDefaultOptions(options);
1847 options[METIS_OPTION_CTYPE] = METIS_CTYPE_SHEM;
1848 options[METIS_OPTION_IPTYPE] = METIS_IPTYPE_GROW;
1849 options[METIS_OPTION_RTYPE] = METIS_RTYPE_GREEDY;
1850 if(metisopt == 4 ) options[METIS_OPTION_MINCONN] = 1;
1851 options[METIS_OPTION_DBGLVL] = 3;
1852 /* if( eg->metisseed ) options[METIS_OPTION_SEED] = eg->metisseed; */
1853
1854 /* Optional weights */
1855 vwgt = NULL;
1856 adjwgt = NULL;
1857
1858 if(metisopt == 2) {
1859 if(info) printf("Starting graph partitioning METIS_PartGraphRecursive.\n");
1860 METIS_PartGraphRecursive(&nn,&ncon,xadj,adjncy,vwgt,&wgtflag,adjwgt,
1861 &nparts,NULL,NULL,options,&edgecut,npart);
1862 }
1863 else if(metisopt == 3 || metisopt == 4) {
1864 if(info) printf("Starting graph partitioning METIS_PartGraphKway.\n");
1865 ncon = 0;
1866 wgtflag = 0;
1867 METIS_PartGraphKway(&nn,&ncon,xadj,adjncy,vwgt,&wgtflag,adjwgt,
1868 &nparts,NULL,NULL,options,&edgecut,npart);
1869 }
1870 else {
1871 printf("Unknown Metis option %d\n",metisopt);
1872 }
1873 if(0) printf("Finished Metis routine for boundary partitioning\n");
1874
1875
1876 free_Ivector(xadj, 0,maxbcelem);
1877 free_Ivector(adjncy, 0,totcon-1);
1878
1879
1880 data->elemconnect = Ivector(1,noelements);
1881 elemconnect = data->elemconnect;
1882 for(i=1;i<=noelements;i++) elemconnect[i] = 0;
1883
1884 /* Set the bulk element partitions by following the boundary elements.
1885 Also set the nodal partitions to apply majority rule later on. */
1886 bcelem = 0;
1887 totpartelems = 0;
1888 for(bc=0;bc<MAXBOUNDARIES;bc++) {
1889 if(bound[bc].created == FALSE) continue;
1890 if(bound[bc].nosides == 0) continue;
1891
1892 for(i=1;i<=bound[bc].nosides;i++) {
1893
1894 GetBoundaryElement(i,&bound[bc],data,sideind,&sideelemtype);
1895 /* GetElementSide(bound[bc].parent[i],bound[bc].side[i],bound[bc].normal[i],
1896 data,sideind,&sideelemtype); */
1897 sidenodes = sideelemtype % 100;
1898 nohits = 0;
1899 for(j=0;j<sidenodes;j++)
1900 if( nodeperm[sideind[j]] ) nohits++;
1901 if( nohits < sidenodes ) continue;
1902
1903 /* Always set the connected element to a negative value as this information is used
1904 in a dirty way later on. */
1905 k = bound[bc].parent[i];
1906 l = bound[bc].parent2[i];
1907
1908 if(0) printf("kl = %d %d %d %d %d\n",k,l,bcelem,npart[bcelem]+1,elemconnect[k]);
1909
1910 /* Have the additition here since npart has indexes [0,nn-1]. */
1911 if(k) {
1912 elemconnect[k] = npart[bcelem]+1;
1913 totpartelems += 1;
1914 }
1915 if(l) {
1916 elemconnect[l] = npart[bcelem]+1;
1917 totpartelems += 1;
1918 }
1919 bcelem++;
1920 }
1921 }
1922
1923 free_Ivector( nodeperm, 1, noknots );
1924 free_Ivector( npart, 0, nn-1 );
1925
1926 data->elemconnectexist = totpartelems;
1927 if(info) printf("Successfully created boundary partitioning for %d bulk elements\n",totpartelems);
1928
1929 return(0);
1930 }
1931 #endif
1932
1933
ExtendBoundaryPartitioning(struct FemType * data,struct BoundaryType * bound,int elemlayers,int info)1934 int ExtendBoundaryPartitioning(struct FemType *data,struct BoundaryType *bound,
1935 int elemlayers,int info)
1936 {
1937 int i,j,k,l,m,n,nonodes,noknots,noelements,totpartelems,nparts,set;
1938 int minpartelems,maxpartelems,refcount,refpart,part;
1939 int *partelems,*elemconnect;
1940 int *invrow,*invcol;
1941
1942 printf("Extending boundary partitioning by majority rule\n");
1943
1944 CreateInverseTopology(data, info);
1945 invrow = data->invtopo.rows;
1946 invcol = data->invtopo.cols;
1947
1948 noelements = data->noelements;
1949 noknots = data->noknots;
1950
1951 if( !data->elemconnectexist ) {
1952 printf("Initial partitioning does not exist!\n");
1953 return(1);
1954 }
1955 elemconnect = data->elemconnect;
1956
1957 nparts = 0;
1958 for(i=1;i<=data->noelements;i++)
1959 nparts = MAX( nparts, elemconnect[i] );
1960
1961 partelems = Ivector(1,nparts);
1962 for(i=1;i<=nparts;i++)
1963 partelems[i] = 0;
1964
1965 totpartelems = 0;
1966 for(i=1;i<=noelements;i++) {
1967 j = elemconnect[i];
1968 if(j) {
1969 partelems[j] += 1;
1970 totpartelems +=1;
1971 }
1972 }
1973
1974 if(info) printf("Initial partitioning with %d elements:\n",totpartelems);
1975 minpartelems = maxpartelems = partelems[1];
1976 for(j=1;j<=nparts;j++) {
1977 minpartelems = MIN( minpartelems, partelems[j] );
1978 maxpartelems = MAX( maxpartelems, partelems[j] );
1979 }
1980 printf("Initial partitioning has %d to %d elements in partition\n",minpartelems,maxpartelems);
1981
1982
1983 /* Counter to find the deminating element */
1984 for(i=1;i<=nparts;i++)
1985 partelems[i] = 0;
1986
1987 /* Apply majority rule to set the remaining bulk elements that
1988 share nodes on the connected boundary. */
1989 for(m=1;m<=elemlayers;m++) {
1990 n = 0;
1991 for(j=1;j<=noelements;j++) {
1992
1993 /* This element is already set */
1994 if( elemconnect[j] > 0 ) continue;
1995
1996 set = FALSE;
1997 nonodes = data->elementtypes[j] % 100;
1998 for(i=0;i<nonodes;i++) {
1999 k = data->topology[j][i];
2000 for(l=invrow[k-1];l<invrow[k];l++) {
2001 part = elemconnect[invcol[l]+1];
2002 if( part > 0 ) {
2003 set = TRUE;
2004 partelems[part] += 1;
2005 }
2006 }
2007 }
2008 /* No hits, cannot set */
2009 if(!set) continue;
2010
2011 refcount = 0;
2012 for(i=0;i<nonodes;i++) {
2013 k = data->topology[j][i];
2014 for(l=invrow[k-1];l<invrow[k];l++) {
2015 part = elemconnect[invcol[l]+1];
2016 if(part > 0 ) {
2017 if( partelems[part] > refcount ) {
2018 refcount = partelems[part];
2019 refpart = part;
2020 }
2021 partelems[part] = 0;
2022 }
2023 }
2024 }
2025
2026 if(0) printf("setting: %d %d %d\n",j,refpart,refcount);
2027
2028 /* Set negative entry first so that this is not used within this loop */
2029 elemconnect[j] = -refpart;
2030 n++;
2031 }
2032
2033 /* Revert to positive indexes at the end */
2034 for(j=1;j<=noelements;j++)
2035 elemconnect[j] = abs( elemconnect[j] );
2036
2037 if(info) printf("Setting partition for %d elements in layer %d\n",n,m);
2038 if( n == 0 ) {
2039 printf("No elements set in the layer, breaking loop %d!\n",m);
2040 break;
2041 }
2042 }
2043
2044
2045 /* Compute the division to partitions */
2046 for(i=1;i<=nparts;i++)
2047 partelems[i] = 0;
2048
2049 totpartelems = 0;
2050 for(i=1;i<=noelements;i++) {
2051 j = elemconnect[i];
2052 if(j) {
2053 partelems[j] += 1;
2054 totpartelems += 1;
2055 }
2056 }
2057
2058 if(info) printf("Extended partitioning of %d elements:\n",totpartelems);
2059 minpartelems = maxpartelems = partelems[1];
2060 for(j=1;j<=nparts;j++) {
2061 minpartelems = MIN( minpartelems, partelems[j] );
2062 maxpartelems = MAX( maxpartelems, partelems[j] );
2063 if(0) printf("Part: %d %d\n",j,partelems[j]);
2064 }
2065 if(info) printf("Extended partitioning has %d to %d elements in partition\n",minpartelems,maxpartelems);
2066
2067 data->elemconnectexist = totpartelems;
2068 free_Ivector( partelems, 1, nparts );
2069
2070 if(info) printf("Successfully extended the partitioning by %d layers\n",elemlayers);
2071
2072 return(0);
2073 }
2074
2075
2076
2077
2078
PartitionSimpleNodes(struct FemType * data,int dimpart[],int dimper[],int partorder,Real corder[],Real parttol,int info)2079 int PartitionSimpleNodes(struct FemType *data,int dimpart[],int dimper[],
2080 int partorder, Real corder[],Real parttol,int info)
2081 /* Do simple partitioning going for the nodes. Then using the the majority rule
2082 define the elemental partitioning. This is suboptimal for Elmer because the
2083 elemental ordering is primary in Elmer. */
2084 {
2085 int i,j,k,k0,kprev,l,ind,minpart,maxpart;
2086 int noknots, noelements,elemsinpart,periodic;
2087 int partitions1,partitions2,partitions3,partitions;
2088 int vpartitions1,vpartitions2,vpartitions3,vpartitions,hit;
2089 int *indx,*nopart,*nodepart;
2090 Real *arrange;
2091 Real x,y,z,cx,cy,cz;
2092
2093 partitions1 = dimpart[0];
2094 partitions2 = dimpart[1];
2095 partitions3 = dimpart[2];
2096 if(data->dim < 3) partitions3 = 1;
2097 partitions = partitions1 * partitions2 * partitions3;
2098
2099 if(partitions1 < 2 && partitions2 < 2 && partitions3 < 2) {
2100 printf("Some of the divisions must be larger than one: %d %d %d\n",
2101 partitions1, partitions2, partitions3 );
2102 bigerror("Partitioning not performed");
2103 }
2104
2105 if(partitions >= data->noelements) {
2106 printf("There must be fever partitions than elements (%d vs %d)!\n",
2107 partitions,data->noelements);
2108 bigerror("Partitioning not performed");
2109 }
2110
2111 if(!data->partitionexist) {
2112 data->partitionexist = TRUE;
2113 data->elempart = Ivector(1,data->noelements);
2114 data->nodepart = Ivector(1,data->noknots);
2115 data->nopartitions = partitions;
2116 }
2117 nodepart = data->nodepart;
2118
2119 vpartitions1 = partitions1;
2120 vpartitions2 = partitions2;
2121 vpartitions3 = partitions3;
2122 periodic = dimper[0] || dimper[1] || dimper[2];
2123 if(periodic) {
2124 if(dimper[0] && partitions1 > 1) vpartitions1 *= 2;
2125 if(dimper[1] && partitions2 > 1) vpartitions2 *= 2;
2126 if(dimper[2] && partitions3 > 1) vpartitions3 *= 2;
2127 }
2128 vpartitions = vpartitions1 * vpartitions2 * vpartitions3;
2129
2130 nopart = Ivector(1,vpartitions);
2131 noelements = data->noelements;
2132 noknots = data->noknots;
2133
2134 if(info) printf("Making a simple partitioning for %d nodes in %d-dimensions.\n",
2135 noknots,data->dim);
2136
2137 arrange = Rvector(1,noknots);
2138 indx = Ivector(1,noknots);
2139
2140 if(partorder) {
2141 cx = corder[0];
2142 cy = corder[1];
2143 cz = corder[2];
2144 }
2145 else {
2146 cx = 1.0;
2147 cy = 0.0001;
2148 cz = cy*cy;
2149 }
2150
2151 z = 0.0;
2152
2153 for(i=1;i<=noknots;i++)
2154 nodepart[i] = 1;
2155
2156 if(vpartitions1 > 1) {
2157 if(info) printf("Ordering 1st direction with (%.3lg*x + %.3lg*y + %.3lg*z)\n",cx,cy,cz);
2158 for(j=1;j<=noknots;j++) {
2159 x = data->x[j];
2160 y = data->y[j];
2161 z = data->z[j];
2162 arrange[j] = cx*x + cy*y + cz*z;
2163 }
2164 SortIndex(noknots,arrange,indx);
2165
2166 for(i=1;i<=noknots;i++) {
2167 ind = indx[i];
2168 k = (i*vpartitions1-1)/noknots+1;
2169 nodepart[ind] = k;
2170 }
2171 }
2172
2173 /* Partition in the 2nd direction taking into account the 1st direction */
2174 if(vpartitions2 > 1) {
2175 if(info) printf("Ordering in the 2nd direction.\n");
2176 for(j=1;j<=noknots;j++) {
2177 x = data->x[j];
2178 y = data->y[j];
2179 z = data->z[j];
2180 arrange[j] = -cy*x + cx*y + cz*z;
2181 }
2182 SortIndex(noknots,arrange,indx);
2183
2184 for(i=1;i<=vpartitions;i++)
2185 nopart[i] = 0;
2186
2187 elemsinpart = noknots / (vpartitions1*vpartitions2);
2188 j = 1;
2189 for(i=1;i<=noknots;i++) {
2190 ind = indx[i];
2191 k0 = (nodepart[ind]-1) * vpartitions2;
2192 for(l=1;l<=vpartitions2;l++) {
2193 hit = FALSE;
2194 if( j < vpartitions ) {
2195 if( nopart[k0+j] >= elemsinpart ) {
2196 j += 1;
2197 hit = TRUE;
2198 }
2199 }
2200 if( j > 1 ) {
2201 if( nopart[k0+j-1] < elemsinpart ) {
2202 j -= 1;
2203 hit = TRUE;
2204 }
2205 }
2206 if( !hit ) break;
2207 }
2208 k = k0 + j;
2209 nopart[k] += 1;
2210 nodepart[ind] = k;
2211 }
2212 }
2213
2214 /* Partition in the 3rd direction taking into account the 1st and 2nd direction */
2215 if(vpartitions3 > 1) {
2216 if(info) printf("Ordering in the 3rd direction.\n");
2217 for(j=1;j<=noknots;j++) {
2218 x = data->x[j];
2219 y = data->y[j];
2220 z = data->z[j];
2221 arrange[j] = -cz*x - cy*y + cx*z;
2222 }
2223 SortIndex(noknots,arrange,indx);
2224
2225 for(i=1;i<=vpartitions;i++)
2226 nopart[i] = 0;
2227
2228 elemsinpart = noknots / (vpartitions1*vpartitions2*vpartitions3);
2229 j = 1;
2230 for(i=1;i<=noknots;i++) {
2231 ind = indx[i];
2232 k0 = (nodepart[ind]-1)*vpartitions3;
2233
2234 for(l=1;l<=vpartitions;l++) {
2235 hit = FALSE;
2236 if( j < vpartitions3 ) {
2237 if( nopart[k0+j] >= elemsinpart ) {
2238 j += 1;
2239 hit = TRUE;
2240 }
2241 }
2242 if( j > 1 ) {
2243 if( nopart[k0+j-1] < elemsinpart ) {
2244 j -= 1;
2245 hit = TRUE;
2246 }
2247 }
2248 if( !hit ) break;
2249 }
2250 k = k0 + j;
2251 nopart[k] += 1;
2252 nodepart[ind] = k;
2253 }
2254 }
2255
2256
2257 /* This piece of code may be used to assign "strides" back to same partition. */
2258 if(parttol > 1.0e-20 ) {
2259 if(0) {
2260 printf("Original partition counts:\n");
2261 for(i=1;i<=partitions;i++)
2262 printf("nopart[%d] = %d\n",i,nopart[i]);
2263 }
2264
2265 for(j=1;j<=noknots;j++) {
2266 x = data->x[j];
2267 y = data->y[j];
2268 z = data->z[j];
2269 arrange[j] = cx*x + cy*y + cz*z;
2270 }
2271 SortIndex(noknots,arrange,indx);
2272
2273 kprev = nodepart[indx[1]];
2274 l = 0;
2275 for(i=2;i<=noknots;i++) {
2276 ind = indx[i];
2277 k = nodepart[ind];
2278 if(k != kprev ) {
2279 if(fabs(arrange[i]-arrange[i-1]) < parttol ) {
2280 nopart[k] -= 1;
2281 nopart[kprev] += 1;
2282 nodepart[ind] = kprev;
2283 k = kprev;
2284 if(0) printf("arrange %d %d %d %d %.6lg %.6lg\n",i,ind,k,kprev,arrange[i],arrange[i-1]);
2285 l++;
2286 }
2287 }
2288 kprev = k;
2289 }
2290 if(0) {
2291 printf("Updated partition counts:\n");
2292 for(i=1;i<=partitions;i++)
2293 printf("nopart[%d] = %d\n",i,nopart[i]);
2294 }
2295 if(info) printf("Number of partition changes due to tolerance: %d\n",l);
2296 }
2297
2298
2299 /* For periodic systems the number of virtual partitions is larger. Now map the mesh so that the
2300 1st and last partition for each direction will be joined */
2301 if(periodic) {
2302 int *partmap;
2303 int p1,p2,p3,q1,q2,q3;
2304 int P,Q;
2305 p1=p2=p3=1;
2306 partmap = Ivector(1,vpartitions);
2307 for(i=1;i<=vpartitions;i++)
2308 partmap[i] = 0;
2309 for(p1=1;p1<=vpartitions1;p1++) {
2310 q1 = p1;
2311 if(dimper[0] && vpartitions1 > 1) {
2312 if(q1==vpartitions1) q1 = 0;
2313 q1 = q1/2 + 1;
2314 }
2315 for(p2=1;p2<=vpartitions2;p2++) {
2316 q2 = p2;
2317 if(dimper[1] && vpartitions2 > 1) {
2318 if(q2==vpartitions2) q2 = 0;
2319 q2 = q2/2 + 1;
2320 }
2321 for(p3=1;p3<=vpartitions3;p3++) {
2322 q3 = p3;
2323 if(dimper[2] && vpartitions3 > 1) {
2324 if(q3==vpartitions3) q3 = 0;
2325 q3 = q3/2 + 1;
2326 }
2327
2328 P = vpartitions3 * vpartitions2 * (p1 - 1) + vpartitions3 * (p2-1) + p3;
2329 Q = partitions3 * partitions2 * (q1 - 1) + partitions3 * (q2-1) + q3;
2330
2331 partmap[P] = Q;
2332 }
2333 }
2334 }
2335 for(i=1;i<=noknots;i++)
2336 nodepart[i] = partmap[nodepart[i]];
2337 free_Ivector(partmap,1,vpartitions);
2338 }
2339
2340
2341 for(i=1;i<=partitions;i++)
2342 nopart[i] = 0;
2343 for(i=1;i<=noknots;i++)
2344 nopart[nodepart[i]] += 1;
2345
2346 minpart = maxpart = nopart[1];
2347 for(i=1;i<=partitions;i++) {
2348 minpart = MIN( nopart[i], minpart );
2349 maxpart = MAX( nopart[i], maxpart );
2350 }
2351
2352 free_Rvector(arrange,1,noelements);
2353 free_Ivector(nopart,1,partitions);
2354 free_Ivector(indx,1,noelements);
2355
2356 PartitionElementsByNodes(data,info);
2357
2358 if(info) printf("Successfully made a partitioning with %d to %d nodes.\n",minpart,maxpart);
2359
2360 return(0);
2361 }
2362
2363
LinearNodes(int elemtype)2364 static int LinearNodes(int elemtype)
2365 {
2366 int elemfamily,nonodes;
2367 int fam2nodemap[] = {0, 1, 2, 3, 4, 4, 5, 6, 8 };
2368
2369 elemfamily = elemtype / 100;
2370 nonodes = fam2nodemap[elemfamily];
2371
2372 return(nonodes);
2373 }
2374
2375
PartitionMetisMesh(struct FemType * data,struct ElmergridType * eg,int partitions,int dual,int info)2376 int PartitionMetisMesh(struct FemType *data,struct ElmergridType *eg,
2377 int partitions,int dual,int info)
2378 /* Perform partitioning using Metis. This uses the elemental routines of Metis that assume that
2379 there exists only one elementtype. If this condition is not met then this routine cannot be
2380 used. If the elements are higher order nodal elements then use only the linear basis. */
2381 {
2382 int i,j,k,periodic, noelements, noknots, sides, highorder;
2383 int nodesd2, etype, numflag,mintype,maxtype,elemtype,minnodes;
2384 int *neededby,*indxper,*inpart;
2385 idx_t *metistopo,*eptr,*npart,*epart;
2386 idx_t ne,nn,ncommon,edgecut,nparts;
2387 idx_t options[METIS_NOPTIONS];
2388
2389 METIS_SetDefaultOptions(options);
2390 options[METIS_OPTION_NUMBERING] = 0;
2391 options[METIS_OPTION_CONTIG] = eg->metiscontig;
2392 options[METIS_OPTION_DBGLVL] = 3;
2393 if( eg->metisseed ) options[METIS_OPTION_SEED] = eg->metisseed;
2394
2395 if(info) printf("Making a Metis partitioning for %d elements in %d-dimensions.\n",
2396 data->noelements,data->dim);
2397
2398 noelements = data->noelements;
2399 noknots = data->noknots;
2400
2401 for(i=1;i<=noelements;i++) {
2402 elemtype = data->elementtypes[i];
2403 nodesd2 = LinearNodes( elemtype );
2404 if(i == 1 ) {
2405 mintype = maxtype = elemtype;
2406 minnodes = nodesd2;
2407 } else {
2408 mintype = MIN( mintype, elemtype );
2409 maxtype = MAX( maxtype, elemtype );
2410 minnodes = MIN( minnodes, nodesd2 );
2411 }
2412 }
2413
2414 if(info) {
2415 if(mintype == maxtype ) {
2416 printf("All elements are of type %d\n",mintype);
2417 }
2418 else {
2419 printf("Minimum element type is %d\n",mintype);
2420 printf("Maximum element type is %d\n",maxtype);
2421 }
2422 }
2423
2424 if( minnodes >= 8) {
2425 ncommon = 4;
2426 }
2427 else if(minnodes > 4 ) {
2428 ncommon = 3;
2429 }
2430 else {
2431 ncommon = 2;
2432 }
2433 if(info) {
2434 printf("Minimum number of linear nodes in elements %d\n",minnodes);
2435 printf("Requiring number of nodes in dual graph %d\n",ncommon);
2436 }
2437
2438 if(!data->partitionexist) {
2439 data->partitionexist = TRUE;
2440 data->elempart = Ivector(1,data->noelements);
2441 data->nodepart = Ivector(1,data->noknots);
2442 data->nopartitions = partitions;
2443 }
2444 inpart = data->elempart;
2445
2446 /* Are there periodic boundaries. This information is used to join the boundaries. */
2447 periodic = data->periodicexist;
2448 if(periodic) {
2449 if(info) printf("There seems to be periodic boundaries\n");
2450 indxper = data->periodic;
2451 }
2452
2453 ne = noelements;
2454 nn = noknots;
2455
2456 neededby = Ivector(1,noknots);
2457
2458 numflag = 0;
2459 nparts = partitions;
2460
2461 /* Mark the nodes that are needed */
2462 k = 0;
2463 for(i=1;i<=noknots;i++)
2464 neededby[i] = 0;
2465 if(periodic) {
2466 for(i=1;i<=noelements;i++) {
2467 nodesd2 = LinearNodes( data->elementtypes[i] );
2468 for(j=0;j<nodesd2;j++) {
2469 neededby[indxper[data->topology[i][j]]] = 1;
2470 k += 1;
2471 }
2472 }
2473 }
2474 else {
2475 for(i=1;i<=noelements;i++) {
2476 nodesd2 = LinearNodes( data->elementtypes[i] );
2477 for(j=0;j<nodesd2;j++) {
2478 neededby[data->topology[i][j]] = 1;
2479 k += 1;
2480 }
2481 }
2482 }
2483
2484 j = 0;
2485 for(i=1;i<=noknots;i++)
2486 if(neededby[i])
2487 neededby[i] = ++j;
2488
2489 nn = j;
2490 if(info) {
2491 if(nn == noknots)
2492 printf("Using all %d possible nodes in the Metis graph\n",nn);
2493 else
2494 printf("Using %d nodes of %d possible nodes in the Metis graph\n",nn,noknots);
2495 printf("Allocating mesh topology of size %d\n",k);
2496 }
2497
2498 eptr = Ivector(0,noelements);
2499 metistopo = Ivector(0,k-1);
2500
2501 k = 0;
2502 eptr[0] = k;
2503 if(periodic) {
2504 for(i=1;i<=noelements;i++) {
2505 nodesd2 = LinearNodes( data->elementtypes[i] );
2506 for(j=0;j<nodesd2;j++) {
2507 metistopo[k] = neededby[indxper[data->topology[i][j]]]-1;
2508 k += 1;
2509 }
2510 eptr[i] = k;
2511 }
2512 }
2513 else if(nn < noknots) {
2514 for(i=1;i<=noelements;i++) {
2515 nodesd2 = LinearNodes( data->elementtypes[i] );
2516 for(j=0;j<nodesd2;j++) {
2517 metistopo[k] = neededby[data->topology[i][j]]-1;
2518 k += 1;
2519 }
2520 eptr[i] = k;
2521 }
2522 }
2523 else {
2524 for(i=1;i<=noelements;i++) {
2525 nodesd2 = LinearNodes( data->elementtypes[i] );
2526 for(j=0;j<nodesd2;j++) {
2527 metistopo[k] = data->topology[i][j]-1;
2528 k += 1;
2529 }
2530 eptr[i] = k;
2531 }
2532 }
2533
2534 npart = Ivector(0,nn-1);
2535 for(i=0;i<nn;i++) npart[i] = 0;
2536
2537 epart = Ivector(0,ne-1);
2538 for(i=0;i<ne;i++) epart[i] = 0;
2539
2540 if(dual) {
2541 if(info) printf("Starting graph partitioning METIS_PartMeshDual.\n");
2542 METIS_PartMeshDual(&ne,&nn,eptr,metistopo,NULL,NULL,&ncommon,
2543 &nparts,NULL,options,&edgecut,epart,npart);
2544 if(info) printf("Finished graph partitioning METIS_PartMeshDual.\n");
2545 }
2546 else {
2547 if(info) printf("Starting graph partitioning METIS_PartMeshNodal.\n");
2548 METIS_PartMeshNodal(&ne,&nn,eptr,metistopo,NULL,NULL,
2549 &nparts,NULL,options,&edgecut,epart,npart);
2550 if(info) printf("Finished graph partitioning METIS_PartMeshNodal.\n");
2551 }
2552
2553 /* Set the partition given by Metis for each element. */
2554 for(i=1;i<=noelements;i++) {
2555 inpart[i] = epart[i-1]+1;
2556 if(inpart[i] < 1 || inpart[i] > partitions)
2557 printf("Invalid partition %d for element %d\n",inpart[i],i);
2558 }
2559
2560 if( highorder ) {
2561 PartitionNodesByElements(data,info);
2562 }
2563 else {
2564 if(info) printf("Set the partition given by Metis for each node\n");
2565 for(i=1;i<=noknots;i++) {
2566 if(periodic)
2567 j = neededby[indxper[i]];
2568 else if(nn < noknots)
2569 j = neededby[i];
2570 else
2571 j = i;
2572 if(!j) printf("Cannot set partitioning for node %d\n",i);
2573 data->nodepart[i] = npart[j-1]+1;
2574 if(data->nodepart[i] < 1 || data->nodepart[i] > partitions)
2575 printf("Invalid partition %d for node %d\n",data->nodepart[i],i);
2576 }
2577 }
2578
2579 free_Ivector(neededby,1,noknots);
2580 free_Ivector(metistopo,0,k-1);
2581 free_Ivector(epart,0,noelements-1);
2582 free_Ivector(npart,0,nn-1);
2583 free_Ivector(eptr,0,noelements+1);
2584
2585 if(info) printf("Successfully made a Metis partition using the element mesh.\n");
2586
2587 return(0);
2588 }
2589
2590
2591
2592 #if USE_METIS
PartitionMetisGraph(struct FemType * data,struct BoundaryType * bound,struct ElmergridType * eg,int partitions,int metisopt,int dual,int info)2593 int PartitionMetisGraph(struct FemType *data,struct BoundaryType *bound,
2594 struct ElmergridType *eg,int partitions,int metisopt,
2595 int dual,int info)
2596 /* Perform partitioning using Metis. One may use either the direct or the dual
2597 graph. The direct graph means that the nodes are first partitioned and the
2598 elements follow. The dual graph means that the elements are partitioned and
2599 the ownership of the nodes will follow. The latter is optimal for Elmer. */
2600 {
2601 int i,j,k,noelements,noknots,errstat;
2602 int nn,ncon,con,maxcon,totcon;
2603 int *xadj,*adjncy,*vwgt,*adjwgt,wgtflag,*npart,**graph;
2604 int numflag,nparts,edgecut,maxconset;
2605 struct CRSType *dualgraph;
2606 idx_t options[METIS_NOPTIONS];
2607
2608 if(info) printf("Making a Metis partitioning for %d nodes in %d-dimensions.\n",
2609 data->noknots,data->dim);
2610 if(partitions < 2 ) {
2611 bigerror("There should be at least two partitions for partitioning!");
2612 }
2613
2614 noknots = data->noknots;
2615 noelements = data->noelements;
2616 maxconset = 0;
2617
2618 if(data->periodicexist && dual ) {
2619 printf("Dual graph not implemented for periodic system!\n");
2620 printf("Enforcing nodal graph for partitioning\n");
2621 dual = FALSE;
2622 }
2623
2624 nparts = partitions;
2625 if( dual ) {
2626 if( eg->partbcz > 1 || eg->partbcr )
2627 PartitionConnectedElements1D(data,bound,eg,info);
2628 else if( eg->partbcmetis > 1 )
2629 PartitionConnectedElementsMetis(data,bound,eg->partbcmetis,metisopt,info);
2630 else if( eg->connect ) {
2631 PartitionConnectedElementsStraight(data,bound,eg,info);
2632 }
2633
2634 if( data->nodeconnectexist ) {
2635 errstat = ExtendBoundaryPartitioning(data,bound,eg->partbclayers,info);
2636 if( errstat ) bigerror("Extend boundary partitioning returned error!");
2637 }
2638
2639 /* Find the number of partitions already used for connected elements */
2640 if( data->elemconnectexist ) {
2641 maxconset = 0;
2642 for(i=1;i<=noelements;i++)
2643 maxconset = MAX( maxconset, data->elemconnect[i] );
2644 nparts -= maxconset;
2645 }
2646
2647 if( nparts == 1 ) {
2648 if(info) printf("Just one partition left, skipping 2nd part of hybrid partitioning!\n");
2649 goto skippart;
2650 }
2651
2652 CreateDualGraph(data,TRUE,info);
2653
2654 /* Create the sparse matrix format graph for Metis */
2655 dualgraph = &data->dualgraph;
2656
2657 nn = dualgraph->rowsize;
2658 xadj = dualgraph->rows;
2659 adjncy = dualgraph->cols;
2660 totcon = dualgraph->colsize;
2661 }
2662 else {
2663 CreateNodalGraph(data,TRUE,info);
2664 maxcon = data->nodalmaxconnections;
2665 nn = noknots;
2666 graph = data->nodalgraph;
2667
2668 /* Compute total number of connections in graph */
2669 totcon = 0;
2670 for(i=1;i<=nn;i++) {
2671 for(j=0;j<maxcon;j++) {
2672 con = graph[j][i];
2673 if(con) totcon++;
2674 }
2675 }
2676
2677 /* Create the sparse matrix format graph for Metis */
2678 xadj = Ivector(0,nn);
2679 adjncy = Ivector(0,totcon-1);
2680 for(i=0;i<totcon;i++)
2681 adjncy[i] = 0;
2682
2683 totcon = 0;
2684 for(i=1;i<=nn;i++) {
2685 xadj[i-1] = totcon;
2686 for(j=0;j<maxcon;j++) {
2687 con = graph[j][i];
2688 if(con) {
2689 adjncy[totcon] = con-1;
2690 totcon++;
2691 }
2692 }
2693 }
2694 xadj[nn] = totcon;
2695 }
2696
2697 if(info) printf("There are %d connections altogether in the graph.\n",totcon);
2698
2699 /* Parameters for Metis */
2700 numflag = 0;
2701 npart = Ivector(0,nn-1);
2702 wgtflag = 0;
2703 ncon = 1;
2704
2705 METIS_SetDefaultOptions(options);
2706
2707 options[METIS_OPTION_NUMBERING] = 0;
2708 options[METIS_OPTION_CONTIG] = eg->metiscontig;
2709 if( eg->metisseed ) options[METIS_OPTION_SEED] = eg->metisseed;
2710
2711 options[METIS_OPTION_CTYPE] = METIS_CTYPE_SHEM;
2712 options[METIS_OPTION_IPTYPE] = METIS_IPTYPE_GROW;
2713 options[METIS_OPTION_RTYPE] = METIS_RTYPE_GREEDY;
2714 if( metisopt == 4 ) options[METIS_OPTION_MINCONN] = 1;
2715 options[METIS_OPTION_DBGLVL] = 3;
2716
2717 /* Optional weights */
2718 vwgt = NULL;
2719 adjwgt = NULL;
2720
2721 if( !dual ) {
2722 /* Create weights if needed */
2723 if(data->periodicexist || eg->connect) {
2724 if(info) printf("Creating weight for %d connections.\n",totcon);
2725 wgtflag = 1;
2726 adjwgt = Ivector(0,totcon-1);
2727 for(i=0;i<totcon;i++)
2728 adjwgt[i] = 1;
2729
2730 if(metisopt != 3) {
2731 printf("For weighted partitioning Metis subroutine METIS_PartGraphKway is enforced\n");
2732 metisopt = 3;
2733 }
2734 }
2735
2736 /* Make the periodic connections stronger */
2737 if(data->periodicexist) {
2738 if(info) printf("Setting periodic connections to dominate %d\n",totcon);
2739 for(i=0;i<noknots;i++) {
2740 j = data->periodic[i+1]-1;
2741 if(j == i) continue;
2742 for(k=xadj[i];k<xadj[i+1];k++)
2743 if(adjncy[k] == j) adjwgt[k] = maxcon;
2744 }
2745 }
2746
2747 /* Make the constraint connections stronger */
2748 if(eg->connect) {
2749 int maxweight;
2750 int con,bc,bctype,sideelemtype,sidenodes;
2751 int j2,ind,ind2;
2752 int sideind[MAXNODESD1];
2753
2754 maxweight = noknots+noelements;
2755 printf("Adding weight of %d for constrained nodes\n",maxweight);
2756
2757 for(con=1;con<=eg->connect;con++) {
2758 bctype = eg->connectbounds[con-1];
2759
2760 for(bc=0;bc<MAXBOUNDARIES;bc++) {
2761 if(bound[bc].created == FALSE) continue;
2762 if(bound[bc].nosides == 0) continue;
2763
2764 for(i=1;i<=bound[bc].nosides;i++) {
2765 if(bound[bc].types[i] != bctype) continue;
2766
2767 GetBoundaryElement(i,&bound[bc],data,sideind,&sideelemtype);
2768 /* GetElementSide(bound[bc].parent[i],bound[bc].side[i],bound[bc].normal[i],
2769 data,sideind,&sideelemtype); */
2770
2771 sidenodes = sideelemtype%100;
2772
2773 for(j=0;j<sidenodes;j++) {
2774 for(j2=0;j2<sidenodes;j2++) {
2775 if(j==j2) continue;
2776
2777 ind = sideind[j]-1;
2778 ind2 = sideind[j2]-1;
2779
2780 for(k=xadj[ind];k<xadj[ind+1];k++)
2781 if(adjncy[k] == ind2) adjwgt[k] = maxweight;
2782 }
2783 }
2784 }
2785 }
2786 }
2787 }
2788 } /* !dual */
2789
2790 if(metisopt == 2) {
2791 if(info) printf("Starting graph partitioning METIS_PartGraphRecursive.\n");
2792 METIS_PartGraphRecursive(&nn,&ncon,xadj,adjncy,vwgt,&wgtflag,adjwgt,
2793 &nparts,NULL,NULL,options,&edgecut,npart);
2794 }
2795 else if( metisopt == 3 || metisopt == 4 ) {
2796 ncon = 1;
2797 wgtflag = 0;
2798 if(info) printf("Starting graph partitioning METIS_PartGraphKway.\n");
2799 METIS_PartGraphKway(&nn, /* number of vertices in the graph */
2800 &ncon, /* number of balancing constraints */
2801 xadj, adjncy, /* the adjacency structure of the graph */
2802 vwgt, /* weights of the vertices */
2803 &wgtflag, /* size of the vectices for computing communication */
2804 adjwgt, /* weight of the edges */
2805 &nparts, /* number of partitions */
2806 NULL, /* weights for each partition and constraint */
2807 NULL, /* allowed load imbalance */
2808 options, /* array of options */
2809 &edgecut, /* the total communication volume */
2810 npart); /* partition vector of the graph */
2811 }
2812 else {
2813 printf("Unknown Metis option %d\n",metisopt);
2814 }
2815
2816 if(info) printf("Finished Metis graph partitioning call.\n");
2817
2818
2819 free_Ivector(adjncy,0,totcon-1);
2820 free_Ivector(xadj,0,nn);
2821 if(wgtflag == 1) free_Ivector(adjwgt,0,totcon-1);
2822
2823 skippart:
2824
2825 if(!data->partitionexist) {
2826 data->partitionexist = TRUE;
2827 data->elempart = Ivector(1,noelements);
2828 data->nodepart = Ivector(1,noknots);
2829 data->nopartitions = partitions;
2830 }
2831
2832 /* Set the partition given by Metis for each node. */
2833 if( dual ) {
2834 if( data->elemconnectexist ) {
2835 for(i=1;i<=noelements;i++) {
2836 j = data->elemconnect[i];
2837 if( nparts == 1 ) {
2838 if(j)
2839 data->elempart[i] = j;
2840 else
2841 data->elempart[i] = maxconset+1;
2842 } else {
2843 if(j < 0)
2844 data->elempart[i] = -j;
2845 else
2846 data->elempart[i] = npart[j-1]+1+maxconset;
2847 }
2848 }
2849 }
2850 else {
2851 for(i=1;i<=nn;i++)
2852 data->elempart[i] = npart[i-1]+1;
2853 }
2854 PartitionNodesByElements(data,info);
2855 }
2856 else {
2857 for(i=1;i<=nn;i++)
2858 data->nodepart[i] = npart[i-1]+1;
2859 PartitionElementsByNodes(data,info);
2860
2861 /* Finally check that the constraint is really honored */
2862 if(!eg->connect) {
2863 int con,bc,bctype,sideelemtype,sidenodes,par;
2864 int ind,sideind[MAXNODESD1];
2865 int *sidehits,sidepartitions;
2866
2867 printf("Checking connection integrity\n");
2868 sidehits = Ivector(1,partitions);
2869
2870 for(con=1;con<=eg->connect;con++) {
2871 bctype = eg->connectbounds[con-1];
2872
2873 for(i=1;i<=partitions;i++)
2874 sidehits[i] = 0;
2875
2876 for(bc=0;bc<MAXBOUNDARIES;bc++) {
2877 if(bound[bc].created == FALSE) continue;
2878 if(bound[bc].nosides == 0) continue;
2879
2880 for(i=1;i<=bound[bc].nosides;i++) {
2881 if(bound[bc].types[i] != bctype) continue;
2882
2883 if(1)
2884 GetBoundaryElement(i,&bound[bc],data,sideind,&sideelemtype);
2885 else
2886 GetElementSide(bound[bc].parent[i],bound[bc].side[i],bound[bc].normal[i],
2887 data,sideind,&sideelemtype);
2888 sidenodes = sideelemtype%100;
2889
2890 for(j=0;j<sidenodes;j++) {
2891 ind = sideind[j];
2892 par = data->nodepart[ind];
2893 sidehits[par] += 1;
2894 }
2895 }
2896 }
2897
2898 sidepartitions = 0;
2899 for(i=1;i<=partitions;i++)
2900 if( sidehits[i] ) sidepartitions += 1;
2901
2902 if(sidepartitions != 1) {
2903 printf("PartitionMetisGraph: side %d belongs to %d partitions\n",bctype,sidepartitions);
2904 bigerror("Parallel constraints might not be set!");
2905 }
2906 }
2907 }
2908 }
2909
2910 if( nparts > 1 ) {
2911 free_Ivector(npart,0,nn-1);
2912 }
2913
2914 if(info) {
2915 if( dual )
2916 printf("Successfully made a Metis partition using the dual graph.\n");
2917 else
2918 printf("Successfully made a Metis partition using the nodal graph.\n");
2919 }
2920
2921 return(0);
2922 }
2923
2924 #endif
2925
2926
CheckPartitioning(struct FemType * data,int info)2927 static void CheckPartitioning(struct FemType *data,int info)
2928 {
2929 int i,j,k,partitions,part,part2,noknots,noelements,mini,maxi,sumi,hit,ind,nodesd2,elemtype;
2930 int *elempart, *nodepart,*elemsinpart,*nodesinpart,*sharedinpart;
2931
2932 noknots = data->noknots;
2933 noelements = data->noelements;
2934 partitions = data->nopartitions;
2935 elemsinpart = Ivector(1,partitions);
2936 nodesinpart = Ivector(1,partitions);
2937 sharedinpart = Ivector(1,partitions);
2938 for(i=1;i<=partitions;i++)
2939 elemsinpart[i] = nodesinpart[i] = sharedinpart[i] = 0;
2940
2941 if(info) printf("Checking for partitioning\n");
2942
2943 /* Check that division of elements */
2944 elempart = data->elempart;
2945 j=0;
2946 for(i=1;i<=data->noelements;i++) {
2947 part = elempart[i];
2948 if(part < 1 || part > partitions)
2949 j++;
2950 else
2951 elemsinpart[part] += 1;
2952 }
2953 if(j) {
2954 printf("Bad initial partitioning: %d elements do not belong anywhere!\n",j);
2955 bigerror("Can't continue with broken partitioning");
2956 }
2957
2958 /* Check the division of nodes */
2959 nodepart = data->nodepart;
2960 j=0;
2961 for(i=1;i<=data->noknots;i++) {
2962 part = nodepart[i];
2963 if(part < 1 || part > partitions)
2964 j++;
2965 else
2966 nodesinpart[part] += 1;
2967 }
2968
2969 if(j) {
2970 printf("Bad initial partitioning: %d nodes do not belong anywhere!\n",j);
2971 bigerror("Can't continue with broken partitioning");
2972 }
2973
2974 if(data->partitiontableexists) {
2975 for(i=1;i<=noknots;i++) {
2976 part = nodepart[i];
2977 for(j=1;j<=data->maxpartitiontable;j++) {
2978 part2 = data->partitiontable[j][i];
2979 if(!part2) break;
2980 if(part != part2) sharedinpart[part2] += 1;
2981 }
2982 }
2983 }
2984
2985 if(info) {
2986 printf("Information on partition bandwidth\n");
2987 if(partitions <= 10) {
2988 printf("Distribution of elements, nodes and shared nodes\n");
2989 printf(" %-10s %-10s %-10s %-10s\n","partition","elements","nodes","shared");
2990 for(i=1;i<=partitions;i++)
2991 printf(" %-10d %-10d %-10d %-10d\n",i,elemsinpart[i],nodesinpart[i],sharedinpart[i]);
2992 }
2993 else {
2994 mini = maxi = elemsinpart[1];
2995 for(i=1;i<=partitions;i++) {
2996 mini = MIN( elemsinpart[i], mini);
2997 maxi = MAX( elemsinpart[i], maxi);
2998 }
2999 printf("Average %d elements with range %d in partition\n",noelements/partitions,maxi-mini);
3000
3001 mini = maxi = nodesinpart[1];
3002 for(i=1;i<=partitions;i++) {
3003 mini = MIN( nodesinpart[i], mini);
3004 maxi = MAX( nodesinpart[i], maxi);
3005 }
3006 printf("Average %d nodes with range %d in partition\n",noknots/partitions,maxi-mini);
3007
3008 sumi = 0;
3009 mini = maxi = sharedinpart[1];
3010 for(i=1;i<=partitions;i++) {
3011 mini = MIN( sharedinpart[i], mini);
3012 maxi = MAX( sharedinpart[i], maxi);
3013 sumi += sharedinpart[i];
3014 }
3015 printf("Average %d shared nodes with range %d in partition\n",sumi/partitions,maxi-mini);
3016 }
3017 }
3018
3019 if(!data->maxpartitiontable) return;
3020
3021 if(0) printf("Checking that each node in elements belongs to nodes\n");
3022 for(i=1;i<=data->noelements;i++) {
3023 part = elempart[i];
3024 elemtype = data->elementtypes[i];
3025 nodesd2 = elemtype % 100;
3026
3027 for(j=0;j < nodesd2;j++) {
3028 ind = data->topology[i][j];
3029
3030 hit = FALSE;
3031 part2 = 0;
3032 for(k=1;k<=data->maxpartitiontable;k++) {
3033 part2 = data->partitiontable[k][ind];
3034 if( part == part2 ) hit = TRUE;
3035 if(hit && !part) break;
3036 }
3037 if(!hit) {
3038 printf("******** Warning *******\n");
3039 printf("Node %d in element %d does not belong to partition %d (%d)\n",ind,i,part,part2);
3040 printf("elemtype = %d nodesd2 = %d\n",elemtype,nodesd2);
3041 for(k=0;k < nodesd2;k++)
3042 printf("ind[%d] = %d\n",k,data->topology[i][k]);
3043 }
3044 }
3045 }
3046
3047 if(0) printf("Checking that each node in partition is shown in partition list\n");
3048 for(i=1;i<=data->noknots;i++) {
3049 part = nodepart[i];
3050
3051 hit = FALSE;
3052 for(j=1;j<=data->maxpartitiontable;j++) {
3053 part2 = data->partitiontable[j][i];
3054 if( part == part2 ) hit = TRUE;
3055 if(hit && !part) break;
3056 }
3057 if(!hit) {
3058 printf("***** Node %d in partition %d is not in partition list\n",i,part);
3059 }
3060 }
3061 }
3062
3063
OptimizePartitioningAtBoundary(struct FemType * data,struct BoundaryType * bound,int info)3064 int OptimizePartitioningAtBoundary(struct FemType *data,struct BoundaryType *bound,int info)
3065 {
3066 int i,j,k,l,boundaryelems,ind,hit1,hit2,fix1,fix2;
3067 int dompart,part1,part2,newmam,mam1,mam2,nodesd2;
3068 int *alteredparent;
3069
3070 if(!data->partitionexist) {
3071 printf("OptimizePartitioningAtBoundary: this should be called only after partitioning\n");
3072 bigerror("Optimization not performed!");
3073 }
3074
3075 if(info) printf("Optimizing the partitioning at boundaries.\n");
3076
3077 /* Memorize the original parent */
3078 alteredparent = Ivector(1,data->noelements);
3079 for(i=1;i<=data->noelements;i++)
3080 alteredparent[i] = data->elempart[i];
3081
3082
3083 /* Set the secondary parent to be a parent also because we want all
3084 internal BCs to be within the same partition.
3085 Also set the nodes of the altered elements to be in the desired partition. */
3086 k = 0;
3087 do {
3088 k++;
3089 boundaryelems = 0;
3090
3091 for(j=0;j < MAXBOUNDARIES;j++) {
3092 if(!bound[j].created) continue;
3093 for(i=1; i <= bound[j].nosides; i++) {
3094 if(bound[j].ediscont)
3095 if(bound[j].discont[i]) continue;
3096
3097 mam1 = bound[j].parent[i];
3098 mam2 = abs(bound[j].parent2[i]);
3099 if(!mam1 || !mam2) continue;
3100
3101 part1 = data->elempart[mam1];
3102 part2 = data->elempart[mam2];
3103 if(part1 == part2) continue;
3104
3105 /* Check if both nodes are enforced to be fixed */
3106 if( data->elemconnectexist ) {
3107 fix1 = ( data->elemconnect[mam1] < 0 );
3108 fix2 = ( data->elemconnect[mam2] < 0 );
3109 if( fix1 && fix2 ) continue;
3110 }
3111 else {
3112 fix1 = fix2 = 0;
3113 }
3114
3115 /* The first iterations check which parents is ruling
3116 thereafter choose pragmatically the other to overcome
3117 oscillating solutions. */
3118 if(fix1 || fix2 ) {
3119 }
3120 if(k < 5) {
3121 hit1 = hit2 = 0;
3122 nodesd2 = data->elementtypes[mam1] % 100;
3123 for(l=0;l < nodesd2;l++) {
3124 ind = data->topology[mam1][l];
3125 if(data->nodepart[ind] == part1) hit1++;
3126 if(data->nodepart[ind] == part2) hit2++;
3127 }
3128 nodesd2 = data->elementtypes[mam2] % 100;
3129 for(l=0;l < nodesd2;l++) {
3130 ind = data->topology[mam2][l];
3131 if(data->nodepart[ind] == part1) hit1++;
3132 if(data->nodepart[ind] == part2) hit2++;
3133 }
3134 fix1 = ( hit1 >= hit2 );
3135 fix2 = !fix1;
3136 }
3137 else {
3138 fix1 = TRUE;
3139 fix2 = FALSE;
3140 }
3141
3142 /* Make the more ruling parent dominate the whole boundary */
3143 if(fix1) {
3144 dompart = part1;
3145 newmam = mam2;
3146 }
3147 else {
3148 dompart = part2;
3149 newmam = mam1;
3150 }
3151
3152 data->elempart[newmam] = dompart;
3153 boundaryelems++;
3154
3155 /* Move the ownership of all nodes to the leading partition */
3156 if(0) {
3157 nodesd2 = data->elementtypes[newmam] % 100;
3158 for(l=0;l < nodesd2;l++) {
3159 ind = data->topology[newmam][l];
3160 data->nodepart[ind] = dompart;
3161 }
3162 }
3163 else if(0) {
3164 nodesd2 = data->elementtypes[mam1] % 100;
3165 for(l=0;l < nodesd2;l++) {
3166 ind = data->topology[mam1][l];
3167 data->nodepart[ind] = dompart;
3168 }
3169 nodesd2 = data->elementtypes[mam2] % 100;
3170 for(l=0;l < nodesd2;l++) {
3171 ind = data->topology[mam2][l];
3172 data->nodepart[ind] = dompart;
3173 }
3174 }
3175
3176 }
3177 }
3178 if(info && boundaryelems)
3179 printf("Round %d: %d bulk elements with BCs removed from interface.\n",
3180 k,boundaryelems);
3181 } while(boundaryelems && k < 10);
3182
3183
3184 j = 0;
3185 for(i=1;i<=data->noelements;i++) {
3186 if(alteredparent[i] == data->elempart[i])
3187 alteredparent[i] = 0;
3188 else
3189 j++;
3190 }
3191 free_Ivector(alteredparent,1,data->noelements);
3192 if(info) printf("Ownership of %d parents was changed at BCs\n",j);
3193
3194
3195
3196 /* Remove the negative secondary parents that were only used to optimize the partitioning */
3197 for(j=0;j < MAXBOUNDARIES;j++) {
3198 if(!bound[j].created) continue;
3199 for(i=1; i <= bound[j].nosides; i++)
3200 bound[j].parent2[i] = MAX(0, bound[j].parent2[i]);
3201 }
3202
3203 if(0) printf("The partitioning at discontinuous gaps was optimized.\n");
3204 return(0);
3205 }
3206
3207
Levelize(int n,int level,int * maxlevel,int * levels,int * rows,int * cols,int * done)3208 static void Levelize(int n,int level,int *maxlevel,int *levels,int *rows,int *cols,int *done)
3209 {
3210 int j,k;
3211
3212 levels[n] = level;
3213 done[n] = TRUE;
3214 *maxlevel = MAX( *maxlevel,level );
3215
3216 for( j=rows[n]; j<rows[n+1]; j++ ) {
3217 k = cols[j];
3218 if ( !done[k] ) Levelize(k,level+1,maxlevel,levels,rows,cols,done);
3219 }
3220 }
3221
3222
RenumberCuthillMckee(int nrows,int * rows,int * cols,int * iperm)3223 static int RenumberCuthillMckee( int nrows, int *rows, int *cols, int *iperm )
3224 {
3225 int i,j,k,n,startn,mindegree,maxlevel,newroot,bw_bef,bw_aft;
3226 int *level,*degree,*done;
3227
3228 done = Ivector(0,nrows-1);
3229 level = Ivector(0,nrows-1);
3230 degree = Ivector(0,nrows-1);
3231
3232 bw_bef = 0;
3233 for(i=0; i<nrows; i++ )
3234 {
3235 for( j=rows[i]; j<rows[i+1]; j++ )
3236 bw_bef = MAX( bw_bef, ABS(cols[j]-i)+1 );
3237 degree[i] = rows[i+1]-rows[i];
3238 }
3239 printf( "RenumberCuthillMckee: Bandwidth before: %d\n", bw_bef );
3240
3241 startn = 0;
3242 mindegree = degree[startn];
3243 for( i=0; i<nrows; i++ ) {
3244 if ( degree[i] < mindegree ) {
3245 startn = i;
3246 mindegree = degree[i];
3247 }
3248 level[i] = 0;
3249 }
3250
3251 maxlevel = 0;
3252 for( i=0; i<nrows; i++ ) done[i]=FALSE;
3253
3254 Levelize( startn,0,&maxlevel,level,rows,cols,done );
3255
3256 newroot = TRUE;
3257 while(newroot) {
3258 newroot = FALSE;
3259 mindegree = degree[startn];
3260 k = startn;
3261
3262 for( i=0; i<nrows; i++ ) {
3263 if ( level[i] == maxlevel ) {
3264 if ( degree[i] < mindegree ) {
3265 k = i;
3266 mindegree = degree[i];
3267 }
3268 }
3269 }
3270
3271 if ( k != startn ) {
3272 j = maxlevel;
3273 maxlevel = 0;
3274 for(i=0; i<nrows; i++ ) done[i]=FALSE;
3275
3276 Levelize( k,0,&maxlevel,level,rows,cols,done );
3277
3278 if ( j > maxlevel ) {
3279 newroot = TRUE;
3280 startn = j;
3281 }
3282 }
3283 }
3284
3285 for(i=0; i<nrows; i++ ) done[i]=-1,iperm[i]=-1;
3286
3287 done[0]=startn;
3288 iperm[startn]=0;
3289 i=1;
3290
3291 for( j=0; j<nrows; j++ ) {
3292 if ( done[j]<0 ) {
3293 for( k=0; k<nrows; k++ ) {
3294 if ( iperm[k]<0 ) {
3295 done[i]=k;
3296 iperm[k]=i;
3297 i++;
3298 break;
3299 }
3300 }
3301 }
3302
3303 for( k=rows[done[j]]; k<rows[done[j]+1]; k++) {
3304 n = cols[k];
3305 if ( iperm[n]<0 ) {
3306 done[i] = n;
3307 iperm[n] = i;
3308 i++;
3309 }
3310 }
3311 }
3312
3313 for( i=0; i<nrows; i++ )
3314 iperm[done[i]] = nrows-1-i;
3315
3316 bw_aft = 0;
3317 for(i=0; i<nrows; i++ )
3318 for( j=rows[i]; j<rows[i+1]; j++ )
3319 bw_aft = MAX( bw_aft, ABS(iperm[cols[j]]-iperm[i])+1 );
3320
3321 printf( "RenumberCuthillMckee: Bandwidth after: %d\n", bw_aft );
3322
3323 free_Ivector(level,0,nrows-1);
3324 free_Ivector(done,0,nrows-1);
3325 free_Ivector(degree,0,nrows-1);
3326
3327 return bw_aft < bw_bef;
3328 }
3329
3330
3331
RenumberPartitions(struct FemType * data,int info)3332 static void RenumberPartitions(struct FemType *data,int info)
3333 /* Minimize bandwidth of partition indexing. This could be favourable if the
3334 communication between neighbouring partitions is faster than between
3335 partitions with larger partition index. Probably there is not much
3336 difference for normal hardware configurations. */
3337 {
3338 int i,j,k,l,m,hit,con,totcon,noelements,noknots,partitions;
3339 int maxneededtimes,totneededtimes;
3340 int part,part1,part2,bw_reduced;
3341 int *nodepart,*elempart;
3342 int *perm;
3343 int *xadj,*adjncy;
3344 int *partparttable[MAXCONNECTIONS];
3345
3346
3347 if(info) printf("Renumbering partitions to minimize bandwidth.\n");
3348
3349 partitions = data->nopartitions;
3350 noelements = data->noelements;
3351 noknots = data->noknots;
3352 maxneededtimes = data->maxpartitiontable;
3353
3354
3355 /* Make the partition-partition list from the node-partition list */
3356 totneededtimes = 0;
3357 totcon = 0;
3358 for(i=1;i<=noknots;i++) {
3359 if(data->partitiontable[2][i] == 0) continue;
3360 for(j=1;j<=maxneededtimes;j++) {
3361 part1 = data->partitiontable[j][i];
3362 if(!part1) break;
3363
3364 for(k=1;k<=maxneededtimes;k++) {
3365 if(k==j) continue;
3366 part2 = data->partitiontable[k][i];
3367 if(!part2) break;
3368
3369 hit = 0;
3370 for(l=1;l<=totneededtimes;l++) {
3371 if(partparttable[l][part1] == part2) {
3372 hit = -1;
3373 break;
3374 }
3375 else if(partparttable[l][part1] == 0) {
3376 totcon++;
3377 partparttable[l][part1] = part2;
3378 hit = 1;
3379 break;
3380 }
3381 }
3382 if(!hit) {
3383 totneededtimes++;
3384 partparttable[totneededtimes] = Ivector(1,partitions);
3385 for(m=1;m<=partitions;m++)
3386 partparttable[totneededtimes][m] = 0;
3387 partparttable[totneededtimes][part1] = part2;
3388 totcon++;
3389 }
3390 }
3391 }
3392 }
3393
3394 if(info) {
3395 printf("There are %d connections altogether\n",totcon);
3396 printf("There are %.3f connections between partitions in average\n",1.0*totcon/partitions);
3397 }
3398
3399
3400 xadj = Ivector(0,partitions);
3401 adjncy = Ivector(0,totcon-1);
3402 for(i=0;i<totcon;i++)
3403 adjncy[i] = 0;
3404
3405 totcon = 0;
3406 for(i=1;i<=partitions;i++) {
3407 xadj[i-1] = totcon;
3408 for(j=1;j<=totneededtimes;j++) {
3409 con = partparttable[j][i];
3410 if(!con) continue;
3411 adjncy[totcon] = con-1;
3412 totcon++;
3413 }
3414 }
3415 xadj[partitions] = totcon;
3416
3417 perm = Ivector(0,partitions-1);
3418 bw_reduced = RenumberCuthillMckee( partitions, xadj, adjncy, perm );
3419
3420
3421 /* Print the new order of partitions */
3422 if(0 && info) {
3423 printf( "Partition order after Cuthill-McKee bandwidth optimization: \n" );
3424 for(i=0;i<partitions;i++)
3425 printf("old=%d new=%d\n",i,perm[i] );
3426 }
3427
3428 /* Use the renumbering or not */
3429 if(bw_reduced) {
3430 if(info) printf("Successful ordering: moving partitions to new positions\n");
3431 nodepart = data->nodepart;
3432 elempart = data->elempart;
3433 for(i=1;i<=noelements;i++)
3434 elempart[i] = perm[elempart[i]-1]+1;
3435 for(i=1;i<=noknots;i++)
3436 nodepart[i] = perm[nodepart[i]-1]+1;
3437 for(i=1;i<=noknots;i++) {
3438 for(j=1;j<=maxneededtimes;j++) {
3439 part = data->partitiontable[j][i];
3440 if(!part) break;
3441 data->partitiontable[j][i] = perm[part-1]+1;
3442 }
3443 }
3444 }
3445
3446 for(i=1;i<=totneededtimes;i++)
3447 free_Ivector(partparttable[i],1,partitions);
3448 free_Ivector(xadj,0,partitions);
3449 free_Ivector(adjncy,0,totcon-1);
3450
3451
3452 }
3453
3454
CheckSharedDeviation(int * neededvector,int partitions,int info)3455 static int CheckSharedDeviation(int *neededvector,int partitions,int info)
3456 {
3457 int i,minshared,maxshared,dshared;
3458 Real sumshared, sumshared2, varshared, needed, det, ave;
3459
3460 sumshared = sumshared2 = 0.0;
3461 minshared = maxshared = neededvector[1];
3462 for(i=1;i<=partitions;i++) {
3463 needed = 1.0 * neededvector[i];
3464 sumshared += needed;
3465 sumshared2 += (needed * needed);
3466 maxshared = MAX(maxshared, neededvector[i]);
3467 minshared = MIN(minshared, neededvector[i]);
3468 }
3469
3470 dshared = maxshared - minshared;
3471
3472 if(info) {
3473 ave = sumshared / partitions;
3474 det = 1.0*sumshared2 / partitions - 1.0*(sumshared/partitions)*(sumshared / partitions);
3475 varshared = sqrt( det );
3476
3477 printf("Average number of elements in partition %.3le\n",ave);
3478 printf("Maximum deviation in ownership %d\n",dshared);
3479 printf("Average deviation in ownership %.3le\n",varshared);
3480 printf("Average relative deviation %.2lf %%\n",100.0*varshared/ave);
3481 }
3482
3483 return(dshared);
3484 }
3485
3486
3487
3488
OptimizePartitioning(struct FemType * data,struct BoundaryType * bound,int noopt,int partbw,int info)3489 int OptimizePartitioning(struct FemType *data,struct BoundaryType *bound,int noopt,
3490 int partbw, int info)
3491 /* Optimize partitioning of elements so that each partition has as closely as possible the
3492 desired amount of elements. Also, ir requested, check that there are no odd couplings within
3493 elelements that are not directly present at the given partition. It is a bit unclear
3494 to which extent these checks are needed in the current Elmer version. */
3495 {
3496 int i,j,k,l,n,m,noelements,partitions,ind,hit;
3497 int noknots,dshared,dshared0;
3498 int *elempart,*nodepart,sharings,maxshared;
3499 int nodesd2,maxneededtimes,*probnodes=NULL,optimize;
3500 int *neededvector;
3501 int somethingdone = 0;
3502 int *elemparts = NULL,*invelemparts = NULL;
3503 int **knows = NULL;
3504
3505 if(!data->partitionexist) {
3506 printf("OptimizePartitioning: this should be called only after partitioning\n");
3507 bigerror("Optimization not performed!");
3508 }
3509
3510 noknots = data->noknots;
3511 noelements = data->noelements;
3512 partitions = data->nopartitions;
3513 if(info) printf("Optimizing for %d partitions\n", partitions);
3514
3515 elempart = data->elempart;
3516 nodepart = data->nodepart;
3517 maxshared = 0;
3518
3519 if( partitions < 2 ) {
3520 printf("OptimizePartitioning: does not make sense for %d partitions\n",partitions);
3521 bigerror("Optimization not performed!");
3522 }
3523
3524 /* Create a table showing to which partitions nodes belong to */
3525 CreatePartitionTable(data,info);
3526 maxneededtimes = data->maxpartitiontable;
3527
3528 /* Renumber the bandwidth of partition-partition connections */
3529 if(partbw) RenumberPartitions(data,info);
3530
3531 /* Check partitioning after table is created for the first time */
3532 printf("Checking partitioning before optimization\n");
3533 CheckPartitioning(data,info);
3534
3535 /* Calculate how many nodes is owned by each partition */
3536 neededvector = Ivector(1,partitions);
3537 for(i=1;i<=partitions;i++)
3538 neededvector[i] = 0;
3539 for(i=1;i<=noknots;i++)
3540 neededvector[nodepart[i]] += 1;
3541
3542 if(!noopt) {
3543 optimize = 1;
3544 probnodes = Ivector(1,noknots);
3545 for(i=1;i<=noknots;i++)
3546 probnodes[i] = 0;
3547 printf("Applying aggressive optimization for load balancing\n");
3548 }
3549 else {
3550 optimize = 1;
3551 }
3552
3553 optimizeownership:
3554
3555 dshared = CheckSharedDeviation(neededvector,partitions,info);
3556
3557 if(!noopt) {
3558
3559 /* Distribute the shared nodes as evenly as possible. */
3560
3561 int nochanges, maxrounds, dtarget;
3562
3563 maxrounds = 5;
3564 dtarget = 3;
3565 nochanges = 0;
3566
3567 for(n=1;n<=maxrounds;n++) {
3568 nochanges = 0;
3569
3570 for(i=1;i<=noknots;i++) {
3571 ind = i;
3572
3573 /* owner partition may only be changed it there a are two of them */
3574 k = data->partitiontable[2][ind];
3575 if(!k) continue;
3576
3577 /* only apply the switch to cases with exactly two partitions
3578 to avoid the nasty multiply coupled nodes. */
3579 if(maxneededtimes > 2)
3580 if(data->partitiontable[3][ind]) continue;
3581
3582 /* Do not change the ownership of nodes that cause topological problems */
3583 if(probnodes[ind]) continue;
3584
3585 j = data->partitiontable[1][ind];
3586
3587 /* Switch the owner to the smaller owner group if possible */
3588 if(abs(neededvector[j] - neededvector[k]) < dtarget ) continue;
3589
3590 if(neededvector[j] < neededvector[k] && nodepart[ind] == k) {
3591 nochanges++;
3592 neededvector[j] += 1;
3593 neededvector[k] -= 1;
3594 nodepart[ind] = j;
3595 }
3596 else if(neededvector[k] < neededvector[j] && nodepart[ind] == j) {
3597 nochanges++;
3598 neededvector[k] += 1;
3599 neededvector[j] -= 1;
3600 nodepart[ind] = k;
3601 }
3602 }
3603 if(info && nochanges) printf("Changed the ownership of %d nodes\n",nochanges);
3604
3605 dshared0 = dshared;
3606 dshared = CheckSharedDeviation(neededvector,partitions,info);
3607
3608 if(dshared >= dshared0) break;
3609 somethingdone = somethingdone + nochanges;
3610 }
3611 }
3612
3613
3614 /* optimizesharing: */
3615
3616 if(info) printf("Checking for problematic sharings\n");
3617 m = 0;
3618
3619 if(partitions > 2) {
3620 if(info) printf("Optimizing sharing for %d partitions\n", partitions);
3621
3622 do {
3623
3624 int i1,i2,e1,e2,owners;
3625
3626 m++;
3627 sharings = 0;
3628 i1 = i2 = 0;
3629 e1 = e2 = 0;
3630
3631 /* There cannot be more partitions sharing a node than there are
3632 node in the elements. Elementtype 827 has 27 nodes. */
3633 maxshared = 27;
3634 if(m == 1 && optimize == 1) {
3635 elemparts = Ivector(1,partitions);
3636 invelemparts = Ivector(1,maxshared);
3637 knows = Imatrix(1,maxshared,1,maxshared);
3638 }
3639
3640 for(j=1;j<=partitions;j++)
3641 elemparts[j] = 0;
3642
3643 for(j=1;j<=maxshared;j++)
3644 invelemparts[j] = 0;
3645
3646 for(j=1;j<=maxshared;j++)
3647 for(k=1;k<=maxshared;k++)
3648 knows[j][k] = 0;
3649
3650 for(i=1;i<=noelements;i++) {
3651 int ownpart;
3652
3653 owners = 0;
3654 nodesd2 = data->elementtypes[i] % 100;
3655 ownpart = FALSE;
3656
3657 /* Check the number of owners in an element */
3658 for(j=0;j < nodesd2;j++) {
3659 ind = data->topology[i][j];
3660 k = nodepart[ind];
3661
3662 /* Mark if the element partition is one of the owners */
3663 if( k == elempart[i]) ownpart = TRUE;
3664 if(!elemparts[k]) {
3665 owners++;
3666 elemparts[k] = owners;
3667 invelemparts[owners] = k;
3668 }
3669 }
3670
3671 /* One strange owner is still ok. */
3672 if(owners - ownpart <= 1) {
3673 /* Nullify the elemparts vector */
3674 for(j=1;j<=owners;j++) {
3675 k = invelemparts[j];
3676 elemparts[k] = 0;
3677 }
3678 continue;
3679 }
3680
3681 /* Check which of the partitions are related by a common node */
3682 for(j=0;j < nodesd2;j++) {
3683 ind = data->topology[i][j];
3684 for(l=1;l<=maxneededtimes;l++) {
3685 e1 = data->partitiontable[l][ind];
3686 if(!e1) break;
3687 e1 = elemparts[e1];
3688 if(!e1) continue;
3689 for(k=l+1;k<=maxneededtimes;k++) {
3690 e2 = data->partitiontable[k][ind];
3691 if(!e2) break;
3692 e2 = elemparts[e2];
3693 if(!e2) continue;
3694 knows[e1][e2] = knows[e2][e1] = TRUE;
3695 }
3696 }
3697 }
3698
3699 /* Check if there are more complex relations:
3700 i.e. two partitions are joined at an element but not at the same node. */
3701 hit = FALSE;
3702 for(j=1;j<=owners;j++)
3703 for(k=j+1;k<=owners;k++) {
3704 if(!knows[j][k]) {
3705 hit += 1;
3706 i1 = invelemparts[j];
3707 i2 = invelemparts[k];
3708 }
3709 }
3710 /* Nullify the elemparts vector */
3711 for(j=1;j<=owners;j++) {
3712 k = invelemparts[j];
3713 elemparts[k] = 0;
3714 }
3715
3716 /* Nullify the knows matrix */
3717 for(j=1;j <= owners;j++)
3718 for(k=1;k <= owners;k++)
3719 knows[j][k] = FALSE;
3720
3721 if(hit) {
3722 e1 = e2 = 0;
3723
3724 if(info) {
3725 if( hit + m > 2 ) printf("Partitions %d and %d in element %d (%d owners) oddly related %d times\n",
3726 i1,i2,i,owners,hit);
3727 }
3728
3729 /* Count the number of nodes with wrong parents */
3730 for(j=0;j < nodesd2;j++) {
3731 ind = data->topology[i][j];
3732 for(l=1;l<=maxneededtimes;l++) {
3733 k = data->partitiontable[l][ind];
3734 if(k == 0) break;
3735 if(k == i1) e1++;
3736 if(k == i2) e2++;
3737 }
3738 }
3739
3740 /* Change the owner of those with less sharings */
3741 for(j=0;j < nodesd2;j++) {
3742 ind = data->topology[i][j];
3743 k = nodepart[ind];
3744 if((k == i1 && e1 < e2) || (k == i2 && e1 >= e2)) {
3745 if(!noopt) probnodes[ind] += 1;
3746 nodepart[ind] = elempart[i];
3747 neededvector[elempart[i]] += 1;
3748 neededvector[k] -= 1;
3749 }
3750 }
3751 sharings++;
3752 }
3753 }
3754
3755 somethingdone += sharings;
3756 if(info && sharings) printf("Changed the ownership of %d nodes\n",sharings);
3757
3758 } while (sharings > 0 && m < 3);
3759
3760 if(info) {
3761 if(sharings)
3762 printf("%d problematic sharings may still exist\n",sharings);
3763 else
3764 printf("There shouldn't be any problematic sharings, knock, knock...\n");
3765 }
3766 }
3767
3768 /* This seems to work also iteratively */
3769 if(!noopt && m+n > 10 && optimize < 50) {
3770 optimize++;
3771 printf("Performing ownership optimization round %d\n",optimize);
3772 goto optimizeownership;
3773 }
3774
3775 free_Ivector(neededvector,1,partitions);
3776
3777 if(!noopt) {
3778 free_Ivector(probnodes,1,noknots);
3779 if(partitions > 2) {
3780 free_Ivector(elemparts,1,partitions);
3781 free_Ivector(invelemparts,1,maxshared);
3782 free_Imatrix(knows,1,maxshared,1,maxshared);
3783 }
3784 }
3785
3786
3787 if(somethingdone) {
3788 if( info ) {
3789 printf("The partitioning was optimized: %d\n",somethingdone);
3790 printf("Checking partitioning after optimization\n");
3791 }
3792 CheckPartitioning(data,info);
3793 }
3794 else {
3795 if(info) printf("Partitioning was not altered\n");
3796 }
3797
3798 return(0);
3799 }
3800
3801
SaveElmerInputPartitioned(struct FemType * data,struct BoundaryType * bound,char * prefix,int decimals,int * parthalo,int indirect,int parthypre,int subparts,int nooverwrite,int info)3802 int SaveElmerInputPartitioned(struct FemType *data,struct BoundaryType *bound,
3803 char *prefix,int decimals,int *parthalo,int indirect,
3804 int parthypre,int subparts,int nooverwrite, int info)
3805 /* Saves the mesh in a form that may be used as input
3806 in Elmer calculations in parallel platforms.
3807 */
3808 {
3809 int noknots,noelements,sumsides,partitions,hit,found,parent,parent2,maxnosides;
3810 int nodesd2,nodesd1,discont,maxelemtype,minelemtype,sidehits,elemsides,side,bctype;
3811 int part,otherpart,part2,part3,elemtype,sideelemtype,*needednodes,*neededtwice;
3812 int **bulktypes,*sidetypes,tottypes,splitsides;
3813 int i,j,k,l,l2,l3,m,n,ind,ind2,sideind[MAXNODESD1],elemhit[MAXNODESD2];
3814 char filename[MAXFILESIZE],outstyle[MAXFILESIZE];
3815 char directoryname[MAXFILESIZE],subdirectoryname[MAXFILESIZE];
3816 int *neededtimes,*elempart,*elementsinpart,*indirectinpart,*sidesinpart;
3817 int maxneededtimes,indirecttype,bcneeded,trueparent,trueparent2,*ownerpart;
3818 int *sharednodes,*ownnodes,reorder,*order=NULL,*invorder=NULL;
3819 int *bcnode,*bcnodedummy,*elementhalo,*neededtimes2;
3820 int partstart,partfin,filesetsize,nofile,nofile2,nobcnodes,hasbcnodes,halomode;
3821 int halobulkelems,halobcs,savethis,fail=0,cdstat,immersed,halocopies,anyparthalo;
3822
3823 FILE *out,*outfiles[MAXPARTITIONS+1];
3824 int sumelementsinpart,sumownnodes,sumsharednodes,sumsidesinpart,sumindirect;
3825
3826
3827
3828 if(info) {
3829 printf("Saving Elmer mesh in partitioned format\n");
3830 if( subparts ) printf("There are %d subpartitions\n",subparts);
3831 }
3832
3833 if(!data->created) {
3834 printf("You tried to save points that were never created.\n");
3835 bigerror("No Elmer mesh files saved!");
3836 }
3837
3838 partitions = data->nopartitions;
3839 if(!partitions) {
3840 printf("Tried to save partiotioned format without partitions!\n");
3841 bigerror("No Elmer mesh files saved!");
3842 }
3843
3844 anyparthalo = FALSE;
3845 for(i=0;i<MAXHALOMODES;i++) {
3846 if( parthalo[i] ) {
3847 printf("Saving halo elements mode active: %d\n",i);
3848 anyparthalo = TRUE;
3849 }
3850 }
3851
3852 if( subparts < 1 && parthalo[3] ) {
3853 printf("There can be no layer halo since there are no layers!\n");
3854 bigerror("No Elmer mesh files saved!");
3855 }
3856
3857
3858 elempart = data->elempart;
3859 ownerpart = data->nodepart;
3860 noelements = data->noelements;
3861 noknots = data->noknots;
3862
3863 minelemtype = 101;
3864 maxelemtype = GetMaxElementType(data);
3865 indirecttype = 0;
3866 halobulkelems = 0;
3867
3868 needednodes = Ivector(1,partitions);
3869 neededtwice = Ivector(1,partitions);
3870 sharednodes = Ivector(1,partitions);
3871 ownnodes = Ivector(1,partitions);
3872 sidetypes = Ivector(minelemtype,maxelemtype);
3873 bulktypes = Imatrix(1,partitions,minelemtype,maxelemtype);
3874 bcnodedummy = Ivector(1,noknots);
3875
3876 /* Order the nodes so that the different partitions have a continuous interval of nodes.
3877 This information is used only just before the saving of node indexes in each instance.
3878 This feature was coded for collaboration with Hypre library that assumes this. */
3879 reorder = parthypre;
3880 if(reorder) {
3881 order = Ivector(1,noknots);
3882 k = 0;
3883 for(j=1;j<=partitions;j++)
3884 for(i=1; i <= noknots; i++)
3885 if(ownerpart[i] == j) {
3886 k++;
3887 order[i] = k;
3888 }
3889 invorder = Ivector(1,noknots);
3890 for(i=1;i<=noknots;i++)
3891 invorder[order[i]] = i;
3892 if(info) printf("Created continuous parallel ordering\n");
3893 }
3894
3895 /* Mark the nodes that are on some boundary in order to create boundary halos */
3896 bcnode = Ivector(1,noknots);
3897 for(i=1;i<=noknots;i++)
3898 bcnode[i] = FALSE;
3899 for(j=0;j < MAXBOUNDARIES;j++) {
3900 for(i=1; i <= bound[j].nosides; i++) {
3901 GetBoundaryElement(i,&bound[j],data,sideind,&sideelemtype);
3902 nodesd1 = sideelemtype%100;
3903 for(l=0;l<nodesd1;l++)
3904 bcnode[sideind[l]] = 1;
3905 }
3906 }
3907 nobcnodes = 0;
3908 for(i=1;i<=noknots;i++) {
3909 if( bcnode[i] == 1) nobcnodes += 1;
3910 }
3911 if(info) printf("Number of boundary nodes at the boundary: %d\n",nobcnodes);
3912
3913
3914 for(i=1;i<=noknots;i++)
3915 bcnodedummy[i] = 0;
3916 for(i=1;i<=noelements;i++) {
3917 k = data->material[i];
3918 elemtype = data->elementtypes[i];
3919 nodesd2 = elemtype%100;
3920
3921 for(j=0;j < nodesd2;j++) {
3922 ind = data->topology[i][j];
3923 if( bcnode[ind] == 1) continue;
3924 if( bcnodedummy[ind] ) {
3925 if( bcnodedummy[ind] != k ) bcnode[ind] = 2;
3926 } else {
3927 bcnodedummy[ind] = k;
3928 }
3929 }
3930 }
3931 nobcnodes = 0;
3932 for(i=1;i<=noknots;i++) {
3933 if( bcnode[i] == 2) nobcnodes += 1;
3934 }
3935 if(info) printf("Number of additional interface nodes: %d\n",nobcnodes);
3936
3937
3938 sprintf(directoryname,"%s",prefix);
3939 sprintf(subdirectoryname,"%s.%d","partitioning",partitions);
3940
3941
3942 #ifdef MINGW32
3943 mkdir(directoryname);
3944 #else
3945 fail = mkdir(directoryname,0750);
3946 #endif
3947 if(info && !fail) printf("Created mesh directory: %s\n",directoryname);
3948 fail = chdir(directoryname);
3949 #ifdef MINGW32
3950 fail = mkdir(subdirectoryname);
3951 #else
3952 fail = mkdir(subdirectoryname,0750);
3953 #endif
3954 if(fail) {
3955 if(info) printf("Reusing existing subdirectory: %s\n",subdirectoryname);
3956 if(nooverwrite) {
3957 printf("Mesh seems to already exist, writing is cancelled!\n");
3958 return(0);
3959 }
3960 }
3961 else {
3962 if(info) printf("Created subdirectory: %s\n",subdirectoryname);
3963 }
3964 cdstat = chdir(subdirectoryname);
3965
3966
3967
3968 if(info) printf("Saving mesh in parallel ElmerSolver format to directory %s/%s.\n",
3969 directoryname,subdirectoryname);
3970
3971 filesetsize = MAXPARTITIONS;
3972 if(partitions > filesetsize)
3973 if(info) printf("Saving %d partitions in maximum sets of %d\n",partitions,filesetsize);
3974
3975 elementsinpart = Ivector(1,partitions);
3976 indirectinpart = Ivector(1,partitions);
3977 sidesinpart = Ivector(1,partitions);
3978 elementhalo = Ivector(1,partitions);
3979 for(i=1;i<=partitions;i++)
3980 elementsinpart[i] = indirectinpart[i] = sidesinpart[i] = elementhalo[i] = 0;
3981
3982 for(j=1;j<=partitions;j++)
3983 for(i=minelemtype;i<=maxelemtype;i++)
3984 bulktypes[j][i] = 0;
3985
3986 /* Compute how many times a node may be needed at maximum */
3987 maxneededtimes = data->maxpartitiontable;
3988 neededtimes = Ivector(1,noknots);
3989 for(i=1;i<=noknots;i++) {
3990 neededtimes[i] = 0;
3991 for(j=1;j<=maxneededtimes;j++)
3992 if(data->partitiontable[j][i]) neededtimes[i] += 1;
3993 }
3994 if(info) printf("Nodes belong to %d partitions in maximum\n",maxneededtimes);
3995
3996
3997 /*********** part.n.elements *********************/
3998 /* Save elements in all partitions and where they are needed */
3999
4000
4001 partstart = 1;
4002 partfin = MIN( partitions, filesetsize );
4003
4004 next_elements_set:
4005
4006 for(part=partstart;part<=partfin;part++) {
4007 sprintf(filename,"%s.%d.%s","part",part,"elements");
4008 nofile = part - partstart + 1;
4009 outfiles[nofile] = fopen(filename,"w");
4010 }
4011
4012 for(i=1;i<=noelements;i++) {
4013 part = elempart[i];
4014
4015 elemtype = data->elementtypes[i];
4016 nodesd2 = elemtype%100;
4017
4018 if(part >= partstart && part <= partfin) {
4019 nofile = part - partstart + 1;
4020 bulktypes[part][elemtype] += 1;
4021 elementsinpart[part] += 1;
4022
4023 fprintf(outfiles[nofile],"%d %d %d ",i,data->material[i],elemtype);
4024
4025 for(j=0;j < nodesd2;j++) {
4026 ind = data->topology[i][j];
4027 if(reorder) ind = order[ind];
4028 fprintf(outfiles[nofile],"%d ",ind);
4029 }
4030 fprintf(outfiles[nofile],"\n");
4031 }
4032
4033
4034 if(!anyparthalo) continue;
4035 halocopies = 0;
4036
4037
4038
4039 for( halomode = 1; halomode <= 4; halomode++) {
4040
4041 if(!parthalo[halomode]) continue;
4042
4043
4044 /* This creates a simple halo when the elements have been partitioned such
4045 that they are in number of "subparts" intervals of z coordinate. */
4046 if( halomode == 3 ) {
4047 if( part <= subparts ) {
4048 int leftright;
4049 for(leftright=-1;leftright <=1;leftright += 2) {
4050 part2 = part+leftright;
4051 if( part2 < 1 || part2 > subparts ) continue;
4052 elementhalo[part2] = i;
4053 halocopies += 1;
4054 }
4055 }
4056 }
4057
4058 /* For the following halos we study the shared and boundary nodes. */
4059 otherpart = 0;
4060 hasbcnodes = FALSE;
4061 for(j=0;j < nodesd2;j++) {
4062 ind = data->topology[i][j];
4063 if(neededtimes[ind] > 1) otherpart++;
4064 if(bcnode[ind]) hasbcnodes = TRUE;
4065 }
4066
4067 /* Classical halo needed by discontinuous Galerkin.
4068 Elements that are attached by joined face element are added to halo */
4069 if( halomode == 1 && otherpart ) {
4070 /* If the saving of halo is requested check it for elements which have at least
4071 two nodes in shared partitions. First make this quick test. */
4072
4073 elemsides = elemtype / 100;
4074 if(elemsides == 8) {
4075 immersed = (otherpart >= 4);
4076 elemsides = 6;
4077 }
4078 else if(elemsides == 7) {
4079 immersed = (otherpart >= 3);
4080 elemsides = 5;
4081 }
4082 else if(elemsides == 6) {
4083 immersed = (otherpart >= 3);
4084 elemsides = 5;
4085 }
4086 else if(elemsides == 5) {
4087 immersed = (otherpart >= 3);
4088 elemsides = 4;
4089 }
4090 else {
4091 immersed = (otherpart >= 2);
4092 }
4093
4094 if( immersed ) {
4095
4096 /* In order for the halo to be present the element should have a boundary
4097 fully immersed in the other partition. This test takes more time. */
4098
4099 for(side=0;side<elemsides;side++) {
4100
4101 GetElementSide(i,side,1,data,&sideind[0],&sideelemtype);
4102
4103 /* Because every node must be on the boundary use the 1st index as the
4104 first test */
4105 for(l=1;l<=neededtimes[sideind[0]];l++) {
4106 part2 = data->partitiontable[l][sideind[0]];
4107
4108 /* We did already save this in partition part */
4109 if(part2 == part) continue;
4110 if(elementhalo[part2] == i) continue;
4111 if(part2 < partstart || part2 > partfin) continue;
4112
4113 /* Now check that also the other nodes on the face are at the interface */
4114 sidehits = 1;
4115 for(k=1;k<sideelemtype%100;k++) {
4116 for(l2=1;l2<=neededtimes[sideind[k]];l2++) {
4117 part3 = data->partitiontable[l2][sideind[k]];
4118 if(part2 == part3) sidehits++;
4119 }
4120 }
4121 if( sidehits < sideelemtype % 100 ) continue;
4122
4123 if(0) printf("Adding halo for partition %d and element %d\n",part2,i);
4124
4125 /* Remember that this element is saved for this partition */
4126 elementhalo[part2] = i;
4127 halocopies += 1;
4128 }
4129 }
4130 }
4131 }
4132
4133 /* This creates a halo that includes all elements attached to boundary or interface nodes */
4134 if( halomode == 2 && hasbcnodes ) {
4135 for(j=0;j < nodesd2;j++) {
4136 ind = data->topology[i][j];
4137
4138 /* If the node is not on the boundary then cycle */
4139 if( !bcnode[ind] ) continue;
4140
4141 /* Check to which partitions this element is associated with */
4142 for(k=1;k<=neededtimes[ind];k++) {
4143 part2 = data->partitiontable[k][ind];
4144
4145 /* This element is already saved to its primary partition */
4146 if( part == part2 ) continue;
4147 if( elementhalo[part2] == i) continue;
4148 if( part2 < partstart || part2 > partfin ) continue;
4149
4150 if(0) printf("Adding bc halo for partition %d and element %d\n",part2,i);
4151
4152 /* Remember that this element is saved for this partition */
4153 elementhalo[part2] = i;
4154 halocopies += 1;
4155 }
4156 }
4157 }
4158
4159 /* This greedy routine makes the halo even if there is just one node in the shared boundary.
4160 Currently not active. */
4161 if( halomode == 4 && otherpart ) {
4162 for(j=0;j < nodesd2;j++) {
4163 ind = data->topology[i][j];
4164 if(neededtimes[ind] == 1) continue;
4165
4166 for(l=1;l<=neededtimes[ind];l++) {
4167 part2 = data->partitiontable[l][ind];
4168
4169 /* We did already save this in partition part */
4170 if(part2 == part) continue;
4171 if( elementhalo[part2] == i) continue;
4172 if(part2 < partstart || part2 > partfin) continue;
4173
4174 /* Remember that this element is saved for this partition */
4175 elementhalo[part2] = i;
4176 halocopies += 1;
4177 }
4178 }
4179 }
4180 }
4181
4182
4183 if( halocopies ) {
4184 halobulkelems += halocopies;
4185
4186 for(part2=partstart; part2 <= partfin; part2++) {
4187 if( part == part2) continue;
4188 if( elementhalo[part2] != i ) continue;
4189
4190 nofile2 = part2 - partstart + 1;
4191 fprintf(outfiles[nofile2],"%d/%d %d %d ",i,part,data->material[i],elemtype);
4192
4193 for(j=0;j < nodesd2;j++) {
4194 ind = data->topology[i][j];
4195 if(reorder) ind = order[ind];
4196 fprintf(outfiles[nofile2],"%d ",ind);
4197 }
4198 fprintf(outfiles[nofile2],"\n");
4199
4200 bulktypes[part2][elemtype] += 1;
4201 elementsinpart[part2] += 1;
4202
4203 /* Add the halo on-the-fly to the partitiontable of the nodes */
4204 for(j=0;j < nodesd2;j++) {
4205 ind = data->topology[i][j];
4206 hit = FALSE;
4207 for(k=1;k<=maxneededtimes;k++) {
4208 part3 = data->partitiontable[k][ind];
4209 if(part3 == part2) hit = TRUE;
4210 if(!part3) break;
4211 }
4212 if(!hit) {
4213 if(k <= maxneededtimes) {
4214 data->partitiontable[k][ind] = part2;
4215 }
4216 else {
4217 maxneededtimes++;
4218 if(0) printf("Allocating new column %d in partitiontable\n",k);
4219 data->partitiontable[k] = Ivector(1,noknots);
4220 for(m=1;m<=noknots;m++)
4221 data->partitiontable[k][m] = 0;
4222 data->partitiontable[k][ind] = part2;
4223 }
4224 }
4225 }
4226 }
4227 } /* halocopies */
4228
4229 } /* noelements */
4230
4231 if( anyparthalo ) {
4232 if(info) printf("There are %d bulk elements in the halo.\n",halobulkelems);
4233 }
4234
4235
4236 for(part=partstart;part<=partfin;part++) {
4237 nofile = part - partstart + 1;
4238 fclose(outfiles[nofile]);
4239 }
4240 if(partfin < partitions) {
4241 partstart = partfin + 1;
4242 partfin = MIN( partfin + filesetsize, partitions);
4243 goto next_elements_set;
4244 }
4245 /* part.n.elements saved */
4246
4247
4248 /* The partitiontable has been changed to include the halo elements. The need for saving the
4249 halo nodes may be checked by looking whether the number of how many partitions needs
4250 the element has changed. */
4251 if(anyparthalo) {
4252 int halonodes;
4253 neededtimes2 = Ivector(1,noknots);
4254 halonodes = 0;
4255
4256 k = 0;
4257 l = 0;
4258 for(i=1;i<=noknots;i++) {
4259 neededtimes2[i] = 0;
4260 for(j=1;j<=maxneededtimes;j++)
4261 if(data->partitiontable[j][i]) {
4262 neededtimes2[i] += 1;
4263 k++;
4264 }
4265 halonodes += neededtimes2[i] - neededtimes[i];
4266 l += neededtimes[i];
4267 }
4268 if(data->maxpartitiontable < maxneededtimes) {
4269 data->maxpartitiontable = maxneededtimes;
4270 if(info) printf("With the halos nodes belong to %d partitions in maximum\n",maxneededtimes);
4271 }
4272 if(info) printf("There are %d additional shared nodes resulting from the halo.\n",halonodes);
4273 }
4274 else {
4275 neededtimes2 = neededtimes;
4276 }
4277
4278 /* Define new BC numbers for indirect connections. These should not be mixed with
4279 existing BCs as they only serve the purpose of automatically creating the matrix structure. */
4280 if(indirect) {
4281 indirecttype = 0;
4282 for(j=0;j < MAXBOUNDARIES;j++)
4283 for(i=1; i <= bound[j].nosides; i++)
4284 if(bound[j].types[i] > indirecttype) indirecttype = bound[j].types[i];
4285 indirecttype++;
4286 if(info) printf("Indirect connections given index %d and elementtype 102.\n",indirecttype);
4287 }
4288
4289 /* The output format is the same for all partitions */
4290 if(data->dim == 2)
4291 sprintf(outstyle,"%%d %%d %%.%dg %%.%dg 0.0\n",decimals,decimals);
4292 else
4293 sprintf(outstyle,"%%d %%d %%.%dg %%.%dg %%.%dg\n",decimals,decimals,decimals);
4294
4295 if(info) printf("Saving mesh for %d partitions\n",partitions);
4296
4297
4298 /*********** part.n.nodes *********************/
4299
4300 partstart = 1;
4301 partfin = MIN( partitions, filesetsize);
4302
4303 for(i=1;i<=partitions;i++) {
4304 needednodes[i] = 0;
4305 neededtwice[i] = 0;
4306 sharednodes[i] = 0;
4307 ownnodes[i] = 0;
4308 }
4309
4310 next_nodes_set:
4311
4312 for(part=partstart;part<=partfin;part++) {
4313 sprintf(filename,"%s.%d.%s","part",part,"nodes");
4314 nofile = part - partstart + 1;
4315 outfiles[nofile] = fopen(filename,"w");
4316 }
4317
4318 for(l=1; l <= noknots; l++) {
4319 i = l;
4320 if(reorder) i=invorder[l];
4321
4322 /* for(j=1;j<=neededtimes2[i];j++) { */
4323 for(j=1;j<=maxneededtimes;j++) {
4324
4325 k = data->partitiontable[j][i];
4326 if(!k) break;
4327
4328 if(k < partstart || k > partfin) continue;
4329 nofile = k - partstart + 1;
4330
4331 ind = i;
4332 if(reorder) ind=order[i];
4333
4334 if(data->dim == 2)
4335 fprintf(outfiles[nofile],outstyle,ind,-1,data->x[i],data->y[i]);
4336 else if(data->dim == 3)
4337 fprintf(outfiles[nofile],outstyle,ind,-1,data->x[i],data->y[i],data->z[i]);
4338
4339 needednodes[k] += 1;
4340 if(k == ownerpart[i])
4341 ownnodes[k] += 1;
4342 else
4343 sharednodes[k] += 1;
4344 }
4345 }
4346
4347 for(part=partstart;part<=partfin;part++) {
4348 nofile = part - partstart + 1;
4349 fclose(outfiles[nofile]);
4350 }
4351 if(partfin < partitions) {
4352 partstart = partfin + 1;
4353 partfin = MIN( partfin + filesetsize, partitions);
4354 goto next_nodes_set;
4355 }
4356 /* part.n.nodes saved */
4357
4358
4359 /*********** part.n.shared *********************/
4360
4361 partstart = 1;
4362 partfin = MIN( partitions, filesetsize );
4363
4364 next_shared_set:
4365
4366 for(part=partstart;part<=partfin;part++) {
4367 sprintf(filename,"%s.%d.%s","part",part,"shared");
4368 nofile = part - partstart + 1;
4369 outfiles[nofile] = fopen(filename,"w");
4370 }
4371
4372 for(l=1; l <= noknots; l++) {
4373 i = l;
4374 if(reorder) i = invorder[l];
4375
4376 if(neededtimes2[i] <= 1) continue;
4377
4378 for(j=1;j<=neededtimes2[i];j++) {
4379 k = data->partitiontable[j][i];
4380
4381 /* if( parthalo[2] && neededtimes[i] == 1 && ownerpart[i] == k ) continue; */
4382
4383 if(k < partstart || k > partfin) continue;
4384 nofile = k - partstart + 1;
4385
4386 ind = i;
4387 if(reorder) ind = order[i];
4388 neededtwice[k] += 1;
4389
4390 /* if( parthalo[2] ) {
4391 fprintf(outfiles[nofile],"%d %d %d",ind,MAX(1,neededtimes[i]),ownerpart[i]);
4392 for(m=1;m<=neededtimes[i];m++)
4393 if(data->partitiontable[m][i] != ownerpart[i]) fprintf(outfiles[nofile]," %d",data->partitiontable[m][i]);
4394 fprintf(outfiles[nofile],"\n");
4395 }
4396 else { */
4397 fprintf(outfiles[nofile],"%d %d %d",ind,neededtimes2[i],ownerpart[i]);
4398 for(m=1;m<=neededtimes2[i];m++)
4399 if(data->partitiontable[m][i] != ownerpart[i]) fprintf(outfiles[nofile]," %d",data->partitiontable[m][i]);
4400 fprintf(outfiles[nofile],"\n");
4401 /* } */
4402 }
4403 }
4404
4405 for(part=partstart;part<=partfin;part++) {
4406 nofile = part - partstart + 1;
4407 fclose(outfiles[nofile]);
4408 }
4409 if(partfin < partitions) {
4410 partstart = partfin + 1;
4411 partfin = MIN( partfin + filesetsize, partitions);
4412 goto next_shared_set;
4413 }
4414 /* part.n.shared saved */
4415
4416
4417
4418
4419
4420 /*********** part.n.boundary *********************/
4421 /* This is still done in partition loop as the subroutines are quite complicated */
4422 discont = FALSE;
4423 splitsides = 0;
4424
4425 maxnosides = 0;
4426 for(j=0;j < MAXBOUNDARIES;j++)
4427 maxnosides = MAX( maxnosides, bound[j].nosides );
4428
4429 halobcs = 0;
4430 for(part=1;part<=partitions;part++) {
4431 int bcneeded2,closeparent,closeparent2,haloelem,saveelem;
4432
4433 sprintf(filename,"%s.%d.%s","part",part,"boundary");
4434 out = fopen(filename,"w");
4435
4436 for(i=minelemtype;i<=maxelemtype;i++)
4437 sidetypes[i] = 0;
4438
4439 sumsides = 0;
4440
4441 /* Loop over boundaries and boundary elements therein. */
4442 for(j=0;j < MAXBOUNDARIES;j++) {
4443
4444 if(bound[j].nosides == 0 ) continue;
4445
4446 /* Normal boundary conditions */
4447 for(i=1; i <= bound[j].nosides; i++) {
4448
4449 GetBoundaryElement(i,&bound[j],data,sideind,&sideelemtype);
4450
4451 bctype = bound[j].types[i];
4452 nodesd1 = sideelemtype%100;
4453
4454 parent = bound[j].parent[i];
4455 parent2 = bound[j].parent2[i];
4456
4457
4458 /* Check whether the side is such that it belongs to the domain.
4459 Then it will be always saved - no matter of the halos. */
4460 trueparent = trueparent2 = FALSE;
4461 if( parent ) trueparent = (elempart[parent] == part);
4462 if( parent2 ) trueparent2 = (elempart[parent2] == part);
4463
4464 saveelem = FALSE;
4465 part2 = 0;
4466
4467 if(trueparent || trueparent2) {
4468 /* Either parent must be associated with this partition, otherwise do not save this (except for halo nodes) */
4469 if( parent && !trueparent ) {
4470 splitsides++;
4471 /* For halo elements also the 2nd true parent exists by construction. */
4472 if(!parthalo[1] && !parthalo[2] && !parthalo[4]) parent = 0;
4473 }
4474 else if( parent2 && !trueparent2 ) {
4475 splitsides++;
4476 if(!parthalo[1] && !parthalo[2] && !parthalo[4]) parent2 = 0;
4477 }
4478 saveelem = TRUE;
4479 }
4480
4481 /* If boundary element is not saved by default then study different halos.
4482 First look for the extruded halo and offset +/-1 in partition index. */
4483 if(!saveelem && parthalo[3] && part <= subparts ) {
4484 int leftright;
4485 for(leftright=-1;leftright <=1;leftright += 2) {
4486 part2 = part+leftright;
4487 if( part2 < 1 || part2 > subparts ) continue;
4488 if( parent ) trueparent = (elempart[parent] == part2);
4489 if( parent2 ) trueparent2 = (elempart[parent2] == part2);
4490
4491 /* If other parent is not true parent then set the parent to zero because it is not
4492 in this partition. */
4493 if(trueparent || trueparent2) {
4494 if( parent && !trueparent )
4495 parent = 0;
4496 else if( parent2 && !trueparent2 )
4497 parent2 = 0;
4498 saveelem = TRUE;
4499 break;
4500 }
4501 }
4502 if(saveelem) halobcs += 1;
4503 }
4504
4505 /* These halo strategies save elements even with just one hit on the original partition
4506 since the aggressive halos ensure that all necessary bulk elements are saved. */
4507 if(!saveelem && (parthalo[2] || parthalo[4]) ) {
4508 bcneeded = 0;
4509 for(l=0;l<nodesd1;l++) {
4510 ind = sideind[l];
4511 for(k=1;k<=neededtimes[ind];k++)
4512 if(part == data->partitiontable[k][ind]) bcneeded++;
4513 }
4514 saveelem = bcneeded;
4515 if(saveelem) halobcs +=1;
4516 }
4517
4518 /* If we have nothing to save then just continue to next boundary element. */
4519 if(!saveelem) continue;
4520
4521 sumsides++;
4522 sidetypes[sideelemtype] += 1;
4523
4524 if( part2 )
4525 fprintf(out,"%d/%d %d %d %d %d",
4526 sumsides,part2,bctype,parent,parent2,sideelemtype);
4527 else
4528 fprintf(out,"%d %d %d %d %d",
4529 sumsides,bctype,parent,parent2,sideelemtype);
4530
4531 if(reorder) {
4532 for(l=0;l<nodesd1;l++)
4533 fprintf(out," %d",order[sideind[l]]);
4534 } else {
4535 for(l=0;l<nodesd1;l++)
4536 fprintf(out," %d",sideind[l]);
4537 }
4538 fprintf(out,"\n");
4539 }
4540 }
4541
4542 /* The discontinuous stuff is more or less obsolete as these things can now be made also
4543 within ElmerSolver. */
4544 /* The second side for discontinuous boundary conditions.
4545 Note that this has not been treated for orphan control. */
4546 for(j=0;j < MAXBOUNDARIES;j++) {
4547 for(i=1; i <= bound[j].nosides; i++) {
4548 if(bound[j].ediscont)
4549 discont = bound[j].discont[i];
4550 else
4551 discont = FALSE;
4552
4553 if(!bound[j].parent2[i] || !discont) continue;
4554
4555 GetElementSide(bound[j].parent2[i],bound[j].side2[i],-bound[j].normal[i],
4556 data,sideind,&sideelemtype);
4557 nodesd1 = sideelemtype%100;
4558
4559 bcneeded = 0;
4560 for(l=0;l<nodesd1;l++) {
4561 ind = sideind[l];
4562 for(k=1;k<=neededtimes[ind];k++)
4563 if(part == data->partitiontable[k][ind]) bcneeded++;
4564 }
4565 if(bcneeded < nodesd1) continue;
4566
4567
4568 trueparent = (elempart[bound[j].parent2[i]] == part);
4569 if(!trueparent) continue;
4570
4571 sumsides++;
4572 fprintf(out,"%d %d %d %d ",
4573 sumsides,bound[j].types[i],bound[j].parent2[i],bound[j].parent[i]);
4574
4575 fprintf(out,"%d ",sideelemtype);
4576 sidetypes[sideelemtype] += 1;
4577 if(reorder) {
4578 for(l=0;l<nodesd1;l++)
4579 fprintf(out,"%d ",order[sideind[l]]);
4580 }
4581 else {
4582 for(l=0;l<nodesd1;l++)
4583 fprintf(out,"%d ",sideind[l]);
4584 }
4585 fprintf(out,"\n");
4586 }
4587 }
4588 sidesinpart[part] = sumsides;
4589
4590
4591 /* Boundary nodes that express indirect couplings between different partitions.
4592 This makes it possible for ElmerSolver to create a matrix connection that
4593 is known to exist. */
4594
4595 if (indirect) {
4596 int maxsides,nodesides,maxnodeconnections,connectednodes,m;
4597 int **nodepairs=NULL,*nodeconnections,**indpairs;
4598
4599 nodeconnections = bcnodedummy;
4600 l = 0;
4601 maxsides = 0;
4602 nodesides = 0;
4603
4604 findindirect:
4605
4606 /* First calculate the maximum number of additional sides */
4607 for(i=1;i<=noelements;i++) {
4608
4609 /* owner partition cannot cause an indirect coupling */
4610 if(elempart[i] == part) continue;
4611
4612 elemtype = data->elementtypes[i];
4613 nodesd2 = elemtype%100;
4614
4615 /* Check how many nodes still belong to this partition,
4616 if more than one there may be indirect coupling. */
4617 for(j=0;j < nodesd2;j++) {
4618 elemhit[j] = FALSE;
4619 ind = data->topology[i][j];
4620 for(k=1;k<=neededtimes[ind];k++)
4621 if(part == data->partitiontable[k][ind]) elemhit[j] = TRUE;
4622 }
4623 bcneeded = 0;
4624 for(j=0;j < nodesd2;j++)
4625 if(elemhit[j]) bcneeded++;
4626 if(bcneeded <= 1) continue;
4627
4628 if(l == 0) {
4629 maxsides += (bcneeded-1)*bcneeded/2;
4630 }
4631 else {
4632 for(j=0;j < nodesd2;j++) {
4633 for(k=j+1;k < nodesd2;k++) {
4634 if(elemhit[j] && elemhit[k]) {
4635 nodesides += 1;
4636
4637 /* The minimum index always first */
4638 if(data->topology[i][j] <= data->topology[i][k]) {
4639 nodepairs[nodesides][1] = data->topology[i][j];
4640 nodepairs[nodesides][2] = data->topology[i][k];
4641 }
4642 else {
4643 nodepairs[nodesides][1] = data->topology[i][k];
4644 nodepairs[nodesides][2] = data->topology[i][j];
4645 }
4646 }
4647 }
4648 }
4649 }
4650 }
4651
4652 /* After first round allocate enough space to memorize all indirect non-element couplings. */
4653 if(l == 0) {
4654 nodepairs = Imatrix(1,maxsides,1,2);
4655 for(i=1;i<=maxsides;i++)
4656 nodepairs[i][1] = nodepairs[i][2] = 0;
4657 l++;
4658 goto findindirect;
4659 }
4660 if(0) printf("Number of non-element connections is %d\n",nodesides);
4661
4662
4663 for(i=1;i<=noknots;i++)
4664 nodeconnections[i] = 0;
4665
4666 for(i=1;i<=nodesides;i++)
4667 nodeconnections[nodepairs[i][1]] += 1;
4668
4669 maxnodeconnections = 0;
4670 for(i=1;i<=noknots;i++)
4671 maxnodeconnections = MAX(maxnodeconnections, nodeconnections[i]);
4672 if(0) printf("Maximum number of node-to-node connections %d\n",maxnodeconnections);
4673
4674 connectednodes = 0;
4675 for(i=1;i<=noknots;i++) {
4676 if(nodeconnections[i] > 0) {
4677 connectednodes++;
4678 nodeconnections[i] = connectednodes;
4679 }
4680 }
4681 if(0) printf("Number of nodes with non-element connections %d\n",connectednodes);
4682
4683 indpairs = Imatrix(1,connectednodes,1,maxnodeconnections);
4684 for(i=1;i<=connectednodes;i++)
4685 for(j=1;j<=maxnodeconnections;j++)
4686 indpairs[i][j] = 0;
4687
4688 for(i=1;i<=nodesides;i++) {
4689 ind = nodeconnections[nodepairs[i][1]];
4690 for(j=1;j<=maxnodeconnections;j++) {
4691 if(indpairs[ind][j] == 0) {
4692 indpairs[ind][j] = i;
4693 break;
4694 }
4695 }
4696 }
4697
4698 /* Remove duplicate connections */
4699 l = 0;
4700 for(i=1;i<=connectednodes;i++) {
4701 for(j=1;j<=maxnodeconnections;j++)
4702 for(k=j+1;k<=maxnodeconnections;k++) {
4703 ind = indpairs[i][j];
4704 ind2 = indpairs[i][k];
4705 if(!ind || !ind2) continue;
4706
4707 if(!nodepairs[ind][1] || !nodepairs[ind][2]) continue;
4708
4709 if(nodepairs[ind][2] == nodepairs[ind2][2]) {
4710 nodepairs[ind2][1] = nodepairs[ind2][2] = 0;
4711 l++;
4712 }
4713 }
4714 }
4715 if(0) printf("Removed %d duplicate connections\n",l);
4716
4717
4718 /* Remove connections that already exist */
4719 m = 0;
4720 for(i=1;i<=noelements;i++) {
4721 if(elempart[i] != part) continue;
4722
4723 elemtype = data->elementtypes[i];
4724 nodesd2 = elemtype%100;
4725
4726 for(j=0;j < nodesd2;j++) {
4727 ind = nodeconnections[data->topology[i][j]];
4728 if(!ind) continue;
4729
4730 for(k=0;k < nodesd2;k++) {
4731 if(j==k) continue;
4732
4733 for(l=1;l<=maxnodeconnections;l++) {
4734 ind2 = indpairs[ind][l];
4735 if(!ind2) break;
4736
4737 if(nodepairs[ind2][1] == data->topology[i][j] && nodepairs[ind2][2] == data->topology[i][k]) {
4738 nodepairs[ind2][1] = nodepairs[ind2][2] = 0;
4739 m++;
4740 }
4741 }
4742 }
4743 }
4744 }
4745 if(0) printf("Removed %d connections that already exists in other elements\n",m);
4746
4747 for(i=1; i <= nodesides; i++) {
4748 ind = nodepairs[i][1];
4749 ind2 = nodepairs[i][2];
4750 if(!ind || !ind2) continue;
4751 sumsides++;
4752
4753 sideelemtype = 102;
4754 if(reorder) {
4755 fprintf(out,"%d %d %d %d %d %d %d\n",
4756 sumsides,indirecttype,0,0,sideelemtype,order[ind],order[ind2]);
4757 } else {
4758 fprintf(out,"%d %d %d %d %d %d %d\n",
4759 sumsides,indirecttype,0,0,sideelemtype,ind,ind2);
4760 }
4761 sidetypes[sideelemtype] += 1;
4762 indirectinpart[part] += 1;
4763 }
4764
4765 /* Finally free some extra space that was allocated */
4766 free_Imatrix(indpairs,1,connectednodes,1,maxnodeconnections);
4767 free_Imatrix(nodepairs,1,maxsides,1,2);
4768 }
4769 /* End of indirect couplings */
4770
4771
4772 fclose(out);
4773 /*********** end of part.n.boundary *********************/
4774
4775
4776
4777
4778 /*********** start of part.n.header *********************/
4779 tottypes = 0;
4780 for(i=minelemtype;i<=maxelemtype;i++) {
4781 if(bulktypes[part][i]) tottypes++;
4782 if(sidetypes[i]) tottypes++;
4783 }
4784
4785 sprintf(filename,"%s.%d.%s","part",part,"header");
4786 out = fopen(filename,"w");
4787 fprintf(out,"%-6d %-6d %-6d\n",
4788 needednodes[part],elementsinpart[part],sumsides);
4789
4790 fprintf(out,"%-6d\n",tottypes);
4791 for(i=minelemtype;i<=maxelemtype;i++)
4792 if(bulktypes[part][i])
4793 fprintf(out,"%-6d %-6d\n",i,bulktypes[part][i]);
4794
4795 for(i=minelemtype;i<=maxelemtype;i++)
4796 if(sidetypes[i])
4797 fprintf(out,"%-6d %-6d\n",i,sidetypes[i]);
4798
4799 fprintf(out,"%-6d %-6d\n",neededtwice[part],0);
4800 fclose(out);
4801
4802 if(info) {
4803 if( indirect ) {
4804 if(part == 1) {
4805 printf(" %-5s %-10s %-10s %-8s %-8s %-8s\n",
4806 "part","elements","nodes","shared","bc elems","indirect");
4807 }
4808 printf(" %-5d %-10d %-10d %-8d %-8d %-8d\n",
4809 part,elementsinpart[part],ownnodes[part],sharednodes[part],sidesinpart[part],
4810 indirectinpart[part]);
4811 }
4812 else {
4813 if( part == 1 ) {
4814 printf(" %-5s %-10s %-10s %-8s %-8s\n",
4815 "part","elements","nodes","shared","bc elems");
4816 }
4817 printf(" %-5d %-10d %-10d %-8d %-8d\n",
4818 part,elementsinpart[part],ownnodes[part],sharednodes[part],sidesinpart[part]);
4819 }
4820 }
4821 }
4822 /*********** end of part.n.header *********************/
4823
4824
4825 sumelementsinpart = sumownnodes = sumsharednodes = sumsidesinpart = sumindirect = 0;
4826 for(i=1;i<=partitions;i++) {
4827 sumelementsinpart += elementsinpart[i];
4828 sumownnodes += ownnodes[i];
4829 sumsharednodes += sharednodes[i];
4830 sumsidesinpart += sidesinpart[i];
4831 sumindirect += indirectinpart[i];
4832 }
4833 n = partitions;
4834 printf("----------------------------------------------------------------------------------------------\n");
4835 printf(" ave %-10.1f %-10.1f %-8.1f %-8.1f %-8.1f\n",
4836 1.0*sumelementsinpart/n,1.0*sumownnodes/n,1.0*sumsharednodes/n,
4837 1.0*sumsidesinpart/n,1.0*sumindirect/n);
4838
4839 if( info && halobcs ) {
4840 printf("Number of boundary elements associated with halo: %d\n",halobcs);
4841 }
4842
4843
4844 #if 0
4845 {
4846 int noparents[2][100];
4847
4848 for(j=0;j<2;j++)
4849 for(i=0;i<100;i++)
4850 noparents[j][i] = 0;
4851
4852 for(j=0;j < MAXBOUNDARIES;j++) {
4853 if(bound[j].nosides == 0 ) continue;
4854 for(i=1; i <= bound[j].nosides; i++) {
4855 GetBoundaryElement(i,&bound[j],data,sideind,&sideelemtype);
4856 parent = bound[j].parent[i];
4857 parent2 = bound[j].parent2[i];
4858 bctype = bound[j].types[i];
4859
4860 k = 0;
4861 if(parent) k++;
4862 if(parent2) k++;
4863
4864 noparents[k][bctype] += 1;
4865 }
4866 }
4867
4868 printf("Number of BC parents\n");
4869 for(i=0;i<100;i++) {
4870 k = 0;
4871 for(j=0;j<2;j++)
4872 k = k + noparents[j][i];
4873 if( k ) {
4874 printf("BC %d: %d %d %d\n",i,noparents[0][i],noparents[1][i],noparents[2][i]);
4875 }
4876 }
4877 }
4878 #endif
4879
4880
4881
4882 if(splitsides && !anyparthalo) {
4883 printf("************************* Warning ****************************\n");
4884 printf("Number or boundary elements split at between parents: %d\n",splitsides);
4885 printf("This could be a problem for internal jump conditions\n");
4886 printf("You could try to use '-halobc' flag as remedy with ElmerSolver.\n");
4887 printf("**************************************************************\n");
4888 }
4889
4890 if(anyparthalo) free_Ivector(neededtimes2,1,noknots);
4891
4892 cdstat = chdir("..");
4893 cdstat = chdir("..");
4894
4895 if(reorder) free_Ivector(order,1,noknots);
4896 free_Ivector(needednodes,1,partitions);
4897 free_Ivector(neededtwice,1,partitions);
4898 free_Ivector(sharednodes,1,partitions);
4899 free_Ivector(ownnodes,1,partitions);
4900 free_Ivector(sidetypes,minelemtype,maxelemtype);
4901 free_Imatrix(bulktypes,1,partitions,minelemtype,maxelemtype);
4902
4903 if(info) printf("Writing of partitioned mesh finished\n");
4904
4905 return(0);
4906 }
4907
4908
4909 #if USE_METIS
ReorderElementsMetis(struct FemType * data,int info)4910 int ReorderElementsMetis(struct FemType *data,int info)
4911 /* Calls the fill reduction ordering algorithm of Metis library. */
4912 {
4913 int i,j,k,nn,totcon,maxcon,con,options[8];
4914 int noelements,noknots,nonodes;
4915 int *xadj,*adjncy,numflag,*iperm=NULL,*perm=NULL,**newtopology;
4916 Real *newx,*newy,*newz = NULL;
4917
4918 noelements = data->noelements;
4919 noknots = data->noknots;
4920
4921 if(info) printf("Reordering %d knots and %d elements using Metis reordering routine.\n",
4922 noknots,noelements);
4923 perm = NULL;
4924 i = CalculateIndexwidth(data,FALSE,perm);
4925 if(info) printf("Indexwidth of the original node order is %d.\n",i);
4926
4927
4928 CreateNodalGraph(data,TRUE,info);
4929 maxcon = data->nodalmaxconnections;
4930
4931 totcon = 0;
4932 for(i=1;i<=noknots;i++) {
4933 for(j=0;j<maxcon;j++) {
4934 con = data->nodalgraph[j][i];
4935 if(con) totcon++;
4936 }
4937 }
4938 if(info) printf("There are %d connections altogether\n",totcon);
4939
4940 xadj = Ivector(0,noknots);
4941 adjncy = Ivector(0,totcon-1);
4942 for(i=0;i<totcon;i++)
4943 adjncy[i] = 0;
4944
4945 totcon = 0;
4946 for(i=1;i<=noknots;i++) {
4947 xadj[i-1] = totcon;
4948 for(j=0;j<maxcon;j++) {
4949 con = data->nodalgraph[j][i];
4950 if(con) {
4951 adjncy[totcon] = con-1;
4952 totcon++;
4953 }
4954 }
4955 }
4956 xadj[noknots] = totcon;
4957
4958 nn = noknots;
4959 numflag = 0;
4960 for(i=0;i<8;i++) options[i] = 0;
4961 perm = Ivector(0,noknots-1);
4962 iperm = Ivector(0,noknots-1);
4963
4964 if(info) printf("Starting Metis reordering routine.\n");
4965
4966 METIS_NodeND(&nn,xadj,adjncy,&numflag,&options[0],perm,iperm);
4967
4968 if(info) printf("Finished Metis reordering routine.\n");
4969
4970 if(info) printf("Moving knots to new positions\n");
4971 newx = Rvector(1,data->noknots);
4972 newy = Rvector(1,data->noknots);
4973 newz = Rvector(1,data->noknots);
4974
4975 for(i=1;i<=data->noknots;i++) {
4976 newx[i] = data->x[perm[i-1]+1];
4977 newy[i] = data->y[perm[i-1]+1];
4978 newz[i] = data->z[perm[i-1]+1];
4979 }
4980
4981 free_Rvector(data->x,1,data->noknots);
4982 free_Rvector(data->y,1,data->noknots);
4983 free_Rvector(data->z,1,data->noknots);
4984
4985 data->x = newx;
4986 data->y = newy;
4987 data->z = newz;
4988
4989
4990 if(info) printf("Chanching the element topology\n");
4991
4992 newtopology = Imatrix(1,noelements,0,data->maxnodes-1);
4993
4994 for(j=1;j<=noelements;j++) {
4995 nonodes = data->elementtypes[j]%100;
4996 for(i=0;i<nonodes;i++) {
4997 k = data->topology[j][i];
4998 newtopology[j][i] = iperm[k-1]+1;
4999 }
5000 }
5001 free_Imatrix(data->topology,1,noelements,0,data->maxnodes-1);
5002 data->topology = newtopology;
5003
5004 i = CalculateIndexwidth(data,FALSE,perm);
5005 if(info) printf("Indexwidth of the new node order is %d.\n",i);
5006
5007 if(0) printf("Deallocating vectors needed for reordering.\n");
5008 free_Ivector(iperm,0,noknots-1);
5009 free_Ivector(perm,0,noknots-1);
5010
5011 return(0);
5012 }
5013 #endif
5014